Actual source code: vecio.c
2: /*
3: This file contains simple binary input routines for vectors. The
4: analogous output routines are within each vector implementation's
5: VecView (with viewer types PETSCVIEWERBINARY)
6: */
8: #include <petscsys.h>
9: #include <petscvec.h>
10: #include <petsc/private/vecimpl.h>
11: #include <petsc/private/viewerimpl.h>
12: #include <petsclayouthdf5.h>
14: PetscErrorCode VecView_Binary(Vec vec,PetscViewer viewer)
15: {
16: PetscErrorCode ierr;
17: PetscBool skipHeader;
18: PetscLayout map;
19: PetscInt tr[2],n,s,N;
20: const PetscScalar *xarray;
24: PetscViewerSetUp(viewer);
25: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
27: VecGetLayout(vec,&map);
28: PetscLayoutGetLocalSize(map,&n);
29: PetscLayoutGetRange(map,&s,NULL);
30: PetscLayoutGetSize(map,&N);
32: tr[0] = VEC_FILE_CLASSID; tr[1] = N;
33: if (!skipHeader) {PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT);}
35: VecGetArrayRead(vec,&xarray);
36: PetscViewerBinaryWriteAll(viewer,xarray,n,s,N,PETSC_SCALAR);
37: VecRestoreArrayRead(vec,&xarray);
39: { /* write to the viewer's .info file */
40: FILE *info;
41: PetscMPIInt rank;
42: PetscViewerFormat format;
43: const char *name,*pre;
45: PetscViewerBinaryGetInfoPointer(viewer,&info);
46: MPI_Comm_rank(PetscObjectComm((PetscObject)vec),&rank);
48: PetscViewerGetFormat(viewer,&format);
49: if (format == PETSC_VIEWER_BINARY_MATLAB) {
50: PetscObjectGetName((PetscObject)vec,&name);
51: if (rank == 0 && info) {
52: PetscFPrintf(PETSC_COMM_SELF,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
53: PetscFPrintf(PETSC_COMM_SELF,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
54: PetscFPrintf(PETSC_COMM_SELF,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
55: }
56: }
58: PetscObjectGetOptionsPrefix((PetscObject)vec,&pre);
59: if (rank == 0 && info) {PetscFPrintf(PETSC_COMM_SELF,info,"-%svecload_block_size %D\n",pre?pre:"",PetscAbs(vec->map->bs));}
60: }
61: return(0);
62: }
64: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
65: {
67: PetscBool skipHeader,flg;
68: PetscInt tr[2],rows,N,n,s,bs;
69: PetscScalar *array;
70: PetscLayout map;
73: PetscViewerSetUp(viewer);
74: PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);
76: VecGetLayout(vec,&map);
77: PetscLayoutGetSize(map,&N);
79: /* read vector header */
80: if (!skipHeader) {
81: PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);
82: if (tr[0] != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Not a vector next in file");
83: if (tr[1] < 0) SETERRQ1(PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Vector size (%D) in file is negative",tr[1]);
84: if (N >= 0 && N != tr[1]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Vector in file different size (%D) than input vector (%D)",tr[1],N);
85: rows = tr[1];
86: } else {
87: if (N < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Vector binary file header was skipped, thus the user must specify the global size of input vector");
88: rows = N;
89: }
91: /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
92: PetscLayoutGetBlockSize(map,&bs);
93: PetscOptionsGetInt(((PetscObject)viewer)->options,((PetscObject)vec)->prefix,"-vecload_block_size",&bs,&flg);
94: if (flg) {VecSetBlockSize(vec,bs);}
95: PetscLayoutGetLocalSize(map,&n);
96: if (N < 0) {VecSetSizes(vec,n,rows);}
97: VecSetUp(vec);
99: /* get vector sizes and check global size */
100: VecGetSize(vec,&N);
101: VecGetLocalSize(vec,&n);
102: VecGetOwnershipRange(vec,&s,NULL);
103: if (N != rows) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Vector in file different size (%D) than input vector (%D)",rows,N);
105: /* read vector values */
106: VecGetArray(vec,&array);
107: PetscViewerBinaryReadAll(viewer,array,n,s,N,PETSC_SCALAR);
108: VecRestoreArray(vec,&array);
109: return(0);
110: }
112: #if defined(PETSC_HAVE_HDF5)
113: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
114: {
115: hid_t scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
116: PetscScalar *x,*arr;
117: const char *vecname;
121: if (!((PetscObject)xin)->name) SETERRQ(PetscObjectComm((PetscObject)xin), PETSC_ERR_SUP, "Vec name must be set with PetscObjectSetName() before VecLoad()");
122: #if defined(PETSC_USE_REAL_SINGLE)
123: scalartype = H5T_NATIVE_FLOAT;
124: #elif defined(PETSC_USE_REAL___FLOAT128)
125: #error "HDF5 output with 128 bit floats not supported."
126: #elif defined(PETSC_USE_REAL___FP16)
127: #error "HDF5 output with 16 bit floats not supported."
128: #else
129: scalartype = H5T_NATIVE_DOUBLE;
130: #endif
131: PetscObjectGetName((PetscObject)xin, &vecname);
132: PetscViewerHDF5Load(viewer,vecname,xin->map,scalartype,(void**)&x);
133: VecSetUp(xin); /* VecSetSizes might have not been called so ensure ops->create has been called */
134: if (!xin->ops->replacearray) {
135: VecGetArray(xin,&arr);
136: PetscArraycpy(arr,x,xin->map->n);
137: PetscFree(x);
138: VecRestoreArray(xin,&arr);
139: } else {
140: VecReplaceArray(xin,x);
141: }
142: return(0);
143: }
144: #endif
146: #if defined(PETSC_HAVE_ADIOS)
147: #include <adios.h>
148: #include <adios_read.h>
149: #include <petsc/private/vieweradiosimpl.h>
150: #include <petsc/private/viewerimpl.h>
152: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
153: {
154: PetscViewer_ADIOS *adios = (PetscViewer_ADIOS*)viewer->data;
155: PetscErrorCode ierr;
156: PetscScalar *x;
157: PetscInt Nfile,N,rstart,n;
158: uint64_t N_t,rstart_t;
159: const char *vecname;
160: ADIOS_VARINFO *v;
161: ADIOS_SELECTION *sel;
164: PetscObjectGetName((PetscObject) xin, &vecname);
166: v = adios_inq_var(adios->adios_fp, vecname);
167: if (v->ndim != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%D)", v->ndim);
168: Nfile = (PetscInt) v->dims[0];
170: /* Set Vec sizes,blocksize,and type if not already set */
171: if ((xin)->map->n < 0 && (xin)->map->N < 0) {
172: VecSetSizes(xin, PETSC_DECIDE, Nfile);
173: }
174: /* If sizes and type already set,check if the vector global size is correct */
175: VecGetSize(xin, &N);
176: VecGetLocalSize(xin, &n);
177: if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);
179: VecGetOwnershipRange(xin,&rstart,NULL);
180: rstart_t = rstart;
181: N_t = n;
182: sel = adios_selection_boundingbox (v->ndim, &rstart_t, &N_t);
183: VecGetArray(xin,&x);
184: adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
185: adios_perform_reads (adios->adios_fp, 1);
186: VecRestoreArray(xin,&x);
187: adios_selection_delete(sel);
189: return(0);
190: }
191: #endif
193: PetscErrorCode VecLoad_Default(Vec newvec, PetscViewer viewer)
194: {
196: PetscBool isbinary;
197: #if defined(PETSC_HAVE_HDF5)
198: PetscBool ishdf5;
199: #endif
200: #if defined(PETSC_HAVE_ADIOS)
201: PetscBool isadios;
202: #endif
205: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
206: #if defined(PETSC_HAVE_HDF5)
207: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
208: #endif
209: #if defined(PETSC_HAVE_ADIOS)
210: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
211: #endif
213: #if defined(PETSC_HAVE_HDF5)
214: if (ishdf5) {
215: if (!((PetscObject)newvec)->name) {
216: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
217: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Since HDF5 format gives ASCII name for each object in file; must use VecLoad() after setting name of Vec with PetscObjectSetName()");
218: }
219: VecLoad_HDF5(newvec, viewer);
220: } else
221: #endif
222: #if defined(PETSC_HAVE_ADIOS)
223: if (isadios) {
224: if (!((PetscObject)newvec)->name) {
225: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
226: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Since ADIOS format gives ASCII name for each object in file; must use VecLoad() after setting name of Vec with PetscObjectSetName()");
227: }
228: VecLoad_ADIOS(newvec, viewer);
229: } else
230: #endif
231: {
232: VecLoad_Binary(newvec, viewer);
233: }
234: return(0);
235: }
237: /*@
238: VecChop - Set all values in the vector with an absolute value less than the tolerance to zero
240: Input Parameters:
241: + v - The vector
242: - tol - The zero tolerance
244: Output Parameters:
245: . v - The chopped vector
247: Level: intermediate
249: .seealso: VecCreate(), VecSet()
250: @*/
251: PetscErrorCode VecChop(Vec v, PetscReal tol)
252: {
253: PetscScalar *a;
254: PetscInt n, i;
258: VecGetLocalSize(v, &n);
259: VecGetArray(v, &a);
260: for (i = 0; i < n; ++i) {
261: if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
262: }
263: VecRestoreArray(v, &a);
264: return(0);
265: }