1 static const char help[] = "Tests MatGetSchurComplement\n";
2
3 #include <petscksp.h>
4
5 /*T
6 Concepts: Mat, Schur Complement
7 T*/
8
Create(MPI_Comm comm,Mat * inA,IS * is0,IS * is1)9 PetscErrorCode Create(MPI_Comm comm,Mat *inA,IS *is0,IS *is1)
10 {
11 PetscErrorCode ierr;
12 Mat A;
13 PetscInt r,rend,M;
14 PetscMPIInt rank;
15
16 PetscFunctionBeginUser;
17 *inA = 0;
18 ierr = MatCreate(comm,&A);CHKERRQ(ierr);
19 ierr = MatSetSizes(A,4,4,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr);
20 ierr = MatSetFromOptions(A);CHKERRQ(ierr);
21 ierr = MatSetUp(A);CHKERRQ(ierr);
22 ierr = MatGetOwnershipRange(A,&r,&rend);CHKERRQ(ierr);
23 ierr = MatGetSize(A,&M,NULL);CHKERRQ(ierr);
24
25 ierr = ISCreateStride(comm,2,r,1,is0);CHKERRQ(ierr);
26 ierr = ISCreateStride(comm,2,r+2,1,is1);CHKERRQ(ierr);
27
28 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
29
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];
33
34 rows[0] = r;
35 rows[1] = r+1;
36 rows[2] = r+2;
37 rows[3] = r+3;
38
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;
44
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;
49
50 cols2[0] = r;
51 cols2[1] = r+2;
52 cols2[2] = (r+4+2)%M;
53
54 cols3[0] = r+1;
55 cols3[1] = r+3;
56 cols3[2] = (r+4+3)%M;
57
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.;
63
64 vals1[0] = RR+6.;
65 vals1[1] = RR+7.;
66 vals1[2] = RR+8.;
67 vals1[3] = RR+9.;
68
69 vals2[0] = RR+10.;
70 vals2[1] = RR+11.;
71 vals2[2] = RR+12.;
72
73 vals3[0] = RR+13.;
74 vals3[1] = RR+14.;
75 vals3[2] = RR+15.;
76 ierr = MatSetValues(A,1,&rows[0],5,cols0,vals0,INSERT_VALUES);CHKERRQ(ierr);
77 ierr = MatSetValues(A,1,&rows[1],4,cols1,vals1,INSERT_VALUES);CHKERRQ(ierr);
78 ierr = MatSetValues(A,1,&rows[2],3,cols2,vals2,INSERT_VALUES);CHKERRQ(ierr);
79 ierr = MatSetValues(A,1,&rows[3],3,cols3,vals3,INSERT_VALUES);CHKERRQ(ierr);
80 }
81 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
82 ierr = MatAssemblyEnd (A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
83 *inA = A;
84 PetscFunctionReturn(0);
85 }
86
Destroy(Mat * A,IS * is0,IS * is1)87 PetscErrorCode Destroy(Mat *A,IS *is0,IS *is1)
88 {
89 PetscErrorCode ierr;
90
91 PetscFunctionBeginUser;
92 ierr = MatDestroy(A);CHKERRQ(ierr);
93 ierr = ISDestroy(is0);CHKERRQ(ierr);
94 ierr = ISDestroy(is1);CHKERRQ(ierr);
95 PetscFunctionReturn(0);
96 }
97
main(int argc,char * argv[])98 int main(int argc,char *argv[])
99 {
100 PetscErrorCode ierr;
101 Mat A,S = NULL,Sexplicit = NULL;
102 MatSchurComplementAinvType ainv_type = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
103 IS is0,is1;
104
105 ierr = PetscInitialize(&argc,&argv,0,help);if (ierr) return ierr;
106 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"ex21","KSP");CHKERRQ(ierr);
107 ierr = 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);CHKERRQ(ierr);
108 ierr = PetscOptionsEnd();CHKERRQ(ierr);
109
110 /* Test the Schur complement one way */
111 ierr = Create(PETSC_COMM_WORLD,&A,&is0,&is1);CHKERRQ(ierr);
112 ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
113 ierr = ISView(is0,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
114 ierr = ISView(is1,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
115 ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_INITIAL_MATRIX,&S,ainv_type,MAT_IGNORE_MATRIX,NULL);CHKERRQ(ierr);
116 ierr = MatSetFromOptions(S);CHKERRQ(ierr);
117 ierr = MatComputeOperator(S,MATAIJ,&Sexplicit);CHKERRQ(ierr);
118 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (0,0) in (1,1)\n");CHKERRQ(ierr);
119 ierr = MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
120 ierr = Destroy(&A,&is0,&is1);CHKERRQ(ierr);
121 ierr = MatDestroy(&S);CHKERRQ(ierr);
122 ierr = MatDestroy(&Sexplicit);CHKERRQ(ierr);
123
124 /* And the other */
125 ierr = Create(PETSC_COMM_WORLD,&A,&is0,&is1);CHKERRQ(ierr);
126 ierr = MatGetSchurComplement(A,is1,is1,is0,is0,MAT_INITIAL_MATRIX,&S,ainv_type,MAT_IGNORE_MATRIX,NULL);CHKERRQ(ierr);
127 ierr = MatSetFromOptions(S);CHKERRQ(ierr);
128 ierr = MatComputeOperator(S,MATAIJ,&Sexplicit);CHKERRQ(ierr);
129 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nExplicit Schur complement of (1,1) in (0,0)\n");CHKERRQ(ierr);
130 ierr = MatView(Sexplicit,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
131 ierr = Destroy(&A,&is0,&is1);CHKERRQ(ierr);
132 ierr = MatDestroy(&S);CHKERRQ(ierr);
133 ierr = MatDestroy(&Sexplicit);CHKERRQ(ierr);
134
135 /* This time just the preconditioning matrix. */
136 ierr = Create(PETSC_COMM_WORLD,&A,&is0,&is1);CHKERRQ(ierr);
137 ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,ainv_type,MAT_INITIAL_MATRIX,&S);CHKERRQ(ierr);
138 ierr = MatSetFromOptions(S);CHKERRQ(ierr);
139 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nPreconditioning Schur complement of (0,0) in (1,1)\n");CHKERRQ(ierr);
140 ierr = MatView(S,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
141 /* Modify and refresh */
142 ierr = MatShift(A,1.);CHKERRQ(ierr);
143 ierr = MatGetSchurComplement(A,is0,is0,is1,is1,MAT_IGNORE_MATRIX,NULL,ainv_type,MAT_REUSE_MATRIX,&S);CHKERRQ(ierr);
144 ierr = PetscPrintf(PETSC_COMM_WORLD,"\nAfter update\n");CHKERRQ(ierr);
145 ierr = MatView(S,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
146 ierr = Destroy(&A,&is0,&is1);CHKERRQ(ierr);
147 ierr = MatDestroy(&S);CHKERRQ(ierr);
148
149 ierr = PetscFinalize();
150 return ierr;
151 }
152
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*/
183