Actual source code: gr2.c

petsc-master 2019-12-13
Report Typos and Errors

  2: /*
  3:    Plots vectors obtained with DMDACreate2d()
  4: */

  6:  #include <petsc/private/dmdaimpl.h>
  7:  #include <petsc/private/glvisvecimpl.h>
  8:  #include <petsc/private/viewerimpl.h>
  9: #include <petsc/private/viewerhdf5impl.h>
 10:  #include <petscdraw.h>

 12: /*
 13:         The data that is passed into the graphics callback
 14: */
 15: typedef struct {
 16:   PetscMPIInt       rank;
 17:   PetscInt          m,n,dof,k;
 18:   PetscReal         xmin,xmax,ymin,ymax,min,max;
 19:   const PetscScalar *xy,*v;
 20:   PetscBool         showaxis,showgrid;
 21:   const char        *name0,*name1;
 22: } ZoomCtx;

 24: /*
 25:        This does the drawing for one particular field
 26:     in one particular set of coordinates. It is a callback
 27:     called from PetscDrawZoom()
 28: */
 29: PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx)
 30: {
 31:   ZoomCtx           *zctx = (ZoomCtx*)ctx;
 32:   PetscErrorCode    ierr;
 33:   PetscInt          m,n,i,j,k,dof,id,c1,c2,c3,c4;
 34:   PetscReal         min,max,x1,x2,x3,x4,y_1,y2,y3,y4;
 35:   const PetscScalar *xy,*v;

 38:   m    = zctx->m;
 39:   n    = zctx->n;
 40:   dof  = zctx->dof;
 41:   k    = zctx->k;
 42:   xy   = zctx->xy;
 43:   v    = zctx->v;
 44:   min  = zctx->min;
 45:   max  = zctx->max;

 47:   /* PetscDraw the contour plot patch */
 48:   PetscDrawCollectiveBegin(draw);
 49:   for (j=0; j<n-1; j++) {
 50:     for (i=0; i<m-1; i++) {
 51:       id   = i+j*m;
 52:       x1   = PetscRealPart(xy[2*id]);
 53:       y_1  = PetscRealPart(xy[2*id+1]);
 54:       c1   = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max);

 56:       id   = i+j*m+1;
 57:       x2   = PetscRealPart(xy[2*id]);
 58:       y2   = PetscRealPart(xy[2*id+1]);
 59:       c2   = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max);

 61:       id   = i+j*m+1+m;
 62:       x3   = PetscRealPart(xy[2*id]);
 63:       y3   = PetscRealPart(xy[2*id+1]);
 64:       c3   = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max);

 66:       id   = i+j*m+m;
 67:       x4   = PetscRealPart(xy[2*id]);
 68:       y4   = PetscRealPart(xy[2*id+1]);
 69:       c4   = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max);

 71:       PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);
 72:       PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);
 73:       if (zctx->showgrid) {
 74:         PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);
 75:         PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);
 76:         PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);
 77:         PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);
 78:       }
 79:     }
 80:   }
 81:   if (zctx->showaxis && !zctx->rank) {
 82:     if (zctx->name0 || zctx->name1) {
 83:       PetscReal xl,yl,xr,yr,x,y;
 84:       PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);
 85:       x  = xl + .30*(xr - xl);
 86:       xl = xl + .01*(xr - xl);
 87:       y  = yr - .30*(yr - yl);
 88:       yl = yl + .01*(yr - yl);
 89:       if (zctx->name0) {PetscDrawString(draw,x,yl,PETSC_DRAW_BLACK,zctx->name0);}
 90:       if (zctx->name1) {PetscDrawStringVertical(draw,xl,y,PETSC_DRAW_BLACK,zctx->name1);}
 91:     }
 92:     /*
 93:        Ideally we would use the PetscDrawAxis object to manage displaying the coordinate limits
 94:        but that may require some refactoring.
 95:     */
 96:     {
 97:       double xmin = (double)zctx->xmin, ymin = (double)zctx->ymin;
 98:       double xmax = (double)zctx->xmax, ymax = (double)zctx->ymax;
 99:       char   value[16]; size_t len; PetscReal w;
100:       PetscSNPrintf(value,16,"%0.2e",xmin);
101:       PetscDrawString(draw,xmin,ymin - .05*(ymax - ymin),PETSC_DRAW_BLACK,value);
102:       PetscSNPrintf(value,16,"%0.2e",xmax);
103:       PetscStrlen(value,&len);
104:       PetscDrawStringGetSize(draw,&w,NULL);
105:       PetscDrawString(draw,xmax - len*w,ymin - .05*(ymax - ymin),PETSC_DRAW_BLACK,value);
106:       PetscSNPrintf(value,16,"%0.2e",ymin);
107:       PetscDrawString(draw,xmin - .05*(xmax - xmin),ymin,PETSC_DRAW_BLACK,value);
108:       PetscSNPrintf(value,16,"%0.2e",ymax);
109:       PetscDrawString(draw,xmin - .05*(xmax - xmin),ymax,PETSC_DRAW_BLACK,value);
110:     }
111:   }
112:   PetscDrawCollectiveEnd(draw);
113:   return(0);
114: }

116: PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer)
117: {
118:   DM                 da,dac,dag;
119:   PetscErrorCode     ierr;
120:   PetscInt           N,s,M,w,ncoors = 4;
121:   const PetscInt     *lx,*ly;
122:   PetscReal          coors[4];
123:   PetscDraw          draw,popup;
124:   PetscBool          isnull,useports = PETSC_FALSE;
125:   MPI_Comm           comm;
126:   Vec                xlocal,xcoor,xcoorl;
127:   DMBoundaryType     bx,by;
128:   DMDAStencilType    st;
129:   ZoomCtx            zctx;
130:   PetscDrawViewPorts *ports = NULL;
131:   PetscViewerFormat  format;
132:   PetscInt           *displayfields;
133:   PetscInt           ndisplayfields,i,nbounds;
134:   const PetscReal    *bounds;

137:   zctx.showgrid = PETSC_FALSE;
138:   zctx.showaxis = PETSC_TRUE;

140:   PetscViewerDrawGetDraw(viewer,0,&draw);
141:   PetscDrawIsNull(draw,&isnull);
142:   if (isnull) return(0);

144:   PetscViewerDrawGetBounds(viewer,&nbounds,&bounds);

146:   VecGetDM(xin,&da);
147:   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");

149:   PetscObjectGetComm((PetscObject)xin,&comm);
150:   MPI_Comm_rank(comm,&zctx.rank);

152:   DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);
153:   DMDAGetOwnershipRanges(da,&lx,&ly,NULL);

155:   /*
156:      Obtain a sequential vector that is going to contain the local values plus ONE layer of
157:      ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will
158:      update the local values pluse ONE layer of ghost values.
159:   */
160:   PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);
161:   if (!xlocal) {
162:     if (bx !=  DM_BOUNDARY_NONE || by !=  DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) {
163:       /*
164:          if original da is not of stencil width one, or periodic or not a box stencil then
165:          create a special DMDA to handle one level of ghost points for graphics
166:       */
167:       DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);
168:       DMSetUp(dac);
169:       PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");
170:     } else {
171:       /* otherwise we can use the da we already have */
172:       dac = da;
173:     }
174:     /* create local vector for holding ghosted values used in graphics */
175:     DMCreateLocalVector(dac,&xlocal);
176:     if (dac != da) {
177:       /* don't keep any public reference of this DMDA, is is only available through xlocal */
178:       PetscObjectDereference((PetscObject)dac);
179:     } else {
180:       /* remove association between xlocal and da, because below we compose in the opposite
181:          direction and if we left this connect we'd get a loop, so the objects could
182:          never be destroyed */
183:       PetscObjectRemoveReference((PetscObject)xlocal,"__PETSc_dm");
184:     }
185:     PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);
186:     PetscObjectDereference((PetscObject)xlocal);
187:   } else {
188:     if (bx !=  DM_BOUNDARY_NONE || by !=  DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) {
189:       VecGetDM(xlocal, &dac);
190:     } else {
191:       dac = da;
192:     }
193:   }

195:   /*
196:       Get local (ghosted) values of vector
197:   */
198:   DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);
199:   DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);
200:   VecGetArrayRead(xlocal,&zctx.v);

202:   /*
203:       Get coordinates of nodes
204:   */
205:   DMGetCoordinates(da,&xcoor);
206:   if (!xcoor) {
207:     DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);
208:     DMGetCoordinates(da,&xcoor);
209:   }

211:   /*
212:       Determine the min and max coordinates in plot
213:   */
214:   VecStrideMin(xcoor,0,NULL,&zctx.xmin);
215:   VecStrideMax(xcoor,0,NULL,&zctx.xmax);
216:   VecStrideMin(xcoor,1,NULL,&zctx.ymin);
217:   VecStrideMax(xcoor,1,NULL,&zctx.ymax);
218:   PetscOptionsGetBool(NULL,NULL,"-draw_contour_axis",&zctx.showaxis,NULL);
219:   if (zctx.showaxis) {
220:     coors[0] = zctx.xmin - .05*(zctx.xmax - zctx.xmin); coors[1] = zctx.ymin - .05*(zctx.ymax - zctx.ymin);
221:     coors[2] = zctx.xmax + .05*(zctx.xmax - zctx.xmin); coors[3] = zctx.ymax + .05*(zctx.ymax - zctx.ymin);
222:   } else {
223:     coors[0] = zctx.xmin; coors[1] = zctx.ymin; coors[2] = zctx.xmax; coors[3] = zctx.ymax;
224:   }
225:   PetscOptionsGetRealArray(NULL,NULL,"-draw_coordinates",coors,&ncoors,NULL);
226:   PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %g %g %g %g\n",(double)coors[0],(double)coors[1],(double)coors[2],(double)coors[3]);

228:   /*
229:       Get local ghosted version of coordinates
230:   */
231:   PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);
232:   if (!xcoorl) {
233:     /* create DMDA to get local version of graphics */
234:     DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);
235:     DMSetUp(dag);
236:     PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");
237:     DMCreateLocalVector(dag,&xcoorl);
238:     PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);
239:     PetscObjectDereference((PetscObject)dag);
240:     PetscObjectDereference((PetscObject)xcoorl);
241:   } else {
242:     VecGetDM(xcoorl,&dag);
243:   }
244:   DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);
245:   DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);
246:   VecGetArrayRead(xcoorl,&zctx.xy);
247:   DMDAGetCoordinateName(da,0,&zctx.name0);
248:   DMDAGetCoordinateName(da,1,&zctx.name1);

250:   /*
251:       Get information about size of area each processor must do graphics for
252:   */
253:   DMDAGetInfo(dac,NULL,&M,&N,NULL,NULL,NULL,NULL,&zctx.dof,NULL,&bx,&by,NULL,NULL);
254:   DMDAGetGhostCorners(dac,NULL,NULL,NULL,&zctx.m,&zctx.n,NULL);
255:   PetscOptionsGetBool(NULL,NULL,"-draw_contour_grid",&zctx.showgrid,NULL);

257:   DMDASelectFields(da,&ndisplayfields,&displayfields);
258:   PetscViewerGetFormat(viewer,&format);
259:   PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);
260:   if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
261:   if (useports) {
262:     PetscViewerDrawGetDraw(viewer,0,&draw);
263:     PetscDrawCheckResizedWindow(draw);
264:     PetscDrawClear(draw);
265:     PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);
266:   }

268:   /*
269:       Loop over each field; drawing each in a different window
270:   */
271:   for (i=0; i<ndisplayfields; i++) {
272:     zctx.k = displayfields[i];

274:     /* determine the min and max value in plot */
275:     VecStrideMin(xin,zctx.k,NULL,&zctx.min);
276:     VecStrideMax(xin,zctx.k,NULL,&zctx.max);
277:     if (zctx.k < nbounds) {
278:       zctx.min = bounds[2*zctx.k];
279:       zctx.max = bounds[2*zctx.k+1];
280:     }
281:     if (zctx.min == zctx.max) {
282:       zctx.min -= 1.e-12;
283:       zctx.max += 1.e-12;
284:     }
285:     PetscInfo2(da,"DMDA 2d contour plot min %g max %g\n",(double)zctx.min,(double)zctx.max);

287:     if (useports) {
288:       PetscDrawViewPortsSet(ports,i);
289:     } else {
290:       const char *title;
291:       PetscViewerDrawGetDraw(viewer,i,&draw);
292:       DMDAGetFieldName(da,zctx.k,&title);
293:       if (title) {PetscDrawSetTitle(draw,title);}
294:     }

296:     PetscDrawGetPopup(draw,&popup);
297:     PetscDrawScalePopup(popup,zctx.min,zctx.max);
298:     PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);
299:     PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);
300:     if (!useports) {PetscDrawSave(draw);}
301:   }
302:   if (useports) {
303:     PetscViewerDrawGetDraw(viewer,0,&draw);
304:     PetscDrawSave(draw);
305:   }

307:   PetscDrawViewPortsDestroy(ports);
308:   PetscFree(displayfields);
309:   VecRestoreArrayRead(xcoorl,&zctx.xy);
310:   VecRestoreArrayRead(xlocal,&zctx.v);
311:   return(0);
312: }

314: #if defined(PETSC_HAVE_HDF5)
315: static PetscErrorCode VecGetHDF5ChunkSize(DM_DA *da, Vec xin, PetscInt dimension, PetscInt timestep, hsize_t *chunkDims)
316: {
317:   PetscMPIInt    comm_size;
319:   hsize_t        chunk_size, target_size, dim;
320:   hsize_t        vec_size = sizeof(PetscScalar)*da->M*da->N*da->P*da->w;
321:   hsize_t        avg_local_vec_size,KiB = 1024,MiB = KiB*KiB,GiB = MiB*KiB,min_size = MiB;
322:   hsize_t        max_chunks = 64*KiB;                                              /* HDF5 internal limitation */
323:   hsize_t        max_chunk_size = 4*GiB;                                           /* HDF5 internal limitation */
324:   PetscInt       zslices=da->p, yslices=da->n, xslices=da->m;

327:   MPI_Comm_size(PetscObjectComm((PetscObject)xin), &comm_size);
328:   avg_local_vec_size = (hsize_t) ceil(vec_size*1.0/comm_size);      /* we will attempt to use this as the chunk size */

330:   target_size = (hsize_t) ceil(PetscMin(vec_size,PetscMin(max_chunk_size,PetscMax(avg_local_vec_size,PetscMax(ceil(vec_size*1.0/max_chunks),min_size)))));
331:   /* following line uses sizeof(PetscReal) instead of sizeof(PetscScalar) because the last dimension of chunkDims[] captures the 2* when complex numbers are being used */
332:   chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(PetscReal);

334:   /*
335:    if size/rank > max_chunk_size, we need radical measures: even going down to
336:    avg_local_vec_size is not enough, so we simply use chunk size of 4 GiB no matter
337:    what, composed in the most efficient way possible.
338:    N.B. this minimises the number of chunks, which may or may not be the optimal
339:    solution. In a BG, for example, the optimal solution is probably to make # chunks = #
340:    IO nodes involved, but this author has no access to a BG to figure out how to
341:    reliably find the right number. And even then it may or may not be enough.
342:    */
343:   if (avg_local_vec_size > max_chunk_size) {
344:     /* check if we can just split local z-axis: is that enough? */
345:     zslices = (PetscInt)ceil(vec_size*1.0/(da->p*max_chunk_size))*zslices;
346:     if (zslices > da->P) {
347:       /* lattice is too large in xy-directions, splitting z only is not enough */
348:       zslices = da->P;
349:       yslices= (PetscInt)ceil(vec_size*1.0/(zslices*da->n*max_chunk_size))*yslices;
350:       if (yslices > da->N) {
351:         /* lattice is too large in x-direction, splitting along z, y is not enough */
352:         yslices = da->N;
353:         xslices= (PetscInt)ceil(vec_size*1.0/(zslices*yslices*da->m*max_chunk_size))*xslices;
354:       }
355:     }
356:     dim = 0;
357:     if (timestep >= 0) {
358:       ++dim;
359:     }
360:     /* prefer to split z-axis, even down to planar slices */
361:     if (dimension == 3) {
362:       chunkDims[dim++] = (hsize_t) da->P/zslices;
363:       chunkDims[dim++] = (hsize_t) da->N/yslices;
364:       chunkDims[dim++] = (hsize_t) da->M/xslices;
365:     } else {
366:       /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */
367:       chunkDims[dim++] = (hsize_t) da->N/yslices;
368:       chunkDims[dim++] = (hsize_t) da->M/xslices;
369:     }
370:     chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(double);
371:   } else {
372:     if (target_size < chunk_size) {
373:       /* only change the defaults if target_size < chunk_size */
374:       dim = 0;
375:       if (timestep >= 0) {
376:         ++dim;
377:       }
378:       /* prefer to split z-axis, even down to planar slices */
379:       if (dimension == 3) {
380:         /* try splitting the z-axis to core-size bits, i.e. divide chunk size by # comm_size in z-direction */
381:         if (target_size >= chunk_size/da->p) {
382:           /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */
383:           chunkDims[dim] = (hsize_t) ceil(da->P*1.0/da->p);
384:         } else {
385:           /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be
386:            radical and let everyone write all they've got */
387:           chunkDims[dim++] = (hsize_t) ceil(da->P*1.0/da->p);
388:           chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n);
389:           chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m);
390:         }
391:       } else {
392:         /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */
393:         if (target_size >= chunk_size/da->n) {
394:           /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */
395:           chunkDims[dim] = (hsize_t) ceil(da->N*1.0/da->n);
396:         } else {
397:           /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be
398:            radical and let everyone write all they've got */
399:           chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n);
400:           chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m);
401:         }

403:       }
404:       chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(double);
405:     } else {
406:       /* precomputed chunks are fine, we don't need to do anything */
407:     }
408:   }
409:   return(0);
410: }
411: #endif

413: #if defined(PETSC_HAVE_HDF5)
414: PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer)
415: {
416:   PetscViewer_HDF5  *hdf5 = (PetscViewer_HDF5*) viewer->data;
417:   DM                dm;
418:   DM_DA             *da;
419:   hid_t             filespace;  /* file dataspace identifier */
420:   hid_t             chunkspace; /* chunk dataset property identifier */
421:   hid_t             dset_id;    /* dataset identifier */
422:   hid_t             memspace;   /* memory dataspace identifier */
423:   hid_t             file_id;
424:   hid_t             group;
425:   hid_t             memscalartype; /* scalar type for mem (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
426:   hid_t             filescalartype; /* scalar type for file (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
427:   hsize_t           dim;
428:   hsize_t           maxDims[6]={0}, dims[6]={0}, chunkDims[6]={0}, count[6]={0}, offset[6]={0}; /* we depend on these being sane later on  */
429:   PetscInt          timestep, dimension;
430:   const PetscScalar *x;
431:   const char        *vecname;
432:   PetscErrorCode    ierr;
433:   PetscBool         dim2;
434:   PetscBool         spoutput;

437:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
438:   PetscViewerHDF5GetTimestep(viewer, &timestep);
439:   PetscViewerHDF5GetBaseDimension2(viewer,&dim2);
440:   PetscViewerHDF5GetSPOutput(viewer,&spoutput);

442:   VecGetDM(xin,&dm);
443:   if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
444:   da = (DM_DA*)dm->data;
445:   DMGetDimension(dm, &dimension);

447:   /* Create the dataspace for the dataset.
448:    *
449:    * dims - holds the current dimensions of the dataset
450:    *
451:    * maxDims - holds the maximum dimensions of the dataset (unlimited
452:    * for the number of time steps with the current dimensions for the
453:    * other dimensions; so only additional time steps can be added).
454:    *
455:    * chunkDims - holds the size of a single time step (required to
456:    * permit extending dataset).
457:    */
458:   dim = 0;
459:   if (timestep >= 0) {
460:     dims[dim]      = timestep+1;
461:     maxDims[dim]   = H5S_UNLIMITED;
462:     chunkDims[dim] = 1;
463:     ++dim;
464:   }
465:   if (dimension == 3) {
466:     PetscHDF5IntCast(da->P,dims+dim);
467:     maxDims[dim]   = dims[dim];
468:     chunkDims[dim] = dims[dim];
469:     ++dim;
470:   }
471:   if (dimension > 1) {
472:     PetscHDF5IntCast(da->N,dims+dim);
473:     maxDims[dim]   = dims[dim];
474:     chunkDims[dim] = dims[dim];
475:     ++dim;
476:   }
477:   PetscHDF5IntCast(da->M,dims+dim);
478:   maxDims[dim]   = dims[dim];
479:   chunkDims[dim] = dims[dim];
480:   ++dim;
481:   if (da->w > 1 || dim2) {
482:     PetscHDF5IntCast(da->w,dims+dim);
483:     maxDims[dim]   = dims[dim];
484:     chunkDims[dim] = dims[dim];
485:     ++dim;
486:   }
487: #if defined(PETSC_USE_COMPLEX)
488:   dims[dim]      = 2;
489:   maxDims[dim]   = dims[dim];
490:   chunkDims[dim] = dims[dim];
491:   ++dim;
492: #endif

494:   VecGetHDF5ChunkSize(da, xin, dimension, timestep, chunkDims);

496:   PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));

498: #if defined(PETSC_USE_REAL_SINGLE)
499:   memscalartype = H5T_NATIVE_FLOAT;
500:   filescalartype = H5T_NATIVE_FLOAT;
501: #elif defined(PETSC_USE_REAL___FLOAT128)
502: #error "HDF5 output with 128 bit floats not supported."
503: #elif defined(PETSC_USE_REAL___FP16)
504: #error "HDF5 output with 16 bit floats not supported."
505: #else
506:   memscalartype = H5T_NATIVE_DOUBLE;
507:   if (spoutput == PETSC_TRUE) filescalartype = H5T_NATIVE_FLOAT;
508:   else filescalartype = H5T_NATIVE_DOUBLE;
509: #endif

511:   /* Create the dataset with default properties and close filespace */
512:   PetscObjectGetName((PetscObject)xin,&vecname);
513:   if (!H5Lexists(group, vecname, H5P_DEFAULT)) {
514:     /* Create chunk */
515:     PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE));
516:     PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims));

518:     PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
519:   } else {
520:     PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT));
521:     PetscStackCallHDF5(H5Dset_extent,(dset_id, dims));
522:   }
523:   PetscStackCallHDF5(H5Sclose,(filespace));

525:   /* Each process defines a dataset and writes it to the hyperslab in the file */
526:   dim = 0;
527:   if (timestep >= 0) {
528:     offset[dim] = timestep;
529:     ++dim;
530:   }
531:   if (dimension == 3) {PetscHDF5IntCast(da->zs,offset + dim++);}
532:   if (dimension > 1)  {PetscHDF5IntCast(da->ys,offset + dim++);}
533:   PetscHDF5IntCast(da->xs/da->w,offset + dim++);
534:   if (da->w > 1 || dim2) offset[dim++] = 0;
535: #if defined(PETSC_USE_COMPLEX)
536:   offset[dim++] = 0;
537: #endif
538:   dim = 0;
539:   if (timestep >= 0) {
540:     count[dim] = 1;
541:     ++dim;
542:   }
543:   if (dimension == 3) {PetscHDF5IntCast(da->ze - da->zs,count + dim++);}
544:   if (dimension > 1)  {PetscHDF5IntCast(da->ye - da->ys,count + dim++);}
545:   PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);
546:   if (da->w > 1 || dim2) {PetscHDF5IntCast(da->w,count + dim++);}
547: #if defined(PETSC_USE_COMPLEX)
548:   count[dim++] = 2;
549: #endif
550:   PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
551:   PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
552:   PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));

554:   VecGetArrayRead(xin, &x);
555:   PetscStackCallHDF5(H5Dwrite,(dset_id, memscalartype, memspace, filespace, hdf5->dxpl_id, x));
556:   PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL));
557:   VecRestoreArrayRead(xin, &x);

559:   #if defined(PETSC_USE_COMPLEX)
560:   {
561:     PetscBool tru = PETSC_TRUE;
562:     PetscViewerHDF5WriteObjectAttribute(viewer,(PetscObject)xin,"complex",PETSC_BOOL,&tru);
563:   }
564:   #endif

566:   /* Close/release resources */
567:   if (group != file_id) {
568:     PetscStackCallHDF5(H5Gclose,(group));
569:   }
570:   PetscStackCallHDF5(H5Sclose,(filespace));
571:   PetscStackCallHDF5(H5Sclose,(memspace));
572:   PetscStackCallHDF5(H5Dclose,(dset_id));
573:   PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);
574:   return(0);
575: }
576: #endif

578: extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer);

580: #if defined(PETSC_HAVE_MPIIO)
581: static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write)
582: {
583:   PetscErrorCode    ierr;
584:   MPI_File          mfdes;
585:   PetscMPIInt       gsizes[4],lsizes[4],lstarts[4],asiz,dof;
586:   MPI_Datatype      view;
587:   const PetscScalar *array;
588:   MPI_Offset        off;
589:   MPI_Aint          ub,ul;
590:   PetscInt          type,rows,vecrows,tr[2];
591:   DM_DA             *dd = (DM_DA*)da->data;
592:   PetscBool         skipheader;

595:   VecGetSize(xin,&vecrows);
596:   PetscViewerBinaryGetSkipHeader(viewer,&skipheader);
597:   if (!write) {
598:     /* Read vector header. */
599:     if (!skipheader) {
600:       PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);
601:       type = tr[0];
602:       rows = tr[1];
603:       if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file");
604:       if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector");
605:     }
606:   } else {
607:     tr[0] = VEC_FILE_CLASSID;
608:     tr[1] = vecrows;
609:     if (!skipheader) {
610:       PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);
611:     }
612:   }

614:   PetscMPIIntCast(dd->w,&dof);
615:   gsizes[0]  = dof;
616:   PetscMPIIntCast(dd->M,gsizes+1);
617:   PetscMPIIntCast(dd->N,gsizes+2);
618:   PetscMPIIntCast(dd->P,gsizes+3);
619:   lsizes[0]  = dof;
620:   PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);
621:   PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);
622:   PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);
623:   lstarts[0] = 0;
624:   PetscMPIIntCast(dd->xs/dof,lstarts+1);
625:   PetscMPIIntCast(dd->ys,lstarts+2);
626:   PetscMPIIntCast(dd->zs,lstarts+3);
627:   MPI_Type_create_subarray(da->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);
628:   MPI_Type_commit(&view);

630:   PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);
631:   PetscViewerBinaryGetMPIIOOffset(viewer,&off);
632:   MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);
633:   VecGetArrayRead(xin,&array);
634:   asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof;
635:   if (write) {
636:     MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);
637:   } else {
638:     MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);
639:   }
640:   MPI_Type_get_extent(view,&ul,&ub);
641:   PetscViewerBinaryAddMPIIOOffset(viewer,ub);
642:   VecRestoreArrayRead(xin,&array);
643:   MPI_Type_free(&view);
644:   return(0);
645: }
646: #endif

648: PetscErrorCode  VecView_MPI_DA(Vec xin,PetscViewer viewer)
649: {
650:   DM                da;
651:   PetscErrorCode    ierr;
652:   PetscInt          dim;
653:   Vec               natural;
654:   PetscBool         isdraw,isvtk,isglvis;
655: #if defined(PETSC_HAVE_HDF5)
656:   PetscBool         ishdf5;
657: #endif
658:   const char        *prefix,*name;
659:   PetscViewerFormat format;

662:   VecGetDM(xin,&da);
663:   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
664:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
665:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);
666: #if defined(PETSC_HAVE_HDF5)
667:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
668: #endif
669:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERGLVIS,&isglvis);
670:   if (isdraw) {
671:     DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);
672:     if (dim == 1) {
673:       VecView_MPI_Draw_DA1d(xin,viewer);
674:     } else if (dim == 2) {
675:       VecView_MPI_Draw_DA2d(xin,viewer);
676:     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim);
677:   } else if (isvtk) {           /* Duplicate the Vec */
678:     Vec Y;
679:     VecDuplicate(xin,&Y);
680:     if (((PetscObject)xin)->name) {
681:       /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */
682:       PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);
683:     }
684:     VecCopy(xin,Y);
685:     {
686:       PetscObject dmvtk;
687:       PetscBool   compatible,compatibleSet;
688:       PetscViewerVTKGetDM(viewer,&dmvtk);
689:       if (dmvtk) {
691:         DMGetCompatibility(da,(DM)dmvtk,&compatible,&compatibleSet);
692:         if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_INCOMP,"Cannot confirm compatibility of DMs associated with Vecs viewed in the same VTK file. Check that grids are the same.");
693:       }
694:       PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_DEFAULT,PETSC_VTK_POINT_FIELD,PETSC_FALSE,(PetscObject)Y);
695:     }
696: #if defined(PETSC_HAVE_HDF5)
697:   } else if (ishdf5) {
698:     VecView_MPI_HDF5_DA(xin,viewer);
699: #endif
700:   } else if (isglvis) {
701:     VecView_GLVis(xin,viewer);
702:   } else {
703: #if defined(PETSC_HAVE_MPIIO)
704:     PetscBool isbinary,isMPIIO;

706:     PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
707:     if (isbinary) {
708:       PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);
709:       if (isMPIIO) {
710:         DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);
711:         return(0);
712:       }
713:     }
714: #endif

716:     /* call viewer on natural ordering */
717:     PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);
718:     DMDACreateNaturalVector(da,&natural);
719:     PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);
720:     DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);
721:     DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);
722:     PetscObjectGetName((PetscObject)xin,&name);
723:     PetscObjectSetName((PetscObject)natural,name);

725:     PetscViewerGetFormat(viewer,&format);
726:     if (format == PETSC_VIEWER_BINARY_MATLAB) {
727:       /* temporarily remove viewer format so it won't trigger in the VecView() */
728:       PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);
729:     }

731:     ((PetscObject)natural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE;
732:     VecView(natural,viewer);
733:     ((PetscObject)natural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE;

735:     if (format == PETSC_VIEWER_BINARY_MATLAB) {
736:       MPI_Comm    comm;
737:       FILE        *info;
738:       const char  *fieldname;
739:       char        fieldbuf[256];
740:       PetscInt    dim,ni,nj,nk,pi,pj,pk,dof,n;

742:       /* set the viewer format back into the viewer */
743:       PetscViewerPopFormat(viewer);
744:       PetscObjectGetComm((PetscObject)viewer,&comm);
745:       PetscViewerBinaryGetInfoPointer(viewer,&info);
746:       DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);
747:       PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");
748:       PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");
749:       if (dim == 1) { PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni); }
750:       if (dim == 2) { PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj); }
751:       if (dim == 3) { PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk); }

753:       for (n=0; n<dof; n++) {
754:         DMDAGetFieldName(da,n,&fieldname);
755:         if (!fieldname || !fieldname[0]) {
756:           PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);
757:           fieldname = fieldbuf;
758:         }
759:         if (dim == 1) { PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1); }
760:         if (dim == 2) { PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1); }
761:         if (dim == 3) { PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);}
762:       }
763:       PetscFPrintf(comm,info,"#$$ clear tmp; \n");
764:       PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");
765:     }

767:     VecDestroy(&natural);
768:   }
769:   return(0);
770: }

772: #if defined(PETSC_HAVE_HDF5)
773: PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer)
774: {
775:   PetscViewer_HDF5 *hdf5 = (PetscViewer_HDF5*) viewer->data;
776:   DM             da;
778:   int            dim,rdim;
779:   hsize_t        dims[6]={0},count[6]={0},offset[6]={0};
780:   PetscInt       dimension,timestep,dofInd;
781:   PetscScalar    *x;
782:   const char     *vecname;
783:   hid_t          filespace; /* file dataspace identifier */
784:   hid_t          dset_id;   /* dataset identifier */
785:   hid_t          memspace;  /* memory dataspace identifier */
786:   hid_t          file_id,group;
787:   hid_t          scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
788:   DM_DA          *dd;
789:   PetscBool      dim2 = PETSC_FALSE;

792: #if defined(PETSC_USE_REAL_SINGLE)
793:   scalartype = H5T_NATIVE_FLOAT;
794: #elif defined(PETSC_USE_REAL___FLOAT128)
795: #error "HDF5 output with 128 bit floats not supported."
796: #elif defined(PETSC_USE_REAL___FP16)
797: #error "HDF5 output with 16 bit floats not supported."
798: #else
799:   scalartype = H5T_NATIVE_DOUBLE;
800: #endif

802:   PetscViewerHDF5OpenGroup(viewer, &file_id, &group);
803:   PetscViewerHDF5GetTimestep(viewer, &timestep);
804:   PetscObjectGetName((PetscObject)xin,&vecname);
805:   VecGetDM(xin,&da);
806:   dd   = (DM_DA*)da->data;
807:   DMGetDimension(da, &dimension);

809:   /* Open dataset */
810:   PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT));

812:   /* Retrieve the dataspace for the dataset */
813:   PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
814:   PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL));

816:   /* Expected dimension for holding the dof's */
817: #if defined(PETSC_USE_COMPLEX)
818:   dofInd = rdim-2;
819: #else
820:   dofInd = rdim-1;
821: #endif

823:   /* The expected number of dimensions, assuming basedimension2 = false */
824:   dim = dimension;
825:   if (dd->w > 1) ++dim;
826:   if (timestep >= 0) ++dim;
827: #if defined(PETSC_USE_COMPLEX)
828:   ++dim;
829: #endif

831:   /* In this case the input dataset have one extra, unexpected dimension. */
832:   if (rdim == dim+1) {
833:     /* In this case the block size unity */
834:     if (dd->w == 1 && dims[dofInd] == 1) dim2 = PETSC_TRUE;

836:     /* Special error message for the case where dof does not match the input file */
837:     else if (dd->w != (PetscInt) dims[dofInd]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Number of dofs in file is %D, not %D as expected",(PetscInt)dims[dofInd],dd->w);

839:   /* Other cases where rdim != dim cannot be handled currently */
840:   } else if (rdim != dim) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file is %d, not %d as expected with dof = %D",rdim,dim,dd->w);

842:   /* Set up the hyperslab size */
843:   dim = 0;
844:   if (timestep >= 0) {
845:     offset[dim] = timestep;
846:     count[dim] = 1;
847:     ++dim;
848:   }
849:   if (dimension == 3) {
850:     PetscHDF5IntCast(dd->zs,offset + dim);
851:     PetscHDF5IntCast(dd->ze - dd->zs,count + dim);
852:     ++dim;
853:   }
854:   if (dimension > 1) {
855:     PetscHDF5IntCast(dd->ys,offset + dim);
856:     PetscHDF5IntCast(dd->ye - dd->ys,count + dim);
857:     ++dim;
858:   }
859:   PetscHDF5IntCast(dd->xs/dd->w,offset + dim);
860:   PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + dim);
861:   ++dim;
862:   if (dd->w > 1 || dim2) {
863:     offset[dim] = 0;
864:     PetscHDF5IntCast(dd->w,count + dim);
865:     ++dim;
866:   }
867: #if defined(PETSC_USE_COMPLEX)
868:   offset[dim] = 0;
869:   count[dim] = 2;
870:   ++dim;
871: #endif

873:   /* Create the memory and filespace */
874:   PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
875:   PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));

877:   VecGetArray(xin, &x);
878:   PetscStackCallHDF5(H5Dread,(dset_id, scalartype, memspace, filespace, hdf5->dxpl_id, x));
879:   VecRestoreArray(xin, &x);

881:   /* Close/release resources */
882:   if (group != file_id) {
883:     PetscStackCallHDF5(H5Gclose,(group));
884:   }
885:   PetscStackCallHDF5(H5Sclose,(filespace));
886:   PetscStackCallHDF5(H5Sclose,(memspace));
887:   PetscStackCallHDF5(H5Dclose,(dset_id));
888:   return(0);
889: }
890: #endif

892: PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer)
893: {
894:   DM             da;
896:   Vec            natural;
897:   const char     *prefix;
898:   PetscInt       bs;
899:   PetscBool      flag;
900:   DM_DA          *dd;
901: #if defined(PETSC_HAVE_MPIIO)
902:   PetscBool isMPIIO;
903: #endif

906:   VecGetDM(xin,&da);
907:   dd   = (DM_DA*)da->data;
908: #if defined(PETSC_HAVE_MPIIO)
909:   PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);
910:   if (isMPIIO) {
911:     DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);
912:     return(0);
913:   }
914: #endif

916:   PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);
917:   DMDACreateNaturalVector(da,&natural);
918:   PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);
919:   PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);
920:   VecLoad(natural,viewer);
921:   DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);
922:   DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);
923:   VecDestroy(&natural);
924:   PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");
925:   PetscOptionsGetInt(NULL,((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);
926:   if (flag && bs != dd->w) {
927:     PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);
928:   }
929:   return(0);
930: }

932: PetscErrorCode  VecLoad_Default_DA(Vec xin, PetscViewer viewer)
933: {
935:   DM             da;
936:   PetscBool      isbinary;
937: #if defined(PETSC_HAVE_HDF5)
938:   PetscBool ishdf5;
939: #endif

942:   VecGetDM(xin,&da);
943:   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");

945: #if defined(PETSC_HAVE_HDF5)
946:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
947: #endif
948:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);

950:   if (isbinary) {
951:     VecLoad_Binary_DA(xin,viewer);
952: #if defined(PETSC_HAVE_HDF5)
953:   } else if (ishdf5) {
954:     VecLoad_HDF5_DA(xin,viewer);
955: #endif
956:   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name);
957:   return(0);
958: }