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";
3
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 */
14
15 #include <petscvec.h>
16 #include <petscviewerhdf5.h>
17
main(int argc,char ** argv)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;
31
32 PetscFunctionBeginUser;
33 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
34
35 /* parse input options to determine I/O format */
36 PetscCall(PetscOptionsGetBool(NULL, NULL, "-hdf5", &usehdf5, NULL));
37
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));
43
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"));
50
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));
62
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));
76
77 /* Create subvector using the index set created above */
78 PetscCall(VecGetSubVector(testvec, loadis, &loadvec));
79 PetscCall(PetscObjectSetName((PetscObject)loadvec, "temp"));
80
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));
91
92 /* Restore subvector to transfer loaded data into parent vector */
93 PetscCall(VecRestoreSubVector(testvec, loadis, &loadvec));
94
95 /* Compute sum of parent vector elements */
96 PetscCall(VecSum(testvec, &sum));
97 PetscCall(VecMean(testvec, &mean));
98
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)));
103
104 /* destroy parent vector, index set and exit */
105 PetscCall(VecDestroy(&testvec));
106 PetscCall(ISDestroy(&loadis));
107 PetscCall(PetscFinalize());
108 return 0;
109 }
110
111 /*TEST
112
113 build:
114 requires: hdf5
115
116 test:
117 nsize: 4
118
119 test:
120 suffix: 2
121 nsize: 4
122 args: -hdf5
123
124 TEST*/
125