Actual source code: ex22.c
petsc-3.13.4 2020-08-01
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: }
127: PetscErrorCode test_solve_matgetvecs(void)
128: {
129: Mat A11, A12,A21, A;
130: KSP ksp;
131: PC pc;
132: Vec b,x, f,h, diag, x1,x2;
133: PetscInt n, np, i,j;
134: Mat tmp[2][2];
135: Vec *tmp_x;
139: PetscPrintf(PETSC_COMM_WORLD, "%s \n", PETSC_FUNCTION_NAME);
141: n = 3;
142: np = 2;
143: /* Create matrices */
144: /* A11 */
145: VecCreate(PETSC_COMM_WORLD, &diag);
146: VecSetSizes(diag, PETSC_DECIDE, n);
147: VecSetFromOptions(diag);
149: VecSet(diag, (1.0/10.0)); /* so inverse = diag(10) */
151: /* As a test, create a diagonal matrix for A11 */
152: MatCreate(PETSC_COMM_WORLD, &A11);
153: MatSetSizes(A11, PETSC_DECIDE, PETSC_DECIDE, n, n);
154: MatSetType(A11, MATAIJ);
155: MatSeqAIJSetPreallocation(A11, n, NULL);
156: MatMPIAIJSetPreallocation(A11, np, NULL,np, NULL);
157: MatDiagonalSet(A11, diag, INSERT_VALUES);
159: VecDestroy(&diag);
161: /* A12 */
162: MatCreate(PETSC_COMM_WORLD, &A12);
163: MatSetSizes(A12, PETSC_DECIDE, PETSC_DECIDE, n, np);
164: MatSetType(A12, MATAIJ);
165: MatSeqAIJSetPreallocation(A12, np, NULL);
166: MatMPIAIJSetPreallocation(A12, np, NULL,np, NULL);
168: for (i=0; i<n; i++) {
169: for (j=0; j<np; j++) {
170: MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
171: }
172: }
173: MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
174: MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
175: MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);
177: /* A21 */
178: MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);
180: /* Create block matrix */
181: tmp[0][0] = A11;
182: tmp[0][1] = A12;
183: tmp[1][0] = A21;
184: tmp[1][1] = NULL;
186: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&tmp[0][0],&A);
187: MatNestSetVecType(A,VECNEST);
188: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
189: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
191: /* Create vectors */
192: MatCreateVecs(A, &b, &x);
193: VecNestGetSubVecs(b,NULL,&tmp_x);
194: f = tmp_x[0];
195: h = tmp_x[1];
197: VecSet(f, 1.0);
198: VecSet(h, 0.0);
200: KSPCreate(PETSC_COMM_WORLD, &ksp);
201: KSPSetOperators(ksp, A, A);
202: KSPGetPC(ksp, &pc);
203: PCSetType(pc, PCNONE);
204: KSPSetFromOptions(ksp);
206: KSPSolve(ksp, b, x);
207: VecNestGetSubVecs(x,NULL,&tmp_x);
208: x1 = tmp_x[0];
209: x2 = tmp_x[1];
211: PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
212: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
213: VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
214: PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
215: VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
216: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
218: KSPDestroy(&ksp);
219: VecDestroy(&x);
220: VecDestroy(&b);
221: MatDestroy(&A11);
222: MatDestroy(&A12);
223: MatDestroy(&A21);
225: MatDestroy(&A);
226: return(0);
227: }
229: int main(int argc, char **args)
230: {
233: PetscInitialize(&argc, &args,(char*)0, help);if (ierr) return ierr;
234: test_solve();
235: test_solve_matgetvecs();
236: PetscFinalize();
237: return ierr;
238: }
240: /*TEST
242: test:
244: test:
245: suffix: 2
246: nsize: 2
248: test:
249: suffix: 3
250: nsize: 2
251: args: -ksp_monitor_short -ksp_type bicg
252: requires: !single
254: TEST*/