xref: /petsc/src/dm/impls/da/gr2.c (revision 1ebf93c6b7d760d592de6ebe6cdc0debaa3caf75) 
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   PetscMPIInt       rank;
16   PetscInt          m,n,dof,k;
17   PetscReal         xmin,xmax,ymin,ymax,min,max;
18   const PetscScalar *xy,*v;
19   PetscBool         showaxis,showgrid;
20   const char        *name0,*name1;
21 } ZoomCtx;
22 
23 /*
24        This does the drawing for one particular field
25     in one particular set of coordinates. It is a callback
26     called from PetscDrawZoom()
27 */
28 #undef __FUNCT__
29 #define __FUNCT__ "VecView_MPI_Draw_DA2d_Zoom"
30 PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx)
31 {
32   ZoomCtx           *zctx = (ZoomCtx*)ctx;
33   PetscErrorCode    ierr;
34   PetscInt          m,n,i,j,k,dof,id,c1,c2,c3,c4;
35   PetscReal         min,max,x1,x2,x3,x4,y_1,y2,y3,y4;
36   const PetscScalar *xy,*v;
37 
38   PetscFunctionBegin;
39   m    = zctx->m;
40   n    = zctx->n;
41   dof  = zctx->dof;
42   k    = zctx->k;
43   xy   = zctx->xy;
44   v    = zctx->v;
45   min  = zctx->min;
46   max  = zctx->max;
47 
48   /* PetscDraw the contour plot patch */
49   ierr = PetscDrawCollectiveBegin(draw);CHKERRQ(ierr);
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   = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max);
56 
57       id   = i+j*m+1;
58       x2   = PetscRealPart(xy[2*id]);
59       y2   = PetscRealPart(xy[2*id+1]);
60       c2   = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max);
61 
62       id   = i+j*m+1+m;
63       x3   = PetscRealPart(xy[2*id]);
64       y3   = PetscRealPart(xy[2*id+1]);
65       c3   = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max);
66 
67       id   = i+j*m+m;
68       x4   = PetscRealPart(xy[2*id]);
69       y4   = PetscRealPart(xy[2*id+1]);
70       c4   = PetscDrawRealToColor(PetscRealPart(v[k+dof*id]),min,max);
71 
72       ierr = PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);CHKERRQ(ierr);
73       ierr = PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);CHKERRQ(ierr);
74       if (zctx->showgrid) {
75         ierr = PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);CHKERRQ(ierr);
76         ierr = PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);CHKERRQ(ierr);
77         ierr = PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);CHKERRQ(ierr);
78         ierr = PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);CHKERRQ(ierr);
79       }
80     }
81   }
82   if (zctx->showaxis && !zctx->rank) {
83     if (zctx->name0 || zctx->name1) {
84       PetscReal xl,yl,xr,yr,x,y;
85       ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr);
86       x  = xl + .30*(xr - xl);
87       xl = xl + .01*(xr - xl);
88       y  = yr - .30*(yr - yl);
89       yl = yl + .01*(yr - yl);
90       if (zctx->name0) {ierr = PetscDrawString(draw,x,yl,PETSC_DRAW_BLACK,zctx->name0);CHKERRQ(ierr);}
91       if (zctx->name1) {ierr = PetscDrawStringVertical(draw,xl,y,PETSC_DRAW_BLACK,zctx->name1);CHKERRQ(ierr);}
92     }
93     /*
94        Ideally we would use the PetscDrawAxis object to manage displaying the coordinate limits
95        but that may require some refactoring.
96     */
97     {
98       double xmin = (double)zctx->xmin, ymin = (double)zctx->ymin;
99       double xmax = (double)zctx->xmax, ymax = (double)zctx->ymax;
100       char   value[16]; size_t len; PetscReal w;
101       ierr = PetscSNPrintf(value,16,"%0.2e",xmin);CHKERRQ(ierr);
102       ierr = PetscDrawString(draw,xmin,ymin - .05*(ymax - ymin),PETSC_DRAW_BLACK,value);CHKERRQ(ierr);
103       ierr = PetscSNPrintf(value,16,"%0.2e",xmax);CHKERRQ(ierr);
104       ierr = PetscStrlen(value,&len);CHKERRQ(ierr);
105       ierr = PetscDrawStringGetSize(draw,&w,NULL);CHKERRQ(ierr);
106       ierr = PetscDrawString(draw,xmax - len*w,ymin - .05*(ymax - ymin),PETSC_DRAW_BLACK,value);CHKERRQ(ierr);
107       ierr = PetscSNPrintf(value,16,"%0.2e",ymin);CHKERRQ(ierr);
108       ierr = PetscDrawString(draw,xmin - .05*(xmax - xmin),ymin,PETSC_DRAW_BLACK,value);CHKERRQ(ierr);
109       ierr = PetscSNPrintf(value,16,"%0.2e",ymax);CHKERRQ(ierr);
110       ierr = PetscDrawString(draw,xmin - .05*(xmax - xmin),ymax,PETSC_DRAW_BLACK,value);CHKERRQ(ierr);
111     }
112   }
113   ierr = PetscDrawCollectiveEnd(draw);CHKERRQ(ierr);
114   PetscFunctionReturn(0);
115 }
116 
117 #undef __FUNCT__
118 #define __FUNCT__ "VecView_MPI_Draw_DA2d"
119 PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer)
120 {
121   DM                 da,dac,dag;
122   PetscErrorCode     ierr;
123   PetscInt           N,s,M,w,ncoors = 4;
124   const PetscInt     *lx,*ly;
125   PetscReal          coors[4];
126   PetscDraw          draw,popup;
127   PetscBool          isnull,useports = PETSC_FALSE;
128   MPI_Comm           comm;
129   Vec                xlocal,xcoor,xcoorl;
130   DMBoundaryType     bx,by;
131   DMDAStencilType    st;
132   ZoomCtx            zctx;
133   PetscDrawViewPorts *ports = NULL;
134   PetscViewerFormat  format;
135   PetscInt           *displayfields;
136   PetscInt           ndisplayfields,i,nbounds;
137   const PetscReal    *bounds;
138 
139   PetscFunctionBegin;
140   zctx.showgrid = PETSC_FALSE;
141   zctx.showaxis = PETSC_TRUE;
142 
143   ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
144   ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr);
145   if (isnull) PetscFunctionReturn(0);
146 
147   ierr = PetscViewerDrawGetBounds(viewer,&nbounds,&bounds);CHKERRQ(ierr);
148 
149   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
150   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
151 
152   ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr);
153   ierr = MPI_Comm_rank(comm,&zctx.rank);CHKERRQ(ierr);
154 
155   ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&bx,&by,0,&st);CHKERRQ(ierr);
156   ierr = DMDAGetOwnershipRanges(da,&lx,&ly,NULL);CHKERRQ(ierr);
157 
158   /*
159      Obtain a sequential vector that is going to contain the local values plus ONE layer of
160      ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will
161      update the local values pluse ONE layer of ghost values.
162   */
163   ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr);
164   if (!xlocal) {
165     if (bx !=  DM_BOUNDARY_NONE || by !=  DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) {
166       /*
167          if original da is not of stencil width one, or periodic or not a box stencil then
168          create a special DMDA to handle one level of ghost points for graphics
169       */
170       ierr = DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);CHKERRQ(ierr);
171       ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr);
172     } else {
173       /* otherwise we can use the da we already have */
174       dac = da;
175     }
176     /* create local vector for holding ghosted values used in graphics */
177     ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr);
178     if (dac != da) {
179       /* don't keep any public reference of this DMDA, is is only available through xlocal */
180       ierr = PetscObjectDereference((PetscObject)dac);CHKERRQ(ierr);
181     } else {
182       /* remove association between xlocal and da, because below we compose in the opposite
183          direction and if we left this connect we'd get a loop, so the objects could
184          never be destroyed */
185       ierr = PetscObjectRemoveReference((PetscObject)xlocal,"__PETSc_dm");CHKERRQ(ierr);
186     }
187     ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr);
188     ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr);
189   } else {
190     if (bx !=  DM_BOUNDARY_NONE || by !=  DM_BOUNDARY_NONE || s != 1 || st != DMDA_STENCIL_BOX) {
191       ierr = VecGetDM(xlocal, &dac);CHKERRQ(ierr);
192     } else {
193       dac = da;
194     }
195   }
196 
197   /*
198       Get local (ghosted) values of vector
199   */
200   ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr);
201   ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr);
202   ierr = VecGetArrayRead(xlocal,&zctx.v);CHKERRQ(ierr);
203 
204   /*
205       Get coordinates of nodes
206   */
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,&zctx.xmin);CHKERRQ(ierr);
217   ierr = VecStrideMax(xcoor,0,NULL,&zctx.xmax);CHKERRQ(ierr);
218   ierr = VecStrideMin(xcoor,1,NULL,&zctx.ymin);CHKERRQ(ierr);
219   ierr = VecStrideMax(xcoor,1,NULL,&zctx.ymax);CHKERRQ(ierr);
220   ierr = PetscOptionsGetBool(NULL,NULL,"-draw_contour_axis",&zctx.showaxis,NULL);CHKERRQ(ierr);
221   if (zctx.showaxis) {
222     coors[0] = zctx.xmin - .05*(zctx.xmax - zctx.xmin); coors[1] = zctx.ymin - .05*(zctx.ymax - zctx.ymin);
223     coors[2] = zctx.xmax + .05*(zctx.xmax - zctx.xmin); coors[3] = zctx.ymax + .05*(zctx.ymax - zctx.ymin);
224   } else {
225     coors[0] = zctx.xmin; coors[1] = zctx.ymin; coors[2] = zctx.xmax; coors[3] = zctx.ymax;
226   }
227   ierr = PetscOptionsGetRealArray(NULL,NULL,"-draw_coordinates",coors,&ncoors,NULL);CHKERRQ(ierr);
228   ierr = PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %g %g %g %g\n",(double)coors[0],(double)coors[1],(double)coors[2],(double)coors[3]);CHKERRQ(ierr);
229 
230   /*
231       Get local ghosted version of coordinates
232   */
233   ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr);
234   if (!xcoorl) {
235     /* create DMDA to get local version of graphics */
236     ierr = DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);CHKERRQ(ierr);
237     ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr);
238     ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr);
239     ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr);
240     ierr = PetscObjectDereference((PetscObject)dag);CHKERRQ(ierr);
241     ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr);
242   } else {
243     ierr = VecGetDM(xcoorl,&dag);CHKERRQ(ierr);
244   }
245   ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr);
246   ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr);
247   ierr = VecGetArrayRead(xcoorl,&zctx.xy);CHKERRQ(ierr);
248   ierr = DMDAGetCoordinateName(da,0,&zctx.name0);CHKERRQ(ierr);
249   ierr = DMDAGetCoordinateName(da,1,&zctx.name1);CHKERRQ(ierr);
250 
251   /*
252       Get information about size of area each processor must do graphics for
253   */
254   ierr = DMDAGetInfo(dac,NULL,&M,&N,NULL,NULL,NULL,NULL,&zctx.dof,NULL,&bx,&by,NULL,NULL);CHKERRQ(ierr);
255   ierr = DMDAGetGhostCorners(dac,NULL,NULL,NULL,&zctx.m,&zctx.n,NULL);CHKERRQ(ierr);
256   ierr = PetscOptionsGetBool(NULL,NULL,"-draw_contour_grid",&zctx.showgrid,NULL);CHKERRQ(ierr);
257 
258   ierr = DMDASelectFields(da,&ndisplayfields,&displayfields);CHKERRQ(ierr);
259   ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
260   ierr = PetscOptionsGetBool(NULL,NULL,"-draw_ports",&useports,NULL);CHKERRQ(ierr);
261   if (format == PETSC_VIEWER_DRAW_PORTS) useports = PETSC_TRUE;
262   if (useports) {
263     ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
264     ierr = PetscDrawCheckResizedWindow(draw);CHKERRQ(ierr);
265     ierr = PetscDrawClear(draw);CHKERRQ(ierr);
266     ierr = PetscDrawViewPortsCreate(draw,ndisplayfields,&ports);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 
275     /* determine the min and max value in plot */
276     ierr = VecStrideMin(xin,zctx.k,NULL,&zctx.min);CHKERRQ(ierr);
277     ierr = VecStrideMax(xin,zctx.k,NULL,&zctx.max);CHKERRQ(ierr);
278     if (zctx.k < nbounds) {
279       zctx.min = bounds[2*zctx.k];
280       zctx.max = bounds[2*zctx.k+1];
281     }
282     if (zctx.min == zctx.max) {
283       zctx.min -= 1.e-12;
284       zctx.max += 1.e-12;
285     }
286     ierr = PetscInfo2(da,"DMDA 2d contour plot min %g max %g\n",(double)zctx.min,(double)zctx.max);CHKERRQ(ierr);
287 
288     if (useports) {
289       ierr = PetscDrawViewPortsSet(ports,i);CHKERRQ(ierr);
290     } else {
291       const char *title;
292       ierr = PetscViewerDrawGetDraw(viewer,i,&draw);CHKERRQ(ierr);
293       ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr);
294       if (title) {ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr);}
295     }
296 
297     ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr);
298     ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);
299     ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr);
300     ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr);
301     if (!useports) {ierr = PetscDrawSave(draw);CHKERRQ(ierr);}
302   }
303   if (useports) {
304     ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr);
305     ierr = PetscDrawSave(draw);CHKERRQ(ierr);
306   }
307 
308   ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr);
309   ierr = PetscFree(displayfields);CHKERRQ(ierr);
310   ierr = VecRestoreArrayRead(xcoorl,&zctx.xy);CHKERRQ(ierr);
311   ierr = VecRestoreArrayRead(xlocal,&zctx.v);CHKERRQ(ierr);
312   PetscFunctionReturn(0);
313 }
314 
315 #if defined(PETSC_HAVE_HDF5)
316 #undef __FUNCT__
317 #define __FUNCT__ "VecGetHDF5ChunkSize"
318 static PetscErrorCode VecGetHDF5ChunkSize(DM_DA *da, Vec xin, PetscInt dimension, PetscInt timestep, hsize_t *chunkDims)
319 {
320   PetscMPIInt    comm_size;
321   PetscErrorCode ierr;
322   hsize_t        chunk_size, target_size, dim;
323   hsize_t        vec_size = sizeof(PetscScalar)*da->M*da->N*da->P*da->w;
324   hsize_t        avg_local_vec_size,KiB = 1024,MiB = KiB*KiB,GiB = MiB*KiB,min_size = MiB;
325   hsize_t        max_chunks = 64*KiB;                                              /* HDF5 internal limitation */
326   hsize_t        max_chunk_size = 4*GiB;                                           /* HDF5 internal limitation */
327   PetscInt       zslices=da->p, yslices=da->n, xslices=da->m;
328 
329   PetscFunctionBegin;
330   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)xin), &comm_size);CHKERRQ(ierr);
331   avg_local_vec_size = (hsize_t) ceil(vec_size*1.0/comm_size);      /* we will attempt to use this as the chunk size */
332 
333   target_size = (hsize_t) ceil(PetscMin(vec_size,PetscMin(max_chunk_size,PetscMax(avg_local_vec_size,PetscMax(ceil(vec_size*1.0/max_chunks),min_size)))));
334   /* following line uses sizeof(PetscReal) instead of sizeof(PetscScalar) because the last dimension of chunkDims[] captures the 2* when complex numbers are being used */
335   chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(PetscReal);
336 
337   /*
338    if size/rank > max_chunk_size, we need radical measures: even going down to
339    avg_local_vec_size is not enough, so we simply use chunk size of 4 GiB no matter
340    what, composed in the most efficient way possible.
341    N.B. this minimises the number of chunks, which may or may not be the optimal
342    solution. In a BG, for example, the optimal solution is probably to make # chunks = #
343    IO nodes involved, but this author has no access to a BG to figure out how to
344    reliably find the right number. And even then it may or may not be enough.
345    */
346   if (avg_local_vec_size > max_chunk_size) {
347     /* check if we can just split local z-axis: is that enough? */
348     zslices = (PetscInt)ceil(vec_size*1.0/(da->p*max_chunk_size))*zslices;
349     if (zslices > da->P) {
350       /* lattice is too large in xy-directions, splitting z only is not enough */
351       zslices = da->P;
352       yslices= (PetscInt)ceil(vec_size*1.0/(zslices*da->n*max_chunk_size))*yslices;
353       if (yslices > da->N) {
354 	/* lattice is too large in x-direction, splitting along z, y is not enough */
355 	yslices = da->N;
356 	xslices= (PetscInt)ceil(vec_size*1.0/(zslices*yslices*da->m*max_chunk_size))*xslices;
357       }
358     }
359     dim = 0;
360     if (timestep >= 0) {
361       ++dim;
362     }
363     /* prefer to split z-axis, even down to planar slices */
364     if (dimension == 3) {
365       chunkDims[dim++] = (hsize_t) da->P/zslices;
366       chunkDims[dim++] = (hsize_t) da->N/yslices;
367       chunkDims[dim++] = (hsize_t) da->M/xslices;
368     } else {
369       /* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */
370       chunkDims[dim++] = (hsize_t) da->N/yslices;
371       chunkDims[dim++] = (hsize_t) da->M/xslices;
372     }
373     chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(double);
374   } else {
375     if (target_size < chunk_size) {
376       /* only change the defaults if target_size < chunk_size */
377       dim = 0;
378       if (timestep >= 0) {
379 	++dim;
380       }
381       /* prefer to split z-axis, even down to planar slices */
382       if (dimension == 3) {
383 	/* try splitting the z-axis to core-size bits, i.e. divide chunk size by # comm_size in z-direction */
384 	if (target_size >= chunk_size/da->p) {
385 	  /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */
386 	  chunkDims[dim] = (hsize_t) ceil(da->P*1.0/da->p);
387 	} else {
388 	  /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be
389            radical and let everyone write all they've got */
390 	  chunkDims[dim++] = (hsize_t) ceil(da->P*1.0/da->p);
391 	  chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n);
392 	  chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m);
393 	}
394       } else {
395 	/* This is a 2D world exceeding 4GiB in size; yes, I've seen them, even used myself */
396 	if (target_size >= chunk_size/da->n) {
397 	  /* just make chunks the size of <local_z>x<whole_world_y>x<whole_world_x>x<dof> */
398 	  chunkDims[dim] = (hsize_t) ceil(da->N*1.0/da->n);
399 	} else {
400 	  /* oops, just splitting the z-axis is NOT ENOUGH, need to split more; let's be
401 	   radical and let everyone write all they've got */
402 	  chunkDims[dim++] = (hsize_t) ceil(da->N*1.0/da->n);
403 	  chunkDims[dim++] = (hsize_t) ceil(da->M*1.0/da->m);
404 	}
405 
406       }
407       chunk_size = (hsize_t) PetscMax(1,chunkDims[0])*PetscMax(1,chunkDims[1])*PetscMax(1,chunkDims[2])*PetscMax(1,chunkDims[3])*PetscMax(1,chunkDims[4])*PetscMax(1,chunkDims[5])*sizeof(double);
408     } else {
409       /* precomputed chunks are fine, we don't need to do anything */
410     }
411   }
412   PetscFunctionReturn(0);
413 }
414 #endif
415 
416 #if defined(PETSC_HAVE_HDF5)
417 #undef __FUNCT__
418 #define __FUNCT__ "VecView_MPI_HDF5_DA"
419 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer)
420 {
421   DM                dm;
422   DM_DA             *da;
423   hid_t             filespace;  /* file dataspace identifier */
424   hid_t             chunkspace; /* chunk dataset property identifier */
425   hid_t             plist_id;   /* property list identifier */
426   hid_t             dset_id;    /* dataset identifier */
427   hid_t             memspace;   /* memory dataspace identifier */
428   hid_t             file_id;
429   hid_t             group;
430   hid_t             memscalartype; /* scalar type for mem (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
431   hid_t             filescalartype; /* scalar type for file (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
432   hsize_t           dim;
433   hsize_t           maxDims[6]={0}, dims[6]={0}, chunkDims[6]={0}, count[6]={0}, offset[6]={0}; /* we depend on these being sane later on  */
434   PetscInt          timestep, dimension;
435   const PetscScalar *x;
436   const char        *vecname;
437   PetscErrorCode    ierr;
438   PetscBool         dim2;
439   PetscBool         spoutput;
440 
441   PetscFunctionBegin;
442   ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr);
443   ierr = PetscViewerHDF5GetTimestep(viewer, &timestep);CHKERRQ(ierr);
444   ierr = PetscViewerHDF5GetBaseDimension2(viewer,&dim2);CHKERRQ(ierr);
445   ierr = PetscViewerHDF5GetSPOutput(viewer,&spoutput);CHKERRQ(ierr);
446 
447   ierr = VecGetDM(xin,&dm);CHKERRQ(ierr);
448   if (!dm) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
449   da = (DM_DA*)dm->data;
450   ierr = DMGetDimension(dm, &dimension);CHKERRQ(ierr);
451 
452   /* Create the dataspace for the dataset.
453    *
454    * dims - holds the current dimensions of the dataset
455    *
456    * maxDims - holds the maximum dimensions of the dataset (unlimited
457    * for the number of time steps with the current dimensions for the
458    * other dimensions; so only additional time steps can be added).
459    *
460    * chunkDims - holds the size of a single time step (required to
461    * permit extending dataset).
462    */
463   dim = 0;
464   if (timestep >= 0) {
465     dims[dim]      = timestep+1;
466     maxDims[dim]   = H5S_UNLIMITED;
467     chunkDims[dim] = 1;
468     ++dim;
469   }
470   if (dimension == 3) {
471     ierr           = PetscHDF5IntCast(da->P,dims+dim);CHKERRQ(ierr);
472     maxDims[dim]   = dims[dim];
473     chunkDims[dim] = dims[dim];
474     ++dim;
475   }
476   if (dimension > 1) {
477     ierr           = PetscHDF5IntCast(da->N,dims+dim);CHKERRQ(ierr);
478     maxDims[dim]   = dims[dim];
479     chunkDims[dim] = dims[dim];
480     ++dim;
481   }
482   ierr           = PetscHDF5IntCast(da->M,dims+dim);CHKERRQ(ierr);
483   maxDims[dim]   = dims[dim];
484   chunkDims[dim] = dims[dim];
485   ++dim;
486   if (da->w > 1 || dim2) {
487     ierr           = PetscHDF5IntCast(da->w,dims+dim);CHKERRQ(ierr);
488     maxDims[dim]   = dims[dim];
489     chunkDims[dim] = dims[dim];
490     ++dim;
491   }
492 #if defined(PETSC_USE_COMPLEX)
493   dims[dim]      = 2;
494   maxDims[dim]   = dims[dim];
495   chunkDims[dim] = dims[dim];
496   ++dim;
497 #endif
498 
499   ierr = VecGetHDF5ChunkSize(da, xin, dimension, timestep, chunkDims); CHKERRQ(ierr);
500 
501   PetscStackCallHDF5Return(filespace,H5Screate_simple,(dim, dims, maxDims));
502 
503 #if defined(PETSC_USE_REAL_SINGLE)
504   memscalartype = H5T_NATIVE_FLOAT;
505   filescalartype = H5T_NATIVE_FLOAT;
506 #elif defined(PETSC_USE_REAL___FLOAT128)
507 #error "HDF5 output with 128 bit floats not supported."
508 #else
509   memscalartype = H5T_NATIVE_DOUBLE;
510   if (spoutput == PETSC_TRUE) filescalartype = H5T_NATIVE_FLOAT;
511   else filescalartype = H5T_NATIVE_DOUBLE;
512 #endif
513 
514   /* Create the dataset with default properties and close filespace */
515   ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr);
516   if (!H5Lexists(group, vecname, H5P_DEFAULT)) {
517     /* Create chunk */
518     PetscStackCallHDF5Return(chunkspace,H5Pcreate,(H5P_DATASET_CREATE));
519     PetscStackCallHDF5(H5Pset_chunk,(chunkspace, dim, chunkDims));
520 
521 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
522     PetscStackCallHDF5Return(dset_id,H5Dcreate2,(group, vecname, filescalartype, filespace, H5P_DEFAULT, chunkspace, H5P_DEFAULT));
523 #else
524     PetscStackCallHDF5Return(dset_id,H5Dcreate,(group, vecname, filescalartype, filespace, H5P_DEFAULT));
525 #endif
526   } else {
527     PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT));
528     PetscStackCallHDF5(H5Dset_extent,(dset_id, dims));
529   }
530   PetscStackCallHDF5(H5Sclose,(filespace));
531 
532   /* Each process defines a dataset and writes it to the hyperslab in the file */
533   dim = 0;
534   if (timestep >= 0) {
535     offset[dim] = timestep;
536     ++dim;
537   }
538   if (dimension == 3) {ierr = PetscHDF5IntCast(da->zs,offset + dim++);CHKERRQ(ierr);}
539   if (dimension > 1)  {ierr = PetscHDF5IntCast(da->ys,offset + dim++);CHKERRQ(ierr);}
540   ierr = PetscHDF5IntCast(da->xs/da->w,offset + dim++);CHKERRQ(ierr);
541   if (da->w > 1 || dim2) offset[dim++] = 0;
542 #if defined(PETSC_USE_COMPLEX)
543   offset[dim++] = 0;
544 #endif
545   dim = 0;
546   if (timestep >= 0) {
547     count[dim] = 1;
548     ++dim;
549   }
550   if (dimension == 3) {ierr = PetscHDF5IntCast(da->ze - da->zs,count + dim++);CHKERRQ(ierr);}
551   if (dimension > 1)  {ierr = PetscHDF5IntCast(da->ye - da->ys,count + dim++);CHKERRQ(ierr);}
552   ierr = PetscHDF5IntCast((da->xe - da->xs)/da->w,count + dim++);CHKERRQ(ierr);
553   if (da->w > 1 || dim2) {ierr = PetscHDF5IntCast(da->w,count + dim++);CHKERRQ(ierr);}
554 #if defined(PETSC_USE_COMPLEX)
555   count[dim++] = 2;
556 #endif
557   PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
558   PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
559   PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
560 
561   /* Create property list for collective dataset write */
562   PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER));
563 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
564   PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE));
565 #endif
566   /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
567 
568   ierr   = VecGetArrayRead(xin, &x);CHKERRQ(ierr);
569   PetscStackCallHDF5(H5Dwrite,(dset_id, memscalartype, memspace, filespace, plist_id, x));
570   PetscStackCallHDF5(H5Fflush,(file_id, H5F_SCOPE_GLOBAL));
571   ierr   = VecRestoreArrayRead(xin, &x);CHKERRQ(ierr);
572 
573   /* Close/release resources */
574   if (group != file_id) {
575     PetscStackCallHDF5(H5Gclose,(group));
576   }
577   PetscStackCallHDF5(H5Pclose,(plist_id));
578   PetscStackCallHDF5(H5Sclose,(filespace));
579   PetscStackCallHDF5(H5Sclose,(memspace));
580   PetscStackCallHDF5(H5Dclose,(dset_id));
581   ierr   = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr);
582   PetscFunctionReturn(0);
583 }
584 #endif
585 
586 extern PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer);
587 
588 #if defined(PETSC_HAVE_MPIIO)
589 #undef __FUNCT__
590 #define __FUNCT__ "DMDAArrayMPIIO"
591 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write)
592 {
593   PetscErrorCode    ierr;
594   MPI_File          mfdes;
595   PetscMPIInt       gsizes[4],lsizes[4],lstarts[4],asiz,dof;
596   MPI_Datatype      view;
597   const PetscScalar *array;
598   MPI_Offset        off;
599   MPI_Aint          ub,ul;
600   PetscInt          type,rows,vecrows,tr[2];
601   DM_DA             *dd = (DM_DA*)da->data;
602   PetscBool         skipheader;
603 
604   PetscFunctionBegin;
605   ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr);
606   ierr = PetscViewerBinaryGetSkipHeader(viewer,&skipheader);CHKERRQ(ierr);
607   if (!write) {
608     /* Read vector header. */
609     if (!skipheader) {
610       ierr = PetscViewerBinaryRead(viewer,tr,2,NULL,PETSC_INT);CHKERRQ(ierr);
611       type = tr[0];
612       rows = tr[1];
613       if (type != VEC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_WRONG,"Not vector next in file");
614       if (rows != vecrows) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector");
615     }
616   } else {
617     tr[0] = VEC_FILE_CLASSID;
618     tr[1] = vecrows;
619     if (!skipheader) {
620       ierr  = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr);
621     }
622   }
623 
624   ierr       = PetscMPIIntCast(dd->w,&dof);CHKERRQ(ierr);
625   gsizes[0]  = dof;
626   ierr       = PetscMPIIntCast(dd->M,gsizes+1);CHKERRQ(ierr);
627   ierr       = PetscMPIIntCast(dd->N,gsizes+2);CHKERRQ(ierr);
628   ierr       = PetscMPIIntCast(dd->P,gsizes+3);CHKERRQ(ierr);
629   lsizes[0]  = dof;
630   ierr       = PetscMPIIntCast((dd->xe-dd->xs)/dof,lsizes+1);CHKERRQ(ierr);
631   ierr       = PetscMPIIntCast(dd->ye-dd->ys,lsizes+2);CHKERRQ(ierr);
632   ierr       = PetscMPIIntCast(dd->ze-dd->zs,lsizes+3);CHKERRQ(ierr);
633   lstarts[0] = 0;
634   ierr       = PetscMPIIntCast(dd->xs/dof,lstarts+1);CHKERRQ(ierr);
635   ierr       = PetscMPIIntCast(dd->ys,lstarts+2);CHKERRQ(ierr);
636   ierr       = PetscMPIIntCast(dd->zs,lstarts+3);CHKERRQ(ierr);
637   ierr       = MPI_Type_create_subarray(da->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr);
638   ierr       = MPI_Type_commit(&view);CHKERRQ(ierr);
639 
640   ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr);
641   ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr);
642   ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char*)"native",MPI_INFO_NULL);CHKERRQ(ierr);
643   ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr);
644   asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof;
645   if (write) {
646     ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr);
647   } else {
648     ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr);
649   }
650   ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr);
651   ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr);
652   ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr);
653   ierr = MPI_Type_free(&view);CHKERRQ(ierr);
654   PetscFunctionReturn(0);
655 }
656 #endif
657 
658 #undef __FUNCT__
659 #define __FUNCT__ "VecView_MPI_DA"
660 PetscErrorCode  VecView_MPI_DA(Vec xin,PetscViewer viewer)
661 {
662   DM                da;
663   PetscErrorCode    ierr;
664   PetscInt          dim;
665   Vec               natural;
666   PetscBool         isdraw,isvtk;
667 #if defined(PETSC_HAVE_HDF5)
668   PetscBool         ishdf5;
669 #endif
670   const char        *prefix,*name;
671   PetscViewerFormat format;
672 
673   PetscFunctionBegin;
674   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
675   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
676   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr);
677   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isvtk);CHKERRQ(ierr);
678 #if defined(PETSC_HAVE_HDF5)
679   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr);
680 #endif
681   if (isdraw) {
682     ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
683     if (dim == 1) {
684       ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr);
685     } else if (dim == 2) {
686       ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr);
687     } else SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim);
688   } else if (isvtk) {           /* Duplicate the Vec and hold a reference to the DM */
689     Vec Y;
690     ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);
691     ierr = VecDuplicate(xin,&Y);CHKERRQ(ierr);
692     if (((PetscObject)xin)->name) {
693       /* If xin was named, copy the name over to Y. The duplicate names are safe because nobody else will ever see Y. */
694       ierr = PetscObjectSetName((PetscObject)Y,((PetscObject)xin)->name);CHKERRQ(ierr);
695     }
696     ierr = VecCopy(xin,Y);CHKERRQ(ierr);
697     ierr = PetscViewerVTKAddField(viewer,(PetscObject)da,DMDAVTKWriteAll,PETSC_VTK_POINT_FIELD,(PetscObject)Y);CHKERRQ(ierr);
698 #if defined(PETSC_HAVE_HDF5)
699   } else if (ishdf5) {
700     ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr);
701 #endif
702   } else {
703 #if defined(PETSC_HAVE_MPIIO)
704     PetscBool isbinary,isMPIIO;
705 
706     ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
707     if (isbinary) {
708       ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr);
709       if (isMPIIO) {
710         ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr);
711         PetscFunctionReturn(0);
712       }
713     }
714 #endif
715 
716     /* call viewer on natural ordering */
717     ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr);
718     ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
719     ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr);
720     ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);
721     ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr);
722     ierr = PetscObjectGetName((PetscObject)xin,&name);CHKERRQ(ierr);
723     ierr = PetscObjectSetName((PetscObject)natural,name);CHKERRQ(ierr);
724 
725     ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr);
726     if (format == PETSC_VIEWER_BINARY_MATLAB) {
727       /* temporarily remove viewer format so it won't trigger in the VecView() */
728       ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);CHKERRQ(ierr);
729     }
730 
731     ierr = VecView(natural,viewer);CHKERRQ(ierr);
732 
733     if (format == PETSC_VIEWER_BINARY_MATLAB) {
734       MPI_Comm    comm;
735       FILE        *info;
736       const char  *fieldname;
737       char        fieldbuf[256];
738       PetscInt    dim,ni,nj,nk,pi,pj,pk,dof,n;
739 
740       /* set the viewer format back into the viewer */
741       ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);
742       ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
743       ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr);
744       ierr = DMDAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,0,0,0,0,0);CHKERRQ(ierr);
745       ierr = PetscFPrintf(comm,info,"#--- begin code written by PetscViewerBinary for MATLAB format ---#\n");CHKERRQ(ierr);
746       ierr = PetscFPrintf(comm,info,"#$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr);
747       if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); }
748       if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); }
749       if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); }
750 
751       for (n=0; n<dof; n++) {
752         ierr = DMDAGetFieldName(da,n,&fieldname);CHKERRQ(ierr);
753         if (!fieldname || !fieldname[0]) {
754           ierr = PetscSNPrintf(fieldbuf,sizeof fieldbuf,"field%D",n);CHKERRQ(ierr);
755           fieldname = fieldbuf;
756         }
757         if (dim == 1) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
758         if (dim == 2) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
759         if (dim == 3) { ierr = PetscFPrintf(comm,info,"#$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);}
760       }
761       ierr = PetscFPrintf(comm,info,"#$$ clear tmp; \n");CHKERRQ(ierr);
762       ierr = PetscFPrintf(comm,info,"#--- end code written by PetscViewerBinary for MATLAB format ---#\n\n");CHKERRQ(ierr);
763     }
764 
765     ierr = VecDestroy(&natural);CHKERRQ(ierr);
766   }
767   PetscFunctionReturn(0);
768 }
769 
770 #if defined(PETSC_HAVE_HDF5)
771 #undef __FUNCT__
772 #define __FUNCT__ "VecLoad_HDF5_DA"
773 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer)
774 {
775   DM             da;
776   PetscErrorCode ierr;
777   hsize_t        dim,rdim;
778   hsize_t        dims[6]={0},count[6]={0},offset[6]={0};
779   PetscInt       dimension,timestep,dofInd;
780   PetscScalar    *x;
781   const char     *vecname;
782   hid_t          filespace; /* file dataspace identifier */
783   hid_t          plist_id;  /* property list identifier */
784   hid_t          dset_id;   /* dataset identifier */
785   hid_t          memspace;  /* memory dataspace identifier */
786   hid_t          file_id,group;
787   hid_t          scalartype; /* scalar type (H5T_NATIVE_FLOAT or H5T_NATIVE_DOUBLE) */
788   DM_DA          *dd;
789   PetscBool      dim2 = PETSC_FALSE;
790 
791   PetscFunctionBegin;
792 #if defined(PETSC_USE_REAL_SINGLE)
793   scalartype = H5T_NATIVE_FLOAT;
794 #elif defined(PETSC_USE_REAL___FLOAT128)
795 #error "HDF5 output with 128 bit floats not supported."
796 #else
797   scalartype = H5T_NATIVE_DOUBLE;
798 #endif
799 
800   ierr = PetscViewerHDF5OpenGroup(viewer, &file_id, &group);CHKERRQ(ierr);
801   ierr = PetscViewerHDF5GetTimestep(viewer, &timestep);CHKERRQ(ierr);
802   ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr);
803   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
804   dd   = (DM_DA*)da->data;
805   ierr = DMGetDimension(da, &dimension);CHKERRQ(ierr);
806 
807   /* Open dataset */
808 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800)
809   PetscStackCallHDF5Return(dset_id,H5Dopen2,(group, vecname, H5P_DEFAULT));
810 #else
811   PetscStackCallHDF5Return(dset_id,H5Dopen,(group, vecname));
812 #endif
813 
814   /* Retrieve the dataspace for the dataset */
815   PetscStackCallHDF5Return(filespace,H5Dget_space,(dset_id));
816   PetscStackCallHDF5Return(rdim,H5Sget_simple_extent_dims,(filespace, dims, NULL));
817 
818   /* Expected dimension for holding the dof's */
819 #if defined(PETSC_USE_COMPLEX)
820   dofInd = rdim-2;
821 #else
822   dofInd = rdim-1;
823 #endif
824 
825   /* The expected number of dimensions, assuming basedimension2 = false */
826   dim = dimension;
827   if (dd->w > 1) ++dim;
828   if (timestep >= 0) ++dim;
829 #if defined(PETSC_USE_COMPLEX)
830   ++dim;
831 #endif
832 
833   /* In this case the input dataset have one extra, unexpected dimension. */
834   if (rdim == dim+1) {
835     /* In this case the block size unity */
836     if (dd->w == 1 && dims[dofInd] == 1) dim2 = PETSC_TRUE;
837 
838     /* Special error message for the case where dof does not match the input file */
839     else if (dd->w != (PetscInt) dims[dofInd]) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Number of dofs in file is %D, not %D as expected",(PetscInt)dims[dofInd],dd->w);
840 
841   /* Other cases where rdim != dim cannot be handled currently */
842   } else if (rdim != dim) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Dimension of array in file is %d, not %d as expected with dof = %D",rdim,dim,dd->w);
843 
844   /* Set up the hyperslab size */
845   dim = 0;
846   if (timestep >= 0) {
847     offset[dim] = timestep;
848     count[dim] = 1;
849     ++dim;
850   }
851   if (dimension == 3) {
852     ierr = PetscHDF5IntCast(dd->zs,offset + dim);CHKERRQ(ierr);
853     ierr = PetscHDF5IntCast(dd->ze - dd->zs,count + dim);CHKERRQ(ierr);
854     ++dim;
855   }
856   if (dimension > 1) {
857     ierr = PetscHDF5IntCast(dd->ys,offset + dim);CHKERRQ(ierr);
858     ierr = PetscHDF5IntCast(dd->ye - dd->ys,count + dim);CHKERRQ(ierr);
859     ++dim;
860   }
861   ierr = PetscHDF5IntCast(dd->xs/dd->w,offset + dim);CHKERRQ(ierr);
862   ierr = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w,count + dim);CHKERRQ(ierr);
863   ++dim;
864   if (dd->w > 1 || dim2) {
865     offset[dim] = 0;
866     ierr = PetscHDF5IntCast(dd->w,count + dim);CHKERRQ(ierr);
867     ++dim;
868   }
869 #if defined(PETSC_USE_COMPLEX)
870   offset[dim] = 0;
871   count[dim] = 2;
872   ++dim;
873 #endif
874 
875   /* Create the memory and filespace */
876   PetscStackCallHDF5Return(memspace,H5Screate_simple,(dim, count, NULL));
877   PetscStackCallHDF5(H5Sselect_hyperslab,(filespace, H5S_SELECT_SET, offset, NULL, count, NULL));
878 
879   /* Create property list for collective dataset write */
880   PetscStackCallHDF5Return(plist_id,H5Pcreate,(H5P_DATASET_XFER));
881 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO)
882   PetscStackCallHDF5(H5Pset_dxpl_mpio,(plist_id, H5FD_MPIO_COLLECTIVE));
883 #endif
884   /* To read dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */
885 
886   ierr   = VecGetArray(xin, &x);CHKERRQ(ierr);
887   PetscStackCallHDF5(H5Dread,(dset_id, scalartype, memspace, filespace, plist_id, x));
888   ierr   = VecRestoreArray(xin, &x);CHKERRQ(ierr);
889 
890   /* Close/release resources */
891   if (group != file_id) {
892     PetscStackCallHDF5(H5Gclose,(group));
893   }
894   PetscStackCallHDF5(H5Pclose,(plist_id));
895   PetscStackCallHDF5(H5Sclose,(filespace));
896   PetscStackCallHDF5(H5Sclose,(memspace));
897   PetscStackCallHDF5(H5Dclose,(dset_id));
898   PetscFunctionReturn(0);
899 }
900 #endif
901 
902 #undef __FUNCT__
903 #define __FUNCT__ "VecLoad_Binary_DA"
904 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer)
905 {
906   DM             da;
907   PetscErrorCode ierr;
908   Vec            natural;
909   const char     *prefix;
910   PetscInt       bs;
911   PetscBool      flag;
912   DM_DA          *dd;
913 #if defined(PETSC_HAVE_MPIIO)
914   PetscBool isMPIIO;
915 #endif
916 
917   PetscFunctionBegin;
918   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
919   dd   = (DM_DA*)da->data;
920 #if defined(PETSC_HAVE_MPIIO)
921   ierr = PetscViewerBinaryGetUseMPIIO(viewer,&isMPIIO);CHKERRQ(ierr);
922   if (isMPIIO) {
923     ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr);
924     PetscFunctionReturn(0);
925   }
926 #endif
927 
928   ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr);
929   ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr);
930   ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr);
931   ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr);
932   ierr = VecLoad(natural,viewer);CHKERRQ(ierr);
933   ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr);
934   ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr);
935   ierr = VecDestroy(&natural);CHKERRQ(ierr);
936   ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr);
937   ierr = PetscOptionsGetInt(NULL,((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr);
938   if (flag && bs != dd->w) {
939     ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr);
940   }
941   PetscFunctionReturn(0);
942 }
943 
944 #undef __FUNCT__
945 #define __FUNCT__ "VecLoad_Default_DA"
946 PetscErrorCode  VecLoad_Default_DA(Vec xin, PetscViewer viewer)
947 {
948   PetscErrorCode ierr;
949   DM             da;
950   PetscBool      isbinary;
951 #if defined(PETSC_HAVE_HDF5)
952   PetscBool ishdf5;
953 #endif
954 
955   PetscFunctionBegin;
956   ierr = VecGetDM(xin,&da);CHKERRQ(ierr);
957   if (!da) SETERRQ(PetscObjectComm((PetscObject)xin),PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
958 
959 #if defined(PETSC_HAVE_HDF5)
960   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr);
961 #endif
962   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr);
963 
964   if (isbinary) {
965     ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr);
966 #if defined(PETSC_HAVE_HDF5)
967   } else if (ishdf5) {
968     ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr);
969 #endif
970   } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name);
971   PetscFunctionReturn(0);
972 }
973