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, (char *)0, 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*/