Actual source code: ex15.c


  2: static char help[] = "Tests VecView() functionality with DMDA objects when using:"
  3:                      "(i) a PetscViewer binary with MPI-IO support; and (ii) when the binary header is skipped.\n\n";

  5: #include <petscdm.h>
  6: #include <petscdmda.h>

  8: #define DMDA_I 5
  9: #define DMDA_J 4
 10: #define DMDA_K 6

 12: const PetscReal dmda_i_val[] = {1.10, 2.3006, 2.32444, 3.44006, 66.9009};
 13: const PetscReal dmda_j_val[] = {0.0, 0.25, 0.5, 0.75};
 14: const PetscReal dmda_k_val[] = {0.0, 1.1, 2.2, 3.3, 4.4, 5.5};

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

 22:   PetscFunctionBeginUser;
 23:   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));

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

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

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

 39:   PetscCall(PetscViewerDestroy(&viewer));
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

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

 49:   PetscFunctionBeginUser;
 50:   PetscCall(PetscObjectGetComm((PetscObject)x, &comm));

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

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

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

 66:   PetscCall(PetscViewerDestroy(&viewer));
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

 70: PetscErrorCode DMDAVecGenerateEntries(DM dm, Vec a)
 71: {
 72:   PetscScalar ****LA_v;
 73:   PetscInt        i, j, k, l, si, sj, sk, ni, nj, nk, M, N, dof;

 75:   PetscFunctionBeginUser;
 76:   PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
 77:   PetscCall(DMDAGetCorners(dm, &si, &sj, &sk, &ni, &nj, &nk));
 78:   PetscCall(DMDAVecGetArrayDOF(dm, a, &LA_v));
 79:   for (k = sk; k < sk + nk; k++) {
 80:     for (j = sj; j < sj + nj; j++) {
 81:       for (i = si; i < si + ni; i++) {
 82:         PetscScalar test_value_s;

 84:         test_value_s = dmda_i_val[i] * ((PetscScalar)i) + dmda_j_val[j] * ((PetscScalar)(i + j * M)) + dmda_k_val[k] * ((PetscScalar)(i + j * M + k * M * N));
 85:         for (l = 0; l < dof; l++) LA_v[k][j][i][l] = (PetscScalar)dof * test_value_s + (PetscScalar)l;
 86:       }
 87:     }
 88:   }
 89:   PetscCall(DMDAVecRestoreArrayDOF(dm, a, &LA_v));
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

 93: PetscErrorCode HeaderlessBinaryReadCheck(DM dm, const char name[])
 94: {
 95:   int         fdes;
 96:   PetscScalar buffer[DMDA_I * DMDA_J * DMDA_K * 10];
 97:   PetscInt    len, d, i, j, k, M, N, dof;
 98:   PetscMPIInt rank;
 99:   PetscBool   dataverified = PETSC_TRUE;

101:   PetscFunctionBeginUser;
102:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
103:   PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
104:   len = DMDA_I * DMDA_J * DMDA_K * dof;
105:   if (rank == 0) {
106:     PetscCall(PetscBinaryOpen(name, FILE_MODE_READ, &fdes));
107:     PetscCall(PetscBinaryRead(fdes, buffer, len, NULL, PETSC_SCALAR));
108:     PetscCall(PetscBinaryClose(fdes));

110:     for (k = 0; k < DMDA_K; k++) {
111:       for (j = 0; j < DMDA_J; j++) {
112:         for (i = 0; i < DMDA_I; i++) {
113:           for (d = 0; d < dof; d++) {
114:             PetscScalar v, test_value_s, test_value;
115:             PetscInt    index;

117:             test_value_s = dmda_i_val[i] * ((PetscScalar)i) + dmda_j_val[j] * ((PetscScalar)(i + j * M)) + dmda_k_val[k] * ((PetscScalar)(i + j * M + k * M * N));
118:             test_value   = (PetscScalar)dof * test_value_s + (PetscScalar)d;

120:             index = dof * (i + j * M + k * M * N) + d;
121:             v     = PetscAbsScalar(test_value - buffer[index]);
122: #if defined(PETSC_USE_COMPLEX)
123:             if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
124:               PetscCall(PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference > 1.0e-10 occurred (delta = (%+1.12e,%+1.12e) [loc %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "(%" PetscInt_FMT ")])\n", (double)PetscRealPart(test_value), (double)PetscImaginaryPart(test_value), i, j, k, d));
125:               dataverified = PETSC_FALSE;
126:             }
127: #else
128:             if (PetscRealPart(v) > 1.0e-10) {
129:               PetscCall(PetscPrintf(PETSC_COMM_SELF, "ERROR: Difference > 1.0e-10 occurred (delta = %+1.12e [loc %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "(%" PetscInt_FMT ")])\n", (double)PetscRealPart(test_value), i, j, k, d));
130:               dataverified = PETSC_FALSE;
131:             }
132: #endif
133:           }
134:         }
135:       }
136:     }
137:     if (dataverified) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Headerless read of data verified for: %s\n", name));
138:   }
139:   PetscFunctionReturn(PETSC_SUCCESS);
140: }

142: PetscErrorCode VecCompare(Vec a, Vec b)
143: {
144:   PetscInt  locmin[2], locmax[2];
145:   PetscReal min[2], max[2];
146:   Vec       ref;

148:   PetscFunctionBeginUser;
149:   PetscCall(VecMin(a, &locmin[0], &min[0]));
150:   PetscCall(VecMax(a, &locmax[0], &max[0]));

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

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

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

162:   PetscCall(VecDuplicate(a, &ref));
163:   PetscCall(VecCopy(a, ref));
164:   PetscCall(VecAXPY(ref, -1.0, b));
165:   PetscCall(VecMin(ref, &locmin[0], &min[0]));
166:   if (PetscAbsReal(min[0]) > 1.0e-10) {
167:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  ERROR: min(a-b) > 1.0e-10\n"));
168:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(a-b) = %+1.10e\n", (double)PetscAbsReal(min[0])));
169:   } else {
170:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "  min(a-b) < 1.0e-10\n"));
171:   }
172:   PetscCall(VecDestroy(&ref));
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

176: PetscErrorCode TestDMDAVec(PetscBool usempiio)
177: {
178:   DM        dm;
179:   Vec       x_ref, x_test;
180:   PetscBool skipheader = PETSC_TRUE;

182:   PetscFunctionBeginUser;
183:   if (!usempiio) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s\n", PETSC_FUNCTION_NAME));
184:   else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s [using mpi-io]\n", PETSC_FUNCTION_NAME));
185:   PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, DMDA_I, DMDA_J, DMDA_K, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 3, 2, NULL, NULL, NULL, &dm));
186:   PetscCall(DMSetFromOptions(dm));
187:   PetscCall(DMSetUp(dm));

189:   PetscCall(DMCreateGlobalVector(dm, &x_ref));
190:   PetscCall(DMDAVecGenerateEntries(dm, x_ref));

192:   if (!usempiio) {
193:     PetscCall(MyVecDump("dmda.pbvec", skipheader, PETSC_FALSE, x_ref));
194:   } else {
195:     PetscCall(MyVecDump("dmda-mpiio.pbvec", skipheader, PETSC_TRUE, x_ref));
196:   }

198:   PetscCall(DMCreateGlobalVector(dm, &x_test));

200:   if (!usempiio) {
201:     PetscCall(MyVecLoad("dmda.pbvec", skipheader, usempiio, x_test));
202:   } else {
203:     PetscCall(MyVecLoad("dmda-mpiio.pbvec", skipheader, usempiio, x_test));
204:   }

206:   PetscCall(VecCompare(x_ref, x_test));

208:   if (!usempiio) {
209:     PetscCall(HeaderlessBinaryReadCheck(dm, "dmda.pbvec"));
210:   } else {
211:     PetscCall(HeaderlessBinaryReadCheck(dm, "dmda-mpiio.pbvec"));
212:   }
213:   PetscCall(VecDestroy(&x_ref));
214:   PetscCall(VecDestroy(&x_test));
215:   PetscCall(DMDestroy(&dm));
216:   PetscFunctionReturn(PETSC_SUCCESS);
217: }

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

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

239: /*TEST

241:    test:

243:    test:
244:       suffix: 2
245:       nsize: 12

247:    test:
248:       suffix: 3
249:       nsize: 12
250:       requires: defined(PETSC_HAVE_MPIIO)
251:       args: -usempiio

253: TEST*/