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*/