Actual source code: ex46.c

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

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

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

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

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

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

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

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

 32:   PetscCall(PetscViewerDestroy(&viewer));
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 34: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

124:     for (i = 0; i < VEC_LEN; i++) {
125:       PetscScalar v;
126:       v = PetscAbsScalar(test_values[i] - buffer[i]);
127: #if defined(PETSC_USE_COMPLEX)
128:       if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
129:         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));
130:         dataverified = PETSC_FALSE;
131:       }
132: #else
133:       if (PetscRealPart(v) > 1.0e-10) {
134:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %" PetscInt_FMT "])\n", (double)PetscRealPart(buffer[i]), i));
135:         dataverified = PETSC_FALSE;
136:       }
137: #endif
138:     }
139:     if (dataverified) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Headerless read of data verified\n"));
140:   }
141:   PetscFunctionReturn(PETSC_SUCCESS);
142: }

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

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

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

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

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

167:   PetscCall(HeaderlessBinaryRead("xH.pbvec"));
168:   PetscFunctionReturn(PETSC_SUCCESS);
169: }

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

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

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

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

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

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

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

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

220: /*TEST

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

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

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

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

250: TEST*/