xref: /petsc/src/dm/impls/da/gr2.c (revision 665c2ded495bb9782a7454dcfef3abf1536c3670)
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;
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   /*
227        get local ghosted version of coordinates
228   */
229   ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr);
230   if (!xcoorl) {
231     /* create DMDA to get local version of graphics */
232     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);
233     ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr);
234     ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr);
235     ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr);
236     ierr = PetscObjectDereference((PetscObject)dag);CHKERRQ(ierr);
237     ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr);
238   } else {
239     ierr = VecGetDM(xcoorl,&dag);CHKERRQ(ierr);
240   }
241   ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr);
242   ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr);
243   ierr = VecGetArray(xcoorl,&zctx.xy);CHKERRQ(ierr);
244 
245   /*
246         Get information about size of area each processor must do graphics for
247   */
248   ierr = DMDAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&bx,&by,0,0);CHKERRQ(ierr);
249   ierr = DMDAGetGhostCorners(dac,0,0,0,&zctx.m,&zctx.n,0);CHKERRQ(ierr);
250 
251   ierr = PetscOptionsGetBool(NULL,"-draw_contour_grid",&zctx.showgrid,NULL);CHKERRQ(ierr);
252 
253   ierr = DMDASelectFields(da,&ndisplayfields,&displayfields);CHKERRQ(ierr);
254 
255   ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
256   ierr = PetscOptionsGetBool(NULL,"-draw_ports",&useports,NULL);CHKERRQ(ierr);
257   if (useports || format == PETSC_VIEWER_DRAW_PORTS) {
258     ierr       = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr);
259     ierr       = PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);CHKERRQ(ierr);
260     zctx.name0 = 0;
261     zctx.name1 = 0;
262   } else {
263     ierr = DMDAGetCoordinateName(da,0,&zctx.name0);CHKERRQ(ierr);
264     ierr = DMDAGetCoordinateName(da,1,&zctx.name1);CHKERRQ(ierr);
265   }
266 
267   /*
268      Loop over each field; drawing each in a different window
269   */
270   for (i=0; i<ndisplayfields; i++) {
271     zctx.k = displayfields[i];
272     if (useports) {
273       ierr = PetscDrawViewPortsSet(ports,i);CHKERRQ(ierr);
274     } else {
275       ierr = PetscViewerDrawGetDraw(viewer,i,&draw);CHKERRQ(ierr);
276     }
277 
278     /*
279         Determine the min and max color in plot
280     */
281     ierr = VecStrideMin(xin,zctx.k,NULL,&zctx.min);CHKERRQ(ierr);
282     ierr = VecStrideMax(xin,zctx.k,NULL,&zctx.max);CHKERRQ(ierr);
283     if (zctx.k < nbounds) {
284       zctx.min = bounds[2*zctx.k];
285       zctx.max = bounds[2*zctx.k+1];
286     }
287     if (zctx.min == zctx.max) {
288       zctx.min -= 1.e-12;
289       zctx.max += 1.e-12;
290     }
291 
292     if (!rank) {
293       const char *title;
294 
295       ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr);
296       if (title) {
297         ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);
298       }
299     }
300     ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr);
301     ierr = PetscInfo2(da,"DMDA 2d contour plot min %G max %G\n",zctx.min,zctx.max);CHKERRQ(ierr);
302 
303     ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr);
304     if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);}
305 
306     zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min);
307 
308     ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr);
309   }
310   ierr = PetscFree(displayfields);CHKERRQ(ierr);
311   ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr);
312 
313   ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr);
314   ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr);
315   PetscFunctionReturn(0);
316 }
317 
318 
319 #if defined(PETSC_HAVE_HDF5)
320 #undef __FUNCT__
321 #define __FUNCT__ "VecView_MPI_HDF5_DA"
322 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer)
323 {
324   DM             dm;
325   DM_DA          *da;
326   hid_t          filespace;  /* file dataspace identifier */
327   hid_t          chunkspace; /* chunk dataset property identifier */
328   hid_t          plist_id;   /* property list identifier */
329   hid_t          dset_id;    /* dataset identifier */
330   hid_t          memspace;   /* memory dataspace identifier */
331   hid_t          file_id;
332   hid_t          group;
333   hid_t          scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
334   herr_t         status;
335   hsize_t        dim;
336   hsize_t        maxDims[6], dims[6], chunkDims[6], count[6], offset[6];
337   PetscInt       timestep;
338   PetscScalar    *x;
339   const char     *vecname;
340   PetscErrorCode ierr;
341 
342   PetscFunctionBegin;
343   ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr);
344   ierr = PetscViewerHDF5GetTimestep(viewer, &timestep);CHKERRQ(ierr);
345 
346   ierr = VecGetDM(xin,&dm);CHKERRQ(ierr);
347   if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
348   da = (DM_DA*)dm->data;
349 
350   /* Create the dataspace for the dataset.
351    *
352    * dims - holds the current dimensions of the dataset
353    *
354    * maxDims - holds the maximum dimensions of the dataset (unlimited
355    * for the number of time steps with the current dimensions for the
356    * other dimensions; so only additional time steps can be added).
357    *
358    * chunkDims - holds the size of a single time step (required to
359    * permit extending dataset).
360    */
361   dim = 0;
362   if (timestep >= 0) {
363     dims[dim]      = timestep+1;
364     maxDims[dim]   = H5S_UNLIMITED;
365     chunkDims[dim] = 1;
366     ++dim;
367   }
368   if (da->dim == 3) {
369     ierr           = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr);
370     maxDims[dim]   = dims[dim];
371     chunkDims[dim] = dims[dim];
372     ++dim;
373   }
374   if (da->dim > 1) {
375     ierr           = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr);
376     maxDims[dim]   = dims[dim];
377     chunkDims[dim] = dims[dim];
378     ++dim;
379   }
380   ierr           = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr);
381   maxDims[dim]   = dims[dim];
382   chunkDims[dim] = dims[dim];
383   ++dim;
384   if (da->w > 1) {
385     ierr           = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr);
386     maxDims[dim]   = dims[dim];
387     chunkDims[dim] = dims[dim];
388     ++dim;
389   }
390 #if defined(PETSC_USE_COMPLEX)
391   dims[dim]      = 2;
392   maxDims[dim]   = dims[dim];
393   chunkDims[dim] = dims[dim];
394   ++dim;
395 #endif
396   filespace = H5Screate_simple(dim, dims, maxDims);
397   if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()");
398 
399 #if defined(PETSC_USE_REAL_SINGLE)
400   scalartype = H5T_NATIVE_FLOAT;
401 #elif defined(PETSC_USE_REAL___FLOAT128)
402 #error "HDF5 output with 128 bit floats not supported."
403 #else
404   scalartype = H5T_NATIVE_DOUBLE;
405 #endif
406 
407   /* Create the dataset with default properties and close filespace */
408   ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr);
409   if (!H5Lexists(group, vecname, H5P_DEFAULT)) {
410     /* Create chunk */
411     chunkspace = H5Pcreate(H5P_DATASET_CREATE);
412     if (chunkspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()");
413     status = H5Pset_chunk(chunkspace, dim, chunkDims);CHKERRQ(status);
414 
415 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
416     dset_id = H5Dcreate2(group, vecname, scalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT);
417 #else
418     dset_id = H5Dcreate(group, vecname, scalartype, filespace, H5P_DEFAULT);
419 #endif
420     if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()");
421   } else {
422     dset_id = H5Dopen2(group, vecname, H5P_DEFAULT);
423     status  = H5Dset_extent(dset_id, dims);CHKERRQ(status);
424   }
425   status = H5Sclose(filespace);CHKERRQ(status);
426 
427   /* Each process defines a dataset and writes it to the hyperslab in the file */
428   dim = 0;
429   if (timestep >= 0) {
430     offset[dim] = timestep;
431     ++dim;
432   }
433   if (da->dim == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);}
434   if (da->dim > 1)  {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);}
435   ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr);
436   if (da->w > 1) offset[dim++] = 0;
437 #if defined(PETSC_USE_COMPLEX)
438   offset[dim++] = 0;
439 #endif
440   dim = 0;
441   if (timestep >= 0) {
442     count[dim] = 1;
443     ++dim;
444   }
445   if (da->dim == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);}
446   if (da->dim > 1)  {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);}
447   ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr);
448   if (da->w > 1) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);}
449 #if defined(PETSC_USE_COMPLEX)
450   count[dim++] = 2;
451 #endif
452   memspace = H5Screate_simple(dim, count, NULL);
453   if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()");
454 
455   filespace = H5Dget_space(dset_id);
456   if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()");
457   status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status);
458 
459   /* Create property list for collective dataset write */
460   plist_id = H5Pcreate(H5P_DATASET_XFER);
461   if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()");
462 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
463   status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status);
464 #endif
465   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
466 
467   ierr   = VecGetArray(xin, &x);CHKERRQ(ierr);
468   status = H5Dwrite(dset_id, scalartype, memspace, filespace, plist_id, x);CHKERRQ(status);
469   status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status);
470   ierr   = VecRestoreArray(xin, &x);CHKERRQ(ierr);
471 
472   /* Close/release resources */
473   if (group != file_id) {
474     status = H5Gclose(group);CHKERRQ(status);
475   }
476   status = H5Pclose(plist_id);CHKERRQ(status);
477   status = H5Sclose(filespace);CHKERRQ(status);
478   status = H5Sclose(memspace);CHKERRQ(status);
479   status = H5Dclose(dset_id);CHKERRQ(status);
480   ierr   = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr);
481   PetscFunctionReturn(0);
482 }
483 #endif
484 
485 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer);
486 
487 #if defined(PETSC_HAVE_MPIIO)
488 #undef __FUNCT__
489 #define __FUNCT__ "DMDAArrayMPIIO"
490 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write)
491 {
492   PetscErrorCode    ierr;
493   MPI_File          mfdes;
494   PetscMPIInt       gsizes[4],lsizes[4],lstarts[4],asiz,dof;
495   MPI_Datatype      view;
496   const PetscScalar *array;
497   MPI_Offset        off;
498   MPI_Aint          ub,ul;
499   PetscInt          type,rows,vecrows,tr[2];
500   DM_DA             *dd = (DM_DA*)da->data;
501 
502   PetscFunctionBegin;
503   ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr);
504   if (!write) {
505     /* Read vector header. */
506     ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr);
507     type = tr[0];
508     rows = tr[1];
509     if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file");
510     if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector");
511   } else {
512     tr[0] = VEC_FILE_CLASSID;
513     tr[1] = vecrows;
514     ierr  = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
515   }
516 
517   ierr       = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr);
518   gsizes[0]  = dof;
519   ierr       = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr);
520   ierr       = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr);
521   ierr       = PetscMPIIntCast(dd->P,gsizes+1);CHKERRQ(ierr);
522   lsizes[0]  = dof;
523   ierr       = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr);
524   ierr       = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr);
525   ierr       = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr);
526   lstarts[0] = 0;
527   ierr       = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr);
528   ierr       = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr);
529   ierr       = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr);
530   ierr       = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr);
531   ierr       = MPI_Type_commit(&view);CHKERRQ(ierr);
532 
533   ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr);
534   ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr);
535   ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr);
536   ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr);
537   asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof;
538   if (write) {
539     ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr);
540   } else {
541     ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr);
542   }
543   ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr);
544   ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr);
545   ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr);
546   ierr = MPI_Type_free(&view);CHKERRQ(ierr);
547   PetscFunctionReturn(0);
548 }
549 #endif
550 
551 EXTERN_C_BEGIN
552 #undef __FUNCT__
553 #define __FUNCT__ "VecView_MPI_DA"
554 PetscErrorCode  VecView_MPI_DA(Vec xin,PetscViewer viewer)
555 {
556   DM                da;
557   PetscErrorCode    ierr;
558   PetscInt          dim;
559   Vec               natural;
560   PetscBool         isdraw,isvtk;
561 #if defined(PETSC_HAVE_HDF5)
562   PetscBool         ishdf5;
563 #endif
564   const char        *prefix,*name;
565   PetscViewerFormat format;
566 
567   PetscFunctionBegin;
568   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
569   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
570   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
571   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr);
572 #if defined(PETSC_HAVE_HDF5)
573   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr);
574 #endif
575   if (isdraw) {
576     ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
577     if (dim == 1) {
578       ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr);
579     } else if (dim == 2) {
580       ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr);
581     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim);
582   } else if (isvtk) {           /* Duplicate the Vec and hold a reference to the DM */
583     Vec Y;
584     ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);
585     ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr);
586     if (((PetscObject)xin)->name) {
587       /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */
588       ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr);
589     }
590     ierr = VecCopy(xin,Y);CHKERRQ(ierr);
591     ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr);
592 #if defined(PETSC_HAVE_HDF5)
593   } else if (ishdf5) {
594     ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr);
595 #endif
596   } else {
597 #if defined(PETSC_HAVE_MPIIO)
598     PetscBool isbinary,isMPIIO;
599 
600     ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
601     if (isbinary) {
602       ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr);
603       if (isMPIIO) {
604         ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr);
605         PetscFunctionReturn(0);
606       }
607     }
608 #endif
609 
610     /* call viewer on natural ordering */
611     ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr);
612     ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
613     ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr);
614     ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);
615     ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);
616     ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr);
617     ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr);
618 
619     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
620     if (format == PETSC_VIEWER_BINARY_MATLAB) {
621       /* temporarily remove viewer format so it won't trigger in the VecView() */
622       ierr = PetscViewerSetFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr);
623     }
624 
625     ierr = VecView(natural,viewer);CHKERRQ(ierr);
626 
627     if (format == PETSC_VIEWER_BINARY_MATLAB) {
628       MPI_Comm    comm;
629       FILE        *info;
630       const char  *fieldname;
631       PetscInt    dim,ni,nj,nk,pi,pj,pk,dof,n;
632       PetscBool   flg;
633 
634       /* set the viewer format back into the viewer */
635       ierr = PetscViewerSetFormat(viewer,format);CHKERRQ(ierr);
636       ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
637       ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr);
638       ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr);
639       ierr = PetscFPrintf(comm,info,"%%--- begin code written by PetscViewerBinary for MATLAB format ---%\n");CHKERRQ(ierr);
640       ierr = PetscFPrintf(comm,info,"%%$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr);
641       if (dim == 1) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); }
642       if (dim == 2) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); }
643       if (dim == 3) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); }
644 
645       for (n=0; n<dof; n++) {
646         ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr);
647         ierr = PetscStrcmp(fieldname,"",&flg);CHKERRQ(ierr);
648         if (!flg) {
649           if (dim == 1) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
650           if (dim == 2) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
651           if (dim == 3) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);}
652         }
653       }
654       ierr = PetscFPrintf(comm,info,"%%$$ clear tmp; \n");CHKERRQ(ierr);
655       ierr = PetscFPrintf(comm,info,"%%--- end code written by PetscViewerBinary for MATLAB format ---%\n\n");CHKERRQ(ierr);
656     }
657 
658     ierr = VecDestroy(&natural);CHKERRQ(ierr);
659   }
660   PetscFunctionReturn(0);
661 }
662 EXTERN_C_END
663 
664 #if defined(PETSC_HAVE_HDF5)
665 #undef __FUNCT__
666 #define __FUNCT__ "VecLoad_HDF5_DA"
667 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer)
668 {
669   DM             da;
670   PetscErrorCode ierr;
671   hsize_t        dim;
672   hsize_t        count[5];
673   hsize_t        offset[5];
674   PetscInt       cnt = 0;
675   PetscScalar    *x;
676   const char     *vecname;
677   hid_t          filespace; /* file dataspace identifier */
678   hid_t          plist_id;  /* property list identifier */
679   hid_t          dset_id;   /* dataset identifier */
680   hid_t          memspace;  /* memory dataspace identifier */
681   hid_t          file_id;
682   herr_t         status;
683   DM_DA          *dd;
684 
685   PetscFunctionBegin;
686   ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr);
687   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
688   dd   = (DM_DA*)da->data;
689 
690   /* Create the dataspace for the dataset */
691   ierr = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1),&dim);CHKERRQ(ierr);
692 #if defined(PETSC_USE_COMPLEX)
693   dim++;
694 #endif
695 
696   /* Create the dataset with default properties and close filespace */
697   ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr);
698 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
699   dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT);
700 #else
701   dset_id = H5Dopen(file_id, vecname);
702 #endif
703   if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname);
704   filespace = H5Dget_space(dset_id);
705 
706   /* Each process defines a dataset and reads it from the hyperslab in the file */
707   cnt = 0;
708   if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->zs,offset + cnt++);CHKERRQ(ierr);}
709   if (dd->dim > 1)  {ierr = PetscHDF5IntCast(dd->ys,offset + cnt++);CHKERRQ(ierr);}
710   ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + cnt++);CHKERRQ(ierr);
711   if (dd->w > 1) offset[cnt++] = 0;
712 #if defined(PETSC_USE_COMPLEX)
713   offset[cnt++] = 0;
714 #endif
715   cnt = 0;
716   if (dd->dim == 3) {ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + cnt++);CHKERRQ(ierr);}
717   if (dd->dim > 1)  {ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + cnt++);CHKERRQ(ierr);}
718   ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + cnt++);CHKERRQ(ierr);
719   if (dd->w > 1) {ierr = PetscHDF5IntCast(dd->w,count + cnt++);CHKERRQ(ierr);}
720 #if defined(PETSC_USE_COMPLEX)
721   count[cnt++] = 2;
722 #endif
723   memspace = H5Screate_simple(dim, count, NULL);
724   if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()");
725 
726   status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status);
727 
728   /* Create property list for collective dataset write */
729   plist_id = H5Pcreate(H5P_DATASET_XFER);
730   if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()");
731 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
732   status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status);
733 #endif
734   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
735 
736   ierr   = VecGetArray(xin, &x);CHKERRQ(ierr);
737   status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status);
738   ierr   = VecRestoreArray(xin, &x);CHKERRQ(ierr);
739 
740   /* Close/release resources */
741   status = H5Pclose(plist_id);CHKERRQ(status);
742   status = H5Sclose(filespace);CHKERRQ(status);
743   status = H5Sclose(memspace);CHKERRQ(status);
744   status = H5Dclose(dset_id);CHKERRQ(status);
745   PetscFunctionReturn(0);
746 }
747 #endif
748 
749 #undef __FUNCT__
750 #define __FUNCT__ "VecLoad_Binary_DA"
751 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer)
752 {
753   DM             da;
754   PetscErrorCode ierr;
755   Vec            natural;
756   const char     *prefix;
757   PetscInt       bs;
758   PetscBool      flag;
759   DM_DA          *dd;
760 #if defined(PETSC_HAVE_MPIIO)
761   PetscBool isMPIIO;
762 #endif
763 
764   PetscFunctionBegin;
765   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
766   dd   = (DM_DA*)da->data;
767 #if defined(PETSC_HAVE_MPIIO)
768   ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr);
769   if (isMPIIO) {
770     ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr);
771     PetscFunctionReturn(0);
772   }
773 #endif
774 
775   ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr);
776   ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
777   ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr);
778   ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr);
779   ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr);
780   ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr);
781   ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr);
782   ierr = VecDestroy(&natural);CHKERRQ(ierr);
783   ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr);
784   ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr);
785   if (flag && bs != dd->w) {
786     ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr);
787   }
788   PetscFunctionReturn(0);
789 }
790 
791 EXTERN_C_BEGIN
792 #undef __FUNCT__
793 #define __FUNCT__ "VecLoad_Default_DA"
794 PetscErrorCode  VecLoad_Default_DA(Vec xin, PetscViewer viewer)
795 {
796   PetscErrorCode ierr;
797   DM             da;
798   PetscBool      isbinary;
799 #if defined(PETSC_HAVE_HDF5)
800   PetscBool ishdf5;
801 #endif
802 
803   PetscFunctionBegin;
804   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
805   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
806 
807 #if defined(PETSC_HAVE_HDF5)
808   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr);
809 #endif
810   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
811 
812   if (isbinary) {
813     ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr);
814 #if defined(PETSC_HAVE_HDF5)
815   } else if (ishdf5) {
816     ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr);
817 #endif
818   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name);
819   PetscFunctionReturn(0);
820 }
821 EXTERN_C_END
822