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 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 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 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 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 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 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 218 int main(int argc, char **args) 219 { 220 PetscBool usempiio = PETSC_FALSE; 221 222 PetscFunctionBeginUser; 223 PetscCall(PetscInitialize(&argc, &args, (char *)0, 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