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: }