1 static const char help[] = "Tests MatGetSchurComplement\n";
2
3 #include <petscksp.h>
4
MatNormDifference(Mat A,Mat B,PetscReal * norm)5 PetscErrorCode MatNormDifference(Mat A, Mat B, PetscReal *norm)
6 {
7 PetscReal bnorm;
8
9 PetscFunctionBegin;
10 PetscCall(MatDuplicate(B, MAT_COPY_VALUES, &B));
11 PetscCall(MatNorm(B, NORM_FROBENIUS, &bnorm));
12 PetscCall(MatAXPY(B, -1.0, A, DIFFERENT_NONZERO_PATTERN));
13 PetscCall(MatNorm(B, NORM_FROBENIUS, norm));
14 PetscCall(MatDestroy(&B));
15 *norm = *norm / bnorm;
16 PetscFunctionReturn(PETSC_SUCCESS);
17 }
18
Create(MPI_Comm comm,Mat * inA,IS * is0,IS * is1)19 PetscErrorCode Create(MPI_Comm comm, Mat *inA, IS *is0, IS *is1)
20 {
21 Mat A;
22 PetscInt r, rend, M;
23 PetscMPIInt rank;
24
25 PetscFunctionBeginUser;
26 *inA = 0;
27 PetscCall(MatCreate(comm, &A));
28 PetscCall(MatSetSizes(A, 4, 4, PETSC_DETERMINE, PETSC_DETERMINE));
29 PetscCall(MatSetFromOptions(A));
30 PetscCall(MatSetUp(A));
31 PetscCall(MatGetOwnershipRange(A, &r, &rend));
32 PetscCall(MatGetSize(A, &M, NULL));
33
34 PetscCall(ISCreateStride(comm, 2, r, 1, is0));
35 PetscCall(ISCreateStride(comm, 2, r + 2, 1, is1));
36
37 PetscCallMPI(MPI_Comm_rank(comm, &rank));
38
39 {
40 PetscInt rows[4], cols0[5], cols1[5], cols2[3], cols3[3];
41 PetscScalar RR = 1000. * rank, vals0[5], vals1[4], vals2[3], vals3[3];
42
43 rows[0] = r;
44 rows[1] = r + 1;
45 rows[2] = r + 2;
46 rows[3] = r + 3;
47
48 cols0[0] = r + 0;
49 cols0[1] = r + 1;
50 cols0[2] = r + 3;
51 cols0[3] = (r + 4) % M;
52 cols0[4] = (r + M - 4) % M;
53
54 cols1[0] = r + 1;
55 cols1[1] = r + 2;
56 cols1[2] = (r + 4 + 1) % M;
57 cols1[3] = (r + M - 4 + 1) % M;
58
59 cols2[0] = r;
60 cols2[1] = r + 2;
61 cols2[2] = (r + 4 + 2) % M;
62
63 cols3[0] = r + 1;
64 cols3[1] = r + 3;
65 cols3[2] = (r + 4 + 3) % M;
66
67 vals0[0] = RR + 1.;
68 vals0[1] = RR + 2.;
69 vals0[2] = RR + 3.;
70 vals0[3] = RR + 4.;
71 vals0[4] = RR + 5.;
72
73 vals1[0] = RR + 6.;
74 vals1[1] = RR + 7.;
75 vals1[2] = RR + 8.;
76 vals1[3] = RR + 9.;
77
78 vals2[0] = RR + 10.;
79 vals2[1] = RR + 11.;
80 vals2[2] = RR + 12.;
81
82 vals3[0] = RR + 13.;
83 vals3[1] = RR + 14.;
84 vals3[2] = RR + 15.;
85 PetscCall(MatSetValues(A, 1, &rows[0], 5, cols0, vals0, INSERT_VALUES));
86 PetscCall(MatSetValues(A, 1, &rows[1], 4, cols1, vals1, INSERT_VALUES));
87 PetscCall(MatSetValues(A, 1, &rows[2], 3, cols2, vals2, INSERT_VALUES));
88 PetscCall(MatSetValues(A, 1, &rows[3], 3, cols3, vals3, INSERT_VALUES));
89 }
90 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
91 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
92 *inA = A;
93 PetscFunctionReturn(PETSC_SUCCESS);
94 }
95
Destroy(Mat * A,IS * is0,IS * is1)96 PetscErrorCode Destroy(Mat *A, IS *is0, IS *is1)
97 {
98 PetscFunctionBeginUser;
99 PetscCall(MatDestroy(A));
100 PetscCall(ISDestroy(is0));
101 PetscCall(ISDestroy(is1));
102 PetscFunctionReturn(PETSC_SUCCESS);
103 }
104
main(int argc,char * argv[])105 int main(int argc, char *argv[])
106 {
107 Mat A, S = NULL, Sexplicit = NULL, Sp, B, C;
108 MatSchurComplementAinvType ainv_type = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
109 IS is0, is1;
110 PetscBool flg;
111 PetscInt m, N = 10, M;
112
113 PetscFunctionBeginUser;
114 PetscCall(PetscInitialize(&argc, &argv, 0, help));
115 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex21", "KSP");
116 PetscCall(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));
117 PetscOptionsEnd();
118
119 /* Test the Schur complement one way */
120 PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
121 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD));
122 PetscCall(ISView(is0, PETSC_VIEWER_STDOUT_WORLD));
123 PetscCall(ISView(is1, PETSC_VIEWER_STDOUT_WORLD));
124 PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_INITIAL_MATRIX, &S, ainv_type, MAT_IGNORE_MATRIX, NULL));
125 PetscCall(MatSetFromOptions(S));
126 PetscCall(MatComputeOperator(S, MATAIJ, &Sexplicit));
127 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nExplicit Schur complement of (0,0) in (1,1)\n"));
128 PetscCall(MatView(Sexplicit, PETSC_VIEWER_STDOUT_WORLD));
129 if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
130 PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_FULL));
131 PetscCall(MatSchurComplementGetPmat(S, MAT_INITIAL_MATRIX, &Sp));
132 PetscCall(MatMultEqual(Sp, Sexplicit, 10, &flg));
133 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Sp != S");
134 PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_DIAG));
135 PetscCall(MatDestroy(&Sp));
136 }
137 PetscCall(Destroy(&A, &is0, &is1));
138 if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
139 PetscCall(MatGetLocalSize(Sexplicit, &m, NULL));
140 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)Sexplicit), m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
141 PetscCall(MatSetRandom(B, NULL));
142 PetscCall(MatMatMult(S, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
143 PetscCall(MatMatMultEqual(Sexplicit, B, C, 10, &flg));
144 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "S*B != C");
145 PetscCall(MatDestroy(&C));
146 PetscCall(MatDestroy(&B));
147 }
148 PetscCall(MatDestroy(&S));
149 PetscCall(MatDestroy(&Sexplicit));
150
151 /* And the other */
152 PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
153 PetscCall(MatGetSchurComplement(A, is1, is1, is0, is0, MAT_INITIAL_MATRIX, &S, ainv_type, MAT_IGNORE_MATRIX, NULL));
154 PetscCall(MatSetFromOptions(S));
155 PetscCall(MatComputeOperator(S, MATAIJ, &Sexplicit));
156 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nExplicit Schur complement of (1,1) in (0,0)\n"));
157 PetscCall(MatView(Sexplicit, PETSC_VIEWER_STDOUT_WORLD));
158
159 /* Test Mat-Mat operations with B AIJ */
160 {
161 Mat B, C, Ce, Cee, Cer;
162 PetscReal err, tol = 10 * PETSC_SMALL;
163 PetscErrorCode (*funcs[])(Mat, Mat, MatReuse, PetscReal, Mat *) = {MatMatMult, MatMatTransposeMult, MatPtAP, MatRARt};
164 const char *names[] = {"MatMatMult", "MatMatTransposeMult", "MatPtAP", "MatRARt"};
165 PetscBool browsacols[] = {PETSC_TRUE, PETSC_FALSE, PETSC_TRUE, PETSC_FALSE};
166
167 for (PetscInt i = 0; i < 4; i++) {
168 PetscCall(MatGetLocalSize(S, NULL, &m));
169 PetscCall(MatGetSize(S, NULL, &M));
170 PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
171 if (browsacols[i]) PetscCall(MatSetSizes(B, m, PETSC_DECIDE, M, 11));
172 else PetscCall(MatSetSizes(B, PETSC_DECIDE, m, 11, M));
173 PetscCall(MatSetType(B, MATAIJ));
174 PetscCall(MatSeqAIJSetPreallocation(B, PETSC_DEFAULT, NULL));
175 PetscCall(MatMPIAIJSetPreallocation(B, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
176 PetscCall(MatSetUp(B));
177 PetscCall(MatSetRandom(B, NULL));
178 PetscCall((*funcs[i])(S, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
179 PetscCall(MatComputeOperator(C, MATAIJ, &Ce));
180 PetscCall(MatMatMult(S, B, MAT_REUSE_MATRIX, PETSC_DETERMINE, &C));
181 PetscCall(MatComputeOperator(C, MATAIJ, &Cer));
182 PetscCall((*funcs[i])(Sexplicit, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &Cee));
183 PetscCall(MatNormDifference(Ce, Cee, &err));
184 PetscCheck(err <= tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in initial %s(): %g", names[i], (double)err);
185 PetscCall(MatNormDifference(Cer, Cee, &err));
186 PetscCheck(err <= tol, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Error in reuse %s(): %g", names[i], (double)err);
187 PetscCall(MatDestroy(&C));
188 PetscCall(MatDestroy(&Ce));
189 PetscCall(MatDestroy(&Cer));
190 PetscCall(MatDestroy(&Cee));
191 PetscCall(MatDestroy(&B));
192 }
193 }
194 if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
195 PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_FULL));
196 PetscCall(MatSchurComplementGetPmat(S, MAT_INITIAL_MATRIX, &Sp));
197 PetscCall(MatMultEqual(Sp, Sexplicit, 10, &flg));
198 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Sp != S");
199 PetscCall(MatSchurComplementSetAinvType(S, MAT_SCHUR_COMPLEMENT_AINV_DIAG));
200 PetscCall(MatDestroy(&Sp));
201 }
202 PetscCall(Destroy(&A, &is0, &is1));
203 if (ainv_type == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
204 PetscCall(MatGetLocalSize(Sexplicit, &m, NULL));
205 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)Sexplicit), m, PETSC_DECIDE, PETSC_DECIDE, N, NULL, &B));
206 PetscCall(MatSetRandom(B, NULL));
207 PetscCall(MatMatMult(S, B, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C));
208 PetscCall(MatMatMultEqual(Sexplicit, B, C, 10, &flg));
209 PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "S*B != C");
210 PetscCall(MatDestroy(&C));
211 PetscCall(MatDestroy(&B));
212 }
213 PetscCall(MatDestroy(&S));
214 PetscCall(MatDestroy(&Sexplicit));
215
216 /* This time just the matrix used to construct the preconditioner. */
217 PetscCall(Create(PETSC_COMM_WORLD, &A, &is0, &is1));
218 PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_IGNORE_MATRIX, NULL, ainv_type, MAT_INITIAL_MATRIX, &S));
219 PetscCall(MatSetFromOptions(S));
220 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nPreconditioning Schur complement of (0,0) in (1,1)\n"));
221 PetscCall(MatView(S, PETSC_VIEWER_STDOUT_WORLD));
222 /* Modify and refresh */
223 PetscCall(MatShift(A, 1.));
224 PetscCall(MatGetSchurComplement(A, is0, is0, is1, is1, MAT_IGNORE_MATRIX, NULL, ainv_type, MAT_REUSE_MATRIX, &S));
225 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nAfter update\n"));
226 PetscCall(MatView(S, PETSC_VIEWER_STDOUT_WORLD));
227 PetscCall(Destroy(&A, &is0, &is1));
228 PetscCall(MatDestroy(&S));
229
230 PetscCall(PetscFinalize());
231 return 0;
232 }
233
234 /*TEST
235 test:
236 suffix: diag_1
237 args: -mat_schur_complement_ainv_type diag
238 nsize: 1
239 test:
240 suffix: blockdiag_1
241 args: -mat_schur_complement_ainv_type blockdiag
242 nsize: 1
243 test:
244 suffix: diag_2
245 args: -mat_schur_complement_ainv_type diag
246 nsize: 2
247 test:
248 suffix: blockdiag_2
249 args: -mat_schur_complement_ainv_type blockdiag
250 nsize: 2
251 test:
252 # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
253 requires: !single
254 suffix: diag_3
255 args: -mat_schur_complement_ainv_type diag -ksp_rtol 1e-12
256 nsize: 3
257 test:
258 # does not work with single because residual norm computed by GMRES recurrence formula becomes invalid
259 requires: !single
260 suffix: blockdiag_3
261 args: -mat_schur_complement_ainv_type blockdiag
262 nsize: 3
263 TEST*/
264