Actual source code: ex46.c

petsc-3.12.2 2019-11-22
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 PetscReal test_values[] = { 0.311256, 88.068, 11.077444, 9953.62, 7.345, 64.8943, 3.1458, 6699.95, 0.00084, 0.0647 };

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

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

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

 31:   VecView(x,viewer);

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

 38:   PetscViewerDestroy(&viewer);
 39:   return(0);
 40: }

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

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

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

 59:   VecLoad(x,viewer);

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

 66:   PetscViewerDestroy(&viewer);
 67:   return(0);
 68: }

 70: PetscErrorCode VecFill(Vec x)
 71: {
 73:   PetscInt       i,s,e;

 76:   VecGetOwnershipRange(x,&s,&e);
 77:   for (i=s; i<e; i++) {
 78:     VecSetValue(x,i,(PetscScalar)test_values[i],INSERT_VALUES);
 79:   }
 80:   VecAssemblyBegin(x);
 81:   VecAssemblyEnd(x);
 82:   return(0);
 83: }

 85: PetscErrorCode VecCompare(Vec a,Vec b)
 86: {
 87:   PetscInt       locmin[2],locmax[2];
 88:   PetscReal      min[2],max[2];
 89:   Vec            ref;

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

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

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

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

106:   VecDuplicate(a,&ref);
107:   VecCopy(a,ref);
108:   VecAXPY(ref,-1.0,b);
109:   VecMin(ref,&locmin[0],&min[0]);
110:   if (PetscAbsReal(min[0]) > 1.0e-10) {
111:     PetscPrintf(PETSC_COMM_WORLD,"  ERROR: min(a-b) > 1.0e-10\n");
112:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) = %+1.10e\n",(double)PetscAbsReal(min[0]));
113:   } else {
114:     PetscPrintf(PETSC_COMM_WORLD,"  min(a-b) < 1.0e-10\n");
115:   }
116:   VecDestroy(&ref);
117:   return(0);
118: }

120: PetscErrorCode HeaderlessBinaryRead(const char name[])
121: {
123:   int            fdes;
124:   PetscScalar    buffer[VEC_LEN];
125:   PetscInt       i;
126:   PetscMPIInt    rank;
127:   PetscBool      dataverified = PETSC_TRUE;

130:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
131:   if (!rank) {
132:     PetscBinaryOpen(name,FILE_MODE_READ,&fdes);
133:     PetscBinaryRead(fdes,buffer,VEC_LEN,NULL,PETSC_SCALAR);
134:     PetscBinaryClose(fdes);

136:     for (i=0; i<VEC_LEN; i++) {
137:       PetscScalar v;
138:       v = PetscAbsScalar(test_values[i]-buffer[i]);
139: #if defined(PETSC_USE_COMPLEX)
140:       if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
141:         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);
142:         dataverified = PETSC_FALSE;
143:       }
144: #else
145:       if (PetscRealPart(v) > 1.0e-10) {
146:         PetscPrintf(PETSC_COMM_SELF,"ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %D])\n",(double)PetscRealPart(buffer[i]),i);
147:         dataverified = PETSC_FALSE;
148:       }
149: #endif
150:     }
151:     if (dataverified) {
152:       PetscPrintf(PETSC_COMM_SELF,"Headerless read of data verified\n");
153:     }
154:   }
155:   return(0);
156: }

158: PetscErrorCode TestBinary(void)
159: {
161:   Vec            x,y;
162:   PetscBool      skipheader = PETSC_TRUE;
163:   PetscBool      usempiio = PETSC_FALSE;
164: 
166:   VecCreate(PETSC_COMM_WORLD,&x);
167:   VecSetSizes(x,PETSC_DECIDE,VEC_LEN);
168:   VecSetFromOptions(x);
169:   VecFill(x);
170:   MyVecDump("xH.pbvec",skipheader,usempiio,x);

172:   VecCreate(PETSC_COMM_WORLD,&y);
173:   VecSetSizes(y,PETSC_DECIDE,VEC_LEN);
174:   VecSetFromOptions(y);

176:   MyVecLoad("xH.pbvec",skipheader,usempiio,y);
177:   VecCompare(x,y);

179:   VecDestroy(&y);
180:   VecDestroy(&x);

182:   HeaderlessBinaryRead("xH.pbvec");
183:   return(0);
184: }

186: #if defined(PETSC_HAVE_MPIIO)
187: PetscErrorCode TestBinaryMPIIO(void)
188: {
190:   Vec            x,y;
191:   PetscBool      skipheader = PETSC_TRUE;
192:   PetscBool      usempiio = PETSC_TRUE;

195:   VecCreate(PETSC_COMM_WORLD,&x);
196:   VecSetSizes(x,PETSC_DECIDE,VEC_LEN);
197:   VecSetFromOptions(x);
198:   VecFill(x);
199:   MyVecDump("xHmpi.pbvec",skipheader,usempiio,x);

201:   VecCreate(PETSC_COMM_WORLD,&y);
202:   VecSetSizes(y,PETSC_DECIDE,VEC_LEN);
203:   VecSetFromOptions(y);

205:   MyVecLoad("xHmpi.pbvec",skipheader,usempiio,y);
206:   VecCompare(x,y);

208:   VecDestroy(&y);
209:   VecDestroy(&x);

211:   HeaderlessBinaryRead("xHmpi.pbvec");
212:   return(0);
213: }
214: #endif

216: int main(int argc,char **args)
217: {
219:   PetscBool      usempiio = PETSC_FALSE;

221:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
222:   PetscOptionsGetBool(NULL,NULL,"-usempiio",&usempiio,NULL);
223:   if (!usempiio) {
224:     TestBinary();
225:   } else {
226: #if defined(PETSC_HAVE_MPIIO)
227:     TestBinaryMPIIO();
228: #else
229:     PetscPrintf(PETSC_COMM_WORLD,"Warning: Executing TestBinaryMPIIO() requires a working MPI-2 implementation\n");
230: #endif
231:   }
232:   PetscFinalize();
233:   return ierr;
234: }


237: /*TEST

239:    test:
240:       output_file: output/ex46_1_p1.out

242:    test:
243:       suffix: 2
244:       nsize: 6
245:       output_file: output/ex46_1_p6.out

247:    test:
248:       suffix: 3
249:       nsize: 12
250:       output_file: output/ex46_1_p12.out

252:    test:
253:       suffix: mpiio
254:       nsize: 6
255:       args: -usempiio
256:       output_file: output/ex46_2_p6.out
257:       requires: mpiio

259: TEST*/