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) PetscCall(TestDMDAVec(PETSC_FALSE));
    226:   else {
    227: #if defined(PETSC_HAVE_MPIIO)
    228:     PetscCall(TestDMDAVec(PETSC_TRUE));
    229: #else
    230:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: Executing TestDMDAVec(PETSC_TRUE) requires a working MPI-2 implementation\n"));
    231: #endif
    232:   }
    233:   PetscCall(PetscFinalize());
    234:   return 0;
    235: }

    237: /*TEST

    239:    test:

    241:    test:
    242:       suffix: 2
    243:       nsize: 12

    245:    test:
    246:       suffix: 3
    247:       nsize: 12
    248:       requires: defined(PETSC_HAVE_MPIIO)
    249:       args: -usempiio

    251: TEST*/