Actual source code: vecio.c

petsc-3.12.2 2019-11-22
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 <petsclayouthdf5.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,NULL,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,NULL,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,*arr;
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); /* VecSetSizes might have not been called so ensure ops->create has been called */
169:   if (!xin->ops->replacearray) {
170:     VecGetArray(xin,&arr);
171:     PetscArraycpy(arr,x,xin->map->n);
172:     PetscFree(x);
173:     VecRestoreArray(xin,&arr);
174:   } else {
175:     VecReplaceArray(xin,x);
176:   }
177:   return(0);
178: }
179: #endif

181: #if defined(PETSC_HAVE_ADIOS)
182: #include <adios.h>
183: #include <adios_read.h>
184:  #include <petsc/private/vieweradiosimpl.h>
185:  #include <petsc/private/viewerimpl.h>

187: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
188: {
189:   PetscViewer_ADIOS *adios = (PetscViewer_ADIOS*)viewer->data;
190:   PetscErrorCode    ierr;
191:   PetscScalar       *x;
192:   PetscInt          Nfile,N,rstart,n;
193:   uint64_t          N_t,rstart_t;
194:   const char        *vecname;
195:   ADIOS_VARINFO     *v;
196:   ADIOS_SELECTION   *sel;

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

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

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

224:   return(0);
225: }
226: #endif

228: #if defined(PETSC_HAVE_ADIOS2)
229: #include <adios2_c.h>
230: #include <petsc/private/vieweradios2impl.h>
231:  #include <petsc/private/viewerimpl.h>

233: PetscErrorCode VecLoad_ADIOS2(Vec xin, PetscViewer viewer)
234: {
235:   PetscViewer_ADIOS2 *adios2 = (PetscViewer_ADIOS2*)viewer->data;
236:   PetscErrorCode     ierr;
237:   PetscScalar        *x;
238:   PetscInt           Nfile,N,rstart,n;
239:   const char         *vecname;

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

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

253:   VecGetOwnershipRange(xin,&rstart,NULL);
254:   VecGetArray(xin,&x);
255:   VecRestoreArray(xin,&x);
256:   return(0);
257: }
258: #endif

260: PetscErrorCode  VecLoad_Default(Vec newvec, PetscViewer viewer)
261: {
263:   PetscBool      isbinary;
264: #if defined(PETSC_HAVE_HDF5)
265:   PetscBool      ishdf5;
266: #endif
267: #if defined(PETSC_HAVE_ADIOS)
268:   PetscBool      isadios;
269: #endif
270: #if defined(PETSC_HAVE_ADIOS2)
271:   PetscBool      isadios2;
272: #endif

275:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
276: #if defined(PETSC_HAVE_HDF5)
277:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
278: #endif
279: #if defined(PETSC_HAVE_ADIOS)
280:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
281: #endif
282: #if defined(PETSC_HAVE_ADIOS2)
283:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS2,&isadios2);
284: #endif

286: #if defined(PETSC_HAVE_HDF5)
287:   if (ishdf5) {
288:     if (!((PetscObject)newvec)->name) {
289:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
290:       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()");
291:     }
292:     VecLoad_HDF5(newvec, viewer);
293:   } else
294: #endif
295: #if defined(PETSC_HAVE_ADIOS)
296:   if (isadios) {
297:     if (!((PetscObject)newvec)->name) {
298:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
299:       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()");
300:     }
301:     VecLoad_ADIOS(newvec, viewer);
302:   } else
303: #endif
304: #if defined(PETSC_HAVE_ADIOS2)
305:   if (isadios2) {
306:     if (!((PetscObject)newvec)->name) {
307:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
308:       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()");
309:     }
310:     VecLoad_ADIOS2(newvec, viewer);
311:   } else
312: #endif
313:   {
314:     VecLoad_Binary(newvec, viewer);
315:   }
316:   return(0);
317: }

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

322:   Input Parameters:
323: + v   - The vector
324: - tol - The zero tolerance

326:   Output Parameters:
327: . v - The chopped vector

329:   Level: intermediate

331: .seealso: VecCreate(), VecSet()
332: @*/
333: PetscErrorCode VecChop(Vec v, PetscReal tol)
334: {
335:   PetscScalar    *a;
336:   PetscInt       n, i;

340:   VecGetLocalSize(v, &n);
341:   VecGetArray(v, &a);
342:   for (i = 0; i < n; ++i) {
343:     if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
344:   }
345:   VecRestoreArray(v, &a);
346:   return(0);
347: }