Actual source code: vecio.c

petsc-3.11.2 2019-05-18
Report Typos and Errors

  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 <petscviewerhdf5.h>

 14: static PetscErrorCode PetscViewerBinaryReadVecHeader_Private(PetscViewer viewer,PetscInt *rows)
 15: {
 17:   MPI_Comm       comm;
 18:   PetscInt       tr[2],type;

 21:   PetscObjectGetComm((PetscObject)viewer,&comm);
 22:   /* Read vector header */
 23:   PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);
 24:   type = tr[0];
 25:   if (type != VEC_FILE_CLASSID) {
 26:     PetscLogEventEnd(VEC_Load,viewer,0,0,0);
 27:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a vector next in file");
 28:   }
 29:   *rows = tr[1];
 30:   return(0);
 31: }

 33: #if defined(PETSC_HAVE_MPIIO)
 34: static PetscErrorCode VecLoad_Binary_MPIIO(Vec vec, PetscViewer viewer)
 35: {
 37:   PetscMPIInt    lsize;
 38:   PetscScalar    *avec;
 39:   MPI_File       mfdes;
 40:   MPI_Offset     off;

 43:   VecGetArray(vec,&avec);
 44:   PetscMPIIntCast(vec->map->n,&lsize);

 46:   PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
 47:   PetscViewerBinaryGetMPIIOOffset(viewer,&off);
 48:   off += vec->map->rstart*sizeof(PetscScalar);
 49:   MPI_File_set_view(mfdes,off,MPIU_SCALAR,MPIU_SCALAR,(char*)"native",MPI_INFO_NULL);
 50:   MPIU_File_read_all(mfdes,avec,lsize,MPIU_SCALAR,MPI_STATUS_IGNORE);
 51:   PetscViewerBinaryAddMPIIOOffset(viewer,vec->map->N*sizeof(PetscScalar));

 53:   VecRestoreArray(vec,&avec);
 54:   VecAssemblyBegin(vec);
 55:   VecAssemblyEnd(vec);
 56:   return(0);
 57: }
 58: #endif

 60: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
 61: {
 62:   PetscMPIInt    size,rank,tag;
 63:   int            fd;
 64:   PetscInt       i,rows = 0,n,*range,N,bs;
 66:   PetscBool      flag,skipheader;
 67:   PetscScalar    *avec,*avecwork;
 68:   MPI_Comm       comm;
 69:   MPI_Request    request;
 70:   MPI_Status     status;
 71: #if defined(PETSC_HAVE_MPIIO)
 72:   PetscBool      useMPIIO;
 73: #endif

 76:   /* force binary viewer to load .info file if it has not yet done so */
 77:   PetscViewerSetUp(viewer);
 78:   PetscObjectGetComm((PetscObject)viewer,&comm);
 79:   MPI_Comm_rank(comm,&rank);
 80:   MPI_Comm_size(comm,&size);

 82:   PetscViewerBinaryGetDescriptor(viewer,&fd);
 83:   PetscViewerBinaryGetSkipHeader(viewer,&skipheader);
 84:   if (!skipheader) {
 85:     PetscViewerBinaryReadVecHeader_Private(viewer,&rows);
 86:   } else {
 87:     VecType vtype;
 88:     VecGetType(vec,&vtype);
 89:     if (!vtype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER, "Vector binary file header was skipped, thus the user must specify the type and length of input vector");
 90:     VecGetSize(vec,&N);
 91:     if (!N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER, "Vector binary file header was skipped, thus the user must specify the length of input vector");
 92:     rows = N;
 93:   }
 94:   /* Set Vec sizes,blocksize,and type if not already set. Block size first so that local sizes will be compatible. */
 95:   PetscOptionsGetInt(NULL,((PetscObject)vec)->prefix, "-vecload_block_size", &bs, &flag);
 96:   if (flag) {
 97:     VecSetBlockSize(vec, bs);
 98:   }
 99:   if (vec->map->n < 0 && vec->map->N < 0) {
100:     VecSetSizes(vec,PETSC_DECIDE,rows);
101:   }

103:   /* If sizes and type already set,check if the vector global size is correct */
104:   VecGetSize(vec, &N);
105:   if (N != rows) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", rows, N);

107: #if defined(PETSC_HAVE_MPIIO)
108:   PetscViewerBinaryGetUseMPIIO(viewer,&useMPIIO);
109:   if (useMPIIO) {
110:     VecLoad_Binary_MPIIO(vec, viewer);
111:     return(0);
112:   }
113: #endif

115:   VecGetLocalSize(vec,&n);
116:   PetscObjectGetNewTag((PetscObject)viewer,&tag);
117:   VecGetArray(vec,&avec);
118:   if (!rank) {
119:     PetscBinaryRead(fd,avec,n,PETSC_SCALAR);

121:     if (size > 1) {
122:       /* read in other chuncks and send to other processors */
123:       /* determine maximum chunck owned by other */
124:       range = vec->map->range;
125:       n = 1;
126:       for (i=1; i<size; i++) n = PetscMax(n,range[i+1] - range[i]);

128:       PetscMalloc1(n,&avecwork);
129:       for (i=1; i<size; i++) {
130:         n    = range[i+1] - range[i];
131:         PetscBinaryRead(fd,avecwork,n,PETSC_SCALAR);
132:         MPI_Isend(avecwork,n,MPIU_SCALAR,i,tag,comm,&request);
133:         MPI_Wait(&request,&status);
134:       }
135:       PetscFree(avecwork);
136:     }
137:   } else {
138:     MPI_Recv(avec,n,MPIU_SCALAR,0,tag,comm,&status);
139:   }

141:   VecRestoreArray(vec,&avec);
142:   VecAssemblyBegin(vec);
143:   VecAssemblyEnd(vec);
144:   return(0);
145: }

147: #if defined(PETSC_HAVE_HDF5)
148: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
149: {
150:   hid_t          scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
151:   PetscScalar    *x;
152:   const char     *vecname;

156:   if (!((PetscObject)xin)->name) SETERRQ(PetscObjectComm((PetscObject)xin), PETSC_ERR_SUP, "Vec name must be set with PetscObjectSetName() before VecLoad()");
157: #if defined(PETSC_USE_REAL_SINGLE)
158:   scalartype = H5T_NATIVE_FLOAT;
159: #elif defined(PETSC_USE_REAL___FLOAT128)
160: #error "HDF5 output with 128 bit floats not supported."
161: #elif defined(PETSC_USE_REAL___FP16)
162: #error "HDF5 output with 16 bit floats not supported."
163: #else
164:   scalartype = H5T_NATIVE_DOUBLE;
165: #endif
166:   PetscObjectGetName((PetscObject)xin, &vecname);
167:   PetscViewerHDF5Load(viewer, vecname, xin->map, scalartype, (void**)&x);
168:   VecSetUp(xin);
169:   VecReplaceArray(xin, x);
170:   return(0);
171: }
172: #endif

174: #if defined(PETSC_HAVE_ADIOS)
175: #include <adios.h>
176: #include <adios_read.h>
177:  #include <petsc/private/vieweradiosimpl.h>
178:  #include <petsc/private/viewerimpl.h>

180: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
181: {
182:   PetscViewer_ADIOS *adios = (PetscViewer_ADIOS*)viewer->data;
183:   PetscErrorCode    ierr;
184:   PetscScalar       *x;
185:   PetscInt          Nfile,N,rstart,n;
186:   uint64_t          N_t,rstart_t;
187:   const char        *vecname;
188:   ADIOS_VARINFO     *v;
189:   ADIOS_SELECTION   *sel;

192:   PetscObjectGetName((PetscObject) xin, &vecname);

194:   v    = adios_inq_var(adios->adios_fp, vecname);
195:   if (v->ndim != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file is not of dimension 1 (%D)", v->ndim);
196:   Nfile = (PetscInt) v->dims[0];

198:   /* Set Vec sizes,blocksize,and type if not already set */
199:   if ((xin)->map->n < 0 && (xin)->map->N < 0) {
200:     VecSetSizes(xin, PETSC_DECIDE, Nfile);
201:   }
202:   /* If sizes and type already set,check if the vector global size is correct */
203:   VecGetSize(xin, &N);
204:   VecGetLocalSize(xin, &n);
205:   if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);
206: 
207:   VecGetOwnershipRange(xin,&rstart,NULL);
208:   rstart_t = rstart;
209:   N_t  = n;
210:   sel  = adios_selection_boundingbox (v->ndim, &rstart_t, &N_t);
211:   VecGetArray(xin,&x);
212:   adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
213:   adios_perform_reads (adios->adios_fp, 1);
214:   VecRestoreArray(xin,&x);
215:   adios_selection_delete(sel);

217:   return(0);
218: }
219: #endif

221: #if defined(PETSC_HAVE_ADIOS2)
222: #include <adios2_c.h>
223: #include <petsc/private/vieweradios2impl.h>
224:  #include <petsc/private/viewerimpl.h>

226: PetscErrorCode VecLoad_ADIOS2(Vec xin, PetscViewer viewer)
227: {
228:   PetscViewer_ADIOS2 *adios2 = (PetscViewer_ADIOS2*)viewer->data;
229:   PetscErrorCode     ierr;
230:   PetscScalar        *x;
231:   PetscInt           Nfile,N,rstart,n;
232:   const char         *vecname;

235:   PetscObjectGetName((PetscObject) xin, &vecname);

237:   /* Set Vec sizes,blocksize,and type if not already set */
238:   if ((xin)->map->n < 0 && (xin)->map->N < 0) {
239:     VecSetSizes(xin, PETSC_DECIDE, Nfile);
240:   }
241:   /* If sizes and type already set,check if the vector global size is correct */
242:   VecGetSize(xin, &N);
243:   VecGetLocalSize(xin, &n);
244:   if (N != Nfile) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Vector in file different length (%D) then input vector (%D)", Nfile, N);

246:   VecGetOwnershipRange(xin,&rstart,NULL);
247:   VecGetArray(xin,&x);
248:   VecRestoreArray(xin,&x);
249:   return(0);
250: }
251: #endif

253: PetscErrorCode  VecLoad_Default(Vec newvec, PetscViewer viewer)
254: {
256:   PetscBool      isbinary;
257: #if defined(PETSC_HAVE_HDF5)
258:   PetscBool      ishdf5;
259: #endif
260: #if defined(PETSC_HAVE_ADIOS)
261:   PetscBool      isadios;
262: #endif
263: #if defined(PETSC_HAVE_ADIOS2)
264:   PetscBool      isadios2;
265: #endif

268:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
269: #if defined(PETSC_HAVE_HDF5)
270:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
271: #endif
272: #if defined(PETSC_HAVE_ADIOS)
273:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
274: #endif
275: #if defined(PETSC_HAVE_ADIOS2)
276:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS2,&isadios2);
277: #endif

279: #if defined(PETSC_HAVE_HDF5)
280:   if (ishdf5) {
281:     if (!((PetscObject)newvec)->name) {
282:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
283:       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()");
284:     }
285:     VecLoad_HDF5(newvec, viewer);
286:   } else
287: #endif
288: #if defined(PETSC_HAVE_ADIOS)
289:   if (isadios) {
290:     if (!((PetscObject)newvec)->name) {
291:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
292:       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()");
293:     }
294:     VecLoad_ADIOS(newvec, viewer);
295:   } else
296: #endif
297: #if defined(PETSC_HAVE_ADIOS2)
298:   if (isadios2) {
299:     if (!((PetscObject)newvec)->name) {
300:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
301:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Since ADIOS2 format gives ASCII name for each object in file; must use VecLoad() after setting name of Vec with PetscObjectSetName()");
302:     }
303:     VecLoad_ADIOS2(newvec, viewer);
304:   } else
305: #endif
306:   {
307:     VecLoad_Binary(newvec, viewer);
308:   }
309:   return(0);
310: }

312: /*@
313:   VecChop - Set all values in the vector with an absolute value less than the tolerance to zero

315:   Input Parameters:
316: + v   - The vector
317: - tol - The zero tolerance

319:   Output Parameters:
320: . v - The chopped vector

322:   Level: intermediate

324: .seealso: VecCreate(), VecSet()
325: @*/
326: PetscErrorCode VecChop(Vec v, PetscReal tol)
327: {
328:   PetscScalar    *a;
329:   PetscInt       n, i;

333:   VecGetLocalSize(v, &n);
334:   VecGetArray(v, &a);
335:   for (i = 0; i < n; ++i) {
336:     if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
337:   }
338:   VecRestoreArray(v, &a);
339:   return(0);
340: }