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;
13: PetscBool flg;
15: PetscFunctionBeginUser;
16: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME));
18: n = 3;
19: np = 2;
20: /* Create matrices */
21: /* A11 */
22: PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
23: PetscCall(VecSetSizes(diag, PETSC_DECIDE, n));
24: PetscCall(VecSetFromOptions(diag));
26: PetscCall(VecSet(diag, (1.0 / 10.0))); /* so inverse = diag(10) */
28: /* As a test, create a diagonal matrix for A11 */
29: PetscCall(MatCreate(PETSC_COMM_WORLD, &A11));
30: PetscCall(MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n));
31: PetscCall(MatSetType(A11, MATAIJ));
32: PetscCall(MatSeqAIJSetPreallocation(A11, n, NULL));
33: PetscCall(MatMPIAIJSetPreallocation(A11, np, NULL, np, NULL));
34: PetscCall(MatDiagonalSet(A11, diag, INSERT_VALUES));
36: PetscCall(VecDestroy(&diag));
38: /* A12 */
39: PetscCall(MatCreate(PETSC_COMM_WORLD, &A12));
40: PetscCall(MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np));
41: PetscCall(MatSetType(A12, MATAIJ));
42: PetscCall(MatSeqAIJSetPreallocation(A12, np, NULL));
43: PetscCall(MatMPIAIJSetPreallocation(A12, np, NULL, np, NULL));
45: for (i = 0; i < n; i++) {
46: for (j = 0; j < np; j++) PetscCall(MatSetValue(A12, i, j, (PetscScalar)(i + j * n), INSERT_VALUES));
47: }
48: PetscCall(MatSetValue(A12, 2, 1, (PetscScalar)(4), INSERT_VALUES));
49: PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
50: PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));
52: /* A21 */
53: PetscCall(MatTranspose(A12, MAT_INITIAL_MATRIX, &A21));
55: A22 = NULL;
57: /* Create block matrix */
58: tmp[0][0] = A11;
59: tmp[0][1] = A12;
60: tmp[1][0] = A21;
61: tmp[1][1] = A22;
63: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A));
64: PetscCall(MatNestSetVecType(A, VECNEST));
65: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
66: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
68: /* Tests MatMissingDiagonal_Nest */
69: PetscCall(MatMissingDiagonal(A, &flg, NULL));
70: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Unexpected %s\n", flg ? "true" : "false"));
72: /* Create vectors */
73: PetscCall(MatCreateVecs(A12, &h, &f));
75: PetscCall(VecSet(f, 1.0));
76: PetscCall(VecSet(h, 0.0));
78: /* Create block vector */
79: tmp_x[0] = f;
80: tmp_x[1] = h;
82: PetscCall(VecCreateNest(PETSC_COMM_WORLD, 2, NULL, tmp_x, &b));
83: PetscCall(VecAssemblyBegin(b));
84: PetscCall(VecAssemblyEnd(b));
85: PetscCall(VecDuplicate(b, &x));
87: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
88: PetscCall(KSPSetOperators(ksp, A, A));
89: PetscCall(KSPSetType(ksp, KSPGMRES));
90: PetscCall(KSPGetPC(ksp, &pc));
91: PetscCall(PCSetType(pc, PCNONE));
92: PetscCall(KSPSetFromOptions(ksp));
94: PetscCall(KSPSolve(ksp, b, x));
96: PetscCall(VecNestGetSubVecs(x, NULL, &_tmp_x));
98: x1 = _tmp_x[0];
99: x2 = _tmp_x[1];
101: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x1 \n"));
102: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
103: PetscCall(VecView(x1, PETSC_VIEWER_STDOUT_WORLD));
104: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x2 \n"));
105: PetscCall(VecView(x2, PETSC_VIEWER_STDOUT_WORLD));
106: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
108: PetscCall(KSPDestroy(&ksp));
109: PetscCall(VecDestroy(&x));
110: PetscCall(VecDestroy(&b));
111: PetscCall(MatDestroy(&A11));
112: PetscCall(MatDestroy(&A12));
113: PetscCall(MatDestroy(&A21));
114: PetscCall(VecDestroy(&f));
115: PetscCall(VecDestroy(&h));
117: PetscCall(MatDestroy(&A));
118: PetscFunctionReturn(PETSC_SUCCESS);
119: }
121: PetscErrorCode test_solve_matgetvecs(void)
122: {
123: Mat A11, A12, A21, A;
124: KSP ksp;
125: PC pc;
126: Vec b, x, f, h, diag, x1, x2;
127: PetscInt n, np, i, j;
128: Mat tmp[2][2];
129: Vec *tmp_x;
131: PetscFunctionBeginUser;
132: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME));
134: n = 3;
135: np = 2;
136: /* Create matrices */
137: /* A11 */
138: PetscCall(VecCreate(PETSC_COMM_WORLD, &diag));
139: PetscCall(VecSetSizes(diag, PETSC_DECIDE, n));
140: PetscCall(VecSetFromOptions(diag));
142: PetscCall(VecSet(diag, (1.0 / 10.0))); /* so inverse = diag(10) */
144: /* As a test, create a diagonal matrix for A11 */
145: PetscCall(MatCreate(PETSC_COMM_WORLD, &A11));
146: PetscCall(MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n));
147: PetscCall(MatSetType(A11, MATAIJ));
148: PetscCall(MatSeqAIJSetPreallocation(A11, n, NULL));
149: PetscCall(MatMPIAIJSetPreallocation(A11, np, NULL, np, NULL));
150: PetscCall(MatDiagonalSet(A11, diag, INSERT_VALUES));
152: PetscCall(VecDestroy(&diag));
154: /* A12 */
155: PetscCall(MatCreate(PETSC_COMM_WORLD, &A12));
156: PetscCall(MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np));
157: PetscCall(MatSetType(A12, MATAIJ));
158: PetscCall(MatSeqAIJSetPreallocation(A12, np, NULL));
159: PetscCall(MatMPIAIJSetPreallocation(A12, np, NULL, np, NULL));
161: for (i = 0; i < n; i++) {
162: for (j = 0; j < np; j++) PetscCall(MatSetValue(A12, i, j, (PetscScalar)(i + j * n), INSERT_VALUES));
163: }
164: PetscCall(MatSetValue(A12, 2, 1, (PetscScalar)(4), INSERT_VALUES));
165: PetscCall(MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY));
166: PetscCall(MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY));
168: /* A21 */
169: PetscCall(MatTranspose(A12, MAT_INITIAL_MATRIX, &A21));
171: /* Create block matrix */
172: tmp[0][0] = A11;
173: tmp[0][1] = A12;
174: tmp[1][0] = A21;
175: tmp[1][1] = NULL;
177: PetscCall(MatCreateNest(PETSC_COMM_WORLD, 2, NULL, 2, NULL, &tmp[0][0], &A));
178: PetscCall(MatNestSetVecType(A, VECNEST));
179: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
180: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
182: /* Create vectors */
183: PetscCall(MatCreateVecs(A, &b, &x));
184: PetscCall(VecNestGetSubVecs(b, NULL, &tmp_x));
185: f = tmp_x[0];
186: h = tmp_x[1];
188: PetscCall(VecSet(f, 1.0));
189: PetscCall(VecSet(h, 0.0));
191: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
192: PetscCall(KSPSetOperators(ksp, A, A));
193: PetscCall(KSPGetPC(ksp, &pc));
194: PetscCall(PCSetType(pc, PCNONE));
195: PetscCall(KSPSetFromOptions(ksp));
197: PetscCall(KSPSolve(ksp, b, x));
198: PetscCall(VecNestGetSubVecs(x, NULL, &tmp_x));
199: x1 = tmp_x[0];
200: x2 = tmp_x[1];
202: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x1 \n"));
203: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL));
204: PetscCall(VecView(x1, PETSC_VIEWER_STDOUT_WORLD));
205: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "x2 \n"));
206: PetscCall(VecView(x2, PETSC_VIEWER_STDOUT_WORLD));
207: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
209: PetscCall(KSPDestroy(&ksp));
210: PetscCall(VecDestroy(&x));
211: PetscCall(VecDestroy(&b));
212: PetscCall(MatDestroy(&A11));
213: PetscCall(MatDestroy(&A12));
214: PetscCall(MatDestroy(&A21));
216: PetscCall(MatDestroy(&A));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: int main(int argc, char **args)
221: {
222: PetscFunctionBeginUser;
223: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
224: PetscCall(test_solve());
225: PetscCall(test_solve_matgetvecs());
226: PetscCall(PetscFinalize());
227: return 0;
228: }
230: /*TEST
232: test:
234: test:
235: suffix: 2
236: nsize: 2
238: test:
239: suffix: 3
240: nsize: 2
241: args: -ksp_monitor_short -ksp_type bicg
242: requires: !single
244: TEST*/