1 static char help[] = "Test VecConcatenate both in serial and parallel.\n";
2
3 #include <petscvec.h>
4
main(int argc,char ** args)5 int main(int argc, char **args)
6 {
7 Vec *x, x_test, y, y_test;
8 IS *x_is;
9 VecScatter y_to_x, x_to_y;
10 PetscInt i, j, nx, shift, x_size, y_size, *idx;
11 PetscScalar *vals;
12 PetscBool flg, x_equal, y_equal;
13
14 PetscFunctionBeginUser;
15 PetscCall(PetscInitialize(&argc, &args, NULL, help));
16 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &nx, &flg));
17 if (!flg) nx = 3;
18
19 y_size = 0;
20 shift = 0;
21 PetscCall(PetscMalloc1(nx, &x));
22 for (i = 0; i < nx; i++) {
23 x_size = 2 * (i + 1);
24 y_size += x_size;
25 PetscCall(VecCreate(PETSC_COMM_WORLD, &x[i]));
26 PetscCall(VecSetSizes(x[i], PETSC_DECIDE, x_size));
27 PetscCall(VecSetFromOptions(x[i]));
28 PetscCall(VecSetUp(x[i]));
29 PetscCall(PetscMalloc2(x_size, &idx, x_size, &vals));
30 for (j = 0; j < x_size; j++) {
31 idx[j] = j;
32 vals[j] = (PetscScalar)(shift + j + 1);
33 }
34 shift += x_size;
35 PetscCall(VecSetValues(x[i], x_size, (const PetscInt *)idx, (const PetscScalar *)vals, INSERT_VALUES));
36 PetscCall(VecAssemblyBegin(x[i]));
37 PetscCall(VecAssemblyEnd(x[i]));
38 PetscCall(PetscFree2(idx, vals));
39 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Original X[%" PetscInt_FMT "] vector\n", i));
40 PetscCall(VecView(x[i], PETSC_VIEWER_STDOUT_WORLD));
41 }
42 PetscCall(VecCreate(PETSC_COMM_WORLD, &y));
43 PetscCall(VecSetSizes(y, PETSC_DECIDE, y_size));
44 PetscCall(VecSetFromOptions(y));
45 PetscCall(VecSetUp(y));
46 PetscCall(PetscMalloc2(y_size, &idx, y_size, &vals));
47 for (j = 0; j < y_size; j++) {
48 idx[j] = j;
49 vals[j] = (PetscScalar)(j + 1);
50 }
51 PetscCall(VecSetValues(y, y_size, (const PetscInt *)idx, (const PetscScalar *)vals, INSERT_VALUES));
52 PetscCall(VecAssemblyBegin(y));
53 PetscCall(VecAssemblyEnd(y));
54 PetscCall(PetscFree2(idx, vals));
55 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Expected Y vector\n"));
56 PetscCall(VecView(y, PETSC_VIEWER_STDOUT_WORLD));
57 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));
58
59 /* ---------- base VecConcatenate() test ----------- */
60 PetscCall(VecConcatenate(nx, (const Vec *)x, &y_test, &x_is));
61 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecConcatenate() for Y = [X[1], X[2], ...]\n"));
62 PetscCall(VecView(y_test, PETSC_VIEWER_STDOUT_WORLD));
63 y_equal = PETSC_FALSE;
64 PetscCall(VecEqual(y_test, y, &y_equal));
65 if (!y_equal) {
66 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " FAIL\n"));
67 } else {
68 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PASS\n"));
69 }
70 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));
71
72 /* ---------- test VecConcatenate() without IS (checks for dangling memory from IS) ----------- */
73 PetscCall(VecDestroy(&y_test));
74 PetscCall(VecConcatenate(nx, (const Vec *)x, &y_test, NULL));
75 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecConcatenate() for Y = [X[1], X[2], ...] w/o IS\n"));
76 PetscCall(VecView(y_test, PETSC_VIEWER_STDOUT_WORLD));
77 y_equal = PETSC_FALSE;
78 PetscCall(VecEqual(y_test, y, &y_equal));
79 if (!y_equal) {
80 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " FAIL\n"));
81 } else {
82 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PASS\n"));
83 }
84 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));
85
86 /* ---------- using index sets on expected Y instead of concatenated Y ----------- */
87 for (i = 0; i < nx; i++) {
88 PetscCall(VecGetSubVector(y, x_is[i], &x_test));
89 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing index set for X[%" PetscInt_FMT "] component\n", i));
90 PetscCall(VecView(x_test, PETSC_VIEWER_STDOUT_WORLD));
91 x_equal = PETSC_FALSE;
92 PetscCall(VecEqual(x_test, x[i], &x_equal));
93 if (!x_equal) {
94 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " FAIL\n"));
95 } else {
96 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PASS\n"));
97 }
98 PetscCall(VecRestoreSubVector(y, x_is[i], &x_test));
99 }
100 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));
101
102 /* ---------- using VecScatter to communicate data from Y to X[i] for all i ----------- */
103 for (i = 0; i < nx; i++) {
104 PetscCall(VecDuplicate(x[i], &x_test));
105 PetscCall(VecZeroEntries(x_test));
106 PetscCall(VecScatterCreate(y, x_is[i], x[i], NULL, &y_to_x));
107 PetscCall(VecScatterBegin(y_to_x, y, x_test, INSERT_VALUES, SCATTER_FORWARD));
108 PetscCall(VecScatterEnd(y_to_x, y, x_test, INSERT_VALUES, SCATTER_FORWARD));
109 PetscCall(VecScatterDestroy(&y_to_x));
110 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecScatter for Y -> X[%" PetscInt_FMT "]\n", i));
111 PetscCall(VecView(x_test, PETSC_VIEWER_STDOUT_WORLD));
112 x_equal = PETSC_FALSE;
113 PetscCall(VecEqual(x_test, x[i], &x_equal));
114 if (!x_equal) {
115 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " FAIL\n"));
116 } else {
117 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PASS\n"));
118 }
119 PetscCall(VecDestroy(&x_test));
120 }
121 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));
122
123 /* ---------- using VecScatter to communicate data from X[i] to Y for all i ----------- */
124 PetscCall(VecZeroEntries(y_test));
125 for (i = 0; i < nx; i++) {
126 PetscCall(VecScatterCreate(x[i], NULL, y, x_is[i], &x_to_y));
127 PetscCall(VecScatterBegin(x_to_y, x[i], y_test, INSERT_VALUES, SCATTER_FORWARD));
128 PetscCall(VecScatterEnd(x_to_y, x[i], y_test, INSERT_VALUES, SCATTER_FORWARD));
129 PetscCall(VecScatterDestroy(&x_to_y));
130 }
131 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecScatter for X[:] -> Y\n"));
132 PetscCall(VecView(y_test, PETSC_VIEWER_STDOUT_WORLD));
133 y_equal = PETSC_FALSE;
134 PetscCall(VecEqual(y_test, y, &y_equal));
135 if (!y_equal) {
136 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " FAIL\n"));
137 } else {
138 PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PASS\n"));
139 }
140 PetscCall(VecDestroy(&y_test));
141
142 PetscCall(VecDestroy(&y));
143 for (i = 0; i < nx; i++) {
144 PetscCall(VecDestroy(&x[i]));
145 PetscCall(ISDestroy(&x_is[i]));
146 }
147 PetscCall(PetscFree(x));
148 PetscCall(PetscFree(x_is));
149 PetscCall(PetscFinalize());
150 return 0;
151 }
152
153 /*TEST
154
155 test:
156 suffix: serial
157
158 test:
159 suffix: parallel
160 nsize: 2
161
162 test:
163 suffix: cuda
164 nsize: 2
165 args: -vec_type cuda
166 requires: cuda
167
168 test:
169 suffix: uneven
170 nsize: 5
171
172 TEST*/
173