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;
23: PetscObjectGetComm((PetscObject)x, &comm);
25: PetscViewerCreate(comm, &viewer);
26: PetscViewerSetType(viewer, PETSCVIEWERBINARY);
27: if (skippheader) PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE);
28: PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
29: if (usempiio) PetscViewerBinarySetUseMPIIO(viewer, PETSC_TRUE);
30: PetscViewerFileSetName(viewer, fname);
32: VecView(x, viewer);
34: PetscViewerBinaryGetUseMPIIO(viewer, &ismpiio);
35: if (ismpiio) PetscPrintf(comm, "*** PetscViewer[write] using MPI-IO ***\n");
36: PetscViewerBinaryGetSkipHeader(viewer, &isskip);
37: if (isskip) PetscPrintf(comm, "*** PetscViewer[write] skipping header ***\n");
39: PetscViewerDestroy(&viewer);
40: return 0;
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;
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 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;
76: DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL);
77: DMDAGetCorners(dm, &si, &sj, &sk, &ni, &nj, &nk);
78: 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: DMDAVecRestoreArrayDOF(dm, a, &LA_v);
90: return 0;
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;
102: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
103: 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: PetscBinaryOpen(name, FILE_MODE_READ, &fdes);
107: PetscBinaryRead(fdes, buffer, len, NULL, PETSC_SCALAR);
108: 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: 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: 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) PetscPrintf(PETSC_COMM_SELF, "Headerless read of data verified for: %s\n", name);
138: }
139: return 0;
140: }
142: PetscErrorCode VecCompare(Vec a, Vec b)
143: {
144: PetscInt locmin[2], locmax[2];
145: PetscReal min[2], max[2];
146: Vec ref;
149: VecMin(a, &locmin[0], &min[0]);
150: VecMax(a, &locmax[0], &max[0]);
152: VecMin(b, &locmin[1], &min[1]);
153: VecMax(b, &locmax[1], &max[1]);
155: PetscPrintf(PETSC_COMM_WORLD, "VecCompare\n");
156: PetscPrintf(PETSC_COMM_WORLD, " min(a) = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[0], locmin[0]);
157: PetscPrintf(PETSC_COMM_WORLD, " max(a) = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[0], locmax[0]);
159: PetscPrintf(PETSC_COMM_WORLD, " min(b) = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[1], locmin[1]);
160: PetscPrintf(PETSC_COMM_WORLD, " max(b) = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[1], locmax[1]);
162: VecDuplicate(a, &ref);
163: VecCopy(a, ref);
164: VecAXPY(ref, -1.0, b);
165: VecMin(ref, &locmin[0], &min[0]);
166: if (PetscAbsReal(min[0]) > 1.0e-10) {
167: PetscPrintf(PETSC_COMM_WORLD, " ERROR: min(a-b) > 1.0e-10\n");
168: PetscPrintf(PETSC_COMM_WORLD, " min(a-b) = %+1.10e\n", (double)PetscAbsReal(min[0]));
169: } else {
170: PetscPrintf(PETSC_COMM_WORLD, " min(a-b) < 1.0e-10\n");
171: }
172: VecDestroy(&ref);
173: return 0;
174: }
176: PetscErrorCode TestDMDAVec(PetscBool usempiio)
177: {
178: DM dm;
179: Vec x_ref, x_test;
180: PetscBool skipheader = PETSC_TRUE;
183: if (!usempiio) PetscPrintf(PETSC_COMM_WORLD, "%s\n", PETSC_FUNCTION_NAME);
184: else PetscPrintf(PETSC_COMM_WORLD, "%s [using mpi-io]\n", PETSC_FUNCTION_NAME);
185: 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: DMSetFromOptions(dm);
187: DMSetUp(dm);
189: DMCreateGlobalVector(dm, &x_ref);
190: DMDAVecGenerateEntries(dm, x_ref);
192: if (!usempiio) {
193: MyVecDump("dmda.pbvec", skipheader, PETSC_FALSE, x_ref);
194: } else {
195: MyVecDump("dmda-mpiio.pbvec", skipheader, PETSC_TRUE, x_ref);
196: }
198: DMCreateGlobalVector(dm, &x_test);
200: if (!usempiio) {
201: MyVecLoad("dmda.pbvec", skipheader, usempiio, x_test);
202: } else {
203: MyVecLoad("dmda-mpiio.pbvec", skipheader, usempiio, x_test);
204: }
206: VecCompare(x_ref, x_test);
208: if (!usempiio) {
209: HeaderlessBinaryReadCheck(dm, "dmda.pbvec");
210: } else {
211: HeaderlessBinaryReadCheck(dm, "dmda-mpiio.pbvec");
212: }
213: VecDestroy(&x_ref);
214: VecDestroy(&x_test);
215: DMDestroy(&dm);
216: return 0;
217: }
219: int main(int argc, char **args)
220: {
221: PetscBool usempiio = PETSC_FALSE;
224: PetscInitialize(&argc, &args, (char *)0, help);
225: PetscOptionsGetBool(NULL, NULL, "-usempiio", &usempiio, NULL);
226: if (!usempiio) {
227: TestDMDAVec(PETSC_FALSE);
228: } else {
229: #if defined(PETSC_HAVE_MPIIO)
230: TestDMDAVec(PETSC_TRUE);
231: #else
232: PetscPrintf(PETSC_COMM_WORLD, "Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n");
233: #endif
234: }
235: 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*/