Actual source code: ex22.c
1: static const char help[] = "Test MatNest solving a linear system\n\n";
3: #include <petscksp.h>
5: PetscErrorCode test_solve(void)
6: {
7: Mat A11, A12, A21, A22, A, tmp[2][2];
8: KSP ksp;
9: PC pc;
10: Vec b, x, f, h, diag, x1, x2;
11: Vec tmp_x[2], *_tmp_x;
12: PetscInt n, np, i, j;
14: PetscFunctionBeginUser;
15: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME));
17: n = 3;
18: np = 2;
19: /* Create matrices */
20: /* A11 */
21: PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
22: PetscCall(VecSetSizes(diag, PETSC_DECIDE, n));
23: PetscCall(VecSetFromOptions(diag));
25: PetscCall(VecSet(diag, 1.0 / 10.0)); /* so inverse = diag(10) */
27: /* As a test, create a diagonal matrix for A11 */
28: PetscCall(MatCreate(PETSC_COMM_WORLD, &A11));
29: PetscCall(MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n));
30: PetscCall(MatSetType(A11, MATAIJ));
31: PetscCall(MatSeqAIJSetPreallocation(A11, n, NULL));
32: PetscCall(MatMPIAIJSetPreallocation(A11, np, NULL, np, NULL));
33: PetscCall(MatDiagonalSet(A11, diag, INSERT_VALUES));
35: PetscCall(VecDestroy(&diag));
37: /* A12 */
38: PetscCall(MatCreate(PETSC_COMM_WORLD, &A12));
39: PetscCall(MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np));
40: PetscCall(MatSetType(A12, MATAIJ));
41: PetscCall(MatSeqAIJSetPreallocation(A12, np, NULL));
42: PetscCall(MatMPIAIJSetPreallocation(A12, np, NULL, np, NULL));
44: for (i = 0; i < n; i++) {
45: for (j = 0; j < np; j++) PetscCall(MatSetValue(A12, i, j, i + j * n, INSERT_VALUES));
46: }
47: PetscCall(MatSetValue(A12, 2, 1, 4, INSERT_VALUES));
48: PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
49: PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));
51: /* A21 */
52: PetscCall(MatTranspose(A12, MAT_INITIAL_MATRIX, &A21));
54: A22 = NULL;
56: /* Create block matrix */
57: tmp[0][0] = A11;
58: tmp[0][1] = A12;
59: tmp[1][0] = A21;
60: tmp[1][1] = A22;
62: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A));
63: PetscCall(MatNestSetVecType(A, VECNEST));
64: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
65: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
67: /* Create vectors */
68: PetscCall(MatCreateVecs(A12, &h, &f));
70: PetscCall(VecSet(f, 1.0));
72: /* Create block vector */
73: tmp_x[0] = f;
74: tmp_x[1] = h;
76: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_x, &b));
77: PetscCall(VecAssemblyBegin(b));
78: PetscCall(VecAssemblyEnd(b));
79: PetscCall(VecDuplicate(b, &x));
81: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
82: PetscCall(KSPSetOperators(ksp, A, A));
83: PetscCall(KSPSetType(ksp, KSPGMRES));
84: PetscCall(KSPGetPC(ksp, &pc));
85: PetscCall(PCSetType(pc, PCNONE));
86: PetscCall(KSPSetFromOptions(ksp));
88: PetscCall(KSPSolve(ksp, b, x));
90: PetscCall(VecNestGetSubVecs(x, NULL, &_tmp_x));
92: x1 = _tmp_x[0];
93: x2 = _tmp_x[1];
95: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x1 \n"));
96: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
97: PetscCall(VecView(x1, PETSC_VIEWER_STDOUT_WORLD));
98: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x2 \n"));
99: PetscCall(VecView(x2, PETSC_VIEWER_STDOUT_WORLD));
100: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
102: PetscCall(KSPDestroy(&ksp));
103: PetscCall(VecDestroy(&x));
104: PetscCall(VecDestroy(&b));
105: PetscCall(MatDestroy(&A11));
106: PetscCall(MatDestroy(&A12));
107: PetscCall(MatDestroy(&A21));
108: PetscCall(VecDestroy(&f));
109: PetscCall(VecDestroy(&h));
111: PetscCall(MatDestroy(&A));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: PetscErrorCode test_solve_matgetvecs(void)
116: {
117: Mat A11, A12, A21, A;
118: KSP ksp;
119: PC pc;
120: Vec b, x, f, h, diag, x1, x2;
121: PetscInt n, np, i, j;
122: Mat tmp[2][2];
123: Vec *tmp_x;
125: PetscFunctionBeginUser;
126: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME));
128: n = 3;
129: np = 2;
130: /* Create matrices */
131: /* A11 */
132: PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
133: PetscCall(VecSetSizes(diag, PETSC_DECIDE, n));
134: PetscCall(VecSetFromOptions(diag));
136: PetscCall(VecSet(diag, 1.0 / 10.0)); /* so inverse = diag(10) */
138: /* As a test, create a diagonal matrix for A11 */
139: PetscCall(MatCreate(PETSC_COMM_WORLD, &A11));
140: PetscCall(MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n));
141: PetscCall(MatSetType(A11, MATAIJ));
142: PetscCall(MatSeqAIJSetPreallocation(A11, n, NULL));
143: PetscCall(MatMPIAIJSetPreallocation(A11, np, NULL, np, NULL));
144: PetscCall(MatDiagonalSet(A11, diag, INSERT_VALUES));
146: PetscCall(VecDestroy(&diag));
148: /* A12 */
149: PetscCall(MatCreate(PETSC_COMM_WORLD, &A12));
150: PetscCall(MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np));
151: PetscCall(MatSetType(A12, MATAIJ));
152: PetscCall(MatSeqAIJSetPreallocation(A12, np, NULL));
153: PetscCall(MatMPIAIJSetPreallocation(A12, np, NULL, np, NULL));
155: for (i = 0; i < n; i++) {
156: for (j = 0; j < np; j++) PetscCall(MatSetValue(A12, i, j, i + j * n, INSERT_VALUES));
157: }
158: PetscCall(MatSetValue(A12, 2, 1, 4, INSERT_VALUES));
159: PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
160: PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));
162: /* A21 */
163: PetscCall(MatTranspose(A12, MAT_INITIAL_MATRIX, &A21));
165: /* Create block matrix */
166: tmp[0][0] = A11;
167: tmp[0][1] = A12;
168: tmp[1][0] = A21;
169: tmp[1][1] = NULL;
171: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A));
172: PetscCall(MatNestSetVecType(A, VECNEST));
173: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
174: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
176: /* Create vectors */
177: PetscCall(MatCreateVecs(A, &b, &x));
178: PetscCall(VecNestGetSubVecs(b, NULL, &tmp_x));
179: f = tmp_x[0];
180: h = tmp_x[1];
182: PetscCall(VecSet(f, 1.0));
183: PetscCall(VecSet(h, 0.0));
185: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
186: PetscCall(KSPSetOperators(ksp, A, A));
187: PetscCall(KSPGetPC(ksp, &pc));
188: PetscCall(PCSetType(pc, PCNONE));
189: PetscCall(KSPSetFromOptions(ksp));
191: PetscCall(KSPSolve(ksp, b, x));
192: PetscCall(VecNestGetSubVecs(x, NULL, &tmp_x));
193: x1 = tmp_x[0];
194: x2 = tmp_x[1];
196: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x1 \n"));
197: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
198: PetscCall(VecView(x1, PETSC_VIEWER_STDOUT_WORLD));
199: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x2 \n"));
200: PetscCall(VecView(x2, PETSC_VIEWER_STDOUT_WORLD));
201: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
203: PetscCall(KSPDestroy(&ksp));
204: PetscCall(VecDestroy(&x));
205: PetscCall(VecDestroy(&b));
206: PetscCall(MatDestroy(&A11));
207: PetscCall(MatDestroy(&A12));
208: PetscCall(MatDestroy(&A21));
210: PetscCall(MatDestroy(&A));
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: int main(int argc, char **args)
215: {
216: PetscFunctionBeginUser;
217: PetscCall(PetscInitialize(&argc, &args, NULL, help));
218: PetscCall(test_solve());
219: PetscCall(test_solve_matgetvecs());
220: PetscCall(PetscFinalize());
221: return 0;
222: }
224: /*TEST
226: test:
228: test:
229: suffix: 2
230: nsize: 2
232: test:
233: suffix: 3
234: nsize: 2
235: args: -ksp_monitor -ksp_type bicg
236: requires: !single
238: TEST*/