1 #define PETSCDM_DLL 2 3 /* 4 Plots vectors obtained with DMDACreate2d() 5 */ 6 7 #include "private/daimpl.h" /*I "petscdm.h" I*/ 8 #include "private/vecimpl.h" 9 10 /* 11 The data that is passed into the graphics callback 12 */ 13 typedef struct { 14 PetscInt m,n,step,k; 15 PetscReal min,max,scale; 16 PetscScalar *xy,*v; 17 PetscBool showgrid; 18 } ZoomCtx; 19 20 /* 21 This does the drawing for one particular field 22 in one particular set of coordinates. It is a callback 23 called from PetscDrawZoom() 24 */ 25 #undef __FUNCT__ 26 #define __FUNCT__ "VecView_MPI_Draw_DA2d_Zoom" 27 PetscErrorCode VecView_MPI_Draw_DA2d_Zoom(PetscDraw draw,void *ctx) 28 { 29 ZoomCtx *zctx = (ZoomCtx*)ctx; 30 PetscErrorCode ierr; 31 PetscInt m,n,i,j,k,step,id,c1,c2,c3,c4; 32 PetscReal s,min,x1,x2,x3,x4,y_1,y2,y3,y4; 33 PetscScalar *v,*xy; 34 35 PetscFunctionBegin; 36 m = zctx->m; 37 n = zctx->n; 38 step = zctx->step; 39 k = zctx->k; 40 v = zctx->v; 41 xy = zctx->xy; 42 s = zctx->scale; 43 min = zctx->min; 44 45 /* PetscDraw the contour plot patch */ 46 for (j=0; j<n-1; j++) { 47 for (i=0; i<m-1; i++) { 48 #if !defined(PETSC_USE_COMPLEX) 49 id = i+j*m; x1 = xy[2*id];y_1 = xy[2*id+1];c1 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 50 id = i+j*m+1; x2 = xy[2*id];y2 = y_1; c2 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 51 id = i+j*m+1+m;x3 = x2; y3 = xy[2*id+1];c3 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 52 id = i+j*m+m; x4 = x1; y4 = y3; c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(v[k+step*id]-min)); 53 #else 54 id = i+j*m; x1 = PetscRealPart(xy[2*id]);y_1 = PetscRealPart(xy[2*id+1]);c1 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 55 id = i+j*m+1; x2 = PetscRealPart(xy[2*id]);y2 = y_1; c2 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 56 id = i+j*m+1+m;x3 = x2; y3 = PetscRealPart(xy[2*id+1]);c3 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 57 id = i+j*m+m; x4 = x1; y4 = y3; c4 = (int)(PETSC_DRAW_BASIC_COLORS+s*(PetscRealPart(v[k+step*id])-min)); 58 #endif 59 ierr = PetscDrawTriangle(draw,x1,y_1,x2,y2,x3,y3,c1,c2,c3);CHKERRQ(ierr); 60 ierr = PetscDrawTriangle(draw,x1,y_1,x3,y3,x4,y4,c1,c3,c4);CHKERRQ(ierr); 61 if (zctx->showgrid) { 62 ierr = PetscDrawLine(draw,x1,y_1,x2,y2,PETSC_DRAW_BLACK);CHKERRQ(ierr); 63 ierr = PetscDrawLine(draw,x2,y2,x3,y3,PETSC_DRAW_BLACK);CHKERRQ(ierr); 64 ierr = PetscDrawLine(draw,x3,y3,x4,y4,PETSC_DRAW_BLACK);CHKERRQ(ierr); 65 ierr = PetscDrawLine(draw,x4,y4,x1,y_1,PETSC_DRAW_BLACK);CHKERRQ(ierr); 66 } 67 } 68 } 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "VecView_MPI_Draw_DA2d" 74 PetscErrorCode VecView_MPI_Draw_DA2d(Vec xin,PetscViewer viewer) 75 { 76 DM da,dac,dag; 77 PetscErrorCode ierr; 78 PetscMPIInt rank; 79 PetscInt igstart,N,s,M,istart,isize,jgstart,w; 80 const PetscInt *lx,*ly; 81 PetscReal coors[4],ymin,ymax,xmin,xmax; 82 PetscDraw draw,popup; 83 PetscBool isnull,useports = PETSC_FALSE; 84 MPI_Comm comm; 85 Vec xlocal,xcoor,xcoorl; 86 DMDAPeriodicType periodic; 87 DMDAStencilType st; 88 ZoomCtx zctx; 89 PetscDrawViewPorts *ports; 90 PetscViewerFormat format; 91 92 PetscFunctionBegin; 93 zctx.showgrid = PETSC_FALSE; 94 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 95 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 96 97 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 98 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 99 100 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 101 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 102 103 ierr = DMDAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&periodic,&st);CHKERRQ(ierr); 104 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,PETSC_NULL);CHKERRQ(ierr); 105 106 /* 107 Obtain a sequential vector that is going to contain the local values plus ONE layer of 108 ghosted values to draw the graphics from. We also need its corresponding DMDA (dac) that will 109 update the local values pluse ONE layer of ghost values. 110 */ 111 ierr = PetscObjectQuery((PetscObject)da,"GraphicsGhosted",(PetscObject*)&xlocal);CHKERRQ(ierr); 112 if (!xlocal) { 113 if (periodic != DMDA_NONPERIODIC || s != 1 || st != DMDA_STENCIL_BOX) { 114 /* 115 if original da is not of stencil width one, or periodic or not a box stencil then 116 create a special DMDA to handle one level of ghost points for graphics 117 */ 118 ierr = DMDACreate2d(comm,DMDA_NONPERIODIC,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);CHKERRQ(ierr); 119 ierr = PetscInfo(da,"Creating auxilary DMDA for managing graphics ghost points\n");CHKERRQ(ierr); 120 } else { 121 /* otherwise we can use the da we already have */ 122 dac = da; 123 } 124 /* create local vector for holding ghosted values used in graphics */ 125 ierr = DMCreateLocalVector(dac,&xlocal);CHKERRQ(ierr); 126 if (dac != da) { 127 /* don't keep any public reference of this DMDA, is is only available through xlocal */ 128 ierr = DMDestroy(dac);CHKERRQ(ierr); 129 } else { 130 /* remove association between xlocal and da, because below we compose in the opposite 131 direction and if we left this connect we'd get a loop, so the objects could 132 never be destroyed */ 133 ierr = PetscObjectCompose((PetscObject)xlocal,"DMDA",0);CHKERRQ(ierr); 134 } 135 ierr = PetscObjectCompose((PetscObject)da,"GraphicsGhosted",(PetscObject)xlocal);CHKERRQ(ierr); 136 ierr = PetscObjectDereference((PetscObject)xlocal);CHKERRQ(ierr); 137 } else { 138 if (periodic == DMDA_NONPERIODIC && s == 1 && st == DMDA_STENCIL_BOX) { 139 dac = da; 140 } else { 141 ierr = PetscObjectQuery((PetscObject)xlocal,"DMDA",(PetscObject*)&dac);CHKERRQ(ierr); 142 } 143 } 144 145 /* 146 Get local (ghosted) values of vector 147 */ 148 ierr = DMGlobalToLocalBegin(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 149 ierr = DMGlobalToLocalEnd(dac,xin,INSERT_VALUES,xlocal);CHKERRQ(ierr); 150 ierr = VecGetArray(xlocal,&zctx.v);CHKERRQ(ierr); 151 152 /* get coordinates of nodes */ 153 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 154 if (!xcoor) { 155 ierr = DMDASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr); 156 ierr = DMDAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 157 } 158 159 /* 160 Determine the min and max coordinates in plot 161 */ 162 ierr = VecStrideMin(xcoor,0,PETSC_NULL,&xmin);CHKERRQ(ierr); 163 ierr = VecStrideMax(xcoor,0,PETSC_NULL,&xmax);CHKERRQ(ierr); 164 ierr = VecStrideMin(xcoor,1,PETSC_NULL,&ymin);CHKERRQ(ierr); 165 ierr = VecStrideMax(xcoor,1,PETSC_NULL,&ymax);CHKERRQ(ierr); 166 coors[0] = xmin - .05*(xmax- xmin); coors[2] = xmax + .05*(xmax - xmin); 167 coors[1] = ymin - .05*(ymax- ymin); coors[3] = ymax + .05*(ymax - ymin); 168 ierr = PetscInfo4(da,"Preparing DMDA 2d contour plot coordinates %G %G %G %G\n",coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 169 170 /* 171 get local ghosted version of coordinates 172 */ 173 ierr = PetscObjectQuery((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject*)&xcoorl);CHKERRQ(ierr); 174 if (!xcoorl) { 175 /* create DMDA to get local version of graphics */ 176 ierr = DMDACreate2d(comm,DMDA_NONPERIODIC,DMDA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);CHKERRQ(ierr); 177 ierr = PetscInfo(dag,"Creating auxilary DMDA for managing graphics coordinates ghost points\n");CHKERRQ(ierr); 178 ierr = DMCreateLocalVector(dag,&xcoorl);CHKERRQ(ierr); 179 ierr = PetscObjectCompose((PetscObject)da,"GraphicsCoordinateGhosted",(PetscObject)xcoorl);CHKERRQ(ierr); 180 ierr = DMDestroy(dag);CHKERRQ(ierr);/* dereference dag */ 181 ierr = PetscObjectDereference((PetscObject)xcoorl);CHKERRQ(ierr); 182 } else { 183 ierr = PetscObjectQuery((PetscObject)xcoorl,"DMDA",(PetscObject*)&dag);CHKERRQ(ierr); 184 } 185 ierr = DMGlobalToLocalBegin(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 186 ierr = DMGlobalToLocalEnd(dag,xcoor,INSERT_VALUES,xcoorl);CHKERRQ(ierr); 187 ierr = VecGetArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 188 189 /* 190 Get information about size of area each processor must do graphics for 191 */ 192 ierr = DMDAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&periodic,0);CHKERRQ(ierr); 193 ierr = DMDAGetGhostCorners(dac,&igstart,&jgstart,0,&zctx.m,&zctx.n,0);CHKERRQ(ierr); 194 ierr = DMDAGetCorners(dac,&istart,0,0,&isize,0,0);CHKERRQ(ierr); 195 196 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_contour_grid",&zctx.showgrid,PETSC_NULL);CHKERRQ(ierr); 197 198 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 199 ierr = PetscOptionsGetBool(PETSC_NULL,"-draw_ports",&useports,PETSC_NULL);CHKERRQ(ierr); 200 if (useports || format == PETSC_VIEWER_DRAW_PORTS){ 201 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 202 ierr = PetscDrawViewPortsCreate(draw,zctx.step,&ports);CHKERRQ(ierr); 203 } 204 /* 205 Loop over each field; drawing each in a different window 206 */ 207 for (zctx.k=0; zctx.k<zctx.step; zctx.k++) { 208 if (useports) { 209 ierr = PetscDrawViewPortsSet(ports,zctx.k);CHKERRQ(ierr); 210 } else { 211 ierr = PetscViewerDrawGetDraw(viewer,zctx.k,&draw);CHKERRQ(ierr); 212 ierr = PetscDrawSynchronizedClear(draw);CHKERRQ(ierr); 213 } 214 215 /* 216 Determine the min and max color in plot 217 */ 218 ierr = VecStrideMin(xin,zctx.k,PETSC_NULL,&zctx.min);CHKERRQ(ierr); 219 ierr = VecStrideMax(xin,zctx.k,PETSC_NULL,&zctx.max);CHKERRQ(ierr); 220 if (zctx.min == zctx.max) { 221 zctx.min -= 1.e-12; 222 zctx.max += 1.e-12; 223 } 224 225 if (!rank) { 226 const char *title; 227 228 ierr = DMDAGetFieldName(da,zctx.k,&title);CHKERRQ(ierr); 229 if (title) { 230 ierr = PetscDrawSetTitle(draw,title);CHKERRQ(ierr); 231 } 232 } 233 ierr = PetscDrawSetCoordinates(draw,coors[0],coors[1],coors[2],coors[3]);CHKERRQ(ierr); 234 ierr = PetscInfo2(da,"DMDA 2d contour plot min %G max %G\n",zctx.min,zctx.max);CHKERRQ(ierr); 235 236 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 237 if (popup) {ierr = PetscDrawScalePopup(popup,zctx.min,zctx.max);CHKERRQ(ierr);} 238 239 zctx.scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/(zctx.max - zctx.min); 240 241 ierr = PetscDrawZoom(draw,VecView_MPI_Draw_DA2d_Zoom,&zctx);CHKERRQ(ierr); 242 } 243 if (useports){ 244 ierr = PetscDrawViewPortsDestroy(ports);CHKERRQ(ierr); 245 } 246 247 ierr = VecRestoreArray(xcoorl,&zctx.xy);CHKERRQ(ierr); 248 ierr = VecRestoreArray(xlocal,&zctx.v);CHKERRQ(ierr); 249 PetscFunctionReturn(0); 250 } 251 252 253 #if defined(PETSC_HAVE_HDF5) 254 #undef __FUNCT__ 255 #define __FUNCT__ "VecView_MPI_HDF5_DA" 256 PetscErrorCode VecView_MPI_HDF5_DA(Vec xin,PetscViewer viewer) 257 { 258 PetscErrorCode ierr; 259 DM dm; 260 DM_DA *da; 261 hsize_t dim,dims[5]; 262 hsize_t count[5]; 263 hsize_t offset[5]; 264 PetscInt cnt = 0; 265 PetscScalar *x; 266 const char *vecname; 267 hid_t filespace; /* file dataspace identifier */ 268 hid_t plist_id; /* property list identifier */ 269 hid_t dset_id; /* dataset identifier */ 270 hid_t memspace; /* memory dataspace identifier */ 271 hid_t file_id; 272 herr_t status; 273 274 PetscFunctionBegin; 275 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 276 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&dm);CHKERRQ(ierr); 277 if (!dm) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 278 da = (DM_DA*)dm->data; 279 280 /* Create the dataspace for the dataset */ 281 dim = PetscHDF5IntCast(da->dim + ((da->w == 1) ? 0 : 1)); 282 if (da->dim == 3) dims[cnt++] = PetscHDF5IntCast(da->P); 283 if (da->dim > 1) dims[cnt++] = PetscHDF5IntCast(da->N); 284 dims[cnt++] = PetscHDF5IntCast(da->M); 285 if (da->w > 1) dims[cnt++] = PetscHDF5IntCast(da->w); 286 #if defined(PETSC_USE_COMPLEX) 287 dim++; 288 dims[cnt++] = 2; 289 #endif 290 filespace = H5Screate_simple(dim, dims, NULL); 291 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 292 293 /* Create the dataset with default properties and close filespace */ 294 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 295 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 296 dset_id = H5Dcreate2(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 297 #else 298 dset_id = H5Dcreate(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT); 299 #endif 300 if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()"); 301 status = H5Sclose(filespace);CHKERRQ(status); 302 303 /* Each process defines a dataset and writes it to the hyperslab in the file */ 304 cnt = 0; 305 if (da->dim == 3) offset[cnt++] = PetscHDF5IntCast(da->zs); 306 if (da->dim > 1) offset[cnt++] = PetscHDF5IntCast(da->ys); 307 offset[cnt++] = PetscHDF5IntCast(da->xs/da->w); 308 if (da->w > 1) offset[cnt++] = 0; 309 #if defined(PETSC_USE_COMPLEX) 310 offset[cnt++] = 0; 311 #endif 312 cnt = 0; 313 if (da->dim == 3) count[cnt++] = PetscHDF5IntCast(da->ze - da->zs); 314 if (da->dim > 1) count[cnt++] = PetscHDF5IntCast(da->ye - da->ys); 315 count[cnt++] = PetscHDF5IntCast((da->xe - da->xs)/da->w); 316 if (da->w > 1) count[cnt++] = PetscHDF5IntCast(da->w); 317 #if defined(PETSC_USE_COMPLEX) 318 count[cnt++] = 2; 319 #endif 320 memspace = H5Screate_simple(dim, count, NULL); 321 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 322 323 324 filespace = H5Dget_space(dset_id); 325 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()"); 326 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 327 328 /* Create property list for collective dataset write */ 329 plist_id = H5Pcreate(H5P_DATASET_XFER); 330 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 331 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 332 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 333 #endif 334 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 335 336 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 337 status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 338 status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status); 339 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 340 341 /* Close/release resources */ 342 status = H5Pclose(plist_id);CHKERRQ(status); 343 status = H5Sclose(filespace);CHKERRQ(status); 344 status = H5Sclose(memspace);CHKERRQ(status); 345 status = H5Dclose(dset_id);CHKERRQ(status); 346 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 347 PetscFunctionReturn(0); 348 } 349 #endif 350 351 EXTERN PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 352 353 #if defined(PETSC_HAVE_MPIIO) 354 #undef __FUNCT__ 355 #define __FUNCT__ "DMDAArrayMPIIO" 356 static PetscErrorCode DMDAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 357 { 358 PetscErrorCode ierr; 359 MPI_File mfdes; 360 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 361 MPI_Datatype view; 362 const PetscScalar *array; 363 MPI_Offset off; 364 MPI_Aint ub,ul; 365 PetscInt type,rows,vecrows,tr[2]; 366 DM_DA *dd = (DM_DA*)da->data; 367 368 PetscFunctionBegin; 369 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 370 if (!write) { 371 /* Read vector header. */ 372 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 373 type = tr[0]; 374 rows = tr[1]; 375 if (type != VEC_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not vector next in file"); 376 if (rows != vecrows) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Vector in file not same size as DMDA vector"); 377 } else { 378 tr[0] = VEC_FILE_CLASSID; 379 tr[1] = vecrows; 380 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 381 } 382 383 dof = PetscMPIIntCast(dd->w); 384 gsizes[0] = dof; gsizes[1] = PetscMPIIntCast(dd->M); gsizes[2] = PetscMPIIntCast(dd->N); gsizes[3] = PetscMPIIntCast(dd->P); 385 lsizes[0] = dof;lsizes[1] = PetscMPIIntCast((dd->xe-dd->xs)/dof); lsizes[2] = PetscMPIIntCast(dd->ye-dd->ys); lsizes[3] = PetscMPIIntCast(dd->ze-dd->zs); 386 lstarts[0] = 0; lstarts[1] = PetscMPIIntCast(dd->xs/dof); lstarts[2] = PetscMPIIntCast(dd->ys); lstarts[3] = PetscMPIIntCast(dd->zs); 387 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 388 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 389 390 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 391 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 392 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char *)"native",MPI_INFO_NULL);CHKERRQ(ierr); 393 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 394 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 395 if (write) { 396 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 397 } else { 398 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 399 } 400 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 401 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 402 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 403 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 404 PetscFunctionReturn(0); 405 } 406 #endif 407 408 EXTERN_C_BEGIN 409 #undef __FUNCT__ 410 #define __FUNCT__ "VecView_MPI_DA" 411 PetscErrorCode PETSCDM_DLLEXPORT VecView_MPI_DA(Vec xin,PetscViewer viewer) 412 { 413 DM da; 414 PetscErrorCode ierr; 415 PetscInt dim; 416 Vec natural; 417 PetscBool isdraw; 418 #if defined(PETSC_HAVE_HDF5) 419 PetscBool ishdf5; 420 #endif 421 const char *prefix; 422 423 PetscFunctionBegin; 424 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 425 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 426 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 427 #if defined(PETSC_HAVE_HDF5) 428 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 429 #endif 430 if (isdraw) { 431 ierr = DMDAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 432 if (dim == 1) { 433 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 434 } else if (dim == 2) { 435 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 436 } else { 437 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DMDA %D",dim); 438 } 439 #if defined(PETSC_HAVE_HDF5) 440 } else if (ishdf5) { 441 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 442 #endif 443 } else { 444 #if defined(PETSC_HAVE_MPIIO) 445 PetscBool isbinary,isMPIIO; 446 447 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 448 if (isbinary) { 449 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 450 if (isMPIIO) { 451 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 452 PetscFunctionReturn(0); 453 } 454 } 455 #endif 456 457 /* call viewer on natural ordering */ 458 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 459 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 460 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 461 ierr = DMDAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 462 ierr = DMDAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 463 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 464 ierr = VecView(natural,viewer);CHKERRQ(ierr); 465 ierr = VecDestroy(natural);CHKERRQ(ierr); 466 } 467 PetscFunctionReturn(0); 468 } 469 EXTERN_C_END 470 471 #if defined(PETSC_HAVE_HDF5) 472 #undef __FUNCT__ 473 #define __FUNCT__ "VecLoad_HDF5_DA" 474 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 475 { 476 DM da; 477 PetscErrorCode ierr; 478 hsize_t dim,dims[5]; 479 hsize_t count[5]; 480 hsize_t offset[5]; 481 PetscInt cnt = 0; 482 PetscScalar *x; 483 const char *vecname; 484 hid_t filespace; /* file dataspace identifier */ 485 hid_t plist_id; /* property list identifier */ 486 hid_t dset_id; /* dataset identifier */ 487 hid_t memspace; /* memory dataspace identifier */ 488 hid_t file_id; 489 herr_t status; 490 DM_DA *dd = (DM_DA*)da->data; 491 492 PetscFunctionBegin; 493 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 494 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 495 496 /* Create the dataspace for the dataset */ 497 dim = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1)); 498 if (dd->dim == 3) dims[cnt++] = PetscHDF5IntCast(dd->P); 499 if (dd->dim > 1) dims[cnt++] = PetscHDF5IntCast(dd->N); 500 dims[cnt++] = PetscHDF5IntCast(dd->M); 501 if (dd->w > 1) PetscHDF5IntCast(dims[cnt++] = dd->w); 502 #if defined(PETSC_USE_COMPLEX) 503 dim++; 504 dims[cnt++] = 2; 505 #endif 506 507 /* Create the dataset with default properties and close filespace */ 508 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 509 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 510 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 511 #else 512 dset_id = H5Dopen(file_id, vecname); 513 #endif 514 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 515 filespace = H5Dget_space(dset_id); 516 517 /* Each process defines a dataset and reads it from the hyperslab in the file */ 518 cnt = 0; 519 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 520 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 521 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 522 if (dd->w > 1) offset[cnt++] = 0; 523 #if defined(PETSC_USE_COMPLEX) 524 offset[cnt++] = 0; 525 #endif 526 cnt = 0; 527 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 528 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 529 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 530 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 531 #if defined(PETSC_USE_COMPLEX) 532 count[cnt++] = 2; 533 #endif 534 memspace = H5Screate_simple(dim, count, NULL); 535 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 536 537 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 538 539 /* Create property list for collective dataset write */ 540 plist_id = H5Pcreate(H5P_DATASET_XFER); 541 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 542 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 543 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 544 #endif 545 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 546 547 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 548 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 549 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 550 551 /* Close/release resources */ 552 status = H5Pclose(plist_id);CHKERRQ(status); 553 status = H5Sclose(filespace);CHKERRQ(status); 554 status = H5Sclose(memspace);CHKERRQ(status); 555 status = H5Dclose(dset_id);CHKERRQ(status); 556 PetscFunctionReturn(0); 557 } 558 #endif 559 560 #undef __FUNCT__ 561 #define __FUNCT__ "VecLoad_Binary_DA" 562 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 563 { 564 DM da; 565 PetscErrorCode ierr; 566 Vec natural; 567 const char *prefix; 568 PetscInt bs; 569 PetscBool flag; 570 DM_DA *dd; 571 #if defined(PETSC_HAVE_MPIIO) 572 PetscBool isMPIIO; 573 #endif 574 575 PetscFunctionBegin; 576 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 577 dd = (DM_DA*)da->data; 578 #if defined(PETSC_HAVE_MPIIO) 579 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 580 if (isMPIIO) { 581 ierr = DMDAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 582 PetscFunctionReturn(0); 583 } 584 #endif 585 586 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 587 ierr = DMDACreateNaturalVector(da,&natural);CHKERRQ(ierr); 588 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 589 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 590 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 591 ierr = DMDANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 592 ierr = DMDANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 593 ierr = VecDestroy(natural);CHKERRQ(ierr); 594 ierr = PetscInfo(xin,"Loading vector from natural ordering into DMDA\n");CHKERRQ(ierr); 595 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 596 if (flag && bs != dd->w) { 597 ierr = PetscInfo2(xin,"Block size in file %D not equal to DMDA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 598 } 599 PetscFunctionReturn(0); 600 } 601 602 EXTERN_C_BEGIN 603 #undef __FUNCT__ 604 #define __FUNCT__ "VecLoad_Default_DA" 605 PetscErrorCode PETSCDM_DLLEXPORT VecLoad_Default_DA(Vec xin, PetscViewer viewer) 606 { 607 PetscErrorCode ierr; 608 DM da; 609 PetscBool isbinary; 610 #if defined(PETSC_HAVE_HDF5) 611 PetscBool ishdf5; 612 #endif 613 614 PetscFunctionBegin; 615 ierr = PetscObjectQuery((PetscObject)xin,"DMDA",(PetscObject*)&da);CHKERRQ(ierr); 616 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA"); 617 618 #if defined(PETSC_HAVE_HDF5) 619 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 620 #endif 621 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 622 623 if (isbinary) { 624 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 625 #if defined(PETSC_HAVE_HDF5) 626 } else if (ishdf5) { 627 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 628 #endif 629 } else { 630 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 631 } 632 PetscFunctionReturn(0); 633 } 634 EXTERN_C_END 635