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