1 #include <../src/ksp/ksp/utils/schurm/schurm.h> /*I "petscksp.h" I*/
2
3 const char *const MatSchurComplementAinvTypes[] = {"DIAG","LUMP","BLOCKDIAG","MatSchurComplementAinvType","MAT_SCHUR_COMPLEMENT_AINV_",NULL};
4
MatCreateVecs_SchurComplement(Mat N,Vec * right,Vec * left)5 PetscErrorCode MatCreateVecs_SchurComplement(Mat N,Vec *right,Vec *left)
6 {
7 Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
8 PetscErrorCode ierr;
9
10 PetscFunctionBegin;
11 if (Na->D) {
12 ierr = MatCreateVecs(Na->D,right,left);CHKERRQ(ierr);
13 PetscFunctionReturn(0);
14 }
15 if (right) {
16 ierr = MatCreateVecs(Na->B,right,NULL);CHKERRQ(ierr);
17 }
18 if (left) {
19 ierr = MatCreateVecs(Na->C,NULL,left);CHKERRQ(ierr);
20 }
21 PetscFunctionReturn(0);
22 }
23
MatView_SchurComplement(Mat N,PetscViewer viewer)24 PetscErrorCode MatView_SchurComplement(Mat N,PetscViewer viewer)
25 {
26 Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
27 PetscErrorCode ierr;
28
29 PetscFunctionBegin;
30 ierr = PetscViewerASCIIPrintf(viewer,"Schur complement A11 - A10 inv(A00) A01\n");CHKERRQ(ierr);
31 if (Na->D) {
32 ierr = PetscViewerASCIIPrintf(viewer,"A11\n");CHKERRQ(ierr);
33 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
34 ierr = MatView(Na->D,viewer);CHKERRQ(ierr);
35 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
36 } else {
37 ierr = PetscViewerASCIIPrintf(viewer,"A11 = 0\n");CHKERRQ(ierr);
38 }
39 ierr = PetscViewerASCIIPrintf(viewer,"A10\n");CHKERRQ(ierr);
40 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
41 ierr = MatView(Na->C,viewer);CHKERRQ(ierr);
42 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
43 ierr = PetscViewerASCIIPrintf(viewer,"KSP of A00\n");CHKERRQ(ierr);
44 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
45 ierr = KSPView(Na->ksp,viewer);CHKERRQ(ierr);
46 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
47 ierr = PetscViewerASCIIPrintf(viewer,"A01\n");CHKERRQ(ierr);
48 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
49 ierr = MatView(Na->B,viewer);CHKERRQ(ierr);
50 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
51 PetscFunctionReturn(0);
52 }
53
54 /*
55 A11^T - A01^T ksptrans(A00,Ap00) A10^T
56 */
MatMultTranspose_SchurComplement(Mat N,Vec x,Vec y)57 PetscErrorCode MatMultTranspose_SchurComplement(Mat N,Vec x,Vec y)
58 {
59 Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
60 PetscErrorCode ierr;
61
62 PetscFunctionBegin;
63 if (!Na->work1) {ierr = MatCreateVecs(Na->A,&Na->work1,NULL);CHKERRQ(ierr);}
64 if (!Na->work2) {ierr = MatCreateVecs(Na->A,&Na->work2,NULL);CHKERRQ(ierr);}
65 ierr = MatMultTranspose(Na->C,x,Na->work1);CHKERRQ(ierr);
66 ierr = KSPSolveTranspose(Na->ksp,Na->work1,Na->work2);CHKERRQ(ierr);
67 ierr = MatMultTranspose(Na->B,Na->work2,y);CHKERRQ(ierr);
68 ierr = VecScale(y,-1.0);CHKERRQ(ierr);
69 if (Na->D) {
70 ierr = MatMultTransposeAdd(Na->D,x,y,y);CHKERRQ(ierr);
71 }
72 PetscFunctionReturn(0);
73 }
74
75 /*
76 A11 - A10 ksp(A00,Ap00) A01
77 */
MatMult_SchurComplement(Mat N,Vec x,Vec y)78 PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
79 {
80 Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
81 PetscErrorCode ierr;
82
83 PetscFunctionBegin;
84 if (!Na->work1) {ierr = MatCreateVecs(Na->A,&Na->work1,NULL);CHKERRQ(ierr);}
85 if (!Na->work2) {ierr = MatCreateVecs(Na->A,&Na->work2,NULL);CHKERRQ(ierr);}
86 ierr = MatMult(Na->B,x,Na->work1);CHKERRQ(ierr);
87 ierr = KSPSolve(Na->ksp,Na->work1,Na->work2);CHKERRQ(ierr);
88 ierr = MatMult(Na->C,Na->work2,y);CHKERRQ(ierr);
89 ierr = VecScale(y,-1.0);CHKERRQ(ierr);
90 if (Na->D) {
91 ierr = MatMultAdd(Na->D,x,y,y);CHKERRQ(ierr);
92 }
93 PetscFunctionReturn(0);
94 }
95
96 /*
97 A11 - A10 ksp(A00,Ap00) A01
98 */
MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)99 PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
100 {
101 Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
102 PetscErrorCode ierr;
103
104 PetscFunctionBegin;
105 if (!Na->work1) {ierr = MatCreateVecs(Na->A,&Na->work1,NULL);CHKERRQ(ierr);}
106 if (!Na->work2) {ierr = MatCreateVecs(Na->A,&Na->work2,NULL);CHKERRQ(ierr);}
107 ierr = MatMult(Na->B,x,Na->work1);CHKERRQ(ierr);
108 ierr = KSPSolve(Na->ksp,Na->work1,Na->work2);CHKERRQ(ierr);
109 if (y == z) {
110 ierr = VecScale(Na->work2,-1.0);CHKERRQ(ierr);
111 ierr = MatMultAdd(Na->C,Na->work2,z,z);CHKERRQ(ierr);
112 } else {
113 ierr = MatMult(Na->C,Na->work2,z);CHKERRQ(ierr);
114 ierr = VecAYPX(z,-1.0,y);CHKERRQ(ierr);
115 }
116 if (Na->D) {
117 ierr = MatMultAdd(Na->D,x,z,z);CHKERRQ(ierr);
118 }
119 PetscFunctionReturn(0);
120 }
121
MatSetFromOptions_SchurComplement(PetscOptionItems * PetscOptionsObject,Mat N)122 PetscErrorCode MatSetFromOptions_SchurComplement(PetscOptionItems *PetscOptionsObject,Mat N)
123 {
124 Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
125 PetscErrorCode ierr;
126
127 PetscFunctionBegin;
128 ierr = PetscOptionsHead(PetscOptionsObject,"MatSchurComplementOptions");CHKERRQ(ierr);
129 Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
130 ierr = PetscOptionsEnum("-mat_schur_complement_ainv_type","Type of approximation for inv(A00) used when assembling Sp = A11 - A10 inv(A00) A01","MatSchurComplementSetAinvType",MatSchurComplementAinvTypes,(PetscEnum)Na->ainvtype,(PetscEnum*)&Na->ainvtype,NULL);CHKERRQ(ierr);
131 ierr = PetscOptionsTail();CHKERRQ(ierr);
132 ierr = KSPSetFromOptions(Na->ksp);CHKERRQ(ierr);
133 PetscFunctionReturn(0);
134 }
135
MatDestroy_SchurComplement(Mat N)136 PetscErrorCode MatDestroy_SchurComplement(Mat N)
137 {
138 Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
139 PetscErrorCode ierr;
140
141 PetscFunctionBegin;
142 ierr = MatDestroy(&Na->A);CHKERRQ(ierr);
143 ierr = MatDestroy(&Na->Ap);CHKERRQ(ierr);
144 ierr = MatDestroy(&Na->B);CHKERRQ(ierr);
145 ierr = MatDestroy(&Na->C);CHKERRQ(ierr);
146 ierr = MatDestroy(&Na->D);CHKERRQ(ierr);
147 ierr = VecDestroy(&Na->work1);CHKERRQ(ierr);
148 ierr = VecDestroy(&Na->work2);CHKERRQ(ierr);
149 ierr = KSPDestroy(&Na->ksp);CHKERRQ(ierr);
150 ierr = PetscFree(N->data);CHKERRQ(ierr);
151 PetscFunctionReturn(0);
152 }
153
154 /*@
155 MatCreateSchurComplement - Creates a new matrix object that behaves like the Schur complement of a matrix
156
157 Collective on A00
158
159 Input Parameters:
160 + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
161 - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
162
163 Output Parameter:
164 . S - the matrix that the Schur complement S = A11 - A10 ksp(A00,Ap00) A01
165
166 Level: intermediate
167
168 Notes:
169 The Schur complement is NOT actually formed! Rather, this
170 object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
171 for Schur complement S and a KSP solver to approximate the action of A^{-1}.
172
173 All four matrices must have the same MPI communicator.
174
175 A00 and A11 must be square matrices.
176
177 MatGetSchurComplement() takes as arguments the index sets for the submatrices and returns both the virtual Schur complement (what this returns) plus
178 a sparse approximation to the true Schur complement (useful for building a preconditioner for the Schur complement).
179
180 MatSchurComplementGetPmat() can be called on the output of this function to generate an explicit approximation to the Schur complement.
181
182 Developer Notes:
183 The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
184 remove redundancy and be clearer and simplier.
185
186
187 .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement(),
188 MatSchurComplementGetPmat()
189
190 @*/
MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat * S)191 PetscErrorCode MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
192 {
193 PetscErrorCode ierr;
194
195 PetscFunctionBegin;
196 ierr = KSPInitializePackage();CHKERRQ(ierr);
197 ierr = MatCreate(PetscObjectComm((PetscObject)A00),S);CHKERRQ(ierr);
198 ierr = MatSetType(*S,MATSCHURCOMPLEMENT);CHKERRQ(ierr);
199 ierr = MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);CHKERRQ(ierr);
200 PetscFunctionReturn(0);
201 }
202
203 /*@
204 MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
205
206 Collective on S
207
208 Input Parameter:
209 + S - matrix obtained with MatCreateSchurComplement (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
210 . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
211 - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
212
213 Level: intermediate
214
215 Notes:
216 The Schur complement is NOT actually formed! Rather, this
217 object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
218 for Schur complement S and a KSP solver to approximate the action of A^{-1}.
219
220 All four matrices must have the same MPI communicator.
221
222 A00 and A11 must be square matrices.
223
224 .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatCreateSchurComplement(), MatGetSchurComplement()
225
226 @*/
MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)227 PetscErrorCode MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
228 {
229 PetscErrorCode ierr;
230 Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
231 PetscBool isschur;
232
233 PetscFunctionBegin;
234 ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
235 if (!isschur) PetscFunctionReturn(0);
236 if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
237 PetscValidHeaderSpecific(A00,MAT_CLASSID,2);
238 PetscValidHeaderSpecific(Ap00,MAT_CLASSID,3);
239 PetscValidHeaderSpecific(A01,MAT_CLASSID,4);
240 PetscValidHeaderSpecific(A10,MAT_CLASSID,5);
241 PetscCheckSameComm(A00,2,Ap00,3);
242 PetscCheckSameComm(A00,2,A01,4);
243 PetscCheckSameComm(A00,2,A10,5);
244 if (A00->rmap->n != A00->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A00 %D do not equal local columns %D",A00->rmap->n,A00->cmap->n);
245 if (A00->rmap->n != Ap00->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A00 %D do not equal local rows of Ap00 %D",A00->rmap->n,Ap00->rmap->n);
246 if (Ap00->rmap->n != Ap00->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of Ap00 %D do not equal local columns %D",Ap00->rmap->n,Ap00->cmap->n);
247 if (A00->cmap->n != A01->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A00 %D do not equal local rows of A01 %D",A00->cmap->n,A01->rmap->n);
248 if (A10->cmap->n != A00->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A10 %D do not equal local rows of A00 %D",A10->cmap->n,A00->rmap->n);
249 if (A11) {
250 PetscValidHeaderSpecific(A11,MAT_CLASSID,6);
251 PetscCheckSameComm(A00,2,A11,6);
252 if (A10->rmap->n != A11->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A10 %D do not equal local rows A11 %D",A10->rmap->n,A11->rmap->n);
253 }
254
255 ierr = MatSetSizes(S,A10->rmap->n,A01->cmap->n,A10->rmap->N,A01->cmap->N);CHKERRQ(ierr);
256 ierr = PetscObjectReference((PetscObject)A00);CHKERRQ(ierr);
257 ierr = PetscObjectReference((PetscObject)Ap00);CHKERRQ(ierr);
258 ierr = PetscObjectReference((PetscObject)A01);CHKERRQ(ierr);
259 ierr = PetscObjectReference((PetscObject)A10);CHKERRQ(ierr);
260 Na->A = A00;
261 Na->Ap = Ap00;
262 Na->B = A01;
263 Na->C = A10;
264 Na->D = A11;
265 if (A11) {
266 ierr = PetscObjectReference((PetscObject)A11);CHKERRQ(ierr);
267 }
268 ierr = MatSetUp(S);CHKERRQ(ierr);
269 ierr = KSPSetOperators(Na->ksp,A00,Ap00);CHKERRQ(ierr);
270 S->assembled = PETSC_TRUE;
271 PetscFunctionReturn(0);
272 }
273
274 /*@
275 MatSchurComplementGetKSP - Gets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
276
277 Not Collective
278
279 Input Parameter:
280 . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
281
282 Output Parameter:
283 . ksp - the linear solver object
284
285 Options Database:
286 . -fieldsplit_<splitname_0>_XXX sets KSP and PC options for the 0-split solver inside the Schur complement used in PCFieldSplit; default <splitname_0> is 0.
287
288 Level: intermediate
289
290 .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
291 @*/
MatSchurComplementGetKSP(Mat S,KSP * ksp)292 PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
293 {
294 Mat_SchurComplement *Na;
295 PetscBool isschur;
296 PetscErrorCode ierr;
297
298 PetscFunctionBegin;
299 PetscValidHeaderSpecific(S,MAT_CLASSID,1);
300 ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
301 if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
302 PetscValidPointer(ksp,2);
303 Na = (Mat_SchurComplement*) S->data;
304 *ksp = Na->ksp;
305 PetscFunctionReturn(0);
306 }
307
308 /*@
309 MatSchurComplementSetKSP - Sets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
310
311 Not Collective
312
313 Input Parameters:
314 + S - matrix created with MatCreateSchurComplement()
315 - ksp - the linear solver object
316
317 Level: developer
318
319 Developer Notes:
320 This is used in PCFieldSplit to reuse the 0-split KSP to implement ksp(A00,Ap00) in S.
321 The KSP operators are overwritten with A00 and Ap00 currently set in S.
322
323 .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
324 @*/
MatSchurComplementSetKSP(Mat S,KSP ksp)325 PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
326 {
327 Mat_SchurComplement *Na;
328 PetscErrorCode ierr;
329 PetscBool isschur;
330
331 PetscFunctionBegin;
332 PetscValidHeaderSpecific(S,MAT_CLASSID,1);
333 ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
334 if (!isschur) PetscFunctionReturn(0);
335 PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
336 Na = (Mat_SchurComplement*) S->data;
337 ierr = PetscObjectReference((PetscObject)ksp);CHKERRQ(ierr);
338 ierr = KSPDestroy(&Na->ksp);CHKERRQ(ierr);
339 Na->ksp = ksp;
340 ierr = KSPSetOperators(Na->ksp, Na->A, Na->Ap);CHKERRQ(ierr);
341 PetscFunctionReturn(0);
342 }
343
344 /*@
345 MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
346
347 Collective on S
348
349 Input Parameters:
350 + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
351 . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
352 - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
353
354 Level: intermediate
355
356 Notes:
357 All four matrices must have the same MPI communicator
358
359 A00 and A11 must be square matrices
360
361 All of the matrices provided must have the same sizes as was used with MatCreateSchurComplement() or MatSchurComplementSetSubMatrices()
362 though they need not be the same matrices.
363
364 .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement()
365
366 @*/
MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)367 PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
368 {
369 PetscErrorCode ierr;
370 Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
371 PetscBool isschur;
372
373 PetscFunctionBegin;
374 PetscValidHeaderSpecific(S,MAT_CLASSID,1);
375 ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
376 if (!isschur) PetscFunctionReturn(0);
377 if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
378 PetscValidHeaderSpecific(A00,MAT_CLASSID,2);
379 PetscValidHeaderSpecific(Ap00,MAT_CLASSID,3);
380 PetscValidHeaderSpecific(A01,MAT_CLASSID,4);
381 PetscValidHeaderSpecific(A10,MAT_CLASSID,5);
382 PetscCheckSameComm(A00,2,Ap00,3);
383 PetscCheckSameComm(A00,2,A01,4);
384 PetscCheckSameComm(A00,2,A10,5);
385 if (A00->rmap->n != A00->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A00 %D do not equal local columns %D",A00->rmap->n,A00->cmap->n);
386 if (A00->rmap->n != Ap00->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A00 %D do not equal local rows of Ap00 %D",A00->rmap->n,Ap00->rmap->n);
387 if (Ap00->rmap->n != Ap00->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of Ap00 %D do not equal local columns %D",Ap00->rmap->n,Ap00->cmap->n);
388 if (A00->cmap->n != A01->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A00 %D do not equal local rows of A01 %D",A00->cmap->n,A01->rmap->n);
389 if (A10->cmap->n != A00->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A10 %D do not equal local rows of A00 %D",A10->cmap->n,A00->rmap->n);
390 if (A11) {
391 PetscValidHeaderSpecific(A11,MAT_CLASSID,6);
392 PetscCheckSameComm(A00,2,A11,6);
393 if (A10->rmap->n != A11->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A10 %D do not equal local rows A11 %D",A10->rmap->n,A11->rmap->n);
394 }
395
396 ierr = PetscObjectReference((PetscObject)A00);CHKERRQ(ierr);
397 ierr = PetscObjectReference((PetscObject)Ap00);CHKERRQ(ierr);
398 ierr = PetscObjectReference((PetscObject)A01);CHKERRQ(ierr);
399 ierr = PetscObjectReference((PetscObject)A10);CHKERRQ(ierr);
400 if (A11) {
401 ierr = PetscObjectReference((PetscObject)A11);CHKERRQ(ierr);
402 }
403
404 ierr = MatDestroy(&Na->A);CHKERRQ(ierr);
405 ierr = MatDestroy(&Na->Ap);CHKERRQ(ierr);
406 ierr = MatDestroy(&Na->B);CHKERRQ(ierr);
407 ierr = MatDestroy(&Na->C);CHKERRQ(ierr);
408 ierr = MatDestroy(&Na->D);CHKERRQ(ierr);
409
410 Na->A = A00;
411 Na->Ap = Ap00;
412 Na->B = A01;
413 Na->C = A10;
414 Na->D = A11;
415
416 ierr = KSPSetOperators(Na->ksp,A00,Ap00);CHKERRQ(ierr);
417 PetscFunctionReturn(0);
418 }
419
420
421 /*@C
422 MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
423
424 Collective on S
425
426 Input Parameter:
427 . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
428
429 Output Parameters:
430 + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
431 - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
432
433 Note: A11 is optional, and thus can be NULL. The submatrices are not increfed before they are returned and should not be modified or destroyed.
434
435 Level: intermediate
436
437 .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
438 @*/
MatSchurComplementGetSubMatrices(Mat S,Mat * A00,Mat * Ap00,Mat * A01,Mat * A10,Mat * A11)439 PetscErrorCode MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
440 {
441 Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
442 PetscErrorCode ierr;
443 PetscBool flg;
444
445 PetscFunctionBegin;
446 PetscValidHeaderSpecific(S,MAT_CLASSID,1);
447 ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);CHKERRQ(ierr);
448 if (!flg) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
449 if (A00) *A00 = Na->A;
450 if (Ap00) *Ap00 = Na->Ap;
451 if (A01) *A01 = Na->B;
452 if (A10) *A10 = Na->C;
453 if (A11) *A11 = Na->D;
454 PetscFunctionReturn(0);
455 }
456
457 #include <petsc/private/kspimpl.h>
458
459 /*@
460 MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
461
462 Collective on M
463
464 Input Parameter:
465 . M - the matrix obtained with MatCreateSchurComplement()
466
467 Output Parameter:
468 . S - the Schur complement matrix
469
470 Note: This can be expensive, so it is mainly for testing
471
472 Level: advanced
473
474 .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
475 @*/
MatSchurComplementComputeExplicitOperator(Mat M,Mat * S)476 PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
477 {
478 Mat B, C, D;
479 KSP ksp;
480 PC pc;
481 PetscBool isLU, isILU;
482 PetscReal fill = 2.0;
483 PetscErrorCode ierr;
484
485 PetscFunctionBegin;
486 ierr = MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);CHKERRQ(ierr);
487 ierr = MatSchurComplementGetKSP(M, &ksp);CHKERRQ(ierr);
488 ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
489 ierr = PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);CHKERRQ(ierr);
490 ierr = PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);CHKERRQ(ierr);
491 if (isLU || isILU) {
492 Mat fact, Bd, AinvB, AinvBd;
493 PetscReal eps = 1.0e-10;
494
495 /* This can be sped up for banded LU */
496 ierr = KSPSetUp(ksp);CHKERRQ(ierr);
497 ierr = PCFactorGetMatrix(pc, &fact);CHKERRQ(ierr);
498 ierr = MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);CHKERRQ(ierr);
499 ierr = MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);CHKERRQ(ierr);
500 ierr = MatMatSolve(fact, Bd, AinvBd);CHKERRQ(ierr);
501 ierr = MatDestroy(&Bd);CHKERRQ(ierr);
502 ierr = MatChop(AinvBd, eps);CHKERRQ(ierr);
503 ierr = MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);CHKERRQ(ierr);
504 ierr = MatDestroy(&AinvBd);CHKERRQ(ierr);
505 ierr = MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
506 ierr = MatDestroy(&AinvB);CHKERRQ(ierr);
507 } else {
508 Mat Ainv;
509
510 ierr = PCComputeOperator(pc, MATAIJ, &Ainv);CHKERRQ(ierr);
511 #if 0
512 /* Symmetric version */
513 ierr = MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
514 #else
515 /* Nonsymmetric version */
516 ierr = MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);CHKERRQ(ierr);
517 #endif
518 ierr = MatDestroy(&Ainv);CHKERRQ(ierr);
519 }
520 if (D) {
521 ierr = MatAXPY(*S, -1.0, D, DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
522 }
523 ierr = MatScale(*S,-1.0);CHKERRQ(ierr);
524 PetscFunctionReturn(0);
525 }
526
527 /* Developer Notes:
528 This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat * newmat,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat * newpmat)529 PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *newpmat)
530 {
531 PetscErrorCode ierr;
532 Mat A=NULL,Ap=NULL,B=NULL,C=NULL,D=NULL;
533 MatReuse reuse;
534
535 PetscFunctionBegin;
536 PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
537 PetscValidType(mat,1);
538 PetscValidHeaderSpecific(isrow0,IS_CLASSID,2);
539 PetscValidHeaderSpecific(iscol0,IS_CLASSID,3);
540 PetscValidHeaderSpecific(isrow1,IS_CLASSID,4);
541 PetscValidHeaderSpecific(iscol1,IS_CLASSID,5);
542 PetscValidLogicalCollectiveEnum(mat,mreuse,6);
543 PetscValidLogicalCollectiveEnum(mat,ainvtype,8);
544 PetscValidLogicalCollectiveEnum(mat,preuse,9);
545 if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(0);
546 if (mreuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*newmat,MAT_CLASSID,7);
547 if (preuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*newpmat,MAT_CLASSID,10);
548
549 if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
550
551 reuse = MAT_INITIAL_MATRIX;
552 if (mreuse == MAT_REUSE_MATRIX) {
553 ierr = MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);CHKERRQ(ierr);
554 if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
555 if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
556 ierr = MatDestroy(&Ap);CHKERRQ(ierr); /* get rid of extra reference */
557 reuse = MAT_REUSE_MATRIX;
558 }
559 ierr = MatCreateSubMatrix(mat,isrow0,iscol0,reuse,&A);CHKERRQ(ierr);
560 ierr = MatCreateSubMatrix(mat,isrow0,iscol1,reuse,&B);CHKERRQ(ierr);
561 ierr = MatCreateSubMatrix(mat,isrow1,iscol0,reuse,&C);CHKERRQ(ierr);
562 ierr = MatCreateSubMatrix(mat,isrow1,iscol1,reuse,&D);CHKERRQ(ierr);
563 switch (mreuse) {
564 case MAT_INITIAL_MATRIX:
565 ierr = MatCreateSchurComplement(A,A,B,C,D,newmat);CHKERRQ(ierr);
566 break;
567 case MAT_REUSE_MATRIX:
568 ierr = MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D);CHKERRQ(ierr);
569 break;
570 default:
571 if (mreuse != MAT_IGNORE_MATRIX) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse %d",(int)mreuse);
572 }
573 if (preuse != MAT_IGNORE_MATRIX) {
574 ierr = MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);CHKERRQ(ierr);
575 }
576 ierr = MatDestroy(&A);CHKERRQ(ierr);
577 ierr = MatDestroy(&B);CHKERRQ(ierr);
578 ierr = MatDestroy(&C);CHKERRQ(ierr);
579 ierr = MatDestroy(&D);CHKERRQ(ierr);
580 PetscFunctionReturn(0);
581 }
582
583 /*@
584 MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
585
586 Collective on A
587
588 Input Parameters:
589 + A - matrix in which the complement is to be taken
590 . isrow0 - rows to eliminate
591 . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
592 . isrow1 - rows in which the Schur complement is formed
593 . iscol1 - columns in which the Schur complement is formed
594 . mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in S
595 . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
596 MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG, or MAT_SCHUR_COMPLEMENT_AINV_LUMP
597 - preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in Sp
598
599 Output Parameters:
600 + S - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
601 - Sp - approximate Schur complement from which a preconditioner can be built
602
603 Note:
604 Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
605 application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
606 the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
607 before forming inv(diag(A00)).
608
609 Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
610 and column index sets. In that case, the user should call PetscObjectComposeFunction() on the *S matrix and pass mreuse of MAT_REUSE_MATRIX to set
611 "MatGetSchurComplement_C" to their function. If their function needs to fall back to the default implementation, it
612 should call MatGetSchurComplement_Basic().
613
614 MatCreateSchurComplement() takes as arguments the four submatrices and returns the virtual Schur complement (what this returns in S).
615
616 MatSchurComplementGetPmat() takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).
617
618 In other words calling MatCreateSchurComplement() followed by MatSchurComplementGetPmat() produces the same output as this function but with slightly different
619 inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.
620
621 Developer Notes:
622 The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
623 remove redundancy and be clearer and simplier.
624
625 Level: advanced
626
627 .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
628 @*/
MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat * S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat * Sp)629 PetscErrorCode MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
630 {
631 PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*) = NULL;
632
633 PetscFunctionBegin;
634 PetscValidHeaderSpecific(A,MAT_CLASSID,1);
635 PetscValidHeaderSpecific(isrow0,IS_CLASSID,2);
636 PetscValidHeaderSpecific(iscol0,IS_CLASSID,3);
637 PetscValidHeaderSpecific(isrow1,IS_CLASSID,4);
638 PetscValidHeaderSpecific(iscol1,IS_CLASSID,5);
639 PetscValidLogicalCollectiveEnum(A,mreuse,6);
640 if (mreuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*S,MAT_CLASSID,7);
641 PetscValidLogicalCollectiveEnum(A,ainvtype,8);
642 PetscValidLogicalCollectiveEnum(A,preuse,9);
643 if (preuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*Sp,MAT_CLASSID,10);
644 PetscValidType(A,1);
645 if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
646 f = NULL;
647 if (mreuse == MAT_REUSE_MATRIX) { /* This is the only situation, in which we can demand that the user pass a non-NULL pointer to non-garbage in S. */
648 ierr = PetscObjectQueryFunction((PetscObject)*S,"MatGetSchurComplement_C",&f);CHKERRQ(ierr);
649 }
650 if (f) {
651 ierr = (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);CHKERRQ(ierr);
652 } else {
653 ierr = MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);CHKERRQ(ierr);
654 }
655 PetscFunctionReturn(0);
656 }
657
658 /*@
659 MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
660
661 Not collective.
662
663 Input Parameters:
664 + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
665 - ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
666 MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG
667
668 Options database:
669 -mat_schur_complement_ainv_type diag | lump | blockdiag
670
671 Note:
672 Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
673 application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
674 the (0,0) block A00 in place of A00^{-1}. This rarely produces a scalable algorithm. Optionally, A00 can be lumped
675 before forming inv(diag(A00)).
676
677 Level: advanced
678
679 .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
680 @*/
MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)681 PetscErrorCode MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
682 {
683 PetscErrorCode ierr;
684 PetscBool isschur;
685 Mat_SchurComplement *schur;
686
687 PetscFunctionBegin;
688 PetscValidHeaderSpecific(S,MAT_CLASSID,1);
689 ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
690 if (!isschur) PetscFunctionReturn(0);
691 PetscValidLogicalCollectiveEnum(S,ainvtype,2);
692 schur = (Mat_SchurComplement*)S->data;
693 if (ainvtype != MAT_SCHUR_COMPLEMENT_AINV_DIAG && ainvtype != MAT_SCHUR_COMPLEMENT_AINV_LUMP && ainvtype != MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %d",(int)ainvtype);
694 schur->ainvtype = ainvtype;
695 PetscFunctionReturn(0);
696 }
697
698 /*@
699 MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
700
701 Not collective.
702
703 Input Parameter:
704 . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
705
706 Output Parameter:
707 . ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
708 MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG
709
710 Note:
711 Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
712 application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
713 the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
714 before forming inv(diag(A00)).
715
716 Level: advanced
717
718 .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
719 @*/
MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType * ainvtype)720 PetscErrorCode MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
721 {
722 PetscErrorCode ierr;
723 PetscBool isschur;
724 Mat_SchurComplement *schur;
725
726 PetscFunctionBegin;
727 PetscValidHeaderSpecific(S,MAT_CLASSID,1);
728 ierr = PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);CHKERRQ(ierr);
729 if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
730 schur = (Mat_SchurComplement*)S->data;
731 if (ainvtype) *ainvtype = schur->ainvtype;
732 PetscFunctionReturn(0);
733 }
734
735 /*@
736 MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
737
738 Collective on A00
739
740 Input Parameters:
741 + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A01,A10, and A11 are optional, implying zero matrices)
742 . ainvtype - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
743 - preuse - MAT_INITIAL_MATRIX for a new Sp, or MAT_REUSE_MATRIX to reuse an existing Sp, or MAT_IGNORE_MATRIX to put nothing in Sp
744
745 Output Parameter:
746 - Spmat - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
747
748 Note:
749 Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
750 application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
751 the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
752 before forming inv(diag(A00)).
753
754 Level: advanced
755
756 .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
757 @*/
MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat * Spmat)758 PetscErrorCode MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
759 {
760 PetscErrorCode ierr;
761 PetscInt N00;
762
763 PetscFunctionBegin;
764 /* Use an appropriate approximate inverse of A00 to form A11 - A10 inv(diag(A00)) A01; a NULL A01, A10 or A11 indicates a zero matrix. */
765 /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
766 if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Spmat: A01, A10 and A11 are all NULL.");
767 PetscValidLogicalCollectiveEnum(A11,preuse,6);
768 if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(0);
769
770 /* A zero size A00 or empty A01 or A10 imply S = A11. */
771 ierr = MatGetSize(A00,&N00,NULL);CHKERRQ(ierr);
772 if (!A01 || !A10 || !N00) {
773 if (preuse == MAT_INITIAL_MATRIX) {
774 ierr = MatDuplicate(A11,MAT_COPY_VALUES,Spmat);CHKERRQ(ierr);
775 } else { /* MAT_REUSE_MATRIX */
776 /* TODO: when can we pass SAME_NONZERO_PATTERN? */
777 ierr = MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
778 }
779 } else {
780 Mat AdB;
781 Vec diag;
782
783 if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
784 ierr = MatDuplicate(A01,MAT_COPY_VALUES,&AdB);CHKERRQ(ierr);
785 ierr = MatCreateVecs(A00,&diag,NULL);CHKERRQ(ierr);
786 if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
787 ierr = MatGetRowSum(A00,diag);CHKERRQ(ierr);
788 } else {
789 ierr = MatGetDiagonal(A00,diag);CHKERRQ(ierr);
790 }
791 ierr = VecReciprocal(diag);CHKERRQ(ierr);
792 ierr = MatDiagonalScale(AdB,diag,NULL);CHKERRQ(ierr);
793 ierr = VecDestroy(&diag);CHKERRQ(ierr);
794 } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
795 Mat A00_inv;
796 MatType type;
797 MPI_Comm comm;
798
799 ierr = PetscObjectGetComm((PetscObject)A00,&comm);CHKERRQ(ierr);
800 ierr = MatGetType(A00,&type);CHKERRQ(ierr);
801 ierr = MatCreate(comm,&A00_inv);CHKERRQ(ierr);
802 ierr = MatSetType(A00_inv,type);CHKERRQ(ierr);
803 ierr = MatInvertBlockDiagonalMat(A00,A00_inv);CHKERRQ(ierr);
804 ierr = MatMatMult(A00_inv,A01,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AdB);CHKERRQ(ierr);
805 ierr = MatDestroy(&A00_inv);CHKERRQ(ierr);
806 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
807 /* Cannot really reuse Spmat in MatMatMult() because of MatAYPX() -->
808 MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult() */
809 ierr = MatDestroy(Spmat);CHKERRQ(ierr);
810 ierr = MatMatMult(A10,AdB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,Spmat);CHKERRQ(ierr);
811 if (!A11) {
812 ierr = MatScale(*Spmat,-1.0);CHKERRQ(ierr);
813 } else {
814 /* TODO: when can we pass SAME_NONZERO_PATTERN? */
815 ierr = MatAYPX(*Spmat,-1,A11,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
816 }
817 ierr = MatDestroy(&AdB);CHKERRQ(ierr);
818 }
819 PetscFunctionReturn(0);
820 }
821
MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat * Spmat)822 PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
823 {
824 Mat A,B,C,D;
825 Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
826 PetscErrorCode ierr;
827
828 PetscFunctionBegin;
829 if (preuse == MAT_IGNORE_MATRIX) PetscFunctionReturn(0);
830 ierr = MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);CHKERRQ(ierr);
831 if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
832 ierr = MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);CHKERRQ(ierr);
833 PetscFunctionReturn(0);
834 }
835
836 /*@
837 MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
838
839 Collective on S
840
841 Input Parameters:
842 + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
843 - preuse - MAT_INITIAL_MATRIX for a new Sp, or MAT_REUSE_MATRIX to reuse an existing Sp, or MAT_IGNORE_MATRIX to put nothing in Sp
844
845 Output Parameter:
846 - Sp - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
847
848 Note:
849 Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
850 application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
851 the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
852 before forming inv(diag(A00)).
853
854 Sometimes users would like to provide problem-specific data in the Schur complement, usually only
855 for special row and column index sets. In that case, the user should call PetscObjectComposeFunction() to set
856 "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
857 it should call MatSchurComplementGetPmat_Basic().
858
859 Developer Notes:
860 The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
861 remove redundancy and be clearer and simplier.
862
863 Level: advanced
864
865 .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
866 @*/
MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat * Sp)867 PetscErrorCode MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
868 {
869 PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);
870
871 PetscFunctionBegin;
872 PetscValidHeaderSpecific(S,MAT_CLASSID,1);
873 PetscValidType(S,1);
874 PetscValidLogicalCollectiveEnum(S,preuse,2);
875 if (preuse != MAT_IGNORE_MATRIX) PetscValidPointer(Sp,3);
876 if (preuse == MAT_REUSE_MATRIX) PetscValidHeaderSpecific(*Sp,MAT_CLASSID,3);
877 if (S->factortype) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
878
879 ierr = PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);CHKERRQ(ierr);
880 if (f) {
881 ierr = (*f)(S,preuse,Sp);CHKERRQ(ierr);
882 } else {
883 ierr = MatSchurComplementGetPmat_Basic(S,preuse,Sp);CHKERRQ(ierr);
884 }
885 PetscFunctionReturn(0);
886 }
887
MatCreate_SchurComplement(Mat N)888 PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
889 {
890 PetscErrorCode ierr;
891 Mat_SchurComplement *Na;
892
893 PetscFunctionBegin;
894 ierr = PetscNewLog(N,&Na);CHKERRQ(ierr);
895 N->data = (void*) Na;
896
897 N->ops->destroy = MatDestroy_SchurComplement;
898 N->ops->getvecs = MatCreateVecs_SchurComplement;
899 N->ops->view = MatView_SchurComplement;
900 N->ops->mult = MatMult_SchurComplement;
901 N->ops->multtranspose = MatMultTranspose_SchurComplement;
902 N->ops->multadd = MatMultAdd_SchurComplement;
903 N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
904 N->assembled = PETSC_FALSE;
905 N->preallocated = PETSC_FALSE;
906
907 ierr = KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);CHKERRQ(ierr);
908 ierr = PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);CHKERRQ(ierr);
909 PetscFunctionReturn(0);
910 }
911