Actual source code: ex22.c
petsc-3.12.2 2019-11-22
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: int n, np, i,j;
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: /* Create vectors */
71: MatCreateVecs(A12, &h, &f);
73: VecSet(f, 1.0);
74: VecSet(h, 0.0);
76: /* Create block vector */
77: tmp_x[0] = f;
78: tmp_x[1] = h;
80: VecCreateNest(PETSC_COMM_WORLD,2,NULL,tmp_x,&b);
81: VecAssemblyBegin(b);
82: VecAssemblyEnd(b);
83: VecDuplicate(b, &x);
85: KSPCreate(PETSC_COMM_WORLD, &ksp);
86: KSPSetOperators(ksp, A, A);
87: KSPSetType(ksp, "gmres");
88: KSPGetPC(ksp, &pc);
89: PCSetType(pc, "none");
90: KSPSetFromOptions(ksp);
92: KSPSolve(ksp, b, x);
94: VecNestGetSubVecs(x,NULL,&_tmp_x);
96: x1 = _tmp_x[0];
97: x2 = _tmp_x[1];
99: PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
100: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
101: VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
102: PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
103: VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
104: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
106: KSPDestroy(&ksp);
107: VecDestroy(&x);
108: VecDestroy(&b);
109: MatDestroy(&A11);
110: MatDestroy(&A12);
111: MatDestroy(&A21);
112: VecDestroy(&f);
113: VecDestroy(&h);
115: MatDestroy(&A);
116: return(0);
117: }
120: PetscErrorCode test_solve_matgetvecs(void)
121: {
122: Mat A11, A12,A21, A;
123: KSP ksp;
124: PC pc;
125: Vec b,x, f,h, diag, x1,x2;
126: int n, np, i,j;
127: Mat tmp[2][2];
128: 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++) {
163: MatSetValue(A12, i,j, (PetscScalar)(i+j*n), INSERT_VALUES);
164: }
165: }
166: MatSetValue(A12, 2,1, (PetscScalar)(4), INSERT_VALUES);
167: MatAssemblyBegin(A12, MAT_FINAL_ASSEMBLY);
168: MatAssemblyEnd(A12, MAT_FINAL_ASSEMBLY);
170: /* A21 */
171: MatTranspose(A12, MAT_INITIAL_MATRIX, &A21);
173: /* Create block matrix */
174: tmp[0][0] = A11;
175: tmp[0][1] = A12;
176: tmp[1][0] = A21;
177: tmp[1][1] = NULL;
179: MatCreateNest(PETSC_COMM_WORLD,2,NULL,2,NULL,&tmp[0][0],&A);
180: MatNestSetVecType(A,VECNEST);
181: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
182: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
184: /* Create vectors */
185: MatCreateVecs(A, &b, &x);
186: VecNestGetSubVecs(b,NULL,&tmp_x);
187: f = tmp_x[0];
188: h = tmp_x[1];
190: VecSet(f, 1.0);
191: VecSet(h, 0.0);
193: KSPCreate(PETSC_COMM_WORLD, &ksp);
194: KSPSetOperators(ksp, A, A);
195: KSPGetPC(ksp, &pc);
196: PCSetType(pc, PCNONE);
197: KSPSetFromOptions(ksp);
199: KSPSolve(ksp, b, x);
200: VecNestGetSubVecs(x,NULL,&tmp_x);
201: x1 = tmp_x[0];
202: x2 = tmp_x[1];
204: PetscPrintf(PETSC_COMM_WORLD, "x1 \n");
205: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO_DETAIL);
206: VecView(x1, PETSC_VIEWER_STDOUT_WORLD);
207: PetscPrintf(PETSC_COMM_WORLD, "x2 \n");
208: VecView(x2, PETSC_VIEWER_STDOUT_WORLD);
209: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
211: KSPDestroy(&ksp);
212: VecDestroy(&x);
213: VecDestroy(&b);
214: MatDestroy(&A11);
215: MatDestroy(&A12);
216: MatDestroy(&A21);
218: MatDestroy(&A);
219: return(0);
220: }
222: int main(int argc, char **args)
223: {
226: PetscInitialize(&argc, &args,(char*)0, help);if (ierr) return ierr;
227: test_solve();
228: test_solve_matgetvecs();
229: PetscFinalize();
230: return ierr;
231: }
233: /*TEST
235: test:
237: test:
238: suffix: 2
239: nsize: 2
241: test:
242: suffix: 3
243: nsize: 2
244: args: -ksp_monitor_short -ksp_type bicg
245: requires: !single
247: TEST*/