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