1 static char help[] = "Nest vector set subvector functionality.\n\n";
2
3 #include <petscvec.h>
4
test_vec_ops(void)5 PetscErrorCode test_vec_ops(void)
6 {
7 Vec X, Y, a, b;
8 Vec c, d, e, f, g, h;
9 PetscScalar val;
10 PetscInt tmp_ind[2];
11 Vec tmp_buf[2];
12
13 PetscFunctionBegin;
14 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "============== %s ==============\n", PETSC_FUNCTION_NAME));
15
16 /* create 4 worker vectors */
17 PetscCall(VecCreate(PETSC_COMM_WORLD, &c));
18 PetscCall(VecSetSizes(c, PETSC_DECIDE, 4));
19 PetscCall(VecSetType(c, VECMPI));
20 PetscCall(VecDuplicate(c, &d));
21 PetscCall(VecDuplicate(c, &e));
22 PetscCall(VecDuplicate(c, &f));
23
24 /* create two more workers of different sizes */
25 PetscCall(VecCreate(PETSC_COMM_WORLD, &g));
26 PetscCall(VecSetSizes(g, PETSC_DECIDE, 6));
27 PetscCall(VecSetType(g, VECMPI));
28 PetscCall(VecCreate(PETSC_COMM_WORLD, &h));
29 PetscCall(VecSetSizes(h, PETSC_DECIDE, 8));
30 PetscCall(VecSetType(h, VECMPI));
31
32 /* set the 6 vectors to some numbers */
33 PetscCall(VecSet(c, 1.0));
34 PetscCall(VecSet(d, 2.0));
35 PetscCall(VecSet(e, 3.0));
36 PetscCall(VecSet(f, 4.0));
37 PetscCall(VecSet(g, 5.0));
38 PetscCall(VecSet(h, 6.0));
39
40 /* assemble a */
41 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "a = [c d] \n"));
42 tmp_buf[0] = c;
43 tmp_buf[1] = d;
44
45 PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &a));
46 PetscCall(VecView(a, PETSC_VIEWER_STDOUT_WORLD));
47 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "a = [d c] \n"));
48 PetscCall(VecNestSetSubVec(a, 1, c));
49 PetscCall(VecNestSetSubVec(a, 0, d));
50 PetscCall(VecAssemblyBegin(a));
51 PetscCall(VecAssemblyEnd(a));
52 PetscCall(VecView(a, PETSC_VIEWER_STDOUT_WORLD));
53
54 /* assemble b */
55 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "b = [e f] \n"));
56 tmp_buf[0] = e;
57 tmp_buf[1] = f;
58
59 PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &b));
60 PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD));
61 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "b = [f e] \n"));
62 PetscCall(VecNestSetSubVec(b, 1, e));
63 PetscCall(VecNestSetSubVec(b, 0, f));
64 PetscCall(VecAssemblyBegin(b));
65 PetscCall(VecAssemblyEnd(b));
66 PetscCall(VecView(b, PETSC_VIEWER_STDOUT_WORLD));
67
68 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X = [a b] \n"));
69 tmp_buf[0] = a;
70 tmp_buf[1] = b;
71
72 PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_buf, &X));
73 PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
74 PetscCall(VecDot(X, X, &val));
75 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %g \n", (double)PetscRealPart(val)));
76
77 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X = [b a] \n"));
78 /* re-order components of X */
79 PetscCall(VecNestSetSubVec(X, 1, a));
80 PetscCall(VecNestSetSubVec(X, 0, b));
81 PetscCall(VecAssemblyBegin(X));
82 PetscCall(VecAssemblyEnd(X));
83 PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
84 PetscCall(VecDot(X, X, &val));
85 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %g \n", (double)PetscRealPart(val)));
86
87 /* re-assemble X */
88 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X = [g h] \n"));
89 PetscCall(VecNestSetSubVec(X, 1, g));
90 PetscCall(VecNestSetSubVec(X, 0, h));
91 PetscCall(VecAssemblyBegin(X));
92 PetscCall(VecAssemblyEnd(X));
93 PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
94 PetscCall(VecDot(X, X, &val));
95 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "X.X = %g \n", (double)PetscRealPart(val)));
96
97 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Y = X \n"));
98 PetscCall(VecDuplicate(X, &Y));
99 PetscCall(VecCopy(X, Y));
100 PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
101 PetscCall(VecDot(Y, Y, &val));
102 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Y.Y = %g \n", (double)PetscRealPart(val)));
103
104 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Y = [a b] \n"));
105 tmp_buf[0] = a;
106 tmp_buf[1] = b;
107 tmp_ind[0] = 0;
108 tmp_ind[1] = 1;
109
110 PetscCall(VecNestSetSubVecs(Y, 2, tmp_ind, tmp_buf));
111 PetscCall(VecView(Y, PETSC_VIEWER_STDOUT_WORLD));
112
113 PetscCall(VecDestroy(&c));
114 PetscCall(VecDestroy(&d));
115 PetscCall(VecDestroy(&e));
116 PetscCall(VecDestroy(&f));
117 PetscCall(VecDestroy(&g));
118 PetscCall(VecDestroy(&h));
119 PetscCall(VecDestroy(&a));
120 PetscCall(VecDestroy(&b));
121 PetscCall(VecDestroy(&X));
122 PetscCall(VecDestroy(&Y));
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125
main(int argc,char ** args)126 int main(int argc, char **args)
127 {
128 PetscFunctionBeginUser;
129 PetscCall(PetscInitialize(&argc, &args, NULL, help));
130 PetscCall(test_vec_ops());
131 PetscCall(PetscFinalize());
132 return 0;
133 }
134
135 /*TEST
136
137 test:
138
139 TEST*/
140