Actual source code: ex15.c
1: static char help[] = "Tests VecView() functionality with DMDA objects when using:"
2: "(i) a PetscViewer binary with MPI-IO support; and (ii) when the binary header is skipped.\n\n";
4: #include <petscdm.h>
5: #include <petscdmda.h>
7: #define DMDA_I 5
8: #define DMDA_J 4
9: #define DMDA_K 6
11: const PetscReal dmda_i_val[] = {1.10, 2.3006, 2.32444, 3.44006, 66.9009};
12: const PetscReal dmda_j_val[] = {0.0, 0.25, 0.5, 0.75};
13: const PetscReal dmda_k_val[] = {0.0, 1.1, 2.2, 3.3, 4.4, 5.5};
15: PetscErrorCode MyVecDump(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x)
16: {
17: MPI_Comm comm;
18: PetscViewer viewer;
19: PetscBool ismpiio, isskip;
21: PetscFunctionBeginUser;
22: PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
24: PetscCall(PetscViewerCreate(comm, &viewer));
25: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
26: if (skippheader) PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE));
27: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
28: if (usempiio) PetscCall(PetscViewerBinarySetUseMPIIO(viewer, PETSC_TRUE));
29: PetscCall(PetscViewerFileSetName(viewer, fname));
31: PetscCall(VecView(x, viewer));
33: PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &ismpiio));
34: if (ismpiio) PetscCall(PetscPrintf(comm, "*** PetscViewer[write] using MPI-IO ***\n"));
35: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &isskip));
36: if (isskip) PetscCall(PetscPrintf(comm, "*** PetscViewer[write] skipping header ***\n"));
38: PetscCall(PetscViewerDestroy(&viewer));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: PetscErrorCode MyVecLoad(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x)
43: {
44: MPI_Comm comm;
45: PetscViewer viewer;
46: PetscBool ismpiio, isskip;
48: PetscFunctionBeginUser;
49: PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
51: PetscCall(PetscViewerCreate(comm, &viewer));
52: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERBINARY));
53: if (skippheader) PetscCall(PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE));
54: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_READ));
55: if (usempiio) PetscCall(PetscViewerBinarySetUseMPIIO(viewer, PETSC_TRUE));
56: PetscCall(PetscViewerFileSetName(viewer, fname));
58: PetscCall(VecLoad(x, viewer));
60: PetscCall(PetscViewerBinaryGetSkipHeader(viewer, &isskip));
61: if (isskip) PetscCall(PetscPrintf(comm, "*** PetscViewer[load] skipping header ***\n"));
62: PetscCall(PetscViewerBinaryGetUseMPIIO(viewer, &ismpiio));
63: if (ismpiio) PetscCall(PetscPrintf(comm, "*** PetscViewer[load] using MPI-IO ***\n"));
65: PetscCall(PetscViewerDestroy(&viewer));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: PetscErrorCode DMDAVecGenerateEntries(DM dm, Vec a)
70: {
71: PetscScalar ****LA_v;
72: PetscInt i, j, k, l, si, sj, sk, ni, nj, nk, M, N, dof;
74: PetscFunctionBeginUser;
75: PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
76: PetscCall(DMDAGetCorners(dm, &si, &sj, &sk, &ni, &nj, &nk));
77: PetscCall(DMDAVecGetArrayDOF(dm, a, &LA_v));
78: for (k = sk; k < sk + nk; k++) {
79: for (j = sj; j < sj + nj; j++) {
80: for (i = si; i < si + ni; i++) {
81: PetscScalar test_value_s;
83: 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));
84: for (l = 0; l < dof; l++) LA_v[k][j][i][l] = (PetscScalar)dof * test_value_s + (PetscScalar)l;
85: }
86: }
87: }
88: PetscCall(DMDAVecRestoreArrayDOF(dm, a, &LA_v));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: PetscErrorCode HeaderlessBinaryReadCheck(DM dm, const char name[])
93: {
94: int fdes;
95: PetscScalar buffer[DMDA_I * DMDA_J * DMDA_K * 10];
96: PetscInt len, d, i, j, k, M, N, dof;
97: PetscMPIInt rank;
98: PetscBool dataverified = PETSC_TRUE;
100: PetscFunctionBeginUser;
101: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
102: PetscCall(DMDAGetInfo(dm, NULL, &M, &N, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
103: len = DMDA_I * DMDA_J * DMDA_K * dof;
104: if (rank == 0) {
105: PetscCall(PetscBinaryOpen(name, FILE_MODE_READ, &fdes));
106: PetscCall(PetscBinaryRead(fdes, buffer, len, NULL, PETSC_SCALAR));
107: PetscCall(PetscBinaryClose(fdes));
109: for (k = 0; k < DMDA_K; k++) {
110: for (j = 0; j < DMDA_J; j++) {
111: for (i = 0; i < DMDA_I; i++) {
112: for (d = 0; d < dof; d++) {
113: PetscScalar v, test_value_s, test_value;
114: PetscInt index;
116: 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));
117: test_value = (PetscScalar)dof * test_value_s + (PetscScalar)d;
119: index = dof * (i + j * M + k * M * N) + d;
120: v = PetscAbsScalar(test_value - buffer[index]);
121: #if defined(PETSC_USE_COMPLEX)
122: if ((PetscRealPart(v) > 1.0e-10) || (PetscImaginaryPart(v) > 1.0e-10)) {
123: 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));
124: dataverified = PETSC_FALSE;
125: }
126: #else
127: if (PetscRealPart(v) > 1.0e-10) {
128: 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));
129: dataverified = PETSC_FALSE;
130: }
131: #endif
132: }
133: }
134: }
135: }
136: if (dataverified) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Headerless read of data verified for: %s\n", name));
137: }
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: PetscErrorCode VecCompare(Vec a, Vec b)
142: {
143: PetscInt locmin[2], locmax[2];
144: PetscReal min[2], max[2];
145: Vec ref;
147: PetscFunctionBeginUser;
148: PetscCall(VecMin(a, &locmin[0], &min[0]));
149: PetscCall(VecMax(a, &locmax[0], &max[0]));
151: PetscCall(VecMin(b, &locmin[1], &min[1]));
152: PetscCall(VecMax(b, &locmax[1], &max[1]));
154: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecCompare\n"));
155: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " min(a) = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[0], locmin[0]));
156: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " max(a) = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[0], locmax[0]));
158: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " min(b) = %+1.2e [loc %" PetscInt_FMT "]\n", (double)min[1], locmin[1]));
159: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " max(b) = %+1.2e [loc %" PetscInt_FMT "]\n", (double)max[1], locmax[1]));
161: PetscCall(VecDuplicate(a, &ref));
162: PetscCall(VecCopy(a, ref));
163: PetscCall(VecAXPY(ref, -1.0, b));
164: PetscCall(VecMin(ref, &locmin[0], &min[0]));
165: if (PetscAbsReal(min[0]) > 1.0e-10) {
166: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " ERROR: min(a-b) > 1.0e-10\n"));
167: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " min(a-b) = %+1.10e\n", (double)PetscAbsReal(min[0])));
168: } else {
169: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " min(a-b) < 1.0e-10\n"));
170: }
171: PetscCall(VecDestroy(&ref));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: PetscErrorCode TestDMDAVec(PetscBool usempiio)
176: {
177: DM dm;
178: Vec x_ref, x_test;
179: PetscBool skipheader = PETSC_TRUE;
181: PetscFunctionBeginUser;
182: if (!usempiio) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s\n", PETSC_FUNCTION_NAME));
183: else PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s [using mpi-io]\n", PETSC_FUNCTION_NAME));
184: 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));
185: PetscCall(DMSetFromOptions(dm));
186: PetscCall(DMSetUp(dm));
188: PetscCall(DMCreateGlobalVector(dm, &x_ref));
189: PetscCall(DMDAVecGenerateEntries(dm, x_ref));
191: if (!usempiio) {
192: PetscCall(MyVecDump("dmda.pbvec", skipheader, PETSC_FALSE, x_ref));
193: } else {
194: PetscCall(MyVecDump("dmda-mpiio.pbvec", skipheader, PETSC_TRUE, x_ref));
195: }
197: PetscCall(DMCreateGlobalVector(dm, &x_test));
199: if (!usempiio) {
200: PetscCall(MyVecLoad("dmda.pbvec", skipheader, usempiio, x_test));
201: } else {
202: PetscCall(MyVecLoad("dmda-mpiio.pbvec", skipheader, usempiio, x_test));
203: }
205: PetscCall(VecCompare(x_ref, x_test));
207: if (!usempiio) {
208: PetscCall(HeaderlessBinaryReadCheck(dm, "dmda.pbvec"));
209: } else {
210: PetscCall(HeaderlessBinaryReadCheck(dm, "dmda-mpiio.pbvec"));
211: }
212: PetscCall(VecDestroy(&x_ref));
213: PetscCall(VecDestroy(&x_test));
214: PetscCall(DMDestroy(&dm));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: int main(int argc, char **args)
219: {
220: PetscBool usempiio = PETSC_FALSE;
222: PetscFunctionBeginUser;
223: PetscCall(PetscInitialize(&argc, &args, NULL, help));
224: PetscCall(PetscOptionsGetBool(NULL, NULL, "-usempiio", &usempiio, NULL));
225: if (!usempiio) {
226: PetscCall(TestDMDAVec(PETSC_FALSE));
227: } else {
228: #if defined(PETSC_HAVE_MPIIO)
229: PetscCall(TestDMDAVec(PETSC_TRUE));
230: #else
231: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n"));
232: #endif
233: }
234: PetscCall(PetscFinalize());
235: return 0;
236: }
238: /*TEST
240: test:
242: test:
243: suffix: 2
244: nsize: 12
246: test:
247: suffix: 3
248: nsize: 12
249: requires: defined(PETSC_HAVE_MPIIO)
250: args: -usempiio
252: TEST*/