Actual source code: schurm.c
petsc-3.9.4 2018-09-11
1: #include <petsc/private/matimpl.h>
2: #include <petscksp.h>
3: const char *const MatSchurComplementAinvTypes[] = {"DIAG","LUMP","BLOCKDIAG","MatSchurComplementAinvType","MAT_SCHUR_COMPLEMENT_AINV_",0};
5: typedef struct {
6: Mat A,Ap,B,C,D;
7: KSP ksp;
8: Vec work1,work2;
9: MatSchurComplementAinvType ainvtype;
10: } Mat_SchurComplement;
12: PetscErrorCode MatCreateVecs_SchurComplement(Mat N,Vec *right,Vec *left)
13: {
14: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
15: PetscErrorCode ierr;
18: if (Na->D) {
19: MatCreateVecs(Na->D,right,left);
20: return(0);
21: }
22: if (right) {
23: MatCreateVecs(Na->B,right,NULL);
24: }
25: if (left) {
26: MatCreateVecs(Na->C,NULL,left);
27: }
28: return(0);
29: }
31: PetscErrorCode MatView_SchurComplement(Mat N,PetscViewer viewer)
32: {
33: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
34: PetscErrorCode ierr;
37: PetscViewerASCIIPrintf(viewer,"Schur complement A11 - A10 inv(A00) A01\n");
38: if (Na->D) {
39: PetscViewerASCIIPrintf(viewer,"A11\n");
40: PetscViewerASCIIPushTab(viewer);
41: MatView(Na->D,viewer);
42: PetscViewerASCIIPopTab(viewer);
43: } else {
44: PetscViewerASCIIPrintf(viewer,"A11 = 0\n");
45: }
46: PetscViewerASCIIPrintf(viewer,"A10\n");
47: PetscViewerASCIIPushTab(viewer);
48: MatView(Na->C,viewer);
49: PetscViewerASCIIPopTab(viewer);
50: PetscViewerASCIIPrintf(viewer,"KSP of A00\n");
51: PetscViewerASCIIPushTab(viewer);
52: KSPView(Na->ksp,viewer);
53: PetscViewerASCIIPopTab(viewer);
54: PetscViewerASCIIPrintf(viewer,"A01\n");
55: PetscViewerASCIIPushTab(viewer);
56: MatView(Na->B,viewer);
57: PetscViewerASCIIPopTab(viewer);
58: return(0);
59: }
61: /*
62: A11^T - A01^T ksptrans(A00,Ap00) A10^T
63: */
64: PetscErrorCode MatMultTranspose_SchurComplement(Mat N,Vec x,Vec y)
65: {
66: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
67: PetscErrorCode ierr;
70: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
71: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
72: MatMultTranspose(Na->C,x,Na->work1);
73: KSPSolveTranspose(Na->ksp,Na->work1,Na->work2);
74: MatMultTranspose(Na->B,Na->work2,y);
75: VecScale(y,-1.0);
76: if (Na->D) {
77: MatMultTransposeAdd(Na->D,x,y,y);
78: }
79: return(0);
80: }
82: /*
83: A11 - A10 ksp(A00,Ap00) A01
84: */
85: PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
86: {
87: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
88: PetscErrorCode ierr;
91: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
92: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
93: MatMult(Na->B,x,Na->work1);
94: KSPSolve(Na->ksp,Na->work1,Na->work2);
95: MatMult(Na->C,Na->work2,y);
96: VecScale(y,-1.0);
97: if (Na->D) {
98: MatMultAdd(Na->D,x,y,y);
99: }
100: return(0);
101: }
103: /*
104: A11 - A10 ksp(A00,Ap00) A01
105: */
106: PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
107: {
108: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
109: PetscErrorCode ierr;
112: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
113: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
114: MatMult(Na->B,x,Na->work1);
115: KSPSolve(Na->ksp,Na->work1,Na->work2);
116: if (y == z) {
117: VecScale(Na->work2,-1.0);
118: MatMultAdd(Na->C,Na->work2,z,z);
119: } else {
120: MatMult(Na->C,Na->work2,z);
121: VecAYPX(z,-1.0,y);
122: }
123: if (Na->D) {
124: MatMultAdd(Na->D,x,z,z);
125: }
126: return(0);
127: }
129: PetscErrorCode MatSetFromOptions_SchurComplement(PetscOptionItems *PetscOptionsObject,Mat N)
130: {
131: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
132: PetscErrorCode ierr;
135: PetscOptionsHead(PetscOptionsObject,"MatSchurComplementOptions");
136: Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
137: 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);
138: PetscOptionsTail();
139: KSPSetFromOptions(Na->ksp);
140: return(0);
141: }
143: PetscErrorCode MatDestroy_SchurComplement(Mat N)
144: {
145: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
146: PetscErrorCode ierr;
149: MatDestroy(&Na->A);
150: MatDestroy(&Na->Ap);
151: MatDestroy(&Na->B);
152: MatDestroy(&Na->C);
153: MatDestroy(&Na->D);
154: VecDestroy(&Na->work1);
155: VecDestroy(&Na->work2);
156: KSPDestroy(&Na->ksp);
157: PetscFree(N->data);
158: return(0);
159: }
161: /*@C
162: MatCreateSchurComplement - Creates a new matrix object that behaves like the Schur complement of a matrix
164: Collective on Mat
166: Input Parameters:
167: + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
168: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
170: Output Parameter:
171: . S - the matrix that the Schur complement S = A11 - A10 ksp(A00,Ap00) A01
173: Level: intermediate
175: Notes: The Schur complement is NOT actually formed! Rather, this
176: object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
177: for Schur complement S and a KSP solver to approximate the action of A^{-1}.
179: All four matrices must have the same MPI communicator.
181: A00 and A11 must be square matrices.
183: MatGetSchurComplement() takes as arguments the index sets for the submatrices and returns both the virtual Schur complement (what this returns) plus
184: a sparse approximation to the true Schur complement (useful for building a preconditioner for the Schur complement).
186: MatSchurComplementGetPmat() can be called on the output of this function to generate an explicit approximation to the Schur complement.
188: Developer Notes: The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
189: remove redundancy and be clearer and simplier.
192: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement(),
193: MatSchurComplementGetPmat()
195: @*/
196: PetscErrorCode MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
197: {
201: KSPInitializePackage();
202: MatCreate(PetscObjectComm((PetscObject)A00),S);
203: MatSetType(*S,MATSCHURCOMPLEMENT);
204: MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);
205: return(0);
206: }
208: /*@
209: MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
211: Collective on Mat
213: Input Parameter:
214: + S - matrix obtained with MatCreateSchurComplement (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
215: . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
216: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
218: Level: intermediate
220: Notes: The Schur complement is NOT actually formed! Rather, this
221: object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
222: for Schur complement S and a KSP solver to approximate the action of A^{-1}.
224: All four matrices must have the same MPI communicator.
226: A00 and A11 must be square matrices.
228: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatCreateSchurComplement(), MatGetSchurComplement()
230: @*/
231: PetscErrorCode MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
232: {
233: PetscErrorCode ierr;
234: Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
235: PetscBool isschur;
238: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
239: if (!isschur) return(0);
240: if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
248: 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);
249: 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);
250: 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);
251: 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);
252: 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);
253: if (A11) {
256: 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);
257: }
259: MatSetSizes(S,A10->rmap->n,A01->cmap->n,A10->rmap->N,A01->cmap->N);
260: PetscObjectReference((PetscObject)A00);
261: PetscObjectReference((PetscObject)Ap00);
262: PetscObjectReference((PetscObject)A01);
263: PetscObjectReference((PetscObject)A10);
264: Na->A = A00;
265: Na->Ap = Ap00;
266: Na->B = A01;
267: Na->C = A10;
268: Na->D = A11;
269: if (A11) {
270: PetscObjectReference((PetscObject)A11);
271: }
272: MatSetUp(S);
273: KSPSetOperators(Na->ksp,A00,Ap00);
274: S->assembled = PETSC_TRUE;
275: return(0);
276: }
278: /*@
279: MatSchurComplementGetKSP - Gets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
281: Not Collective
283: Input Parameter:
284: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
286: Output Parameter:
287: . ksp - the linear solver object
289: Options Database:
290: . -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.
292: Level: intermediate
294: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
295: @*/
296: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
297: {
298: Mat_SchurComplement *Na;
299: PetscBool isschur;
300: PetscErrorCode ierr;
304: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
305: if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
307: Na = (Mat_SchurComplement*) S->data;
308: *ksp = Na->ksp;
309: return(0);
310: }
312: /*@
313: MatSchurComplementSetKSP - Sets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
315: Not Collective
317: Input Parameters:
318: + S - matrix created with MatCreateSchurComplement()
319: - ksp - the linear solver object
321: Level: developer
323: Developer Notes:
324: This is used in PCFieldSplit to reuse the 0-split KSP to implement ksp(A00,Ap00) in S.
325: The KSP operators are overwritten with A00 and Ap00 currently set in S.
327: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
328: @*/
329: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
330: {
331: Mat_SchurComplement *Na;
332: PetscErrorCode ierr;
333: PetscBool isschur;
337: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
338: if (!isschur) return(0);
340: Na = (Mat_SchurComplement*) S->data;
341: PetscObjectReference((PetscObject)ksp);
342: KSPDestroy(&Na->ksp);
343: Na->ksp = ksp;
344: KSPSetOperators(Na->ksp, Na->A, Na->Ap);
345: return(0);
346: }
348: /*@
349: MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
351: Collective on Mat
353: Input Parameters:
354: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
355: . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
356: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
358: Level: intermediate
360: Notes: 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: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement()
369: @*/
370: PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
371: {
372: PetscErrorCode ierr;
373: Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
374: PetscBool isschur;
378: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
379: if (!isschur) return(0);
380: if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
388: 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);
389: 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);
390: 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);
391: 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);
392: 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);
393: if (A11) {
396: 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);
397: }
399: PetscObjectReference((PetscObject)A00);
400: PetscObjectReference((PetscObject)Ap00);
401: PetscObjectReference((PetscObject)A01);
402: PetscObjectReference((PetscObject)A10);
403: if (A11) {
404: PetscObjectReference((PetscObject)A11);
405: }
407: MatDestroy(&Na->A);
408: MatDestroy(&Na->Ap);
409: MatDestroy(&Na->B);
410: MatDestroy(&Na->C);
411: MatDestroy(&Na->D);
413: Na->A = A00;
414: Na->Ap = Ap00;
415: Na->B = A01;
416: Na->C = A10;
417: Na->D = A11;
419: KSPSetOperators(Na->ksp,A00,Ap00);
420: return(0);
421: }
424: /*@C
425: MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
427: Collective on Mat
429: Input Parameter:
430: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
432: Output Paramters:
433: + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
434: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
436: 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.
438: Level: intermediate
440: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
441: @*/
442: PetscErrorCode MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
443: {
444: Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
445: PetscErrorCode ierr;
446: PetscBool flg;
450: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);
451: if (!flg) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
452: if (A00) *A00 = Na->A;
453: if (Ap00) *Ap00 = Na->Ap;
454: if (A01) *A01 = Na->B;
455: if (A10) *A10 = Na->C;
456: if (A11) *A11 = Na->D;
457: return(0);
458: }
460: #include <petsc/private/kspimpl.h>
462: /*@
463: MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
465: Collective on Mat
467: Input Parameter:
468: . M - the matrix obtained with MatCreateSchurComplement()
470: Output Parameter:
471: . S - the Schur complement matrix
473: Note: This can be expensive, so it is mainly for testing
475: Level: advanced
477: .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
478: @*/
479: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
480: {
481: Mat B, C, D;
482: KSP ksp;
483: PC pc;
484: PetscBool isLU, isILU;
485: PetscReal fill = 2.0;
489: MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);
490: MatSchurComplementGetKSP(M, &ksp);
491: KSPGetPC(ksp, &pc);
492: PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);
493: PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);
494: if (isLU || isILU) {
495: Mat fact, Bd, AinvB, AinvBd;
496: PetscReal eps = 1.0e-10;
498: /* This can be sped up for banded LU */
499: KSPSetUp(ksp);
500: PCFactorGetMatrix(pc, &fact);
501: MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);
502: MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);
503: MatMatSolve(fact, Bd, AinvBd);
504: MatDestroy(&Bd);
505: MatChop(AinvBd, eps);
506: MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);
507: MatDestroy(&AinvBd);
508: MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);
509: MatDestroy(&AinvB);
510: } else {
511: Mat Ainvd, Ainv;
513: PCComputeExplicitOperator(pc, &Ainvd);
514: MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);
515: MatDestroy(&Ainvd);
516: #if 0
517: /* Symmetric version */
518: MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);
519: #else
520: /* Nonsymmetric version */
521: MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);
522: #endif
523: MatDestroy(&Ainv);
524: }
525: if (D) {
526: MatAXPY(*S, -1.0, D, DIFFERENT_NONZERO_PATTERN);
527: }
528: MatScale(*S,-1.0);
529: return(0);
530: }
532: /* Developer Notes: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
533: PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *newpmat)
534: {
536: Mat A=0,Ap=0,B=0,C=0,D=0;
537: MatReuse reuse;
549: if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) return(0);
553: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
555: reuse = MAT_INITIAL_MATRIX;
556: if (mreuse == MAT_REUSE_MATRIX) {
557: MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);
558: if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
559: if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
560: MatDestroy(&Ap); /* get rid of extra reference */
561: reuse = MAT_REUSE_MATRIX;
562: }
563: MatCreateSubMatrix(mat,isrow0,iscol0,reuse,&A);
564: MatCreateSubMatrix(mat,isrow0,iscol1,reuse,&B);
565: MatCreateSubMatrix(mat,isrow1,iscol0,reuse,&C);
566: MatCreateSubMatrix(mat,isrow1,iscol1,reuse,&D);
567: switch (mreuse) {
568: case MAT_INITIAL_MATRIX:
569: MatCreateSchurComplement(A,A,B,C,D,newmat);
570: break;
571: case MAT_REUSE_MATRIX:
572: MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D);
573: break;
574: default:
575: if (mreuse != MAT_IGNORE_MATRIX) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse %d",(int)mreuse);
576: }
577: if (preuse != MAT_IGNORE_MATRIX) {
578: MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);
579: }
580: MatDestroy(&A);
581: MatDestroy(&B);
582: MatDestroy(&C);
583: MatDestroy(&D);
584: return(0);
585: }
587: /*@
588: MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
590: Collective on Mat
592: Input Parameters:
593: + A - matrix in which the complement is to be taken
594: . isrow0 - rows to eliminate
595: . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
596: . isrow1 - rows in which the Schur complement is formed
597: . iscol1 - columns in which the Schur complement is formed
598: . mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in S
599: . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
600: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG, or MAT_SCHUR_COMPLEMENT_AINV_LUMP
601: - preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in Sp
603: Output Parameters:
604: + S - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
605: - Sp - approximate Schur complement from which a preconditioner can be built
607: Note:
608: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
609: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
610: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
611: before forming inv(diag(A00)).
613: Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
614: 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
615: "MatGetSchurComplement_C" to their function. If their function needs to fall back to the default implementation, it
616: should call MatGetSchurComplement_Basic().
618: MatCreateSchurComplement() takes as arguments the four submatrices and returns the virtual Schur complement (what this returns in S).
620: MatSchurComplementGetPmat() takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).
622: In other words calling MatCreateSchurComplement() followed by MatSchurComplementGetPmat() produces the same output as this function but with slightly different
623: inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.
625: Developer Notes: The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
626: remove redundancy and be clearer and simplier.
628: Level: advanced
630: Concepts: matrices^submatrices
632: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
633: @*/
634: PetscErrorCode MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
635: {
636: PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*) = NULL;
650: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
651: f = NULL;
652: 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. */
653: PetscObjectQueryFunction((PetscObject)*S,"MatGetSchurComplement_C",&f);
654: }
655: if (f) {
656: (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
657: } else {
658: MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
659: }
660: return(0);
661: }
663: /*@
664: MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
666: Not collective.
668: Input Parameters:
669: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
670: - ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
671: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG
673: Options database:
674: -mat_schur_complement_ainv_type diag | lump | blockdiag
676: Note:
677: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
678: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
679: the (0,0) block A00 in place of A00^{-1}. This rarely produces a scalable algorithm. Optionally, A00 can be lumped
680: before forming inv(diag(A00)).
682: Level: advanced
684: Concepts: matrices^submatrices
686: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
687: @*/
688: PetscErrorCode MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
689: {
690: PetscErrorCode ierr;
691: PetscBool isschur;
692: Mat_SchurComplement *schur;
696: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
697: if (!isschur) return(0);
699: schur = (Mat_SchurComplement*)S->data;
700: 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);
701: schur->ainvtype = ainvtype;
702: return(0);
703: }
705: /*@
706: MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
708: Not collective.
710: Input Parameter:
711: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
713: Output Parameter:
714: . ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
715: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG
717: Note:
718: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
719: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
720: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
721: before forming inv(diag(A00)).
723: Level: advanced
725: Concepts: matrices^submatrices
727: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
728: @*/
729: PetscErrorCode MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
730: {
731: PetscErrorCode ierr;
732: PetscBool isschur;
733: Mat_SchurComplement *schur;
737: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
738: if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
739: schur = (Mat_SchurComplement*)S->data;
740: if (ainvtype) *ainvtype = schur->ainvtype;
741: return(0);
742: }
744: /*@
745: MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
747: Collective on Mat
749: Input Parameters:
750: + 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)
751: . ainvtype - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
752: - 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
754: Output Parameter:
755: - Spmat - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
757: Note:
758: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
759: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
760: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
761: before forming inv(diag(A00)).
763: Level: advanced
765: Concepts: matrices^submatrices
767: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
768: @*/
769: PetscErrorCode MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
770: {
772: PetscInt N00;
775: /* 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. */
776: /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
777: if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Spmat: A01, A10 and A11 are all NULL.");
779: if (preuse == MAT_IGNORE_MATRIX) return(0);
781: /* A zero size A00 or empty A01 or A10 imply S = A11. */
782: MatGetSize(A00,&N00,NULL);
783: if (!A01 || !A10 || !N00) {
784: if (preuse == MAT_INITIAL_MATRIX) {
785: MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
786: } else { /* MAT_REUSE_MATRIX */
787: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
788: MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
789: }
790: } else {
791: Mat AdB;
792: Vec diag;
794: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
795: MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
796: MatCreateVecs(A00,&diag,NULL);
797: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
798: MatGetRowSum(A00,diag);
799: } else {
800: MatGetDiagonal(A00,diag);
801: }
802: VecReciprocal(diag);
803: MatDiagonalScale(AdB,diag,NULL);
804: VecDestroy(&diag);
805: } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
806: Mat A00_inv;
807: MatType type;
808: MPI_Comm comm;
810: PetscObjectGetComm((PetscObject)A00,&comm);
811: MatGetType(A00,&type);
812: MatCreate(comm,&A00_inv);
813: MatSetType(A00_inv,type);
814: MatInvertBlockDiagonalMat(A00,A00_inv);
815: MatMatMult(A00_inv,A01,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AdB);
816: MatDestroy(&A00_inv);
817: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
818: /* Cannot really reuse Spmat in MatMatMult() because of MatAYPX() -->
819: MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult() */
820: MatDestroy(Spmat);
821: MatMatMult(A10,AdB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,Spmat);
822: if (!A11) {
823: MatScale(*Spmat,-1.0);
824: } else {
825: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
826: MatAYPX(*Spmat,-1,A11,DIFFERENT_NONZERO_PATTERN);
827: }
828: MatDestroy(&AdB);
829: }
830: return(0);
831: }
833: PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
834: {
835: Mat A,B,C,D;
836: Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
837: PetscErrorCode ierr;
840: if (preuse == MAT_IGNORE_MATRIX) return(0);
841: MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
842: if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
843: MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);
844: return(0);
845: }
847: /*@
848: MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
850: Collective on Mat
852: Input Parameters:
853: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
854: - 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
856: Output Parameter:
857: - Sp - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
859: Note:
860: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
861: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
862: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
863: before forming inv(diag(A00)).
865: Sometimes users would like to provide problem-specific data in the Schur complement, usually only
866: for special row and column index sets. In that case, the user should call PetscObjectComposeFunction() to set
867: "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
868: it should call MatSchurComplementGetPmat_Basic().
870: Developer Notes: The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
871: remove redundancy and be clearer and simplier.
873: Level: advanced
875: Concepts: matrices^submatrices
877: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
878: @*/
879: PetscErrorCode MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
880: {
881: PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);
889: if (S->factortype) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
891: PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);
892: if (f) {
893: (*f)(S,preuse,Sp);
894: } else {
895: MatSchurComplementGetPmat_Basic(S,preuse,Sp);
896: }
897: return(0);
898: }
900: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
901: {
902: PetscErrorCode ierr;
903: Mat_SchurComplement *Na;
906: PetscNewLog(N,&Na);
907: N->data = (void*) Na;
909: N->ops->destroy = MatDestroy_SchurComplement;
910: N->ops->getvecs = MatCreateVecs_SchurComplement;
911: N->ops->view = MatView_SchurComplement;
912: N->ops->mult = MatMult_SchurComplement;
913: N->ops->multtranspose = MatMultTranspose_SchurComplement;
914: N->ops->multadd = MatMultAdd_SchurComplement;
915: N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
916: N->assembled = PETSC_FALSE;
917: N->preallocated = PETSC_FALSE;
919: KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
920: PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
921: return(0);
922: }
924: static PetscBool KSPMatRegisterAllCalled;
926: /*@C
927: KSPMatRegisterAll - Registers all matrix implementations in the KSP package.
929: Not Collective
931: Level: advanced
933: .keywords: Mat, KSP, register, all
935: .seealso: MatRegisterAll(), KSPInitializePackage()
936: @*/
937: PetscErrorCode KSPMatRegisterAll(void)
938: {
942: if (KSPMatRegisterAllCalled) return(0);
943: KSPMatRegisterAllCalled = PETSC_TRUE;
944: MatRegister(MATSCHURCOMPLEMENT,MatCreate_SchurComplement);
945: return(0);
946: }