Actual source code: ex44.c
1: static char help[] = "Test VecConcatenate both in serial and parallel.\n";
3: #include <petscvec.h>
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;
14: PetscFunctionBeginUser;
15: PetscCall(PetscInitialize(&argc, &args, NULL, help));
16: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nx", &nx, &flg));
17: if (!flg) nx = 3;
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"));
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"));
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"));
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"));
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(VecScatterCreate(y, x_is[i], x[i], NULL, &y_to_x));
106: PetscCall(VecScatterBegin(y_to_x, y, x_test, INSERT_VALUES, SCATTER_FORWARD));
107: PetscCall(VecScatterEnd(y_to_x, y, x_test, INSERT_VALUES, SCATTER_FORWARD));
108: PetscCall(VecScatterDestroy(&y_to_x));
109: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecScatter for Y -> X[%" PetscInt_FMT "]\n", i));
110: PetscCall(VecView(x_test, PETSC_VIEWER_STDOUT_WORLD));
111: x_equal = PETSC_FALSE;
112: PetscCall(VecEqual(x_test, x[i], &x_equal));
113: if (!x_equal) {
114: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " FAIL\n"));
115: } else {
116: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PASS\n"));
117: }
118: PetscCall(VecDestroy(&x_test));
119: }
120: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "------------------------------------------------------------\n"));
122: /* ---------- using VecScatter to communicate data from X[i] to Y for all i ----------- */
123: PetscCall(VecZeroEntries(y_test));
124: for (i = 0; i < nx; i++) {
125: PetscCall(VecScatterCreate(x[i], NULL, y, x_is[i], &x_to_y));
126: PetscCall(VecScatterBegin(x_to_y, x[i], y_test, INSERT_VALUES, SCATTER_FORWARD));
127: PetscCall(VecScatterEnd(x_to_y, x[i], y_test, INSERT_VALUES, SCATTER_FORWARD));
128: PetscCall(VecScatterDestroy(&x_to_y));
129: }
130: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Testing VecScatter for X[:] -> Y\n"));
131: PetscCall(VecView(y_test, PETSC_VIEWER_STDOUT_WORLD));
132: y_equal = PETSC_FALSE;
133: PetscCall(VecEqual(y_test, y, &y_equal));
134: if (!y_equal) {
135: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " FAIL\n"));
136: } else {
137: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " PASS\n"));
138: }
139: PetscCall(VecDestroy(&y_test));
141: PetscCall(VecDestroy(&y));
142: for (i = 0; i < nx; i++) {
143: PetscCall(VecDestroy(&x[i]));
144: PetscCall(ISDestroy(&x_is[i]));
145: }
146: PetscCall(PetscFree(x));
147: PetscCall(PetscFree(x_is));
148: PetscCall(PetscFinalize());
149: return 0;
150: }
152: /*TEST
154: test:
155: suffix: serial
157: test:
158: suffix: parallel
159: nsize: 2
161: test:
162: suffix: cuda
163: nsize: 2
164: args: -vec_type cuda
165: requires: cuda
167: test:
168: suffix: uneven
169: nsize: 5
171: TEST*/