Actual source code: ex51.c

  1: static char help[] = "Test integrity of subvector data, use \n\
  2: use -hdf5 to specify HDF5 viewer format for subvector I/O \n\n";

  4: /*
  5:    Tests for transfer of data from subvectors to parent vectors after
  6:    loading data into subvector. This routine does the following : creates
  7:    a vector of size 50, sets it to 2 and saves it to disk. Creates a
  8:    vector of size 100, set it to 1 and extracts the last 50 elements
  9:    as a subvector. Loads the saved vector from disk into the subvector
 10:    and restores the subvector. To verify that the data has been loaded
 11:    into the parent vector, the sum of its elements is calculated.
 12:    The arithmetic mean is also calculated in order to test VecMean().
 13: */

 15: #include <petscvec.h>
 16: #include <petscviewerhdf5.h>

 18: int main(int argc, char **argv)
 19: {
 20:   Vec         testvec;         /* parent vector of size 100 */
 21:   Vec         loadvec;         /* subvector extracted from the parent vector */
 22:   Vec         writevec;        /* vector used to save data to be loaded by loadvec */
 23:   IS          loadis;          /* index set to extract last 50 elements of testvec */
 24:   PetscInt    low, high;       /* used to store vecownership output */
 25:   PetscInt    issize, isstart; /* index set params */
 26:   PetscInt    skipuntil = 50;  /* parameter to slice the last N elements of parent vec */
 27:   PetscViewer viewer;          /* viewer for I/O */
 28:   PetscScalar sum;             /* used to test sum of parent vector elements */
 29:   PetscScalar mean;            /* used to test mean of parent vector elements */
 30:   PetscBool   usehdf5 = PETSC_FALSE;

 32:   PetscFunctionBeginUser;
 33:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));

 35:   /* parse input options to determine I/O format */
 36:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-hdf5", &usehdf5, NULL));

 38:   /* Create parent vector with 100 elements, set it to 1 */
 39:   PetscCall(VecCreate(PETSC_COMM_WORLD, &testvec));
 40:   PetscCall(VecSetSizes(testvec, PETSC_DECIDE, 100));
 41:   PetscCall(VecSetUp(testvec));
 42:   PetscCall(VecSet(testvec, (PetscScalar)1));

 44:   /* Create a vector with 50 elements, set it to 2. */
 45:   PetscCall(VecCreate(PETSC_COMM_WORLD, &writevec));
 46:   PetscCall(VecSetSizes(writevec, PETSC_DECIDE, 50));
 47:   PetscCall(VecSetUp(writevec));
 48:   PetscCall(VecSet(writevec, (PetscScalar)2));
 49:   PetscCall(PetscObjectSetName((PetscObject)writevec, "temp"));

 51:   /* Save to disk in specified format, destroy vector & viewer */
 52:   if (usehdf5) {
 53:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "writing vector in hdf5 to vector.dat ...\n"));
 54:     PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE, &viewer));
 55:   } else {
 56:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "writing vector in binary to vector.dat ...\n"));
 57:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_WRITE, &viewer));
 58:   }
 59:   PetscCall(VecView(writevec, viewer));
 60:   PetscCall(VecDestroy(&writevec));
 61:   PetscCall(PetscViewerDestroy(&viewer));

 63:   /* Create index sets on each mpi rank to select the last 50 elements of parent vec */
 64:   PetscCall(VecGetOwnershipRange(testvec, &low, &high));
 65:   if (low >= skipuntil) {
 66:     isstart = low;
 67:     issize  = high - low;
 68:   } else if (low <= skipuntil && high >= skipuntil) {
 69:     isstart = skipuntil;
 70:     issize  = high - skipuntil;
 71:   } else {
 72:     isstart = low;
 73:     issize  = 0;
 74:   }
 75:   PetscCall(ISCreateStride(PETSC_COMM_WORLD, issize, isstart, 1, &loadis));

 77:   /* Create subvector using the index set created above */
 78:   PetscCall(VecGetSubVector(testvec, loadis, &loadvec));
 79:   PetscCall(PetscObjectSetName((PetscObject)loadvec, "temp"));

 81:   /* Load the previously saved vector into the subvector, destroy viewer */
 82:   if (usehdf5) {
 83:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "reading vector in hdf5 from vector.dat ...\n"));
 84:     PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_READ, &viewer));
 85:   } else {
 86:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "reading vector in binary from vector.dat ...\n"));
 87:     PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "vector.dat", FILE_MODE_READ, &viewer));
 88:   }
 89:   PetscCall(VecLoad(loadvec, viewer));
 90:   PetscCall(PetscViewerDestroy(&viewer));

 92:   /* Restore subvector to transfer loaded data into parent vector */
 93:   PetscCall(VecRestoreSubVector(testvec, loadis, &loadvec));

 95:   /* Compute sum of parent vector elements */
 96:   PetscCall(VecSum(testvec, &sum));
 97:   PetscCall(VecMean(testvec, &mean));

 99:   /* to verify that the loaded data has been transferred */
100:   PetscCheck(sum == 150, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Data has not been transferred from subvector to parent vector");
101:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecSum on parent vec is : %e\n", (double)PetscAbsScalar(sum)));
102:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecMean on parent vec is : %e\n", (double)PetscAbsScalar(mean)));

104:   /* destroy parent vector, index set and exit */
105:   PetscCall(VecDestroy(&testvec));
106:   PetscCall(ISDestroy(&loadis));
107:   PetscCall(PetscFinalize());
108:   return 0;
109: }

111: /*TEST

113:   build:
114:     requires: hdf5

116:   test:
117:     nsize: 4

119:   test:
120:     suffix: 2
121:     nsize: 4
122:     args: -hdf5

124: TEST*/