Actual source code: schurm.c
1: #include <../src/ksp/ksp/utils/schurm/schurm.h>
3: const char *const MatSchurComplementAinvTypes[] = {"DIAG","LUMP","BLOCKDIAG","MatSchurComplementAinvType","MAT_SCHUR_COMPLEMENT_AINV_",NULL};
5: PetscErrorCode MatCreateVecs_SchurComplement(Mat N,Vec *right,Vec *left)
6: {
7: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
8: PetscErrorCode ierr;
11: if (Na->D) {
12: MatCreateVecs(Na->D,right,left);
13: return(0);
14: }
15: if (right) {
16: MatCreateVecs(Na->B,right,NULL);
17: }
18: if (left) {
19: MatCreateVecs(Na->C,NULL,left);
20: }
21: return(0);
22: }
24: PetscErrorCode MatView_SchurComplement(Mat N,PetscViewer viewer)
25: {
26: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
27: PetscErrorCode ierr;
30: PetscViewerASCIIPrintf(viewer,"Schur complement A11 - A10 inv(A00) A01\n");
31: if (Na->D) {
32: PetscViewerASCIIPrintf(viewer,"A11\n");
33: PetscViewerASCIIPushTab(viewer);
34: MatView(Na->D,viewer);
35: PetscViewerASCIIPopTab(viewer);
36: } else {
37: PetscViewerASCIIPrintf(viewer,"A11 = 0\n");
38: }
39: PetscViewerASCIIPrintf(viewer,"A10\n");
40: PetscViewerASCIIPushTab(viewer);
41: MatView(Na->C,viewer);
42: PetscViewerASCIIPopTab(viewer);
43: PetscViewerASCIIPrintf(viewer,"KSP of A00\n");
44: PetscViewerASCIIPushTab(viewer);
45: KSPView(Na->ksp,viewer);
46: PetscViewerASCIIPopTab(viewer);
47: PetscViewerASCIIPrintf(viewer,"A01\n");
48: PetscViewerASCIIPushTab(viewer);
49: MatView(Na->B,viewer);
50: PetscViewerASCIIPopTab(viewer);
51: return(0);
52: }
54: /*
55: A11^T - A01^T ksptrans(A00,Ap00) A10^T
56: */
57: PetscErrorCode MatMultTranspose_SchurComplement(Mat N,Vec x,Vec y)
58: {
59: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
60: PetscErrorCode ierr;
63: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
64: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
65: MatMultTranspose(Na->C,x,Na->work1);
66: KSPSolveTranspose(Na->ksp,Na->work1,Na->work2);
67: MatMultTranspose(Na->B,Na->work2,y);
68: VecScale(y,-1.0);
69: if (Na->D) {
70: MatMultTransposeAdd(Na->D,x,y,y);
71: }
72: return(0);
73: }
75: /*
76: A11 - A10 ksp(A00,Ap00) A01
77: */
78: PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
79: {
80: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
81: PetscErrorCode ierr;
84: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
85: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
86: MatMult(Na->B,x,Na->work1);
87: KSPSolve(Na->ksp,Na->work1,Na->work2);
88: MatMult(Na->C,Na->work2,y);
89: VecScale(y,-1.0);
90: if (Na->D) {
91: MatMultAdd(Na->D,x,y,y);
92: }
93: return(0);
94: }
96: /*
97: A11 - A10 ksp(A00,Ap00) A01
98: */
99: PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
100: {
101: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
102: PetscErrorCode ierr;
105: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
106: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
107: MatMult(Na->B,x,Na->work1);
108: KSPSolve(Na->ksp,Na->work1,Na->work2);
109: if (y == z) {
110: VecScale(Na->work2,-1.0);
111: MatMultAdd(Na->C,Na->work2,z,z);
112: } else {
113: MatMult(Na->C,Na->work2,z);
114: VecAYPX(z,-1.0,y);
115: }
116: if (Na->D) {
117: MatMultAdd(Na->D,x,z,z);
118: }
119: return(0);
120: }
122: PetscErrorCode MatSetFromOptions_SchurComplement(PetscOptionItems *PetscOptionsObject,Mat N)
123: {
124: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
125: PetscErrorCode ierr;
128: PetscOptionsHead(PetscOptionsObject,"MatSchurComplementOptions");
129: Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
130: PetscOptionsEnum("-mat_schur_complement_ainv_type","Type of approximation for DIAGFORM(A00) used when assembling Sp = A11 - A10 inv(DIAGFORM(A00)) A01","MatSchurComplementSetAinvType",MatSchurComplementAinvTypes,(PetscEnum)Na->ainvtype,(PetscEnum*)&Na->ainvtype,NULL);
131: PetscOptionsTail();
132: KSPSetFromOptions(Na->ksp);
133: return(0);
134: }
136: PetscErrorCode MatDestroy_SchurComplement(Mat N)
137: {
138: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
139: PetscErrorCode ierr;
142: MatDestroy(&Na->A);
143: MatDestroy(&Na->Ap);
144: MatDestroy(&Na->B);
145: MatDestroy(&Na->C);
146: MatDestroy(&Na->D);
147: VecDestroy(&Na->work1);
148: VecDestroy(&Na->work2);
149: KSPDestroy(&Na->ksp);
150: PetscFree(N->data);
151: return(0);
152: }
154: /*@
155: MatCreateSchurComplement - Creates a new Mat that behaves like the Schur complement of a matrix
157: Collective on A00
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}
163: Output Parameter:
164: . S - the matrix that the Schur complement S = A11 - A10 ksp(A00,Ap00) A01
166: Level: intermediate
168: Notes:
169: The Schur complement is NOT explicitly formed! Rather, this function returns a virtual Schur complement
170: that can compute 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}.
173: All four matrices must have the same MPI communicator.
175: A00 and A11 must be square matrices.
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 Schur complement (useful for building a preconditioner for the Schur complement) which can be obtained from this
179: matrix with MatSchurComplementGetPmat()
181: Developer Notes:
182: The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
183: remove redundancy and be clearer and simpler.
185: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement(),
186: MatSchurComplementGetPmat(), MatSchurComplementSetSubMatrices()
188: @*/
189: PetscErrorCode MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
190: {
194: KSPInitializePackage();
195: MatCreate(PetscObjectComm((PetscObject)A00),S);
196: MatSetType(*S,MATSCHURCOMPLEMENT);
197: MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);
198: return(0);
199: }
201: /*@
202: MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
204: Collective on S
206: Input Parameters:
207: + S - matrix obtained with MatSetType(S,MATSCHURCOMPLEMENT)
208: . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
209: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
211: Level: intermediate
213: Notes:
214: The Schur complement is NOT explicitly formed! Rather, this
215: object performs the matrix-vector product of the Schur complement by using formula S = A11 - A10 ksp(A00,Ap00) A01
217: All four matrices must have the same MPI communicator.
219: A00 and A11 must be square matrices.
221: This is to be used in the context of code such as
222: $ MatSetType(S,MATSCHURCOMPLEMENT);
223: $ MatSchurComplementSetSubMatrices(S,...);
225: while MatSchurComplementUpdateSubMatrices() should only be called after MatCreateSchurComplement() or MatSchurComplementSetSubMatrices()
227: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatCreateSchurComplement(), MatGetSchurComplement()
229: @*/
230: PetscErrorCode MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
231: {
232: PetscErrorCode ierr;
233: Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
234: PetscBool isschur;
237: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
238: if (!isschur) return(0);
239: if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
247: 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);
248: 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);
249: 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);
250: 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);
251: 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);
252: if (A11) {
255: 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);
256: }
258: MatSetSizes(S,A10->rmap->n,A01->cmap->n,A10->rmap->N,A01->cmap->N);
259: PetscObjectReference((PetscObject)A00);
260: PetscObjectReference((PetscObject)Ap00);
261: PetscObjectReference((PetscObject)A01);
262: PetscObjectReference((PetscObject)A10);
263: Na->A = A00;
264: Na->Ap = Ap00;
265: Na->B = A01;
266: Na->C = A10;
267: Na->D = A11;
268: if (A11) {
269: PetscObjectReference((PetscObject)A11);
270: }
271: MatSetUp(S);
272: KSPSetOperators(Na->ksp,A00,Ap00);
273: S->assembled = PETSC_TRUE;
274: return(0);
275: }
277: /*@
278: MatSchurComplementGetKSP - Gets the KSP object that is used to solve with A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
280: Not Collective
282: Input Parameter:
283: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
285: Output Parameter:
286: . ksp - the linear solver object
288: Options Database:
289: . -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.
291: Level: intermediate
293: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
294: @*/
295: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
296: {
297: Mat_SchurComplement *Na;
298: PetscBool isschur;
299: PetscErrorCode ierr;
303: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
304: if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
306: Na = (Mat_SchurComplement*) S->data;
307: *ksp = Na->ksp;
308: return(0);
309: }
311: /*@
312: MatSchurComplementSetKSP - Sets the KSP object that is used to solve with A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
314: Not Collective
316: Input Parameters:
317: + S - matrix created with MatCreateSchurComplement()
318: - ksp - the linear solver object
320: Level: developer
322: Developer Notes:
323: This is used in PCFieldSplit to reuse the 0-split KSP to implement ksp(A00,Ap00) in S.
324: The KSP operators are overwritten with A00 and Ap00 currently set in S.
326: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
327: @*/
328: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
329: {
330: Mat_SchurComplement *Na;
331: PetscErrorCode ierr;
332: PetscBool isschur;
336: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
337: if (!isschur) return(0);
339: Na = (Mat_SchurComplement*) S->data;
340: PetscObjectReference((PetscObject)ksp);
341: KSPDestroy(&Na->ksp);
342: Na->ksp = ksp;
343: KSPSetOperators(Na->ksp, Na->A, Na->Ap);
344: return(0);
345: }
347: /*@
348: MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
350: Collective on S
352: Input Parameters:
353: + S - matrix obtained with MatCreateSchurComplement() (or MatSchurSetSubMatrices()) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
354: . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
355: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
357: Level: intermediate
359: Notes:
360: All four matrices must have the same MPI communicator
362: A00 and A11 must be square matrices
364: All of the matrices provided must have the same sizes as was used with MatCreateSchurComplement() or MatSchurComplementSetSubMatrices()
365: though they need not be the same matrices.
367: This can only be called after MatCreateSchurComplement() or MatSchurComplementSetSubMatrices(), it cannot be called immediately after MatSetType(S,MATSCHURCOMPLEMENT);
369: Developer Notes:
370: This code is almost identical to MatSchurComplementSetSubMatrices(). The API should be refactored.
372: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement()
374: @*/
375: PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
376: {
377: PetscErrorCode ierr;
378: Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
379: PetscBool isschur;
383: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
384: if (!isschur) return(0);
385: if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
393: 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);
394: 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);
395: 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);
396: 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);
397: 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);
398: if (A11) {
401: 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);
402: }
404: PetscObjectReference((PetscObject)A00);
405: PetscObjectReference((PetscObject)Ap00);
406: PetscObjectReference((PetscObject)A01);
407: PetscObjectReference((PetscObject)A10);
408: if (A11) {
409: PetscObjectReference((PetscObject)A11);
410: }
412: MatDestroy(&Na->A);
413: MatDestroy(&Na->Ap);
414: MatDestroy(&Na->B);
415: MatDestroy(&Na->C);
416: MatDestroy(&Na->D);
418: Na->A = A00;
419: Na->Ap = Ap00;
420: Na->B = A01;
421: Na->C = A10;
422: Na->D = A11;
424: KSPSetOperators(Na->ksp,A00,Ap00);
425: return(0);
426: }
428: /*@C
429: MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
431: Collective on S
433: Input Parameter:
434: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
436: Output Parameters:
437: + A00 - the upper-left block of the original matrix A = [A00 A01; A10 A11]
438: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
439: . A01 - the upper-right block of the original matrix A = [A00 A01; A10 A11]
440: . A10 - the lower-left block of the original matrix A = [A00 A01; A10 A11]
441: - A11 - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]
443: Note: A11 is optional, and thus can be NULL. The reference counts of the submatrices are not increased before they are returned and the matrices should not be modified or destroyed.
445: Level: intermediate
447: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
448: @*/
449: PetscErrorCode MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
450: {
451: Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
452: PetscErrorCode ierr;
453: PetscBool flg;
457: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);
458: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
459: if (A00) *A00 = Na->A;
460: if (Ap00) *Ap00 = Na->Ap;
461: if (A01) *A01 = Na->B;
462: if (A10) *A10 = Na->C;
463: if (A11) *A11 = Na->D;
464: return(0);
465: }
467: #include <petsc/private/kspimpl.h>
469: /*@
470: MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
472: Collective on M
474: Input Parameter:
475: . M - the matrix obtained with MatCreateSchurComplement()
477: Output Parameter:
478: . S - the Schur complement matrix
480: Notes:
481: This can be expensive, so it is mainly for testing
483: Use MatSchurComplementGetPmat() to get a sparse approximation for the Schur complement suitable for use in building a preconditioner
485: Level: advanced
487: .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate(), MatSchurComplementGetPmat()
488: @*/
489: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat A, Mat *S)
490: {
491: Mat B, C, D, Bd, AinvBd;
492: KSP ksp;
493: PetscInt n,N,m,M;
497: MatSchurComplementGetSubMatrices(A, NULL, NULL, &B, &C, &D);
498: MatSchurComplementGetKSP(A, &ksp);
499: KSPSetUp(ksp);
500: MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);
501: MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);
502: KSPMatSolve(ksp, Bd, AinvBd);
503: MatDestroy(&Bd);
504: MatChop(AinvBd, PETSC_SMALL);
505: if (D) {
506: MatGetLocalSize(D, &m, &n);
507: MatGetSize(D, &M, &N);
508: MatCreateDense(PetscObjectComm((PetscObject)A), m, n, M, N, NULL, S);
509: }
510: MatMatMult(C, AinvBd, D ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DEFAULT, S);
511: MatDestroy(&AinvBd);
512: if (D) {
513: MatAXPY(*S, -1.0, D, DIFFERENT_NONZERO_PATTERN);
514: }
515: MatConvert(*S, MATAIJ, MAT_INPLACE_MATRIX, S);
516: MatScale(*S, -1.0);
517: return(0);
518: }
520: /* Developer Notes:
521: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
522: PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *Sp)
523: {
525: Mat A=NULL,Ap=NULL,B=NULL,C=NULL,D=NULL;
526: MatReuse reuse;
538: if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) return(0);
542: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
544: reuse = MAT_INITIAL_MATRIX;
545: if (mreuse == MAT_REUSE_MATRIX) {
546: MatSchurComplementGetSubMatrices(*S,&A,&Ap,&B,&C,&D);
547: if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
548: if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
549: MatDestroy(&Ap); /* get rid of extra reference */
550: reuse = MAT_REUSE_MATRIX;
551: }
552: MatCreateSubMatrix(mat,isrow0,iscol0,reuse,&A);
553: MatCreateSubMatrix(mat,isrow0,iscol1,reuse,&B);
554: MatCreateSubMatrix(mat,isrow1,iscol0,reuse,&C);
555: MatCreateSubMatrix(mat,isrow1,iscol1,reuse,&D);
556: switch (mreuse) {
557: case MAT_INITIAL_MATRIX:
558: MatCreateSchurComplement(A,A,B,C,D,S);
559: break;
560: case MAT_REUSE_MATRIX:
561: MatSchurComplementUpdateSubMatrices(*S,A,A,B,C,D);
562: break;
563: default:
564: if (mreuse != MAT_IGNORE_MATRIX) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse %d",(int)mreuse);
565: }
566: if (preuse != MAT_IGNORE_MATRIX) {
567: MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,Sp);
568: }
569: MatDestroy(&A);
570: MatDestroy(&B);
571: MatDestroy(&C);
572: MatDestroy(&D);
573: return(0);
574: }
576: /*@
577: MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
579: Collective on A
581: Input Parameters:
582: + A - matrix in which the complement is to be taken
583: . isrow0 - rows to eliminate
584: . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
585: . isrow1 - rows in which the Schur complement is formed
586: . iscol1 - columns in which the Schur complement is formed
587: . mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in S
588: . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
589: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG, or MAT_SCHUR_COMPLEMENT_AINV_LUMP
590: - preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in Sp
592: Output Parameters:
593: + S - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
594: - Sp - approximate Schur complement from which a preconditioner can be built A11 - A10 inv(DIAGFORM(A00)) A01
596: Note:
597: Since the real Schur complement is usually dense, providing a good approximation to Sp usually requires
598: application-specific information.
600: Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
601: 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
602: "MatGetSchurComplement_C" to their function. If their function needs to fall back to the default implementation, it
603: should call MatGetSchurComplement_Basic().
605: MatCreateSchurComplement() takes as arguments the four submatrices and returns the virtual Schur complement (what this function returns in S).
607: MatSchurComplementGetPmat() takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).
609: In other words calling MatCreateSchurComplement() followed by MatSchurComplementGetPmat() produces the same output as this function but with slightly different
610: inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.
612: Developer Notes:
613: The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
614: remove redundancy and be clearer and simpler.
616: Level: advanced
618: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
619: @*/
620: PetscErrorCode MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
621: {
622: PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*) = NULL;
636: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
637: f = NULL;
638: 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. */
639: PetscObjectQueryFunction((PetscObject)*S,"MatGetSchurComplement_C",&f);
640: }
641: if (f) {
642: (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
643: } else {
644: MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
645: }
646: return(0);
647: }
649: /*@
650: MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
652: Not collective.
654: Input Parameters:
655: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
656: - ainvtype - type of approximation to be used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
657: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG
659: Options database:
660: -mat_schur_complement_ainv_type diag | lump | blockdiag
662: Level: advanced
664: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
665: @*/
666: PetscErrorCode MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
667: {
668: PetscErrorCode ierr;
669: PetscBool isschur;
670: Mat_SchurComplement *schur;
674: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
675: if (!isschur) return(0);
677: schur = (Mat_SchurComplement*)S->data;
678: 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);
679: schur->ainvtype = ainvtype;
680: return(0);
681: }
683: /*@
684: MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
686: Not collective.
688: Input Parameter:
689: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
691: Output Parameter:
692: . ainvtype - type of approximation used to form approximate Schur complement Sp = A11 - A10 inv(DIAGFORM(A00)) A01:
693: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG
695: Level: advanced
697: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
698: @*/
699: PetscErrorCode MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
700: {
701: PetscErrorCode ierr;
702: PetscBool isschur;
703: Mat_SchurComplement *schur;
707: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
708: if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
709: schur = (Mat_SchurComplement*)S->data;
710: if (ainvtype) *ainvtype = schur->ainvtype;
711: return(0);
712: }
714: /*@
715: MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by explicitly assembling the sparse matrix
716: Sp = A11 - A10 inv(DIAGFORM(A00)) A01
718: Collective on A00
720: Input Parameters:
721: + A00 - the upper-left part of the original matrix A = [A00 A01; A10 A11]
722: . A01 - (optional) the upper-right part of the original matrix A = [A00 A01; A10 A11]
723: . A10 - (optional) the lower-left part of the original matrix A = [A00 A01; A10 A11]
724: . A11 - (optional) the lower-right part of the original matrix A = [A00 A01; A10 A11]
725: . ainvtype - type of approximation for DIAGFORM(A00) used when forming Sp = A11 - A10 inv(DIAGFORM(A00)) A01. See MatSchurComplementAinvType.
726: - 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
728: Output Parameter:
729: - Sp - approximate Schur complement suitable for preconditioning the true Schur complement S = A11 - A10 inv(A00) A01
731: Level: advanced
733: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
734: @*/
735: PetscErrorCode MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
736: {
738: PetscInt N00;
741: /* Use an appropriate approximate inverse of A00 to form A11 - A10 inv(DIAGFORM(A00)) A01; a NULL A01, A10 or A11 indicates a zero matrix. */
742: /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
743: if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Sp: A01, A10 and A11 are all NULL.");
745: if (preuse == MAT_IGNORE_MATRIX) return(0);
747: /* A zero size A00 or empty A01 or A10 imply S = A11. */
748: MatGetSize(A00,&N00,NULL);
749: if (!A01 || !A10 || !N00) {
750: if (preuse == MAT_INITIAL_MATRIX) {
751: MatDuplicate(A11,MAT_COPY_VALUES,Sp);
752: } else { /* MAT_REUSE_MATRIX */
753: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
754: MatCopy(A11,*Sp,DIFFERENT_NONZERO_PATTERN);
755: }
756: } else {
757: Mat AdB;
758: Vec diag;
760: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
761: MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
762: MatCreateVecs(A00,&diag,NULL);
763: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
764: MatGetRowSum(A00,diag);
765: } else {
766: MatGetDiagonal(A00,diag);
767: }
768: VecReciprocal(diag);
769: MatDiagonalScale(AdB,diag,NULL);
770: VecDestroy(&diag);
771: } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
772: Mat A00_inv;
773: MatType type;
774: MPI_Comm comm;
776: PetscObjectGetComm((PetscObject)A00,&comm);
777: MatGetType(A00,&type);
778: MatCreate(comm,&A00_inv);
779: MatSetType(A00_inv,type);
780: MatInvertBlockDiagonalMat(A00,A00_inv);
781: MatMatMult(A00_inv,A01,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AdB);
782: MatDestroy(&A00_inv);
783: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
784: /* Cannot really reuse Sp in MatMatMult() because of MatAYPX() -->
785: MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult() */
786: MatDestroy(Sp);
787: MatMatMult(A10,AdB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,Sp);
788: if (!A11) {
789: MatScale(*Sp,-1.0);
790: } else {
791: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
792: MatAYPX(*Sp,-1,A11,DIFFERENT_NONZERO_PATTERN);
793: }
794: MatDestroy(&AdB);
795: }
796: return(0);
797: }
799: PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Sp)
800: {
801: Mat A,B,C,D;
802: Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
803: PetscErrorCode ierr;
806: if (preuse == MAT_IGNORE_MATRIX) return(0);
807: MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
808: if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
809: MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Sp);
810: return(0);
811: }
813: /*@
814: MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(DIAGFORM(A00)) A01
816: Collective on S
818: Input Parameters:
819: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) that implements the action of A11 - A10 ksp(A00,Ap00) A01
820: - 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
822: Output Parameter:
823: - Sp - approximate Schur complement suitable for preconditioning the exact Schur complement S = A11 - A10 inv(A00) A01
825: Note:
826: The approximation of Sp depends on the the argument passed to to MatSchurComplementSetAinvType()
827: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG or
828: -mat_schur_complement_ainv_type <diag,lump,blockdiag>
830: Sometimes users would like to provide problem-specific data in the Schur complement, usually only
831: for special row and column index sets. In that case, the user should call PetscObjectComposeFunction() to set
832: "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
833: it should call MatSchurComplementGetPmat_Basic().
835: Developer Notes:
836: The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
837: remove redundancy and be clearer and simpler.
839: This routine should be called MatSchurComplementCreatePmat()
841: Level: advanced
843: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
844: @*/
845: PetscErrorCode MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
846: {
847: PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);
855: if (S->factortype) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
857: PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);
858: if (f) {
859: (*f)(S,preuse,Sp);
860: } else {
861: MatSchurComplementGetPmat_Basic(S,preuse,Sp);
862: }
863: return(0);
864: }
866: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
867: {
868: PetscErrorCode ierr;
869: Mat_SchurComplement *Na;
872: PetscNewLog(N,&Na);
873: N->data = (void*) Na;
875: N->ops->destroy = MatDestroy_SchurComplement;
876: N->ops->getvecs = MatCreateVecs_SchurComplement;
877: N->ops->view = MatView_SchurComplement;
878: N->ops->mult = MatMult_SchurComplement;
879: N->ops->multtranspose = MatMultTranspose_SchurComplement;
880: N->ops->multadd = MatMultAdd_SchurComplement;
881: N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
882: N->assembled = PETSC_FALSE;
883: N->preallocated = PETSC_FALSE;
885: KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
886: PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
887: return(0);
888: }