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));
71: PetscCall(VecSet(h, 0.0));
73: /* Create block vector */
74: tmp_x[0] = f;
75: tmp_x[1] = h;
77: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_x, &b));
78: PetscCall(VecAssemblyBegin(b));
79: PetscCall(VecAssemblyEnd(b));
80: PetscCall(VecDuplicate(b, &x));
82: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
83: PetscCall(KSPSetOperators(ksp, A, A));
84: PetscCall(KSPSetType(ksp, KSPGMRES));
85: PetscCall(KSPGetPC(ksp, &pc));
86: PetscCall(PCSetType(pc, PCNONE));
87: PetscCall(KSPSetFromOptions(ksp));
89: PetscCall(KSPSolve(ksp, b, x));
91: PetscCall(VecNestGetSubVecs(x, NULL, &_tmp_x));
93: x1 = _tmp_x[0];
94: x2 = _tmp_x[1];
96: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x1 \n"));
97: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
98: PetscCall(VecView(x1, PETSC_VIEWER_STDOUT_WORLD));
99: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x2 \n"));
100: PetscCall(VecView(x2, PETSC_VIEWER_STDOUT_WORLD));
101: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
103: PetscCall(KSPDestroy(&ksp));
104: PetscCall(VecDestroy(&x));
105: PetscCall(VecDestroy(&b));
106: PetscCall(MatDestroy(&A11));
107: PetscCall(MatDestroy(&A12));
108: PetscCall(MatDestroy(&A21));
109: PetscCall(VecDestroy(&f));
110: PetscCall(VecDestroy(&h));
112: PetscCall(MatDestroy(&A));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: PetscErrorCode test_solve_matgetvecs(void)
117: {
118: Mat A11, A12, A21, A;
119: KSP ksp;
120: PC pc;
121: Vec b, x, f, h, diag, x1, x2;
122: PetscInt n, np, i, j;
123: Mat tmp[2][2];
124: Vec *tmp_x;
126: PetscFunctionBeginUser;
127: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME));
129: n = 3;
130: np = 2;
131: /* Create matrices */
132: /* A11 */
133: PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
134: PetscCall(VecSetSizes(diag, PETSC_DECIDE, n));
135: PetscCall(VecSetFromOptions(diag));
137: PetscCall(VecSet(diag, 1.0 / 10.0)); /* so inverse = diag(10) */
139: /* As a test, create a diagonal matrix for A11 */
140: PetscCall(MatCreate(PETSC_COMM_WORLD, &A11));
141: PetscCall(MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n));
142: PetscCall(MatSetType(A11, MATAIJ));
143: PetscCall(MatSeqAIJSetPreallocation(A11, n, NULL));
144: PetscCall(MatMPIAIJSetPreallocation(A11, np, NULL, np, NULL));
145: PetscCall(MatDiagonalSet(A11, diag, INSERT_VALUES));
147: PetscCall(VecDestroy(&diag));
149: /* A12 */
150: PetscCall(MatCreate(PETSC_COMM_WORLD, &A12));
151: PetscCall(MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np));
152: PetscCall(MatSetType(A12, MATAIJ));
153: PetscCall(MatSeqAIJSetPreallocation(A12, np, NULL));
154: PetscCall(MatMPIAIJSetPreallocation(A12, np, NULL, np, NULL));
156: for (i = 0; i < n; i++) {
157: for (j = 0; j < np; j++) PetscCall(MatSetValue(A12, i, j, i + j * n, INSERT_VALUES));
158: }
159: PetscCall(MatSetValue(A12, 2, 1, 4, INSERT_VALUES));
160: PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
161: PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));
163: /* A21 */
164: PetscCall(MatTranspose(A12, MAT_INITIAL_MATRIX, &A21));
166: /* Create block matrix */
167: tmp[0][0] = A11;
168: tmp[0][1] = A12;
169: tmp[1][0] = A21;
170: tmp[1][1] = NULL;
172: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A));
173: PetscCall(MatNestSetVecType(A, VECNEST));
174: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
175: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
177: /* Create vectors */
178: PetscCall(MatCreateVecs(A, &b, &x));
179: PetscCall(VecNestGetSubVecs(b, NULL, &tmp_x));
180: f = tmp_x[0];
181: h = tmp_x[1];
183: PetscCall(VecSet(f, 1.0));
184: PetscCall(VecSet(h, 0.0));
186: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
187: PetscCall(KSPSetOperators(ksp, A, A));
188: PetscCall(KSPGetPC(ksp, &pc));
189: PetscCall(PCSetType(pc, PCNONE));
190: PetscCall(KSPSetFromOptions(ksp));
192: PetscCall(KSPSolve(ksp, b, x));
193: PetscCall(VecNestGetSubVecs(x, NULL, &tmp_x));
194: x1 = tmp_x[0];
195: x2 = tmp_x[1];
197: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x1 \n"));
198: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
199: PetscCall(VecView(x1, PETSC_VIEWER_STDOUT_WORLD));
200: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x2 \n"));
201: PetscCall(VecView(x2, PETSC_VIEWER_STDOUT_WORLD));
202: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
204: PetscCall(KSPDestroy(&ksp));
205: PetscCall(VecDestroy(&x));
206: PetscCall(VecDestroy(&b));
207: PetscCall(MatDestroy(&A11));
208: PetscCall(MatDestroy(&A12));
209: PetscCall(MatDestroy(&A21));
211: PetscCall(MatDestroy(&A));
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: int main(int argc, char **args)
216: {
217: PetscFunctionBeginUser;
218: PetscCall(PetscInitialize(&argc, &args, NULL, help));
219: PetscCall(test_solve());
220: PetscCall(test_solve_matgetvecs());
221: PetscCall(PetscFinalize());
222: return 0;
223: }
225: /*TEST
227: test:
229: test:
230: suffix: 2
231: nsize: 2
233: test:
234: suffix: 3
235: nsize: 2
236: args: -ksp_monitor_short -ksp_type bicg
237: requires: !single
239: TEST*/