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:   PetscBool         skipHeader;
 17:   PetscLayout       map;
 18:   PetscInt          tr[2],n,s,N;
 19:   const PetscScalar *xarray;

 22:   PetscViewerSetUp(viewer);
 23:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);

 25:   VecGetLayout(vec,&map);
 26:   PetscLayoutGetLocalSize(map,&n);
 27:   PetscLayoutGetRange(map,&s,NULL);
 28:   PetscLayoutGetSize(map,&N);

 30:   tr[0] = VEC_FILE_CLASSID; tr[1] = N;
 31:   if (!skipHeader) PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT);

 33:   VecGetArrayRead(vec,&xarray);
 34:   PetscViewerBinaryWriteAll(viewer,xarray,n,s,N,PETSC_SCALAR);
 35:   VecRestoreArrayRead(vec,&xarray);

 37:   { /* write to the viewer's .info file */
 38:     FILE             *info;
 39:     PetscMPIInt       rank;
 40:     PetscViewerFormat format;
 41:     const char        *name,*pre;

 43:     PetscViewerBinaryGetInfoPointer(viewer,&info);
 44:     MPI_Comm_rank(PetscObjectComm((PetscObject)vec),&rank);

 46:     PetscViewerGetFormat(viewer,&format);
 47:     if (format == PETSC_VIEWER_BINARY_MATLAB) {
 48:       PetscObjectGetName((PetscObject)vec,&name);
 49:       if (rank == 0 && info) {
 50:         PetscFPrintf(PETSC_COMM_SELF,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
 51:         PetscFPrintf(PETSC_COMM_SELF,info,"#$$ Set.%s = PetscBinaryRead(fd);\n",name);
 52:         PetscFPrintf(PETSC_COMM_SELF,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
 53:       }
 54:     }

 56:     PetscObjectGetOptionsPrefix((PetscObject)vec,&pre);
 57:     if (rank == 0 && info) PetscFPrintf(PETSC_COMM_SELF,info,"-%svecload_block_size %" PetscInt_FMT "\n",pre?pre:"",PetscAbs(vec->map->bs));
 58:   }
 59:   return 0;
 60: }

 62: PetscErrorCode VecLoad_Binary(Vec vec, PetscViewer viewer)
 63: {
 64:   PetscBool      skipHeader,flg;
 65:   PetscInt       tr[2],rows,N,n,s,bs;
 66:   PetscScalar    *array;
 67:   PetscLayout    map;

 69:   PetscViewerSetUp(viewer);
 70:   PetscViewerBinaryGetSkipHeader(viewer,&skipHeader);

 72:   VecGetLayout(vec,&map);
 73:   PetscLayoutGetSize(map,&N);

 75:   /* read vector header */
 76:   if (!skipHeader) {
 77:     PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);
 81:     rows = tr[1];
 82:   } else {
 84:     rows = N;
 85:   }

 87:   /* set vector size, blocksize, and type if not already set; block size first so that local sizes will be compatible. */
 88:   PetscLayoutGetBlockSize(map,&bs);
 89:   PetscOptionsGetInt(((PetscObject)viewer)->options,((PetscObject)vec)->prefix,"-vecload_block_size",&bs,&flg);
 90:   if (flg) VecSetBlockSize(vec,bs);
 91:   PetscLayoutGetLocalSize(map,&n);
 92:   if (N < 0) VecSetSizes(vec,n,rows);
 93:   VecSetUp(vec);

 95:   /* get vector sizes and check global size */
 96:   VecGetSize(vec,&N);
 97:   VecGetLocalSize(vec,&n);
 98:   VecGetOwnershipRange(vec,&s,NULL);

101:   /* read vector values */
102:   VecGetArray(vec,&array);
103:   PetscViewerBinaryReadAll(viewer,array,n,s,N,PETSC_SCALAR);
104:   VecRestoreArray(vec,&array);
105:   return 0;
106: }

108: #if defined(PETSC_HAVE_HDF5)
109: PetscErrorCode VecLoad_HDF5(Vec xin, PetscViewer viewer)
110: {
111:   hid_t          scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
112:   PetscScalar    *x,*arr;
113:   const char     *vecname;

116: #if defined(PETSC_USE_REAL_SINGLE)
117:   scalartype = H5T_NATIVE_FLOAT;
118: #elif defined(PETSC_USE_REAL___FLOAT128)
119: #error "HDF5 output with 128 bit floats not supported."
120: #elif defined(PETSC_USE_REAL___FP16)
121: #error "HDF5 output with 16 bit floats not supported."
122: #else
123:   scalartype = H5T_NATIVE_DOUBLE;
124: #endif
125:   PetscObjectGetName((PetscObject)xin, &vecname);
126:   PetscViewerHDF5Load(viewer,vecname,xin->map,scalartype,(void**)&x);
127:   VecSetUp(xin); /* VecSetSizes might have not been called so ensure ops->create has been called */
128:   if (!xin->ops->replacearray) {
129:     VecGetArray(xin,&arr);
130:     PetscArraycpy(arr,x,xin->map->n);
131:     PetscFree(x);
132:     VecRestoreArray(xin,&arr);
133:   } else {
134:     VecReplaceArray(xin,x);
135:   }
136:   return 0;
137: }
138: #endif

140: #if defined(PETSC_HAVE_ADIOS)
141: #include <adios.h>
142: #include <adios_read.h>
143: #include <petsc/private/vieweradiosimpl.h>
144: #include <petsc/private/viewerimpl.h>

146: PetscErrorCode VecLoad_ADIOS(Vec xin, PetscViewer viewer)
147: {
148:   PetscViewer_ADIOS *adios = (PetscViewer_ADIOS*)viewer->data;
149:   PetscScalar       *x;
150:   PetscInt          Nfile,N,rstart,n;
151:   uint64_t          N_t,rstart_t;
152:   const char        *vecname;
153:   ADIOS_VARINFO     *v;
154:   ADIOS_SELECTION   *sel;

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

158:   v    = adios_inq_var(adios->adios_fp, vecname);
160:   Nfile = (PetscInt) v->dims[0];

162:   /* Set Vec sizes,blocksize,and type if not already set */
163:   if ((xin)->map->n < 0 && (xin)->map->N < 0) VecSetSizes(xin, PETSC_DECIDE, Nfile);
164:   /* If sizes and type already set,check if the vector global size is correct */
165:   VecGetSize(xin, &N);
166:   VecGetLocalSize(xin, &n);

169:   VecGetOwnershipRange(xin,&rstart,NULL);
170:   rstart_t = rstart;
171:   N_t  = n;
172:   sel  = adios_selection_boundingbox (v->ndim, &rstart_t, &N_t);
173:   VecGetArray(xin,&x);
174:   adios_schedule_read(adios->adios_fp, sel, vecname, 0, 1, x);
175:   adios_perform_reads (adios->adios_fp, 1);
176:   VecRestoreArray(xin,&x);
177:   adios_selection_delete(sel);

179:   return 0;
180: }
181: #endif

183: PetscErrorCode  VecLoad_Default(Vec newvec, PetscViewer viewer)
184: {
185:   PetscBool      isbinary;
186: #if defined(PETSC_HAVE_HDF5)
187:   PetscBool      ishdf5;
188: #endif
189: #if defined(PETSC_HAVE_ADIOS)
190:   PetscBool      isadios;
191: #endif

193:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
194: #if defined(PETSC_HAVE_HDF5)
195:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
196: #endif
197: #if defined(PETSC_HAVE_ADIOS)
198:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
199: #endif

201: #if defined(PETSC_HAVE_HDF5)
202:   if (ishdf5) {
203:     if (!((PetscObject)newvec)->name) {
204:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
205:       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()");
206:     }
207:     VecLoad_HDF5(newvec, viewer);
208:   } else
209: #endif
210: #if defined(PETSC_HAVE_ADIOS)
211:   if (isadios) {
212:     if (!((PetscObject)newvec)->name) {
213:       PetscLogEventEnd(VEC_Load,viewer,0,0,0);
214:       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()");
215:     }
216:     VecLoad_ADIOS(newvec, viewer);
217:   } else
218: #endif
219:   {
220:     VecLoad_Binary(newvec, viewer);
221:   }
222:   return 0;
223: }

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

228:   Input Parameters:
229: + v   - The vector
230: - tol - The zero tolerance

232:   Output Parameters:
233: . v - The chopped vector

235:   Level: intermediate

237: .seealso: VecCreate(), VecSet()
238: @*/
239: PetscErrorCode VecChop(Vec v, PetscReal tol)
240: {
241:   PetscScalar    *a;
242:   PetscInt       n, i;

244:   VecGetLocalSize(v, &n);
245:   VecGetArray(v, &a);
246:   for (i = 0; i < n; ++i) {
247:     if (PetscAbsScalar(a[i]) < tol) a[i] = 0.0;
248:   }
249:   VecRestoreArray(v, &a);
250:   return 0;
251: }