xref: /petsc/src/dm/impls/da/gr2.c (revision 00de8ff0695ff394d09a2c60082aeaab5870b6e2)
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 EXTERN_C_BEGIN
554 #undef __FUNCT__
555 #define __FUNCT__ "VecView_MPI_DA"
556 PetscErrorCode  VecView_MPI_DA(Vec xin,PetscViewer viewer)
557 {
558   DM                da;
559   PetscErrorCode    ierr;
560   PetscInt          dim;
561   Vec               natural;
562   PetscBool         isdraw,isvtk;
563 #if defined(PETSC_HAVE_HDF5)
564   PetscBool         ishdf5;
565 #endif
566   const char        *prefix,*name;
567   PetscViewerFormat format;
568 
569   PetscFunctionBegin;
570   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
571   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
572   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
573   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr);
574 #if defined(PETSC_HAVE_HDF5)
575   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr);
576 #endif
577   if (isdraw) {
578     ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
579     if (dim == 1) {
580       ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr);
581     } else if (dim == 2) {
582       ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr);
583     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim);
584   } else if (isvtk) {           /* Duplicate the Vec and hold a reference to the DM */
585     Vec Y;
586     ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);
587     ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr);
588     if (((PetscObject)xin)->name) {
589       /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */
590       ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr);
591     }
592     ierr = VecCopy(xin,Y);CHKERRQ(ierr);
593     ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr);
594 #if defined(PETSC_HAVE_HDF5)
595   } else if (ishdf5) {
596     ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr);
597 #endif
598   } else {
599 #if defined(PETSC_HAVE_MPIIO)
600     PetscBool isbinary,isMPIIO;
601 
602     ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
603     if (isbinary) {
604       ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr);
605       if (isMPIIO) {
606         ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr);
607         PetscFunctionReturn(0);
608       }
609     }
610 #endif
611 
612     /* call viewer on natural ordering */
613     ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr);
614     ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
615     ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr);
616     ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);
617     ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);
618     ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr);
619     ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr);
620 
621     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
622     if (format == PETSC_VIEWER_BINARY_MATLAB) {
623       /* temporarily remove viewer format so it won't trigger in the VecView() */
624       ierr = PetscViewerSetFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr);
625     }
626 
627     ierr = VecView(natural,viewer);CHKERRQ(ierr);
628 
629     if (format == PETSC_VIEWER_BINARY_MATLAB) {
630       MPI_Comm    comm;
631       FILE        *info;
632       const char  *fieldname;
633       PetscInt    dim,ni,nj,nk,pi,pj,pk,dof,n;
634       PetscBool   flg;
635 
636       /* set the viewer format back into the viewer */
637       ierr = PetscViewerSetFormat(viewer,format);CHKERRQ(ierr);
638       ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
639       ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr);
640       ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr);
641       ierr = PetscFPrintf(comm,info,"%%--- begin code written by PetscViewerBinary for MATLAB format ---%\n");CHKERRQ(ierr);
642       ierr = PetscFPrintf(comm,info,"%%$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr);
643       if (dim == 1) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); }
644       if (dim == 2) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); }
645       if (dim == 3) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); }
646 
647       for (n=0; n<dof; n++) {
648         ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr);
649         ierr = PetscStrcmp(fieldname,"",&flg);CHKERRQ(ierr);
650         if (!flg) {
651           if (dim == 1) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
652           if (dim == 2) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
653           if (dim == 3) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);}
654         }
655       }
656       ierr = PetscFPrintf(comm,info,"%%$$ clear tmp; \n");CHKERRQ(ierr);
657       ierr = PetscFPrintf(comm,info,"%%--- end code written by PetscViewerBinary for MATLAB format ---%\n\n");CHKERRQ(ierr);
658     }
659 
660     ierr = VecDestroy(&natural);CHKERRQ(ierr);
661   }
662   PetscFunctionReturn(0);
663 }
664 EXTERN_C_END
665 
666 #if defined(PETSC_HAVE_HDF5)
667 #undef __FUNCT__
668 #define __FUNCT__ "VecLoad_HDF5_DA"
669 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer)
670 {
671   DM             da;
672   PetscErrorCode ierr;
673   hsize_t        dim;
674   hsize_t        count[5];
675   hsize_t        offset[5];
676   PetscInt       cnt = 0;
677   PetscScalar    *x;
678   const char     *vecname;
679   hid_t          filespace; /* file dataspace identifier */
680   hid_t          plist_id;  /* property list identifier */
681   hid_t          dset_id;   /* dataset identifier */
682   hid_t          memspace;  /* memory dataspace identifier */
683   hid_t          file_id;
684   herr_t         status;
685   DM_DA          *dd;
686 
687   PetscFunctionBegin;
688   ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr);
689   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
690   dd   = (DM_DA*)da->data;
691 
692   /* Create the dataspace for the dataset */
693   ierr = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1),&dim);CHKERRQ(ierr);
694 #if defined(PETSC_USE_COMPLEX)
695   dim++;
696 #endif
697 
698   /* Create the dataset with default properties and close filespace */
699   ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr);
700 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
701   dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT);
702 #else
703   dset_id = H5Dopen(file_id, vecname);
704 #endif
705   if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname);
706   filespace = H5Dget_space(dset_id);
707 
708   /* Each process defines a dataset and reads it from the hyperslab in the file */
709   cnt = 0;
710   if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->zs,offset + cnt++);CHKERRQ(ierr);}
711   if (dd->dim > 1)  {ierr = PetscHDF5IntCast(dd->ys,offset + cnt++);CHKERRQ(ierr);}
712   ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + cnt++);CHKERRQ(ierr);
713   if (dd->w > 1) offset[cnt++] = 0;
714 #if defined(PETSC_USE_COMPLEX)
715   offset[cnt++] = 0;
716 #endif
717   cnt = 0;
718   if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + cnt++);CHKERRQ(ierr);}
719   if (dd->dim > 1)  {ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + cnt++);CHKERRQ(ierr);}
720   ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + cnt++);CHKERRQ(ierr);
721   if (dd->w > 1) {ierr = PetscHDF5IntCast(dd->w,count + cnt++);CHKERRQ(ierr);}
722 #if defined(PETSC_USE_COMPLEX)
723   count[cnt++] = 2;
724 #endif
725   memspace = H5Screate_simple(dim, count, NULL);
726   if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()");
727 
728   status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status);
729 
730   /* Create property list for collective dataset write */
731   plist_id = H5Pcreate(H5P_DATASET_XFER);
732   if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()");
733 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
734   status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status);
735 #endif
736   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
737 
738   ierr   = VecGetArray(xin, &x);CHKERRQ(ierr);
739   status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status);
740   ierr   = VecRestoreArray(xin, &x);CHKERRQ(ierr);
741 
742   /* Close/release resources */
743   status = H5Pclose(plist_id);CHKERRQ(status);
744   status = H5Sclose(filespace);CHKERRQ(status);
745   status = H5Sclose(memspace);CHKERRQ(status);
746   status = H5Dclose(dset_id);CHKERRQ(status);
747   PetscFunctionReturn(0);
748 }
749 #endif
750 
751 #undef __FUNCT__
752 #define __FUNCT__ "VecLoad_Binary_DA"
753 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer)
754 {
755   DM             da;
756   PetscErrorCode ierr;
757   Vec            natural;
758   const char     *prefix;
759   PetscInt       bs;
760   PetscBool      flag;
761   DM_DA          *dd;
762 #if defined(PETSC_HAVE_MPIIO)
763   PetscBool isMPIIO;
764 #endif
765 
766   PetscFunctionBegin;
767   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
768   dd   = (DM_DA*)da->data;
769 #if defined(PETSC_HAVE_MPIIO)
770   ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr);
771   if (isMPIIO) {
772     ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr);
773     PetscFunctionReturn(0);
774   }
775 #endif
776 
777   ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr);
778   ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
779   ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr);
780   ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr);
781   ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr);
782   ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr);
783   ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr);
784   ierr = VecDestroy(&natural);CHKERRQ(ierr);
785   ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr);
786   ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr);
787   if (flag && bs != dd->w) {
788     ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr);
789   }
790   PetscFunctionReturn(0);
791 }
792 
793 EXTERN_C_BEGIN
794 #undef __FUNCT__
795 #define __FUNCT__ "VecLoad_Default_DA"
796 PetscErrorCode  VecLoad_Default_DA(Vec xin, PetscViewer viewer)
797 {
798   PetscErrorCode ierr;
799   DM             da;
800   PetscBool      isbinary;
801 #if defined(PETSC_HAVE_HDF5)
802   PetscBool ishdf5;
803 #endif
804 
805   PetscFunctionBegin;
806   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
807   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
808 
809 #if defined(PETSC_HAVE_HDF5)
810   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr);
811 #endif
812   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
813 
814   if (isbinary) {
815     ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr);
816 #if defined(PETSC_HAVE_HDF5)
817   } else if (ishdf5) {
818     ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr);
819 #endif
820   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name);
821   PetscFunctionReturn(0);
822 }
823 EXTERN_C_END
824