xref: /petsc/src/dm/impls/da/gr2.c (revision b90c6cbe8a28dbbfeb8633979a88a1769a41a59e)
1 
2 /*
3    Plots vectors obtained with DMDACreate2d()
4 */
5 
6 #include <petsc-private/dmdaimpl.h>      /*I  "petscdmda.h"   I*/
7 #include <petsc-private/vecimpl.h>
8 #include <petscdraw.h>
9 #include <petscviewerhdf5.h>
10 
11 /*
12         The data that is passed into the graphics callback
13 */
14 typedef struct {
15   PetscInt    m,n,step,k;
16   PetscReal   min,max,scale;
17   PetscScalar *xy,*v;
18   PetscBool   showgrid;
19   const char  *name0,*name1;
20 } ZoomCtx;
21 
22 /*
23        This does the drawing for one particular field
24     in one particular set of coordinates. It is a callback
25     called from PetscDrawZoom()
26 */
27 #undef __FUNCT__
28 #define __FUNCT__ "VecView_MPI_Draw_DA2d_Zoom"
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,step,id,c1,c2,c3,c4;
34   PetscReal      s,min,max,x1,x2,x3,x4,y_1,y2,y3,y4,xmin = PETSC_MAX_REAL,xmax = PETSC_MIN_REAL,ymin = PETSC_MAX_REAL,ymax = PETSC_MIN_REAL;
35   PetscReal      xminf,xmaxf,yminf,ymaxf,w;
36   PetscScalar    *v,*xy;
37   char           value[16];
38   size_t         len;
39 
40   PetscFunctionBegin;
41   m    = zctx->m;
42   n    = zctx->n;
43   step = zctx->step;
44   k    = zctx->k;
45   v    = zctx->v;
46   xy   = zctx->xy;
47   s    = zctx->scale;
48   min  = zctx->min;
49   max  = zctx->max;
50 
51   /* PetscDraw the contour plot patch */
52   for (j=0; j<n-1; j++) {
53     for (i=0; i<m-1; i++) {
54       id   = i+j*m;
55       x1   = PetscRealPart(xy[2*id]);
56       y_1  = PetscRealPart(xy[2*id+1]);
57       c1   = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscClipInterval(PetscRealPart(v[k+step*id]),min,max)-min));
58       xmin = PetscMin(xmin,x1);
59       ymin = PetscMin(ymin,y_1);
60       xmax = PetscMax(xmax,x1);
61       ymax = PetscMax(ymax,y_1);
62 
63       id   = i+j*m+1;
64       x2   = PetscRealPart(xy[2*id]);
65       y2   = y_1;
66       c2   = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscClipInterval(PetscRealPart(v[k+step*id]),min,max)-min));
67       xmin = PetscMin(xmin,x2);
68       xmax = PetscMax(xmax,x2);
69 
70       id   = i+j*m+1+m;
71       x3   = x2;
72       y3   = PetscRealPart(xy[2*id+1]);
73       c3   = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscClipInterval(PetscRealPart(v[k+step*id]),min,max)-min));
74       ymin = PetscMin(ymin,y3);
75       ymax = PetscMax(ymax,y3);
76 
77       id = i+j*m+m;
78       x4 = x1;
79       y4 = y3;
80       c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscClipInterval(PetscRealPart(v[k+step*id]),min,max)-min));
81 
82       ierr = PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);CHKERRQ(ierr);
83       ierr = PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);CHKERRQ(ierr);
84       if (zctx->showgrid) {
85         ierr = PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);CHKERRQ(ierr);
86         ierr = PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);CHKERRQ(ierr);
87         ierr = PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);CHKERRQ(ierr);
88         ierr = PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);CHKERRQ(ierr);
89       }
90     }
91   }
92   if (zctx->name0) {
93     PetscReal xl,yl,xr,yr,x,y;
94     ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr);
95     x    = xl + .3*(xr - xl);
96     xl   = xl + .01*(xr - xl);
97     y    = yr - .3*(yr - yl);
98     yl   = yl + .01*(yr - yl);
99     ierr = PetscDrawString(draw,x,yl,PETSC_DRAW_BLACK,zctx->name0);CHKERRQ(ierr);
100     ierr = PetscDrawStringVertical(draw,xl,y,PETSC_DRAW_BLACK,zctx->name1);CHKERRQ(ierr);
101   }
102   /*
103      Ideally we would use the PetscDrawAxis object to manage displaying the coordinate limits
104      but that may require some refactoring.
105   */
106   ierr = MPI_Allreduce(&xmin,&xminf,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)draw));CHKERRQ(ierr);
107   ierr = MPI_Allreduce(&xmax,&xmaxf,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)draw));CHKERRQ(ierr);
108   ierr = MPI_Allreduce(&ymin,&yminf,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)draw));CHKERRQ(ierr);
109   ierr = MPI_Allreduce(&ymax,&ymaxf,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)draw));CHKERRQ(ierr);
110   ierr = PetscSNPrintf(value,16,"%f",xminf);CHKERRQ(ierr);
111   ierr = PetscDrawString(draw,xminf,yminf - .05*(ymaxf - yminf),PETSC_DRAW_BLACK,value);CHKERRQ(ierr);
112   ierr = PetscSNPrintf(value,16,"%f",xmaxf);CHKERRQ(ierr);
113   ierr = PetscStrlen(value,&len);CHKERRQ(ierr);
114   ierr = PetscDrawStringGetSize(draw,&w,NULL);CHKERRQ(ierr);
115   ierr = PetscDrawString(draw,xmaxf - len*w,yminf - .05*(ymaxf - yminf),PETSC_DRAW_BLACK,value);CHKERRQ(ierr);
116   ierr = PetscSNPrintf(value,16,"%f",yminf);CHKERRQ(ierr);
117   ierr = PetscDrawString(draw,xminf - .05*(xmaxf - xminf),yminf,PETSC_DRAW_BLACK,value);CHKERRQ(ierr);
118   ierr = PetscSNPrintf(value,16,"%f",ymaxf);CHKERRQ(ierr);
119   ierr = PetscDrawString(draw,xminf - .05*(xmaxf - xminf),ymaxf,PETSC_DRAW_BLACK,value);CHKERRQ(ierr);
120   PetscFunctionReturn(0);
121 }
122 
123 #undef __FUNCT__
124 #define __FUNCT__ "VecView_MPI_Draw_DA2d"
125 PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer)
126 {
127   DM                 da,dac,dag;
128   PetscErrorCode     ierr;
129   PetscMPIInt        rank;
130   PetscInt           N,s,M,w,ncoors = 4;
131   const PetscInt     *lx,*ly;
132   PetscReal          coors[4],ymin,ymax,xmin,xmax;
133   PetscDraw          draw,popup;
134   PetscBool          isnull,useports = PETSC_FALSE;
135   MPI_Comm           comm;
136   Vec                xlocal,xcoor,xcoorl;
137   DMDABoundaryType   bx,by;
138   DMDAStencilType    st;
139   ZoomCtx            zctx;
140   PetscDrawViewPorts *ports = NULL;
141   PetscViewerFormat  format;
142   PetscInt           *displayfields;
143   PetscInt           ndisplayfields,i,nbounds;
144   const PetscReal    *bounds;
145 
146   PetscFunctionBegin;
147   zctx.showgrid = PETSC_FALSE;
148 
149   ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
150   ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0);
151   ierr = PetscViewerDrawGetBounds(viewer,&nbounds,&bounds);CHKERRQ(ierr);
152 
153   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
154   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
155 
156   ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr);
157   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
158 
159   ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);CHKERRQ(ierr);
160   ierr = DMDAGetOwnershipRanges(da,&lx,&ly,NULL);CHKERRQ(ierr);
161 
162   /*
163         Obtain a sequential vector that is going to contain the local values plus ONE layer of
164      ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will
165      update the local values pluse ONE layer of ghost values.
166   */
167   ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr);
168   if (!xlocal) {
169     if (bx !=  DMDA_BOUNDARY_NONE || by !=  DMDA_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) {
170       /*
171          if original da is not of stencil width one, or periodic or not a box stencil then
172          create a special DMDA to handle one level of ghost points for graphics
173       */
174       ierr = DMDACreate2d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);CHKERRQ(ierr);
175       ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr);
176     } else {
177       /* otherwise we can use the da we already have */
178       dac = da;
179     }
180     /* create local vector for holding ghosted values used in graphics */
181     ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr);
182     if (dac != da) {
183       /* don't keep any public reference of this DMDA, is is only available through xlocal */
184       ierr = PetscObjectDereference((PetscObject)dac);CHKERRQ(ierr);
185     } else {
186       /* remove association between xlocal and da, because below we compose in the opposite
187          direction and if we left this connect we'd get a loop, so the objects could
188          never be destroyed */
189       ierr = PetscObjectRemoveReference((PetscObject)xlocal,"__PETSc_dm");CHKERRQ(ierr);
190     }
191     ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr);
192     ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr);
193   } else {
194     if (bx !=  DMDA_BOUNDARY_NONE || by !=  DMDA_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) {
195       ierr = VecGetDM(xlocal, &dac);CHKERRQ(ierr);
196     } else {
197       dac = da;
198     }
199   }
200 
201   /*
202       Get local (ghosted) values of vector
203   */
204   ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr);
205   ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr);
206   ierr = VecGetArray(xlocal,&zctx.v);CHKERRQ(ierr);
207 
208   /* get coordinates of nodes */
209   ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr);
210   if (!xcoor) {
211     ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr);
212     ierr = DMGetCoordinates(da,&xcoor);CHKERRQ(ierr);
213   }
214 
215   /*
216       Determine the min and max  coordinates in plot
217   */
218   ierr     = VecStrideMin(xcoor,0,NULL,&xmin);CHKERRQ(ierr);
219   ierr     = VecStrideMax(xcoor,0,NULL,&xmax);CHKERRQ(ierr);
220   ierr     = VecStrideMin(xcoor,1,NULL,&ymin);CHKERRQ(ierr);
221   ierr     = VecStrideMax(xcoor,1,NULL,&ymax);CHKERRQ(ierr);
222   coors[0] = xmin - .05*(xmax- xmin); coors[2] = xmax + .05*(xmax - xmin);
223   coors[1] = ymin - .05*(ymax- ymin); coors[3] = ymax + .05*(ymax - ymin);
224   ierr     = PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %G %G %G %G\n",coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr);
225 
226   ierr = PetscOptionsGetRealArray(NULL,"-draw_coordinates",coors,&ncoors,NULL);CHKERRQ(ierr);
227 
228   /*
229        get local ghosted version of coordinates
230   */
231   ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr);
232   if (!xcoorl) {
233     /* create DMDA to get local version of graphics */
234     ierr = DMDACreate2d(comm,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);CHKERRQ(ierr);
235     ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr);
236     ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr);
237     ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr);
238     ierr = PetscObjectDereference((PetscObject)dag);CHKERRQ(ierr);
239     ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr);
240   } else {
241     ierr = VecGetDM(xcoorl,&dag);CHKERRQ(ierr);
242   }
243   ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr);
244   ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr);
245   ierr = VecGetArray(xcoorl,&zctx.xy);CHKERRQ(ierr);
246 
247   /*
248         Get information about size of area each processor must do graphics for
249   */
250   ierr = DMDAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&bx,&by,0,0);CHKERRQ(ierr);
251   ierr = DMDAGetGhostCorners(dac,0,0,0,&zctx.m,&zctx.n,0);CHKERRQ(ierr);
252 
253   ierr = PetscOptionsGetBool(NULL,"-draw_contour_grid",&zctx.showgrid,NULL);CHKERRQ(ierr);
254 
255   ierr = DMDASelectFields(da,&ndisplayfields,&displayfields);CHKERRQ(ierr);
256 
257   ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
258   ierr = PetscOptionsGetBool(NULL,"-draw_ports",&useports,NULL);CHKERRQ(ierr);
259   if (useports || format == PETSC_VIEWER_DRAW_PORTS) {
260     ierr       = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr);
261     ierr       = PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);CHKERRQ(ierr);
262     zctx.name0 = 0;
263     zctx.name1 = 0;
264   } else {
265     ierr = DMDAGetCoordinateName(da,0,&zctx.name0);CHKERRQ(ierr);
266     ierr = DMDAGetCoordinateName(da,1,&zctx.name1);CHKERRQ(ierr);
267   }
268 
269   /*
270      Loop over each field; drawing each in a different window
271   */
272   for (i=0; i<ndisplayfields; i++) {
273     zctx.k = displayfields[i];
274     if (useports) {
275       ierr = PetscDrawViewPortsSet(ports,i);CHKERRQ(ierr);
276     } else {
277       ierr = PetscViewerDrawGetDraw(viewer,i,&draw);CHKERRQ(ierr);
278     }
279 
280     /*
281         Determine the min and max color in plot
282     */
283     ierr = VecStrideMin(xin,zctx.k,NULL,&zctx.min);CHKERRQ(ierr);
284     ierr = VecStrideMax(xin,zctx.k,NULL,&zctx.max);CHKERRQ(ierr);
285     if (zctx.k < nbounds) {
286       zctx.min = bounds[2*zctx.k];
287       zctx.max = bounds[2*zctx.k+1];
288     }
289     if (zctx.min == zctx.max) {
290       zctx.min -= 1.e-12;
291       zctx.max += 1.e-12;
292     }
293 
294     if (!rank) {
295       const char *title;
296 
297       ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr);
298       if (title) {
299         ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);
300       }
301     }
302     ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr);
303     ierr = PetscInfo2(da,"DMDA 2d contour plot min %G max %G\n",zctx.min,zctx.max);CHKERRQ(ierr);
304 
305     ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr);
306     if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);}
307 
308     zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min);
309 
310     ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr);
311   }
312   ierr = PetscFree(displayfields);CHKERRQ(ierr);
313   ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr);
314 
315   ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr);
316   ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr);
317   PetscFunctionReturn(0);
318 }
319 
320 
321 #if defined(PETSC_HAVE_HDF5)
322 #undef __FUNCT__
323 #define __FUNCT__ "VecView_MPI_HDF5_DA"
324 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer)
325 {
326   DM             dm;
327   DM_DA          *da;
328   hid_t          filespace;  /* file dataspace identifier */
329   hid_t          chunkspace; /* chunk dataset property identifier */
330   hid_t          plist_id;   /* property list identifier */
331   hid_t          dset_id;    /* dataset identifier */
332   hid_t          memspace;   /* memory dataspace identifier */
333   hid_t          file_id;
334   hid_t          group;
335   hid_t          scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
336   herr_t         status;
337   hsize_t        dim;
338   hsize_t        maxDims[6], dims[6], chunkDims[6], count[6], offset[6];
339   PetscInt       timestep;
340   PetscScalar    *x;
341   const char     *vecname;
342   PetscErrorCode ierr;
343 
344   PetscFunctionBegin;
345   ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr);
346   ierr = PetscViewerHDF5GetTimestep(viewer, &timestep);CHKERRQ(ierr);
347 
348   ierr = VecGetDM(xin,&dm);CHKERRQ(ierr);
349   if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
350   da = (DM_DA*)dm->data;
351 
352   /* Create the dataspace for the dataset.
353    *
354    * dims - holds the current dimensions of the dataset
355    *
356    * maxDims - holds the maximum dimensions of the dataset (unlimited
357    * for the number of time steps with the current dimensions for the
358    * other dimensions; so only additional time steps can be added).
359    *
360    * chunkDims - holds the size of a single time step (required to
361    * permit extending dataset).
362    */
363   dim = 0;
364   if (timestep >= 0) {
365     dims[dim]      = timestep+1;
366     maxDims[dim]   = H5S_UNLIMITED;
367     chunkDims[dim] = 1;
368     ++dim;
369   }
370   if (da->dim == 3) {
371     ierr           = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr);
372     maxDims[dim]   = dims[dim];
373     chunkDims[dim] = dims[dim];
374     ++dim;
375   }
376   if (da->dim > 1) {
377     ierr           = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr);
378     maxDims[dim]   = dims[dim];
379     chunkDims[dim] = dims[dim];
380     ++dim;
381   }
382   ierr           = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr);
383   maxDims[dim]   = dims[dim];
384   chunkDims[dim] = dims[dim];
385   ++dim;
386   if (da->w > 1) {
387     ierr           = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr);
388     maxDims[dim]   = dims[dim];
389     chunkDims[dim] = dims[dim];
390     ++dim;
391   }
392 #if defined(PETSC_USE_COMPLEX)
393   dims[dim]      = 2;
394   maxDims[dim]   = dims[dim];
395   chunkDims[dim] = dims[dim];
396   ++dim;
397 #endif
398   filespace = H5Screate_simple(dim, dims, maxDims);
399   if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()");
400 
401 #if defined(PETSC_USE_REAL_SINGLE)
402   scalartype = H5T_NATIVE_FLOAT;
403 #elif defined(PETSC_USE_REAL___FLOAT128)
404 #error "HDF5 output with 128 bit floats not supported."
405 #else
406   scalartype = H5T_NATIVE_DOUBLE;
407 #endif
408 
409   /* Create the dataset with default properties and close filespace */
410   ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr);
411   if (!H5Lexists(group, vecname, H5P_DEFAULT)) {
412     /* Create chunk */
413     chunkspace = H5Pcreate(H5P_DATASET_CREATE);
414     if (chunkspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()");
415     status = H5Pset_chunk(chunkspace, dim, chunkDims);CHKERRQ(status);
416 
417 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
418     dset_id = H5Dcreate2(group, vecname, scalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT);
419 #else
420     dset_id = H5Dcreate(group, vecname, scalartype, filespace, H5P_DEFAULT);
421 #endif
422     if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()");
423   } else {
424     dset_id = H5Dopen2(group, vecname, H5P_DEFAULT);
425     status  = H5Dset_extent(dset_id, dims);CHKERRQ(status);
426   }
427   status = H5Sclose(filespace);CHKERRQ(status);
428 
429   /* Each process defines a dataset and writes it to the hyperslab in the file */
430   dim = 0;
431   if (timestep >= 0) {
432     offset[dim] = timestep;
433     ++dim;
434   }
435   if (da->dim == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);}
436   if (da->dim > 1)  {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);}
437   ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr);
438   if (da->w > 1) offset[dim++] = 0;
439 #if defined(PETSC_USE_COMPLEX)
440   offset[dim++] = 0;
441 #endif
442   dim = 0;
443   if (timestep >= 0) {
444     count[dim] = 1;
445     ++dim;
446   }
447   if (da->dim == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);}
448   if (da->dim > 1)  {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);}
449   ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr);
450   if (da->w > 1) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);}
451 #if defined(PETSC_USE_COMPLEX)
452   count[dim++] = 2;
453 #endif
454   memspace = H5Screate_simple(dim, count, NULL);
455   if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()");
456 
457   filespace = H5Dget_space(dset_id);
458   if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()");
459   status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status);
460 
461   /* Create property list for collective dataset write */
462   plist_id = H5Pcreate(H5P_DATASET_XFER);
463   if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()");
464 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
465   status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status);
466 #endif
467   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
468 
469   ierr   = VecGetArray(xin, &x);CHKERRQ(ierr);
470   status = H5Dwrite(dset_id, scalartype, memspace, filespace, plist_id, x);CHKERRQ(status);
471   status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status);
472   ierr   = VecRestoreArray(xin, &x);CHKERRQ(ierr);
473 
474   /* Close/release resources */
475   if (group != file_id) {
476     status = H5Gclose(group);CHKERRQ(status);
477   }
478   status = H5Pclose(plist_id);CHKERRQ(status);
479   status = H5Sclose(filespace);CHKERRQ(status);
480   status = H5Sclose(memspace);CHKERRQ(status);
481   status = H5Dclose(dset_id);CHKERRQ(status);
482   ierr   = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr);
483   PetscFunctionReturn(0);
484 }
485 #endif
486 
487 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer);
488 
489 #if defined(PETSC_HAVE_MPIIO)
490 #undef __FUNCT__
491 #define __FUNCT__ "DMDAArrayMPIIO"
492 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write)
493 {
494   PetscErrorCode    ierr;
495   MPI_File          mfdes;
496   PetscMPIInt       gsizes[4],lsizes[4],lstarts[4],asiz,dof;
497   MPI_Datatype      view;
498   const PetscScalar *array;
499   MPI_Offset        off;
500   MPI_Aint          ub,ul;
501   PetscInt          type,rows,vecrows,tr[2];
502   DM_DA             *dd = (DM_DA*)da->data;
503 
504   PetscFunctionBegin;
505   ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr);
506   if (!write) {
507     /* Read vector header. */
508     ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr);
509     type = tr[0];
510     rows = tr[1];
511     if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file");
512     if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector");
513   } else {
514     tr[0] = VEC_FILE_CLASSID;
515     tr[1] = vecrows;
516     ierr  = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
517   }
518 
519   ierr       = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr);
520   gsizes[0]  = dof;
521   ierr       = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr);
522   ierr       = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr);
523   ierr       = PetscMPIIntCast(dd->P,gsizes+1);CHKERRQ(ierr);
524   lsizes[0]  = dof;
525   ierr       = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr);
526   ierr       = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr);
527   ierr       = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr);
528   lstarts[0] = 0;
529   ierr       = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr);
530   ierr       = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr);
531   ierr       = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr);
532   ierr       = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr);
533   ierr       = MPI_Type_commit(&view);CHKERRQ(ierr);
534 
535   ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr);
536   ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr);
537   ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr);
538   ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr);
539   asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof;
540   if (write) {
541     ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr);
542   } else {
543     ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr);
544   }
545   ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr);
546   ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr);
547   ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr);
548   ierr = MPI_Type_free(&view);CHKERRQ(ierr);
549   PetscFunctionReturn(0);
550 }
551 #endif
552 
553 #undef __FUNCT__
554 #define __FUNCT__ "VecView_MPI_DA"
555 PetscErrorCode  VecView_MPI_DA(Vec xin,PetscViewer viewer)
556 {
557   DM                da;
558   PetscErrorCode    ierr;
559   PetscInt          dim;
560   Vec               natural;
561   PetscBool         isdraw,isvtk;
562 #if defined(PETSC_HAVE_HDF5)
563   PetscBool         ishdf5;
564 #endif
565   const char        *prefix,*name;
566   PetscViewerFormat format;
567 
568   PetscFunctionBegin;
569   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
570   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
571   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
572   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr);
573 #if defined(PETSC_HAVE_HDF5)
574   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr);
575 #endif
576   if (isdraw) {
577     ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
578     if (dim == 1) {
579       ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr);
580     } else if (dim == 2) {
581       ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr);
582     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim);
583   } else if (isvtk) {           /* Duplicate the Vec and hold a reference to the DM */
584     Vec Y;
585     ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);
586     ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr);
587     if (((PetscObject)xin)->name) {
588       /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */
589       ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr);
590     }
591     ierr = VecCopy(xin,Y);CHKERRQ(ierr);
592     ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr);
593 #if defined(PETSC_HAVE_HDF5)
594   } else if (ishdf5) {
595     ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr);
596 #endif
597   } else {
598 #if defined(PETSC_HAVE_MPIIO)
599     PetscBool isbinary,isMPIIO;
600 
601     ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
602     if (isbinary) {
603       ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr);
604       if (isMPIIO) {
605         ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr);
606         PetscFunctionReturn(0);
607       }
608     }
609 #endif
610 
611     /* call viewer on natural ordering */
612     ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr);
613     ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
614     ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr);
615     ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);
616     ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);
617     ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr);
618     ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr);
619 
620     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
621     if (format == PETSC_VIEWER_BINARY_MATLAB) {
622       /* temporarily remove viewer format so it won't trigger in the VecView() */
623       ierr = PetscViewerSetFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr);
624     }
625 
626     ierr = VecView(natural,viewer);CHKERRQ(ierr);
627 
628     if (format == PETSC_VIEWER_BINARY_MATLAB) {
629       MPI_Comm    comm;
630       FILE        *info;
631       const char  *fieldname;
632       PetscInt    dim,ni,nj,nk,pi,pj,pk,dof,n;
633       PetscBool   flg;
634 
635       /* set the viewer format back into the viewer */
636       ierr = PetscViewerSetFormat(viewer,format);CHKERRQ(ierr);
637       ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
638       ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr);
639       ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr);
640       ierr = PetscFPrintf(comm,info,"%%--- begin code written by PetscViewerBinary for MATLAB format ---%\n");CHKERRQ(ierr);
641       ierr = PetscFPrintf(comm,info,"%%$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr);
642       if (dim == 1) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); }
643       if (dim == 2) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); }
644       if (dim == 3) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); }
645 
646       for (n=0; n<dof; n++) {
647         ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr);
648         ierr = PetscStrcmp(fieldname,"",&flg);CHKERRQ(ierr);
649         if (!flg) {
650           if (dim == 1) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
651           if (dim == 2) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
652           if (dim == 3) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);}
653         }
654       }
655       ierr = PetscFPrintf(comm,info,"%%$$ clear tmp; \n");CHKERRQ(ierr);
656       ierr = PetscFPrintf(comm,info,"%%--- end code written by PetscViewerBinary for MATLAB format ---%\n\n");CHKERRQ(ierr);
657     }
658 
659     ierr = VecDestroy(&natural);CHKERRQ(ierr);
660   }
661   PetscFunctionReturn(0);
662 }
663 
664 #if defined(PETSC_HAVE_HDF5)
665 #undef __FUNCT__
666 #define __FUNCT__ "VecLoad_HDF5_DA"
667 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer)
668 {
669   DM             da;
670   PetscErrorCode ierr;
671   hsize_t        dim;
672   hsize_t        count[5];
673   hsize_t        offset[5];
674   PetscInt       cnt = 0;
675   PetscScalar    *x;
676   const char     *vecname;
677   hid_t          filespace; /* file dataspace identifier */
678   hid_t          plist_id;  /* property list identifier */
679   hid_t          dset_id;   /* dataset identifier */
680   hid_t          memspace;  /* memory dataspace identifier */
681   hid_t          file_id;
682   herr_t         status;
683   DM_DA          *dd;
684 
685   PetscFunctionBegin;
686   ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr);
687   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
688   dd   = (DM_DA*)da->data;
689 
690   /* Create the dataspace for the dataset */
691   ierr = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1),&dim);CHKERRQ(ierr);
692 #if defined(PETSC_USE_COMPLEX)
693   dim++;
694 #endif
695 
696   /* Create the dataset with default properties and close filespace */
697   ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr);
698 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
699   dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT);
700 #else
701   dset_id = H5Dopen(file_id, vecname);
702 #endif
703   if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname);
704   filespace = H5Dget_space(dset_id);
705 
706   /* Each process defines a dataset and reads it from the hyperslab in the file */
707   cnt = 0;
708   if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->zs,offset + cnt++);CHKERRQ(ierr);}
709   if (dd->dim > 1)  {ierr = PetscHDF5IntCast(dd->ys,offset + cnt++);CHKERRQ(ierr);}
710   ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + cnt++);CHKERRQ(ierr);
711   if (dd->w > 1) offset[cnt++] = 0;
712 #if defined(PETSC_USE_COMPLEX)
713   offset[cnt++] = 0;
714 #endif
715   cnt = 0;
716   if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + cnt++);CHKERRQ(ierr);}
717   if (dd->dim > 1)  {ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + cnt++);CHKERRQ(ierr);}
718   ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + cnt++);CHKERRQ(ierr);
719   if (dd->w > 1) {ierr = PetscHDF5IntCast(dd->w,count + cnt++);CHKERRQ(ierr);}
720 #if defined(PETSC_USE_COMPLEX)
721   count[cnt++] = 2;
722 #endif
723   memspace = H5Screate_simple(dim, count, NULL);
724   if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()");
725 
726   status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status);
727 
728   /* Create property list for collective dataset write */
729   plist_id = H5Pcreate(H5P_DATASET_XFER);
730   if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()");
731 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
732   status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status);
733 #endif
734   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
735 
736   ierr   = VecGetArray(xin, &x);CHKERRQ(ierr);
737   status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status);
738   ierr   = VecRestoreArray(xin, &x);CHKERRQ(ierr);
739 
740   /* Close/release resources */
741   status = H5Pclose(plist_id);CHKERRQ(status);
742   status = H5Sclose(filespace);CHKERRQ(status);
743   status = H5Sclose(memspace);CHKERRQ(status);
744   status = H5Dclose(dset_id);CHKERRQ(status);
745   PetscFunctionReturn(0);
746 }
747 #endif
748 
749 #undef __FUNCT__
750 #define __FUNCT__ "VecLoad_Binary_DA"
751 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer)
752 {
753   DM             da;
754   PetscErrorCode ierr;
755   Vec            natural;
756   const char     *prefix;
757   PetscInt       bs;
758   PetscBool      flag;
759   DM_DA          *dd;
760 #if defined(PETSC_HAVE_MPIIO)
761   PetscBool isMPIIO;
762 #endif
763 
764   PetscFunctionBegin;
765   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
766   dd   = (DM_DA*)da->data;
767 #if defined(PETSC_HAVE_MPIIO)
768   ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr);
769   if (isMPIIO) {
770     ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr);
771     PetscFunctionReturn(0);
772   }
773 #endif
774 
775   ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr);
776   ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
777   ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr);
778   ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr);
779   ierr = VecLoad(natural,viewer);CHKERRQ(ierr);
780   ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr);
781   ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr);
782   ierr = VecDestroy(&natural);CHKERRQ(ierr);
783   ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr);
784   ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr);
785   if (flag && bs != dd->w) {
786     ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr);
787   }
788   PetscFunctionReturn(0);
789 }
790 
791 #undef __FUNCT__
792 #define __FUNCT__ "VecLoad_Default_DA"
793 PetscErrorCode  VecLoad_Default_DA(Vec xin, PetscViewer viewer)
794 {
795   PetscErrorCode ierr;
796   DM             da;
797   PetscBool      isbinary;
798 #if defined(PETSC_HAVE_HDF5)
799   PetscBool ishdf5;
800 #endif
801 
802   PetscFunctionBegin;
803   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
804   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
805 
806 #if defined(PETSC_HAVE_HDF5)
807   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr);
808 #endif
809   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
810 
811   if (isbinary) {
812     ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr);
813 #if defined(PETSC_HAVE_HDF5)
814   } else if (ishdf5) {
815     ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr);
816 #endif
817   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name);
818   PetscFunctionReturn(0);
819 }
820