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;

 16:   PetscFunctionBeginUser;
 17:   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));

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

 26:   PetscCall(VecView(x, viewer));

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

 33:   PetscCall(PetscViewerDestroy(&viewer));
 34:   PetscFunctionReturn(PETSC_SUCCESS);
 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;

 43:   PetscFunctionBeginUser;
 44:   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));

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

 53:   PetscCall(VecLoad(x, viewer));

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

 60:   PetscCall(PetscViewerDestroy(&viewer));
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

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

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

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

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

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

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

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

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

118:   PetscFunctionBeginUser;
119:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
120:   if (rank == 0) {
121:     PetscCall(PetscBinaryOpen(name, FILE_MODE_READ, &fdes));
122:     PetscCall(PetscBinaryRead(fdes, buffer, VEC_LEN, NULL, PETSC_SCALAR));
123:     PetscCall(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:         PetscCall(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:         PetscCall(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) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Headerless read of data verified\n"));
141:   }
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }

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

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

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

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

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

168:   PetscCall(HeaderlessBinaryRead("xH.pbvec"));
169:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

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

196:   PetscCall(HeaderlessBinaryRead("xHmpi.pbvec"));
197:   PetscFunctionReturn(PETSC_SUCCESS);
198: }
199: #endif

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

205:   PetscFunctionBeginUser;
206:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
207:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-usempiio", &usempiio, NULL));
208:   if (!usempiio) {
209:     PetscCall(TestBinary());
210:   } else {
211: #if defined(PETSC_HAVE_MPIIO)
212:     PetscCall(TestBinaryMPIIO());
213: #else
214:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Executing TestBinaryMPIIO() requires a working MPI-2 implementation\n"));
215: #endif
216:   }
217:   PetscCall(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*/