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++) {
47: MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
48: }
49: }
50: MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
51: MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
52: MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);
54: /* A21 */
55: MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);
57: A22 = NULL;
59: /* Create block matrix */
60: tmp[0][0] = A11;
61: tmp[0][1] = A12;
62: tmp[1][0] = A21;
63: tmp[1][1] = A22;
65: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&tmp[0][0],&A);
66: MatNestSetVecType(A,VECNEST);
67: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
70: /* Tests MatMissingDiagonal_Nest */
71: MatMissingDiagonal(A,&flg,NULL);
72: if (!flg) {
73: PetscPrintf(PETSC_COMM_WORLD,"Unexpected %s\n",flg ? "true" : "false");
74: }
76: /* Create vectors */
77: MatCreateVecs(A12, &h, &f);
79: VecSet(f, 1.0);
80: VecSet(h, 0.0);
82: /* Create block vector */
83: tmp_x[0] = f;
84: tmp_x[1] = h;
86: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_x,&b);
87: VecAssemblyBegin(b);
88: VecAssemblyEnd(b);
89: VecDuplicate(b, &x);
91: KSPCreate(PETSC_COMM_WORLD, &ksp);
92: KSPSetOperators(ksp, A, A);
93: KSPSetType(ksp, KSPGMRES);
94: KSPGetPC(ksp, &pc);
95: PCSetType(pc, PCNONE);
96: KSPSetFromOptions(ksp);
98: KSPSolve(ksp, b, x);
100: VecNestGetSubVecs(x,NULL,&_tmp_x);
102: x1 = _tmp_x[0];
103: x2 = _tmp_x[1];
105: PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
106: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
107: VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
108: PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
109: VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
110: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
112: KSPDestroy(&ksp);
113: VecDestroy(&x);
114: VecDestroy(&b);
115: MatDestroy(&A11);
116: MatDestroy(&A12);
117: MatDestroy(&A21);
118: VecDestroy(&f);
119: VecDestroy(&h);
121: MatDestroy(&A);
122: return 0;
123: }
125: PetscErrorCode test_solve_matgetvecs(void)
126: {
127: Mat A11, A12,A21, A;
128: KSP ksp;
129: PC pc;
130: Vec b,x, f,h, diag, x1,x2;
131: PetscInt n, np, i,j;
132: Mat tmp[2][2];
133: Vec *tmp_x;
136: PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME);
138: n = 3;
139: np = 2;
140: /* Create matrices */
141: /* A11 */
142: VecCreate(PETSC_COMM_WORLD, &diag);
143: VecSetSizes(diag, PETSC_DECIDE, n);
144: VecSetFromOptions(diag);
146: VecSet(diag, (1.0/10.0)); /* so inverse = diag(10) */
148: /* As a test, create a diagonal matrix for A11 */
149: MatCreate(PETSC_COMM_WORLD, &A11);
150: MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n);
151: MatSetType(A11, MATAIJ);
152: MatSeqAIJSetPreallocation(A11, n, NULL);
153: MatMPIAIJSetPreallocation(A11, np, NULL,np, NULL);
154: MatDiagonalSet(A11, diag, INSERT_VALUES);
156: VecDestroy(&diag);
158: /* A12 */
159: MatCreate(PETSC_COMM_WORLD, &A12);
160: MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np);
161: MatSetType(A12, MATAIJ);
162: MatSeqAIJSetPreallocation(A12, np, NULL);
163: MatMPIAIJSetPreallocation(A12, np, NULL,np, NULL);
165: for (i=0; i<n; i++) {
166: for (j=0; j<np; j++) {
167: MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
168: }
169: }
170: MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
171: MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
172: MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);
174: /* A21 */
175: MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);
177: /* Create block matrix */
178: tmp[0][0] = A11;
179: tmp[0][1] = A12;
180: tmp[1][0] = A21;
181: tmp[1][1] = NULL;
183: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&tmp[0][0],&A);
184: MatNestSetVecType(A,VECNEST);
185: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
186: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
188: /* Create vectors */
189: MatCreateVecs(A, &b, &x);
190: VecNestGetSubVecs(b,NULL,&tmp_x);
191: f = tmp_x[0];
192: h = tmp_x[1];
194: VecSet(f, 1.0);
195: VecSet(h, 0.0);
197: KSPCreate(PETSC_COMM_WORLD, &ksp);
198: KSPSetOperators(ksp, A, A);
199: KSPGetPC(ksp, &pc);
200: PCSetType(pc, PCNONE);
201: KSPSetFromOptions(ksp);
203: KSPSolve(ksp, b, x);
204: VecNestGetSubVecs(x,NULL,&tmp_x);
205: x1 = tmp_x[0];
206: x2 = tmp_x[1];
208: PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
209: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
210: VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
211: PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
212: VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
213: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
215: KSPDestroy(&ksp);
216: VecDestroy(&x);
217: VecDestroy(&b);
218: MatDestroy(&A11);
219: MatDestroy(&A12);
220: MatDestroy(&A21);
222: MatDestroy(&A);
223: return 0;
224: }
226: int main(int argc, char **args)
227: {
229: PetscInitialize(&argc, &args,(char*)0, help);
230: test_solve();
231: test_solve_matgetvecs();
232: PetscFinalize();
233: return 0;
234: }
236: /*TEST
238: test:
240: test:
241: suffix: 2
242: nsize: 2
244: test:
245: suffix: 3
246: nsize: 2
247: args: -ksp_monitor_short -ksp_type bicg
248: requires: !single
250: TEST*/