Actual source code: ex21.c
petsc-3.14.6 2021-03-30
1: static const char help[] = "Tests MatGetSchurComplement\n";
3: #include <petscksp.h>
5: /*T
6: Concepts: Mat, Schur Complement
7: T*/
9: PetscErrorCode Create(MPI_Comm comm,Mat *inA,IS *is0,IS *is1)
10: {
12: Mat A;
13: PetscInt r,rend,M;
14: PetscMPIInt rank;
17: *inA = 0;
18: MatCreate(comm,&A);
19: MatSetSizes(A,4,4,PETSC_DETERMINE,PETSC_DETERMINE);
20: MatSetFromOptions(A);
21: MatSetUp(A);
22: MatGetOwnershipRange(A,&r,&rend);
23: MatGetSize(A,&M,NULL);
25: ISCreateStride(comm,2,r,1,is0);
26: ISCreateStride(comm,2,r+2,1,is1);
28: MPI_Comm_rank(comm,&rank);
30: {
31: PetscInt rows[4],cols0[5],cols1[5],cols2[3],cols3[3];
32: PetscScalar RR = 1000.*rank,vals0[5],vals1[4],vals2[3],vals3[3];
34: rows[0] = r;
35: rows[1] = r+1;
36: rows[2] = r+2;
37: rows[3] = r+3;
39: cols0[0] = r+0;
40: cols0[1] = r+1;
41: cols0[2] = r+3;
42: cols0[3] = (r+4)%M;
43: cols0[4] = (r+M-4)%M;
45: cols1[0] = r+1;
46: cols1[1] = r+2;
47: cols1[2] = (r+4+1)%M;
48: cols1[3] = (r+M-4+1)%M;
50: cols2[0] = r;
51: cols2[1] = r+2;
52: cols2[2] = (r+4+2)%M;
54: cols3[0] = r+1;
55: cols3[1] = r+3;
56: cols3[2] = (r+4+3)%M;
58: vals0[0] = RR+1.;
59: vals0[1] = RR+2.;
60: vals0[2] = RR+3.;
61: vals0[3] = RR+4.;
62: vals0[4] = RR+5.;
64: vals1[0] = RR+6.;
65: vals1[1] = RR+7.;
66: vals1[2] = RR+8.;
67: vals1[3] = RR+9.;
69: vals2[0] = RR+10.;
70: vals2[1] = RR+11.;
71: vals2[2] = RR+12.;
73: vals3[0] = RR+13.;
74: vals3[1] = RR+14.;
75: vals3[2] = RR+15.;
76: MatSetValues(A,1,&rows[0],5,cols0,vals0,INSERT_VALUES);
77: MatSetValues(A,1,&rows[1],4,cols1,vals1,INSERT_VALUES);
78: MatSetValues(A,1,&rows[2],3,cols2,vals2,INSERT_VALUES);
79: MatSetValues(A,1,&rows[3],3,cols3,vals3,INSERT_VALUES);
80: }
81: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
82: MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);
83: *inA = A;
84: return(0);
85: }
87: PetscErrorCode Destroy(Mat *A,IS *is0,IS *is1)
88: {
92: MatDestroy(A);
93: ISDestroy(is0);
94: ISDestroy(is1);
95: return(0);
96: }
98: int main(int argc,char *argv[])
99: {
101: Mat A,S = NULL,Sexplicit = NULL;
102: MatSchurComplementAinvType ainv_type = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
103: IS is0,is1;
105: PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
106: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21","KSP");
107: PetscOptionsEnum("-mat_schur_complement_ainv_type","Type of approximation for inv(A00) used when assembling Sp = A11 - A10 inv(A00) A01","MatSchurComplementAinvType",MatSchurComplementAinvTypes,(PetscEnum)ainv_type,(PetscEnum*)&ainv_type,NULL);
108: PetscOptionsEnd();
110: /* Test the Schur complement one way */
111: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
112: MatView(A,PETSC_VIEWER_STDOUT_WORLD);
113: ISView(is0,PETSC_VIEWER_STDOUT_WORLD);
114: ISView(is1,PETSC_VIEWER_STDOUT_WORLD);
115: MatGetSchurComplement(A,is0,is0,is1,is1,MAT_INITIAL_MATRIX,&S,ainv_type,MAT_IGNORE_MATRIX,NULL);
116: MatSetFromOptions(S);
117: MatComputeOperator(S,MATAIJ,&Sexplicit);
118: PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (0,0) in (1,1)\n");
119: MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);
120: Destroy(&A,&is0,&is1);
121: MatDestroy(&S);
122: MatDestroy(&Sexplicit);
124: /* And the other */
125: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
126: MatGetSchurComplement(A,is1,is1,is0,is0,MAT_INITIAL_MATRIX,&S,ainv_type,MAT_IGNORE_MATRIX,NULL);
127: MatSetFromOptions(S);
128: MatComputeOperator(S,MATAIJ,&Sexplicit);
129: PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (1,1) in (0,0)\n");
130: MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);
131: Destroy(&A,&is0,&is1);
132: MatDestroy(&S);
133: MatDestroy(&Sexplicit);
135: /* This time just the preconditioning matrix. */
136: Create(PETSC_COMM_WORLD,&A,&is0,&is1);
137: MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,ainv_type,MAT_INITIAL_MATRIX,&S);
138: MatSetFromOptions(S);
139: PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioning Schur complement of (0,0) in (1,1)\n");
140: MatView(S,PETSC_VIEWER_STDOUT_WORLD);
141: /* Modify and refresh */
142: MatShift(A,1.);
143: MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,ainv_type,MAT_REUSE_MATRIX,&S);
144: PetscPrintf(PETSC_COMM_WORLD,"\nAfter update\n");
145: MatView(S,PETSC_VIEWER_STDOUT_WORLD);
146: Destroy(&A,&is0,&is1);
147: MatDestroy(&S);
149: PetscFinalize();
150: return ierr;
151: }
153: /*TEST
154: test:
155: suffix: diag_1
156: args: -mat_schur_complement_ainv_type diag
157: nsize: 1
158: test:
159: suffix: blockdiag_1
160: args: -mat_schur_complement_ainv_type blockdiag
161: nsize: 1
162: test:
163: suffix: diag_2
164: args: -mat_schur_complement_ainv_type diag
165: nsize: 2
166: test:
167: suffix: blockdiag_2
168: args: -mat_schur_complement_ainv_type blockdiag
169: nsize: 2
170: test:
171: # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
172: requires: !single
173: suffix: diag_3
174: args: -mat_schur_complement_ainv_type diag
175: nsize: 3
176: test:
177: # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
178: requires: !single
179: suffix: blockdiag_3
180: args: -mat_schur_complement_ainv_type blockdiag
181: nsize: 3
182: TEST*/