Actual source code: ex15.c

petsc-3.8.4 2018-03-24
Report Typos and Errors

  2: static char help[] = "Tests VecView() functionality with DMDA objects when using:"\
  3: "(i) a PetscViewer binary with MPI-IO support; and (ii) when the binary header is skipped.\n\n";

  5:  #include <petscdm.h>
  6:  #include <petscdmda.h>

  8: #define DMDA_I 5
  9: #define DMDA_J 4
 10: #define DMDA_K 6

 12: const PetscReal dmda_i_val[] = { 1.10, 2.3006, 2.32444, 3.44006, 66.9009 };
 13: const PetscReal dmda_j_val[] = { 0.0, 0.25, 0.5, 0.75 };
 14: const PetscReal dmda_k_val[] = { 0.0, 1.1, 2.2, 3.3, 4.4, 5.5 };

 16: PetscErrorCode MyVecDump(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
 17: {
 18:   MPI_Comm       comm;
 19:   PetscViewer    viewer;
 20:   PetscBool      ismpiio,isskip;

 24:   PetscObjectGetComm((PetscObject)x,&comm);

 26:   PetscViewerCreate(comm,&viewer);
 27:   PetscViewerSetType(viewer,PETSCVIEWERBINARY);
 28:   if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
 29:   PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
 30:   if (usempiio) { PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE); }
 31:   PetscViewerFileSetName(viewer,fname);

 33:   VecView(x,viewer);

 35:   PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
 36:   if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[write] using MPI-IO ***\n"); }
 37:   PetscViewerBinaryGetSkipHeader(viewer,&isskip);
 38:   if (isskip) { PetscPrintf(comm,"*** PetscViewer[write] skipping header ***\n"); }

 40:   PetscViewerDestroy(&viewer);
 41:   return(0);
 42: }

 44: PetscErrorCode MyVecLoad(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)
 45: {
 46:   MPI_Comm       comm;
 47:   PetscViewer    viewer;
 48:   PetscBool      ismpiio,isskip;

 52:   PetscObjectGetComm((PetscObject)x,&comm);

 54:   PetscViewerCreate(comm,&viewer);
 55:   PetscViewerSetType(viewer,PETSCVIEWERBINARY);
 56:   if (skippheader) { PetscViewerBinarySetSkipHeader(viewer,PETSC_TRUE); }
 57:   PetscViewerFileSetMode(viewer,FILE_MODE_READ);
 58:   if (usempiio) { PetscViewerBinarySetUseMPIIO(viewer,PETSC_TRUE); }
 59:   PetscViewerFileSetName(viewer,fname);

 61:   VecLoad(x,viewer);

 63:   PetscViewerBinaryGetSkipHeader(viewer,&isskip);
 64:   if (isskip) { PetscPrintf(comm,"*** PetscViewer[load] skipping header ***\n"); }
 65:   PetscViewerBinaryGetUseMPIIO(viewer,&ismpiio);
 66:   if (ismpiio) { PetscPrintf(comm,"*** PetscViewer[load] using MPI-IO ***\n"); }

 68:   PetscViewerDestroy(&viewer);
 69:   return(0);
 70: }

 72: PetscErrorCode DMDAVecGenerateEntries(DM dm,Vec a)
 73: {
 74:   PetscScalar    ****LA_v;
 75:   PetscInt       i,j,k,si,sj,sk,ni,nj,nk,M,N;

 79:   DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
 80:   DMDAGetCorners(dm,&si,&sj,&sk,&ni,&nj,&nk);
 81:   DMDAVecGetArrayDOF(dm,a,&LA_v);
 82:   for (k=sk; k<sk+nk; k++) {
 83:     for (j=sj; j<sj+nj; j++) {
 84:       for (i=si; i<si+ni; i++) {
 85:         PetscScalar test_value_s;
 86: 
 87:         test_value_s = dmda_i_val[i]*((PetscScalar)i) + dmda_j_val[j]*((PetscScalar)(i+j*M)) + dmda_k_val[k]*((PetscScalar)(i + j*M + k*M*N));
 88:         LA_v[k][j][i][0] = 3.0 * test_value_s;
 89:         LA_v[k][j][i][1] = 3.0 * test_value_s + 1.0;
 90:         LA_v[k][j][i][2] = 3.0 * test_value_s + 2.0;
 91:       }
 92:     }
 93:   }
 94:   DMDAVecRestoreArrayDOF(dm,a,&LA_v);
 95:   return(0);
 96: }

 98: PetscErrorCode HeaderlessBinaryReadCheck(DM dm,const char name[])
 99: {
101:   int            fdes;
102:   PetscScalar    buffer[DMDA_I*DMDA_J*DMDA_K*3];
103:   PetscInt       len,d,i,j,k,M,N;
104:   PetscMPIInt    rank;
105:   PetscBool      dataverified = PETSC_TRUE;

108:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
109:   DMDAGetInfo(dm,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);
110:   len = DMDA_I*DMDA_J*DMDA_K*3;
111:   if (!rank) {
112:     PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
113:     PetscBinaryRead(fdes,buffer,len,PETSC_SCALAR);
114:     PetscBinaryClose(fdes);

116:     for (k=0; k<DMDA_K; k++) {
117:       for (j=0; j<DMDA_J; j++) {
118:         for (i=0; i<DMDA_I; i++) {
119:           for (d=0; d<3; d++) {
120:             PetscScalar v,test_value_s,test_value;
121:             PetscInt    index;

123:             test_value_s = dmda_i_val[i]*((PetscScalar)i) + dmda_j_val[j]*((PetscScalar)(i+j*M)) + dmda_k_val[k]*((PetscScalar)(i + j*M + k*M*N));
124:             test_value = 3.0 * test_value_s + (PetscScalar)d;

126:             index = 3*(i + j*M + k*M*N) + d;
127:             v = PetscAbsScalar(test_value-buffer[index]);
128: #if defined(PETSC_USE_COMPLEX)
129:             if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
130:               PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = (%+1.12e,%+1.12e) [loc %D,%D,%D(%D)])\n",(double)PetscRealPart(test_value),(double)PetscImaginaryPart(test_value),i,j,k,d);
131:               dataverified = PETSC_FALSE;
132:             }
133: #else
134:             if (PetscRealPart(v) > 1.0e-10) {
135:               PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %D,%D,%D(%D)])\n",(double)PetscRealPart(test_value),i,j,k,d);
136:               dataverified = PETSC_FALSE;
137:             }
138: #endif
139:           }
140:         }
141:       }
142:     }
143:     if (dataverified) {
144:       PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified for: %s\n",name);
145:     }
146:   }
147:   return(0);
148: }

150: PetscErrorCode VecCompare(Vec a,Vec b)
151: {
152:   PetscInt       locmin[2],locmax[2];
153:   PetscReal      min[2],max[2];
154:   Vec            ref;

158:   VecMin(a,&locmin[0],&min[0]);
159:   VecMax(a,&locmax[0],&max[0]);

161:   VecMin(b,&locmin[1],&min[1]);
162:   VecMax(b,&locmax[1],&max[1]);

164:   PetscPrintf(PETSC_COMM_WORLD,"VecCompare\n");
165:   PetscPrintf(PETSC_COMM_WORLD,"  min(a)   = %+1.2e [loc %D]\n",(double)min[0],locmin[0]);
166:   PetscPrintf(PETSC_COMM_WORLD,"  max(a)   = %+1.2e [loc %D]\n",(double)max[0],locmax[0]);

168:   PetscPrintf(PETSC_COMM_WORLD,"  min(b)   = %+1.2e [loc %D]\n",(double)min[1],locmin[1]);
169:   PetscPrintf(PETSC_COMM_WORLD,"  max(b)   = %+1.2e [loc %D]\n",(double)max[1],locmax[1]);

171:   VecDuplicate(a,&ref);
172:   VecCopy(a,ref);
173:   VecAXPY(ref,-1.0,b);
174:   VecMin(ref,&locmin[0],&min[0]);
175:   if (PetscAbsReal(min[0]) > 1.0e-10) {
176:     PetscPrintf(PETSC_COMM_WORLD,"  ERROR: min(a-b) > 1.0e-10\n");
177:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
178:   } else {
179:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) < 1.0e-10\n");
180:   }
181:   VecDestroy(&ref);
182:   return(0);
183: }

185: PetscErrorCode TestDMDAVec(PetscBool usempiio)
186: {
187:   DM             dm;
188:   Vec            x_ref,x_test;
189:   PetscBool      skipheader = PETSC_TRUE;

193:   if (!usempiio) { PetscPrintf(PETSC_COMM_WORLD,"%s\n",PETSC_FUNCTION_NAME); }
194:   else { PetscPrintf(PETSC_COMM_WORLD,"%s [using mpi-io]\n",PETSC_FUNCTION_NAME); }
195:   DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,DMDA_I,DMDA_J,DMDA_K,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
196:                         3,2,NULL,NULL,NULL,&dm);
197:   DMSetFromOptions(dm);
198:   DMSetUp(dm);

200:   DMCreateGlobalVector(dm,&x_ref);
201:   DMDAVecGenerateEntries(dm,x_ref);

203:   if (!usempiio) {
204:     MyVecDump("dmda.pbvec",skipheader,PETSC_FALSE,x_ref);
205:   } else {
206:     MyVecDump("dmda-mpiio.pbvec",skipheader,PETSC_TRUE,x_ref);
207:   }

209:   DMCreateGlobalVector(dm,&x_test);

211:   if (!usempiio) {
212:     MyVecLoad("dmda.pbvec",skipheader,usempiio,x_test);
213:   } else {
214:     MyVecLoad("dmda-mpiio.pbvec",skipheader,usempiio,x_test);
215:   }

217:   VecCompare(x_ref,x_test);

219:   if (!usempiio) {
220:     HeaderlessBinaryReadCheck(dm,"dmda.pbvec");
221:   } else {
222:     HeaderlessBinaryReadCheck(dm,"dmda-mpiio.pbvec");
223:   }
224:   VecDestroy(&x_ref);
225:   VecDestroy(&x_test);
226:   DMDestroy(&dm);
227:   return(0);
228: }

230: int main(int argc,char **args)
231: {
233:   PetscBool      usempiio = PETSC_FALSE;

235:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
236:   PetscOptionsGetBool(NULL,NULL,"-usempiio",&usempiio,NULL);
237:   if (!usempiio) {
238:     TestDMDAVec(PETSC_FALSE);
239:   } else {
240: #if defined(PETSC_HAVE_MPIIO)
241:     TestDMDAVec(PETSC_TRUE);
242: #else
243:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n");
244: #endif
245:   }
246:   PetscFinalize();
247:   return ierr;
248: }

250: /*TEST

252:    test:
253:       suffix: 3
254:       nsize: 12
255:       requires: define(PETSC_HAVE_MPIIO) !define(PETSC_HAVE_LIBMSMPI)
256:       args: -usempiio

258: TEST*/