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