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;
17: PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME);
19: n = 3;
20: np = 2;
21: /* Create matrices */
22: /* A11 */
23: VecCreate(PETSC_COMM_WORLD, &diag);
24: VecSetSizes(diag, PETSC_DECIDE, n);
25: VecSetFromOptions(diag);
27: VecSet(diag, (1.0/10.0)); /* so inverse = diag(10) */
29: /* As a test, create a diagonal matrix for A11 */
30: MatCreate(PETSC_COMM_WORLD, &A11);
31: MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n);
32: MatSetType(A11, MATAIJ);
33: MatSeqAIJSetPreallocation(A11, n, NULL);
34: MatMPIAIJSetPreallocation(A11, np, NULL,np, NULL);
35: MatDiagonalSet(A11, diag, INSERT_VALUES);
37: VecDestroy(&diag);
39: /* A12 */
40: MatCreate(PETSC_COMM_WORLD, &A12);
41: MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np);
42: MatSetType(A12, MATAIJ);
43: MatSeqAIJSetPreallocation(A12, np, NULL);
44: MatMPIAIJSetPreallocation(A12, np, NULL,np, NULL);
46: for (i=0; i<n; i++) {
47: for (j=0; j<np; j++) {
48: MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
49: }
50: }
51: MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
52: MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
53: MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);
55: /* A21 */
56: MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);
58: A22 = NULL;
60: /* Create block matrix */
61: tmp[0][0] = A11;
62: tmp[0][1] = A12;
63: tmp[1][0] = A21;
64: tmp[1][1] = A22;
66: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&tmp[0][0],&A);
67: MatNestSetVecType(A,VECNEST);
68: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
71: /* Tests MatMissingDiagonal_Nest */
72: MatMissingDiagonal(A,&flg,NULL);
73: if (!flg) {
74: PetscPrintf(PETSC_COMM_WORLD,"Unexpected %s\n",flg ? "true" : "false");
75: }
77: /* Create vectors */
78: MatCreateVecs(A12, &h, &f);
80: VecSet(f, 1.0);
81: VecSet(h, 0.0);
83: /* Create block vector */
84: tmp_x[0] = f;
85: tmp_x[1] = h;
87: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_x,&b);
88: VecAssemblyBegin(b);
89: VecAssemblyEnd(b);
90: VecDuplicate(b, &x);
92: KSPCreate(PETSC_COMM_WORLD, &ksp);
93: KSPSetOperators(ksp, A, A);
94: KSPSetType(ksp, KSPGMRES);
95: KSPGetPC(ksp, &pc);
96: PCSetType(pc, PCNONE);
97: KSPSetFromOptions(ksp);
99: KSPSolve(ksp, b, x);
101: VecNestGetSubVecs(x,NULL,&_tmp_x);
103: x1 = _tmp_x[0];
104: x2 = _tmp_x[1];
106: PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
107: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
108: VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
109: PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
110: VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
111: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
113: KSPDestroy(&ksp);
114: VecDestroy(&x);
115: VecDestroy(&b);
116: MatDestroy(&A11);
117: MatDestroy(&A12);
118: MatDestroy(&A21);
119: VecDestroy(&f);
120: VecDestroy(&h);
122: MatDestroy(&A);
123: return(0);
124: }
126: PetscErrorCode test_solve_matgetvecs(void)
127: {
128: Mat A11, A12,A21, A;
129: KSP ksp;
130: PC pc;
131: Vec b,x, f,h, diag, x1,x2;
132: PetscInt n, np, i,j;
133: Mat tmp[2][2];
134: Vec *tmp_x;
138: PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME);
140: n = 3;
141: np = 2;
142: /* Create matrices */
143: /* A11 */
144: VecCreate(PETSC_COMM_WORLD, &diag);
145: VecSetSizes(diag, PETSC_DECIDE, n);
146: VecSetFromOptions(diag);
148: VecSet(diag, (1.0/10.0)); /* so inverse = diag(10) */
150: /* As a test, create a diagonal matrix for A11 */
151: MatCreate(PETSC_COMM_WORLD, &A11);
152: MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n);
153: MatSetType(A11, MATAIJ);
154: MatSeqAIJSetPreallocation(A11, n, NULL);
155: MatMPIAIJSetPreallocation(A11, np, NULL,np, NULL);
156: MatDiagonalSet(A11, diag, INSERT_VALUES);
158: VecDestroy(&diag);
160: /* A12 */
161: MatCreate(PETSC_COMM_WORLD, &A12);
162: MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np);
163: MatSetType(A12, MATAIJ);
164: MatSeqAIJSetPreallocation(A12, np, NULL);
165: MatMPIAIJSetPreallocation(A12, np, NULL,np, NULL);
167: for (i=0; i<n; i++) {
168: for (j=0; j<np; j++) {
169: MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
170: }
171: }
172: MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
173: MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
174: MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);
176: /* A21 */
177: MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);
179: /* Create block matrix */
180: tmp[0][0] = A11;
181: tmp[0][1] = A12;
182: tmp[1][0] = A21;
183: tmp[1][1] = NULL;
185: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&tmp[0][0],&A);
186: MatNestSetVecType(A,VECNEST);
187: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
188: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
190: /* Create vectors */
191: MatCreateVecs(A, &b, &x);
192: VecNestGetSubVecs(b,NULL,&tmp_x);
193: f = tmp_x[0];
194: h = tmp_x[1];
196: VecSet(f, 1.0);
197: VecSet(h, 0.0);
199: KSPCreate(PETSC_COMM_WORLD, &ksp);
200: KSPSetOperators(ksp, A, A);
201: KSPGetPC(ksp, &pc);
202: PCSetType(pc, PCNONE);
203: KSPSetFromOptions(ksp);
205: KSPSolve(ksp, b, x);
206: VecNestGetSubVecs(x,NULL,&tmp_x);
207: x1 = tmp_x[0];
208: x2 = tmp_x[1];
210: PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
211: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
212: VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
213: PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
214: VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
215: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
217: KSPDestroy(&ksp);
218: VecDestroy(&x);
219: VecDestroy(&b);
220: MatDestroy(&A11);
221: MatDestroy(&A12);
222: MatDestroy(&A21);
224: MatDestroy(&A);
225: return(0);
226: }
228: int main(int argc, char **args)
229: {
232: PetscInitialize(&argc, &args,(char*)0, help);if (ierr) return ierr;
233: test_solve();
234: test_solve_matgetvecs();
235: PetscFinalize();
236: return ierr;
237: }
239: /*TEST
241: test:
243: test:
244: suffix: 2
245: nsize: 2
247: test:
248: suffix: 3
249: nsize: 2
250: args: -ksp_monitor_short -ksp_type bicg
251: requires: !single
253: TEST*/