Actual source code: ex3.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  2: static char help[] = "Tests PetscViewerBinary VecView()/VecLoad() function correctly when binary header is skipped.\n\n";

  4: /*T
  5:  Concepts: viewers^skipheader
  6: T*/

  8: #include <petscviewer.h>
  9: #include <petscvec.h>

 11: #define VEC_LEN 10
 12: const PetscScalar test_values[] = { 0.311256, 88.068, 11.077444, 9953.62, 7.345, 64.8943, 3.1458, 6699.95, 0.00084, 0.0647 };

 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) { PetscViewerBinarySetMPIIO(viewer); }
 31:   PetscViewerFileSetName(viewer,fname);

 33:   VecView(x,viewer);

 35:   PetscViewerBinaryGetMPIIO(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: }

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

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

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

 63:   VecLoad(x,viewer);

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

 70:   PetscViewerDestroy(&viewer);
 71:   return(0);
 72: }

 76: PetscErrorCode VecFill(Vec x)
 77: {
 79:   PetscInt       i,s,e;

 82:   VecGetOwnershipRange(x,&s,&e);
 83:   for (i=s; i<e; i++) {
 84:     VecSetValue(x,i,test_values[i],INSERT_VALUES);
 85:   }
 86:   VecAssemblyBegin(x);
 87:   VecAssemblyEnd(x);
 88:   return(0);
 89: }

 93: PetscErrorCode VecCompare(Vec a,Vec b)
 94: {
 95:   PetscInt       locmin[2],locmax[2];
 96:   PetscReal      min[2],max[2];
 97:   Vec            ref;

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

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

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

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

114:   VecDuplicate(a,&ref);
115:   VecCopy(a,ref);
116:   VecAXPY(ref,-1.0,b);
117:   VecMin(ref,&locmin[0],&min[0]);
118:   if (PetscAbsReal(min[0]) > 1.0e-10) {
119:     PetscPrintf(PETSC_COMM_WORLD,"  ERROR: min(a-b) > 1.0e-10\n");
120:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
121:   } else {
122:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) < 1.0e-10\n");
123:   }
124:   VecDestroy(&ref);
125:   return(0);
126: }

130: PetscErrorCode HeaderlessBinaryRead(const char name[])
131: {
133:   int            fdes;
134:   PetscScalar    buffer[VEC_LEN];
135:   PetscInt       i;
136:   PetscMPIInt    rank;
137:   PetscBool      dataverified = PETSC_TRUE;

140:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
141:   if (!rank) {
142:     PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
143:     PetscBinaryRead(fdes,buffer,VEC_LEN,PETSC_SCALAR);
144:     PetscBinaryClose(fdes);

146:     for (i=0; i<VEC_LEN; i++) {
147:       PetscScalar v;
148:       v = PetscAbsScalar(test_values[i]-buffer[i]);
149: #if defined(PETSC_USE_COMPLEX)
150:       if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
151:         PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = (%+1.12e,%+1.12e) [loc %D])\n",(double)PetscRealPart(buffer[i]),(double)PetscImaginaryPart(buffer[i]),i);
152:         dataverified = PETSC_FALSE;
153:       }
154: #else
155:       if (PetscRealPart(v) > 1.0e-10) {
156:         PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %D])\n",(double)PetscRealPart(buffer[i]),i);
157:         dataverified = PETSC_FALSE;
158:       }
159: #endif
160:     }
161:     if (dataverified) {
162:       PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified\n");
163:     }
164:   }
165:   return(0);
166: }

170: PetscErrorCode TestBinary(void)
171: {
173:   Vec            x,y;
174:   PetscBool      skipheader = PETSC_TRUE;
175:   PetscBool      usempiio = PETSC_FALSE;
176: 
178:   VecCreate(PETSC_COMM_WORLD,&x);
179:   VecSetSizes(x,PETSC_DECIDE,VEC_LEN);
180:   VecSetFromOptions(x);
181:   VecFill(x);
182:   MyVecDump("xH.pbvec",skipheader,usempiio,x);

184:   VecCreate(PETSC_COMM_WORLD,&y);
185:   VecSetSizes(y,PETSC_DECIDE,VEC_LEN);
186:   VecSetFromOptions(y);

188:   MyVecLoad("xH.pbvec",skipheader,usempiio,y);
189:   VecCompare(x,y);

191:   VecDestroy(&y);
192:   VecDestroy(&x);

194:   HeaderlessBinaryRead("xH.pbvec");
195:   return(0);
196: }

198: #if defined(PETSC_HAVE_MPIIO)
201: PetscErrorCode TestBinaryMPIIO(void)
202: {
204:   Vec            x,y;
205:   PetscBool      skipheader = PETSC_TRUE;
206:   PetscBool      usempiio = PETSC_TRUE;

209:   VecCreate(PETSC_COMM_WORLD,&x);
210:   VecSetSizes(x,PETSC_DECIDE,VEC_LEN);
211:   VecSetFromOptions(x);
212:   VecFill(x);
213:   MyVecDump("xHmpi.pbvec",skipheader,usempiio,x);

215:   VecCreate(PETSC_COMM_WORLD,&y);
216:   VecSetSizes(y,PETSC_DECIDE,VEC_LEN);
217:   VecSetFromOptions(y);

219:   MyVecLoad("xHmpi.pbvec",skipheader,usempiio,y);
220:   VecCompare(x,y);

222:   VecDestroy(&y);
223:   VecDestroy(&x);

225:   HeaderlessBinaryRead("xHmpi.pbvec");
226:   return(0);
227: }
228: #endif

232: int main(int argc,char **args)
233: {
235:   PetscBool      usempiio = PETSC_FALSE;
236: 
237:   PetscInitialize(&argc,&args,(char*)0,help);

239:   PetscOptionsGetBool(NULL,"-usempiio",&usempiio,NULL);
240:   if (!usempiio) {
241:     TestBinary();
242:   } else {
243: #if defined(PETSC_HAVE_MPIIO)
244:     TestBinaryMPIIO();
245: #else
246:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestBinaryMPIIO() requires a working MPI-2 implementation\n");
247: #endif
248:   }
249:   PetscFinalize();
250:   return 0;
251: }