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