Actual source code: ex46.c


  2: static char help[] = "Tests PetscViewerBinary VecView()/VecLoad() function correctly when binary header is skipped.\n\n";

  4: #include <petscviewer.h>
  5: #include <petscvec.h>

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

 10: PetscErrorCode MyVecDump(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x)
 11: {
 12:   MPI_Comm    comm;
 13:   PetscViewer viewer;
 14:   PetscBool   ismpiio, isskip;

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

 19:   PetscViewerCreate(comm, &viewer);
 20:   PetscViewerSetType(viewer, PETSCVIEWERBINARY);
 21:   if (skippheader) PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE);
 22:   PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
 23:   if (usempiio) PetscViewerBinarySetUseMPIIO(viewer, PETSC_TRUE);
 24:   PetscViewerFileSetName(viewer, fname);

 26:   VecView(x, viewer);

 28:   PetscViewerBinaryGetUseMPIIO(viewer, &ismpiio);
 29:   if (ismpiio) PetscPrintf(comm, "*** PetscViewer[write] using MPI-IO ***\n");
 30:   PetscViewerBinaryGetSkipHeader(viewer, &isskip);
 31:   if (isskip) PetscPrintf(comm, "*** PetscViewer[write] skipping header ***\n");

 33:   PetscViewerDestroy(&viewer);
 34:   return 0;
 35: }

 37: PetscErrorCode MyVecLoad(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x)
 38: {
 39:   MPI_Comm    comm;
 40:   PetscViewer viewer;
 41:   PetscBool   ismpiio, isskip;

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

 46:   PetscViewerCreate(comm, &viewer);
 47:   PetscViewerSetType(viewer, PETSCVIEWERBINARY);
 48:   if (skippheader) PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE);
 49:   PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 50:   if (usempiio) PetscViewerBinarySetUseMPIIO(viewer, PETSC_TRUE);
 51:   PetscViewerFileSetName(viewer, fname);

 53:   VecLoad(x, viewer);

 55:   PetscViewerBinaryGetSkipHeader(viewer, &isskip);
 56:   if (isskip) PetscPrintf(comm, "*** PetscViewer[load] skipping header ***\n");
 57:   PetscViewerBinaryGetUseMPIIO(viewer, &ismpiio);
 58:   if (ismpiio) PetscPrintf(comm, "*** PetscViewer[load] using MPI-IO ***\n");

 60:   PetscViewerDestroy(&viewer);
 61:   return 0;
 62: }

 64: PetscErrorCode VecFill(Vec x)
 65: {
 66:   PetscInt i, s, e;

 69:   VecGetOwnershipRange(x, &s, &e);
 70:   for (i = s; i < e; i++) VecSetValue(x, i, (PetscScalar)test_values[i], INSERT_VALUES);
 71:   VecAssemblyBegin(x);
 72:   VecAssemblyEnd(x);
 73:   return 0;
 74: }

 76: PetscErrorCode VecCompare(Vec a, Vec b)
 77: {
 78:   PetscInt  locmin[2], locmax[2];
 79:   PetscReal min[2], max[2];
 80:   Vec       ref;

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

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

 89:   PetscPrintf(PETSC_COMM_WORLD, "VecCompare\n");
 90:   PetscPrintf(PETSC_COMM_WORLD, "  min(a)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[0], locmin[0]);
 91:   PetscPrintf(PETSC_COMM_WORLD, "  max(a)   = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[0], locmax[0]);

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

 96:   VecDuplicate(a, &ref);
 97:   VecCopy(a, ref);
 98:   VecAXPY(ref, -1.0, b);
 99:   VecMin(ref, &locmin[0], &min[0]);
100:   if (PetscAbsReal(min[0]) > 1.0e-10) {
101:     PetscPrintf(PETSC_COMM_WORLD, "  ERROR: min(a-b) > 1.0e-10\n");
102:     PetscPrintf(PETSC_COMM_WORLD, "  min(a-b) = %+1.10e\n", (double)PetscAbsReal(min[0]));
103:   } else {
104:     PetscPrintf(PETSC_COMM_WORLD, "  min(a-b) < 1.0e-10\n");
105:   }
106:   VecDestroy(&ref);
107:   return 0;
108: }

110: PetscErrorCode HeaderlessBinaryRead(const char name[])
111: {
112:   int         fdes;
113:   PetscScalar buffer[VEC_LEN];
114:   PetscInt    i;
115:   PetscMPIInt rank;
116:   PetscBool   dataverified = PETSC_TRUE;

119:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
120:   if (rank == 0) {
121:     PetscBinaryOpen(name, FILE_MODE_READ, &fdes);
122:     PetscBinaryRead(fdes, buffer, VEC_LEN, NULL, PETSC_SCALAR);
123:     PetscBinaryClose(fdes);

125:     for (i = 0; i < VEC_LEN; i++) {
126:       PetscScalar v;
127:       v = PetscAbsScalar(test_values[i] - buffer[i]);
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 %" PetscInt_FMT "])\n", (double)PetscRealPart(buffer[i]), (double)PetscImaginaryPart(buffer[i]), i);
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 %" PetscInt_FMT "])\n", (double)PetscRealPart(buffer[i]), i);
136:         dataverified = PETSC_FALSE;
137:       }
138: #endif
139:     }
140:     if (dataverified) PetscPrintf(PETSC_COMM_SELF, "Headerless read of data verified\n");
141:   }
142:   return 0;
143: }

145: PetscErrorCode TestBinary(void)
146: {
147:   Vec       x, y;
148:   PetscBool skipheader = PETSC_TRUE;
149:   PetscBool usempiio   = PETSC_FALSE;

152:   VecCreate(PETSC_COMM_WORLD, &x);
153:   VecSetSizes(x, PETSC_DECIDE, VEC_LEN);
154:   VecSetFromOptions(x);
155:   VecFill(x);
156:   MyVecDump("xH.pbvec", skipheader, usempiio, x);

158:   VecCreate(PETSC_COMM_WORLD, &y);
159:   VecSetSizes(y, PETSC_DECIDE, VEC_LEN);
160:   VecSetFromOptions(y);

162:   MyVecLoad("xH.pbvec", skipheader, usempiio, y);
163:   VecCompare(x, y);

165:   VecDestroy(&y);
166:   VecDestroy(&x);

168:   HeaderlessBinaryRead("xH.pbvec");
169:   return 0;
170: }

172: #if defined(PETSC_HAVE_MPIIO)
173: PetscErrorCode TestBinaryMPIIO(void)
174: {
175:   Vec       x, y;
176:   PetscBool skipheader = PETSC_TRUE;
177:   PetscBool usempiio   = PETSC_TRUE;

180:   VecCreate(PETSC_COMM_WORLD, &x);
181:   VecSetSizes(x, PETSC_DECIDE, VEC_LEN);
182:   VecSetFromOptions(x);
183:   VecFill(x);
184:   MyVecDump("xHmpi.pbvec", skipheader, usempiio, x);

186:   VecCreate(PETSC_COMM_WORLD, &y);
187:   VecSetSizes(y, PETSC_DECIDE, VEC_LEN);
188:   VecSetFromOptions(y);

190:   MyVecLoad("xHmpi.pbvec", skipheader, usempiio, y);
191:   VecCompare(x, y);

193:   VecDestroy(&y);
194:   VecDestroy(&x);

196:   HeaderlessBinaryRead("xHmpi.pbvec");
197:   return 0;
198: }
199: #endif

201: int main(int argc, char **args)
202: {
203:   PetscBool usempiio = PETSC_FALSE;

206:   PetscInitialize(&argc, &args, (char *)0, help);
207:   PetscOptionsGetBool(NULL, NULL, "-usempiio", &usempiio, NULL);
208:   if (!usempiio) {
209:     TestBinary();
210:   } else {
211: #if defined(PETSC_HAVE_MPIIO)
212:     TestBinaryMPIIO();
213: #else
214:     PetscPrintf(PETSC_COMM_WORLD, "Warning: Executing TestBinaryMPIIO() requires a working MPI-2 implementation\n");
215: #endif
216:   }
217:   PetscFinalize();
218:   return 0;
219: }

221: /*TEST

223:    test:
224:       output_file: output/ex46_1_p1.out

226:    test:
227:       suffix: 2
228:       nsize: 6
229:       output_file: output/ex46_1_p6.out

231:    test:
232:       suffix: 3
233:       nsize: 12
234:       output_file: output/ex46_1_p12.out

236:    testset:
237:       requires: mpiio
238:       args: -usempiio
239:       test:
240:          suffix: mpiio_1
241:          output_file: output/ex46_2_p1.out
242:       test:
243:          suffix: mpiio_2
244:          nsize: 6
245:          output_file: output/ex46_2_p6.out
246:       test:
247:          suffix: mpiio_3
248:          nsize: 12
249:          output_file: output/ex46_2_p12.out

251: TEST*/