Actual source code: ex51.c

petsc-3.14.3 2021-01-09
Report Typos and Errors
  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 it's elements is calculated.
 12: */

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

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

 31:   PetscInitialize(&argc, &argv, (char*) 0, help);if (ierr) return ierr;

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

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

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

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

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

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

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

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

 93:   /* Compute sum of parent vector elements */
 94:   VecSum(testvec, &sum);

 96:   /* to verify that the loaded data has been transferred */
 97:   if (sum != 150) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_PLIB,"Data has not been transferred from subvector to parent vector");
 98:   PetscPrintf(PETSC_COMM_WORLD,"VecSum on parent vec is : %e\n",sum);

100:   /* destroy parent vector, index set and exit */
101:   VecDestroy(&testvec);
102:   ISDestroy(&loadis);
103:   PetscFinalize();
104:   return ierr;
105: }

107: /*TEST

109:   build:
110:     requires: hdf5

112:   test:
113:     nsize: 4

115:   test:
116:     suffix: 2
117:     nsize: 4
118:     args: -hdf5

120: TEST*/