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, (char *)0, 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(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"));
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));
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: }
153: /*TEST
155: test:
156: suffix: serial
158: test:
159: suffix: parallel
160: nsize: 2
162: test:
163: suffix: cuda
164: nsize: 2
165: args: -vec_type cuda
166: requires: cuda
168: test:
169: suffix: uneven
170: nsize: 5
172: TEST*/