1 #define PETSCDM_DLL 2 3 /* 4 Plots vectors obtained with DACreate2d() 5 */ 6 7 #include "private/daimpl.h" /*I "petscda.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 DAPeriodicType periodic; 87 DAStencilType 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,"DA",(PetscObject*)&da);CHKERRQ(ierr); 98 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DA"); 99 100 ierr = PetscObjectGetComm((PetscObject)xin,&comm);CHKERRQ(ierr); 101 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 102 103 ierr = DAGetInfo(da,0,&M,&N,0,&zctx.m,&zctx.n,0,&w,&s,&periodic,&st);CHKERRQ(ierr); 104 ierr = DAGetOwnershipRanges(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 DA (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 != DA_NONPERIODIC || s != 1 || st != DA_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 DA to handle one level of ghost points for graphics 117 */ 118 ierr = DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,M,N,zctx.m,zctx.n,w,1,lx,ly,&dac);CHKERRQ(ierr); 119 ierr = PetscInfo(da,"Creating auxilary DA 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 = DACreateLocalVector(dac,&xlocal);CHKERRQ(ierr); 126 if (dac != da) { 127 /* don't keep any public reference of this DA, 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,"DA",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 == DA_NONPERIODIC && s == 1 && st == DA_STENCIL_BOX) { 139 dac = da; 140 } else { 141 ierr = PetscObjectQuery((PetscObject)xlocal,"DA",(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 = DAGetCoordinates(da,&xcoor);CHKERRQ(ierr); 154 if (!xcoor) { 155 ierr = DASetUniformCoordinates(da,0.0,1.0,0.0,1.0,0.0,0.0);CHKERRQ(ierr); 156 ierr = DAGetCoordinates(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 DA 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 DA to get local version of graphics */ 176 ierr = DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_BOX,M,N,zctx.m,zctx.n,2,1,lx,ly,&dag);CHKERRQ(ierr); 177 ierr = PetscInfo(dag,"Creating auxilary DA for managing graphics coordinates ghost points\n");CHKERRQ(ierr); 178 ierr = DACreateLocalVector(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,"DA",(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 = DAGetInfo(dac,0,&M,&N,0,0,0,0,&zctx.step,0,&periodic,0);CHKERRQ(ierr); 193 ierr = DAGetGhostCorners(dac,&igstart,&jgstart,0,&zctx.m,&zctx.n,0);CHKERRQ(ierr); 194 ierr = DAGetCorners(dac,&istart,0,0,&isize,0,0);CHKERRQ(ierr); 195 196 ierr = PetscOptionsGetTruth(PETSC_NULL,"-draw_contour_grid",&zctx.showgrid,PETSC_NULL);CHKERRQ(ierr); 197 198 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 199 ierr = PetscOptionsGetTruth(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 = DAGetFieldName(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,"DA 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 da; 260 hsize_t dim,dims[5]; 261 hsize_t count[5]; 262 hsize_t offset[5]; 263 PetscInt cnt = 0; 264 PetscScalar *x; 265 const char *vecname; 266 hid_t filespace; /* file dataspace identifier */ 267 hid_t plist_id; /* property list identifier */ 268 hid_t dset_id; /* dataset identifier */ 269 hid_t memspace; /* memory dataspace identifier */ 270 hid_t file_id; 271 herr_t status; 272 273 PetscFunctionBegin; 274 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 275 ierr = PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);CHKERRQ(ierr); 276 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DA"); 277 278 /* Create the dataspace for the dataset */ 279 dim = PetscHDF5IntCast(da->dim + ((da->w == 1) ? 0 : 1)); 280 if (da->dim == 3) dims[cnt++] = PetscHDF5IntCast(da->P); 281 if (da->dim > 1) dims[cnt++] = PetscHDF5IntCast(da->N); 282 dims[cnt++] = PetscHDF5IntCast(da->M); 283 if (da->w > 1) dims[cnt++] = PetscHDF5IntCast(da->w); 284 #if defined(PETSC_USE_COMPLEX) 285 dim++; 286 dims[cnt++] = 2; 287 #endif 288 filespace = H5Screate_simple(dim, dims, NULL); 289 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 290 291 /* Create the dataset with default properties and close filespace */ 292 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 293 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 294 dset_id = H5Dcreate2(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 295 #else 296 dset_id = H5Dcreate(file_id, vecname, H5T_NATIVE_DOUBLE, filespace, H5P_DEFAULT); 297 #endif 298 if (dset_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dcreate2()"); 299 status = H5Sclose(filespace);CHKERRQ(status); 300 301 /* Each process defines a dataset and writes it to the hyperslab in the file */ 302 cnt = 0; 303 if (da->dim == 3) offset[cnt++] = PetscHDF5IntCast(da->zs); 304 if (da->dim > 1) offset[cnt++] = PetscHDF5IntCast(da->ys); 305 offset[cnt++] = PetscHDF5IntCast(da->xs/da->w); 306 if (da->w > 1) offset[cnt++] = 0; 307 #if defined(PETSC_USE_COMPLEX) 308 offset[cnt++] = 0; 309 #endif 310 cnt = 0; 311 if (da->dim == 3) count[cnt++] = PetscHDF5IntCast(da->ze - da->zs); 312 if (da->dim > 1) count[cnt++] = PetscHDF5IntCast(da->ye - da->ys); 313 count[cnt++] = PetscHDF5IntCast((da->xe - da->xs)/da->w); 314 if (da->w > 1) count[cnt++] = PetscHDF5IntCast(da->w); 315 #if defined(PETSC_USE_COMPLEX) 316 count[cnt++] = 2; 317 #endif 318 memspace = H5Screate_simple(dim, count, NULL); 319 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 320 321 322 filespace = H5Dget_space(dset_id); 323 if (filespace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dget_space()"); 324 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 325 326 /* Create property list for collective dataset write */ 327 plist_id = H5Pcreate(H5P_DATASET_XFER); 328 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 329 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 330 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 331 #endif 332 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 333 334 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 335 status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 336 status = H5Fflush(file_id, H5F_SCOPE_GLOBAL);CHKERRQ(status); 337 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 338 339 /* Close/release resources */ 340 status = H5Pclose(plist_id);CHKERRQ(status); 341 status = H5Sclose(filespace);CHKERRQ(status); 342 status = H5Sclose(memspace);CHKERRQ(status); 343 status = H5Dclose(dset_id);CHKERRQ(status); 344 ierr = PetscInfo1(xin,"Wrote Vec object with name %s\n",vecname);CHKERRQ(ierr); 345 PetscFunctionReturn(0); 346 } 347 #endif 348 349 EXTERN PetscErrorCode VecView_MPI_Draw_DA1d(Vec,PetscViewer); 350 351 #if defined(PETSC_HAVE_MPIIO) 352 #undef __FUNCT__ 353 #define __FUNCT__ "DAArrayMPIIO" 354 static PetscErrorCode DAArrayMPIIO(DM da,PetscViewer viewer,Vec xin,PetscBool write) 355 { 356 PetscErrorCode ierr; 357 MPI_File mfdes; 358 PetscMPIInt gsizes[4],lsizes[4],lstarts[4],asiz,dof; 359 MPI_Datatype view; 360 const PetscScalar *array; 361 MPI_Offset off; 362 MPI_Aint ub,ul; 363 PetscInt type,rows,vecrows,tr[2]; 364 DM_DA *dd = (DM_DA*)da->data; 365 366 PetscFunctionBegin; 367 ierr = VecGetSize(xin,&vecrows);CHKERRQ(ierr); 368 if (!write) { 369 /* Read vector header. */ 370 ierr = PetscViewerBinaryRead(viewer,tr,2,PETSC_INT);CHKERRQ(ierr); 371 type = tr[0]; 372 rows = tr[1]; 373 if (type != VEC_FILE_CLASSID) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_WRONG,"Not vector next in file"); 374 if (rows != vecrows) SETERRQ(((PetscObject)da)->comm,PETSC_ERR_ARG_SIZ,"Vector in file not same size as DA vector"); 375 } else { 376 tr[0] = VEC_FILE_CLASSID; 377 tr[1] = vecrows; 378 ierr = PetscViewerBinaryWrite(viewer,tr,2,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 379 } 380 381 dof = PetscMPIIntCast(dd->w); 382 gsizes[0] = dof; gsizes[1] = PetscMPIIntCast(dd->M); gsizes[2] = PetscMPIIntCast(dd->N); gsizes[3] = PetscMPIIntCast(dd->P); 383 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); 384 lstarts[0] = 0; lstarts[1] = PetscMPIIntCast(dd->xs/dof); lstarts[2] = PetscMPIIntCast(dd->ys); lstarts[3] = PetscMPIIntCast(dd->zs); 385 ierr = MPI_Type_create_subarray(dd->dim+1,gsizes,lsizes,lstarts,MPI_ORDER_FORTRAN,MPIU_SCALAR,&view);CHKERRQ(ierr); 386 ierr = MPI_Type_commit(&view);CHKERRQ(ierr); 387 388 ierr = PetscViewerBinaryGetMPIIODescriptor(viewer,&mfdes);CHKERRQ(ierr); 389 ierr = PetscViewerBinaryGetMPIIOOffset(viewer,&off);CHKERRQ(ierr); 390 ierr = MPI_File_set_view(mfdes,off,MPIU_SCALAR,view,(char *)"native",MPI_INFO_NULL);CHKERRQ(ierr); 391 ierr = VecGetArrayRead(xin,&array);CHKERRQ(ierr); 392 asiz = lsizes[1]*(lsizes[2] > 0 ? lsizes[2] : 1)*(lsizes[3] > 0 ? lsizes[3] : 1)*dof; 393 if (write) { 394 ierr = MPIU_File_write_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 395 } else { 396 ierr = MPIU_File_read_all(mfdes,(PetscScalar*)array,asiz,MPIU_SCALAR,MPI_STATUS_IGNORE);CHKERRQ(ierr); 397 } 398 ierr = MPI_Type_get_extent(view,&ul,&ub);CHKERRQ(ierr); 399 ierr = PetscViewerBinaryAddMPIIOOffset(viewer,ub);CHKERRQ(ierr); 400 ierr = VecRestoreArrayRead(xin,&array);CHKERRQ(ierr); 401 ierr = MPI_Type_free(&view);CHKERRQ(ierr); 402 PetscFunctionReturn(0); 403 } 404 #endif 405 406 EXTERN_C_BEGIN 407 #undef __FUNCT__ 408 #define __FUNCT__ "VecView_MPI_DA" 409 PetscErrorCode PETSCDM_DLLEXPORT VecView_MPI_DA(Vec xin,PetscViewer viewer) 410 { 411 DM da; 412 PetscErrorCode ierr; 413 PetscInt dim; 414 Vec natural; 415 PetscBool isdraw; 416 #if defined(PETSC_HAVE_HDF5) 417 PetscBool ishdf5; 418 #endif 419 const char *prefix; 420 421 PetscFunctionBegin; 422 ierr = PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);CHKERRQ(ierr); 423 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DA"); 424 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 425 #if defined(PETSC_HAVE_HDF5) 426 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 427 #endif 428 if (isdraw) { 429 ierr = DAGetInfo(da,&dim,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); 430 if (dim == 1) { 431 ierr = VecView_MPI_Draw_DA1d(xin,viewer);CHKERRQ(ierr); 432 } else if (dim == 2) { 433 ierr = VecView_MPI_Draw_DA2d(xin,viewer);CHKERRQ(ierr); 434 } else { 435 SETERRQ1(((PetscObject)da)->comm,PETSC_ERR_SUP,"Cannot graphically view vector associated with this dimensional DA %D",dim); 436 } 437 #if defined(PETSC_HAVE_HDF5) 438 } else if (ishdf5) { 439 ierr = VecView_MPI_HDF5_DA(xin,viewer);CHKERRQ(ierr); 440 #endif 441 } else { 442 #if defined(PETSC_HAVE_MPIIO) 443 PetscBool isbinary,isMPIIO; 444 445 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 446 if (isbinary) { 447 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 448 if (isMPIIO) { 449 ierr = DAArrayMPIIO(da,viewer,xin,PETSC_TRUE);CHKERRQ(ierr); 450 PetscFunctionReturn(0); 451 } 452 } 453 #endif 454 455 /* call viewer on natural ordering */ 456 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 457 ierr = DACreateNaturalVector(da,&natural);CHKERRQ(ierr); 458 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 459 ierr = DAGlobalToNaturalBegin(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 460 ierr = DAGlobalToNaturalEnd(da,xin,INSERT_VALUES,natural);CHKERRQ(ierr); 461 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 462 ierr = VecView(natural,viewer);CHKERRQ(ierr); 463 ierr = VecDestroy(natural);CHKERRQ(ierr); 464 } 465 PetscFunctionReturn(0); 466 } 467 EXTERN_C_END 468 469 #if defined(PETSC_HAVE_HDF5) 470 #undef __FUNCT__ 471 #define __FUNCT__ "VecLoad_HDF5_DA" 472 PetscErrorCode VecLoad_HDF5_DA(Vec xin, PetscViewer viewer) 473 { 474 DM da; 475 PetscErrorCode ierr; 476 hsize_t dim,dims[5]; 477 hsize_t count[5]; 478 hsize_t offset[5]; 479 PetscInt cnt = 0; 480 PetscScalar *x; 481 const char *vecname; 482 hid_t filespace; /* file dataspace identifier */ 483 hid_t plist_id; /* property list identifier */ 484 hid_t dset_id; /* dataset identifier */ 485 hid_t memspace; /* memory dataspace identifier */ 486 hid_t file_id; 487 herr_t status; 488 DM_DA *dd = (DM_DA*)da->data; 489 490 PetscFunctionBegin; 491 ierr = PetscViewerHDF5GetFileId(viewer, &file_id);CHKERRQ(ierr); 492 ierr = PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);CHKERRQ(ierr); 493 494 /* Create the dataspace for the dataset */ 495 dim = PetscHDF5IntCast(dd->dim + ((dd->w == 1) ? 0 : 1)); 496 if (dd->dim == 3) dims[cnt++] = PetscHDF5IntCast(dd->P); 497 if (dd->dim > 1) dims[cnt++] = PetscHDF5IntCast(dd->N); 498 dims[cnt++] = PetscHDF5IntCast(dd->M); 499 if (dd->w > 1) PetscHDF5IntCast(dims[cnt++] = dd->w); 500 #if defined(PETSC_USE_COMPLEX) 501 dim++; 502 dims[cnt++] = 2; 503 #endif 504 505 /* Create the dataset with default properties and close filespace */ 506 ierr = PetscObjectGetName((PetscObject)xin,&vecname);CHKERRQ(ierr); 507 #if (H5_VERS_MAJOR * 10000 + H5_VERS_MINOR * 100 + H5_VERS_RELEASE >= 10800) 508 dset_id = H5Dopen2(file_id, vecname, H5P_DEFAULT); 509 #else 510 dset_id = H5Dopen(file_id, vecname); 511 #endif 512 if (dset_id == -1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Dopen2() with Vec named %s",vecname); 513 filespace = H5Dget_space(dset_id); 514 515 /* Each process defines a dataset and reads it from the hyperslab in the file */ 516 cnt = 0; 517 if (dd->dim == 3) offset[cnt++] = PetscHDF5IntCast(dd->zs); 518 if (dd->dim > 1) offset[cnt++] = PetscHDF5IntCast(dd->ys); 519 offset[cnt++] = PetscHDF5IntCast(dd->xs/dd->w); 520 if (dd->w > 1) offset[cnt++] = 0; 521 #if defined(PETSC_USE_COMPLEX) 522 offset[cnt++] = 0; 523 #endif 524 cnt = 0; 525 if (dd->dim == 3) count[cnt++] = PetscHDF5IntCast(dd->ze - dd->zs); 526 if (dd->dim > 1) count[cnt++] = PetscHDF5IntCast(dd->ye - dd->ys); 527 count[cnt++] = PetscHDF5IntCast((dd->xe - dd->xs)/dd->w); 528 if (dd->w > 1) count[cnt++] = PetscHDF5IntCast(dd->w); 529 #if defined(PETSC_USE_COMPLEX) 530 count[cnt++] = 2; 531 #endif 532 memspace = H5Screate_simple(dim, count, NULL); 533 if (memspace == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Screate_simple()"); 534 535 status = H5Sselect_hyperslab(filespace, H5S_SELECT_SET, offset, NULL, count, NULL);CHKERRQ(status); 536 537 /* Create property list for collective dataset write */ 538 plist_id = H5Pcreate(H5P_DATASET_XFER); 539 if (plist_id == -1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Cannot H5Pcreate()"); 540 #if defined(PETSC_HAVE_H5PSET_FAPL_MPIO) 541 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);CHKERRQ(status); 542 #endif 543 /* To write dataset independently use H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_INDEPENDENT) */ 544 545 ierr = VecGetArray(xin, &x);CHKERRQ(ierr); 546 status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, x);CHKERRQ(status); 547 ierr = VecRestoreArray(xin, &x);CHKERRQ(ierr); 548 549 /* Close/release resources */ 550 status = H5Pclose(plist_id);CHKERRQ(status); 551 status = H5Sclose(filespace);CHKERRQ(status); 552 status = H5Sclose(memspace);CHKERRQ(status); 553 status = H5Dclose(dset_id);CHKERRQ(status); 554 PetscFunctionReturn(0); 555 } 556 #endif 557 558 #undef __FUNCT__ 559 #define __FUNCT__ "VecLoad_Binary_DA" 560 PetscErrorCode VecLoad_Binary_DA(Vec xin, PetscViewer viewer) 561 { 562 DM da; 563 PetscErrorCode ierr; 564 Vec natural; 565 const char *prefix; 566 PetscInt bs; 567 PetscBool flag; 568 DM_DA *dd; 569 #if defined(PETSC_HAVE_MPIIO) 570 PetscBool isMPIIO; 571 #endif 572 573 PetscFunctionBegin; 574 ierr = PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);CHKERRQ(ierr); 575 dd = (DM_DA*)da->data; 576 #if defined(PETSC_HAVE_MPIIO) 577 ierr = PetscViewerBinaryGetMPIIO(viewer,&isMPIIO);CHKERRQ(ierr); 578 if (isMPIIO) { 579 ierr = DAArrayMPIIO(da,viewer,xin,PETSC_FALSE);CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 #endif 583 584 ierr = PetscObjectGetOptionsPrefix((PetscObject)xin,&prefix);CHKERRQ(ierr); 585 ierr = DACreateNaturalVector(da,&natural);CHKERRQ(ierr); 586 ierr = PetscObjectSetName((PetscObject)natural,((PetscObject)xin)->name);CHKERRQ(ierr); 587 ierr = PetscObjectSetOptionsPrefix((PetscObject)natural,prefix);CHKERRQ(ierr); 588 ierr = VecLoad_Binary(natural,viewer);CHKERRQ(ierr); 589 ierr = DANaturalToGlobalBegin(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 590 ierr = DANaturalToGlobalEnd(da,natural,INSERT_VALUES,xin);CHKERRQ(ierr); 591 ierr = VecDestroy(natural);CHKERRQ(ierr); 592 ierr = PetscInfo(xin,"Loading vector from natural ordering into DA\n");CHKERRQ(ierr); 593 ierr = PetscOptionsGetInt(((PetscObject)xin)->prefix,"-vecload_block_size",&bs,&flag);CHKERRQ(ierr); 594 if (flag && bs != dd->w) { 595 ierr = PetscInfo2(xin,"Block size in file %D not equal to DA's dof %D\n",bs,dd->w);CHKERRQ(ierr); 596 } 597 PetscFunctionReturn(0); 598 } 599 600 EXTERN_C_BEGIN 601 #undef __FUNCT__ 602 #define __FUNCT__ "VecLoad_Default_DA" 603 PetscErrorCode PETSCDM_DLLEXPORT VecLoad_Default_DA(Vec xin, PetscViewer viewer) 604 { 605 PetscErrorCode ierr; 606 DM da; 607 PetscBool isbinary; 608 #if defined(PETSC_HAVE_HDF5) 609 PetscBool ishdf5; 610 #endif 611 612 PetscFunctionBegin; 613 ierr = PetscObjectQuery((PetscObject)xin,"DA",(PetscObject*)&da);CHKERRQ(ierr); 614 if (!da) SETERRQ(((PetscObject)xin)->comm,PETSC_ERR_ARG_WRONG,"Vector not generated from a DA"); 615 616 #if defined(PETSC_HAVE_HDF5) 617 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);CHKERRQ(ierr); 618 #endif 619 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 620 621 if (isbinary) { 622 ierr = VecLoad_Binary_DA(xin,viewer);CHKERRQ(ierr); 623 #if defined(PETSC_HAVE_HDF5) 624 } else if (ishdf5) { 625 ierr = VecLoad_HDF5_DA(xin,viewer);CHKERRQ(ierr); 626 #endif 627 } else { 628 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for vector loading", ((PetscObject)viewer)->type_name); 629 } 630 PetscFunctionReturn(0); 631 } 632 EXTERN_C_END 633