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