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";
3
4 #include <petscdm.h>
5 #include <petscdmda.h>
6
7 #define DMDA_I 5
8 #define DMDA_J 4
9 #define DMDA_K 6
10
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};
14
MyVecDump(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)15 PetscErrorCode MyVecDump(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x)
16 {
17 MPI_Comm comm;
18 PetscViewer viewer;
19 PetscBool ismpiio, isskip;
20
21 PetscFunctionBeginUser;
22 PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
23
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));
30
31 PetscCall(VecView(x, viewer));
32
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"));
37
38 PetscCall(PetscViewerDestroy(&viewer));
39 PetscFunctionReturn(PETSC_SUCCESS);
40 }
41
MyVecLoad(const char fname[],PetscBool skippheader,PetscBool usempiio,Vec x)42 PetscErrorCode MyVecLoad(const char fname[], PetscBool skippheader, PetscBool usempiio, Vec x)
43 {
44 MPI_Comm comm;
45 PetscViewer viewer;
46 PetscBool ismpiio, isskip;
47
48 PetscFunctionBeginUser;
49 PetscCall(PetscObjectGetComm((PetscObject)x, &comm));
50
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));
57
58 PetscCall(VecLoad(x, viewer));
59
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"));
64
65 PetscCall(PetscViewerDestroy(&viewer));
66 PetscFunctionReturn(PETSC_SUCCESS);
67 }
68
DMDAVecGenerateEntries(DM dm,Vec a)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;
73
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;
82
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 }
91
HeaderlessBinaryReadCheck(DM dm,const char name[])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;
99
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));
108
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;
115
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;
118
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 }
140
VecCompare(Vec a,Vec b)141 PetscErrorCode VecCompare(Vec a, Vec b)
142 {
143 PetscInt locmin[2], locmax[2];
144 PetscReal min[2], max[2];
145 Vec ref;
146
147 PetscFunctionBeginUser;
148 PetscCall(VecMin(a, &locmin[0], &min[0]));
149 PetscCall(VecMax(a, &locmax[0], &max[0]));
150
151 PetscCall(VecMin(b, &locmin[1], &min[1]));
152 PetscCall(VecMax(b, &locmax[1], &max[1]));
153
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]));
157
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]));
160
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 }
174
TestDMDAVec(PetscBool usempiio)175 PetscErrorCode TestDMDAVec(PetscBool usempiio)
176 {
177 DM dm;
178 Vec x_ref, x_test;
179 PetscBool skipheader = PETSC_TRUE;
180
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));
187
188 PetscCall(DMCreateGlobalVector(dm, &x_ref));
189 PetscCall(DMDAVecGenerateEntries(dm, x_ref));
190
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 }
196
197 PetscCall(DMCreateGlobalVector(dm, &x_test));
198
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 }
204
205 PetscCall(VecCompare(x_ref, x_test));
206
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 }
217
main(int argc,char ** args)218 int main(int argc, char **args)
219 {
220 PetscBool usempiio = PETSC_FALSE;
221
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 }
237
238 /*TEST
239
240 test:
241
242 test:
243 suffix: 2
244 nsize: 12
245
246 test:
247 suffix: 3
248 nsize: 12
249 requires: defined(PETSC_HAVE_MPIIO)
250 args: -usempiio
251
252 TEST*/
253