Actual source code: schurm.c
petsc-3.10.5 2019-03-28
1: #include <../src/ksp/ksp/utils/schurm/schurm.h>
3: const char *const MatSchurComplementAinvTypes[] = {"DIAG","LUMP","BLOCKDIAG","MatSchurComplementAinvType","MAT_SCHUR_COMPLEMENT_AINV_",0};
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 inv(A00) used when assembling Sp = A11 - A10 inv(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 matrix object that behaves like the Schur complement of a matrix
157: Collective on Mat
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 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}.
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 true Schur complement (useful for building a preconditioner for the Schur complement).
180: MatSchurComplementGetPmat() can be called on the output of this function to generate an explicit approximation to the Schur complement.
182: Developer Notes:
183: The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
184: remove redundancy and be clearer and simplier.
187: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement(),
188: MatSchurComplementGetPmat()
190: @*/
191: PetscErrorCode MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
192: {
196: KSPInitializePackage();
197: MatCreate(PetscObjectComm((PetscObject)A00),S);
198: MatSetType(*S,MATSCHURCOMPLEMENT);
199: MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);
200: return(0);
201: }
203: /*@
204: MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
206: Collective on Mat
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}.
213: Level: intermediate
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}.
220: All four matrices must have the same MPI communicator.
222: A00 and A11 must be square matrices.
224: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatCreateSchurComplement(), MatGetSchurComplement()
226: @*/
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;
234: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
235: if (!isschur) return(0);
236: if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
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) {
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: }
255: MatSetSizes(S,A10->rmap->n,A01->cmap->n,A10->rmap->N,A01->cmap->N);
256: PetscObjectReference((PetscObject)A00);
257: PetscObjectReference((PetscObject)Ap00);
258: PetscObjectReference((PetscObject)A01);
259: PetscObjectReference((PetscObject)A10);
260: Na->A = A00;
261: Na->Ap = Ap00;
262: Na->B = A01;
263: Na->C = A10;
264: Na->D = A11;
265: if (A11) {
266: PetscObjectReference((PetscObject)A11);
267: }
268: MatSetUp(S);
269: KSPSetOperators(Na->ksp,A00,Ap00);
270: S->assembled = PETSC_TRUE;
271: return(0);
272: }
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
277: Not Collective
279: Input Parameter:
280: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
282: Output Parameter:
283: . ksp - the linear solver object
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.
288: Level: intermediate
290: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
291: @*/
292: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
293: {
294: Mat_SchurComplement *Na;
295: PetscBool isschur;
296: PetscErrorCode ierr;
300: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
301: if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
303: Na = (Mat_SchurComplement*) S->data;
304: *ksp = Na->ksp;
305: return(0);
306: }
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
311: Not Collective
313: Input Parameters:
314: + S - matrix created with MatCreateSchurComplement()
315: - ksp - the linear solver object
317: Level: developer
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.
323: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
324: @*/
325: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
326: {
327: Mat_SchurComplement *Na;
328: PetscErrorCode ierr;
329: PetscBool isschur;
333: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
334: if (!isschur) return(0);
336: Na = (Mat_SchurComplement*) S->data;
337: PetscObjectReference((PetscObject)ksp);
338: KSPDestroy(&Na->ksp);
339: Na->ksp = ksp;
340: KSPSetOperators(Na->ksp, Na->A, Na->Ap);
341: return(0);
342: }
344: /*@
345: MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
347: Collective on Mat
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}.
354: Level: intermediate
356: Notes:
357: All four matrices must have the same MPI communicator
359: A00 and A11 must be square matrices
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.
364: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement()
366: @*/
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;
375: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
376: if (!isschur) return(0);
377: if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
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) {
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: }
396: PetscObjectReference((PetscObject)A00);
397: PetscObjectReference((PetscObject)Ap00);
398: PetscObjectReference((PetscObject)A01);
399: PetscObjectReference((PetscObject)A10);
400: if (A11) {
401: PetscObjectReference((PetscObject)A11);
402: }
404: MatDestroy(&Na->A);
405: MatDestroy(&Na->Ap);
406: MatDestroy(&Na->B);
407: MatDestroy(&Na->C);
408: MatDestroy(&Na->D);
410: Na->A = A00;
411: Na->Ap = Ap00;
412: Na->B = A01;
413: Na->C = A10;
414: Na->D = A11;
416: KSPSetOperators(Na->ksp,A00,Ap00);
417: return(0);
418: }
421: /*@C
422: MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
424: Collective on Mat
426: Input Parameter:
427: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
429: Output Paramters:
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}.
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.
435: Level: intermediate
437: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
438: @*/
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;
447: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);
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: return(0);
455: }
457: #include <petsc/private/kspimpl.h>
459: /*@
460: MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
462: Collective on Mat
464: Input Parameter:
465: . M - the matrix obtained with MatCreateSchurComplement()
467: Output Parameter:
468: . S - the Schur complement matrix
470: Note: This can be expensive, so it is mainly for testing
472: Level: advanced
474: .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
475: @*/
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;
486: MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);
487: MatSchurComplementGetKSP(M, &ksp);
488: KSPGetPC(ksp, &pc);
489: PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);
490: PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);
491: if (isLU || isILU) {
492: Mat fact, Bd, AinvB, AinvBd;
493: PetscReal eps = 1.0e-10;
495: /* This can be sped up for banded LU */
496: KSPSetUp(ksp);
497: PCFactorGetMatrix(pc, &fact);
498: MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);
499: MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);
500: MatMatSolve(fact, Bd, AinvBd);
501: MatDestroy(&Bd);
502: MatChop(AinvBd, eps);
503: MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);
504: MatDestroy(&AinvBd);
505: MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);
506: MatDestroy(&AinvB);
507: } else {
508: Mat Ainvd, Ainv;
510: PCComputeExplicitOperator(pc, &Ainvd);
511: MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);
512: MatDestroy(&Ainvd);
513: #if 0
514: /* Symmetric version */
515: MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);
516: #else
517: /* Nonsymmetric version */
518: MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);
519: #endif
520: MatDestroy(&Ainv);
521: }
522: if (D) {
523: MatAXPY(*S, -1.0, D, DIFFERENT_NONZERO_PATTERN);
524: }
525: MatScale(*S,-1.0);
526: return(0);
527: }
529: /* Developer Notes:
530: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
531: PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *newpmat)
532: {
534: Mat A=0,Ap=0,B=0,C=0,D=0;
535: MatReuse reuse;
547: if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) return(0);
551: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
553: reuse = MAT_INITIAL_MATRIX;
554: if (mreuse == MAT_REUSE_MATRIX) {
555: MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);
556: if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
557: if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
558: MatDestroy(&Ap); /* get rid of extra reference */
559: reuse = MAT_REUSE_MATRIX;
560: }
561: MatCreateSubMatrix(mat,isrow0,iscol0,reuse,&A);
562: MatCreateSubMatrix(mat,isrow0,iscol1,reuse,&B);
563: MatCreateSubMatrix(mat,isrow1,iscol0,reuse,&C);
564: MatCreateSubMatrix(mat,isrow1,iscol1,reuse,&D);
565: switch (mreuse) {
566: case MAT_INITIAL_MATRIX:
567: MatCreateSchurComplement(A,A,B,C,D,newmat);
568: break;
569: case MAT_REUSE_MATRIX:
570: MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D);
571: break;
572: default:
573: if (mreuse != MAT_IGNORE_MATRIX) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse %d",(int)mreuse);
574: }
575: if (preuse != MAT_IGNORE_MATRIX) {
576: MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);
577: }
578: MatDestroy(&A);
579: MatDestroy(&B);
580: MatDestroy(&C);
581: MatDestroy(&D);
582: return(0);
583: }
585: /*@
586: MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
588: Collective on Mat
590: Input Parameters:
591: + A - matrix in which the complement is to be taken
592: . isrow0 - rows to eliminate
593: . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
594: . isrow1 - rows in which the Schur complement is formed
595: . iscol1 - columns in which the Schur complement is formed
596: . mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in S
597: . ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
598: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG, or MAT_SCHUR_COMPLEMENT_AINV_LUMP
599: - preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in Sp
601: Output Parameters:
602: + S - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
603: - Sp - approximate Schur complement from which a preconditioner can be built
605: Note:
606: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
607: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
608: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
609: before forming inv(diag(A00)).
611: Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
612: 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
613: "MatGetSchurComplement_C" to their function. If their function needs to fall back to the default implementation, it
614: should call MatGetSchurComplement_Basic().
616: MatCreateSchurComplement() takes as arguments the four submatrices and returns the virtual Schur complement (what this returns in S).
618: MatSchurComplementGetPmat() takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).
620: In other words calling MatCreateSchurComplement() followed by MatSchurComplementGetPmat() produces the same output as this function but with slightly different
621: inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.
623: Developer Notes:
624: The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
625: remove redundancy and be clearer and simplier.
627: Level: advanced
629: Concepts: matrices^submatrices
631: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
632: @*/
633: PetscErrorCode MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
634: {
635: PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*) = NULL;
649: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
650: f = NULL;
651: 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. */
652: PetscObjectQueryFunction((PetscObject)*S,"MatGetSchurComplement_C",&f);
653: }
654: if (f) {
655: (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
656: } else {
657: MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
658: }
659: return(0);
660: }
662: /*@
663: MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
665: Not collective.
667: Input Parameters:
668: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
669: - ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
670: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG
672: Options database:
673: -mat_schur_complement_ainv_type diag | lump | blockdiag
675: Note:
676: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
677: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
678: the (0,0) block A00 in place of A00^{-1}. This rarely produces a scalable algorithm. Optionally, A00 can be lumped
679: before forming inv(diag(A00)).
681: Level: advanced
683: Concepts: matrices^submatrices
685: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
686: @*/
687: PetscErrorCode MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
688: {
689: PetscErrorCode ierr;
690: PetscBool isschur;
691: Mat_SchurComplement *schur;
695: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
696: if (!isschur) return(0);
698: schur = (Mat_SchurComplement*)S->data;
699: 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);
700: schur->ainvtype = ainvtype;
701: return(0);
702: }
704: /*@
705: MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
707: Not collective.
709: Input Parameter:
710: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
712: Output Parameter:
713: . ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
714: MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG
716: Note:
717: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
718: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
719: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
720: before forming inv(diag(A00)).
722: Level: advanced
724: Concepts: matrices^submatrices
726: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
727: @*/
728: PetscErrorCode MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
729: {
730: PetscErrorCode ierr;
731: PetscBool isschur;
732: Mat_SchurComplement *schur;
736: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
737: if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
738: schur = (Mat_SchurComplement*)S->data;
739: if (ainvtype) *ainvtype = schur->ainvtype;
740: return(0);
741: }
743: /*@
744: MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
746: Collective on Mat
748: Input Parameters:
749: + 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)
750: . ainvtype - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
751: - 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
753: Output Parameter:
754: - Spmat - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
756: Note:
757: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
758: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
759: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
760: before forming inv(diag(A00)).
762: Level: advanced
764: Concepts: matrices^submatrices
766: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
767: @*/
768: PetscErrorCode MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
769: {
771: PetscInt N00;
774: /* 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. */
775: /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
776: if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Spmat: A01, A10 and A11 are all NULL.");
778: if (preuse == MAT_IGNORE_MATRIX) return(0);
780: /* A zero size A00 or empty A01 or A10 imply S = A11. */
781: MatGetSize(A00,&N00,NULL);
782: if (!A01 || !A10 || !N00) {
783: if (preuse == MAT_INITIAL_MATRIX) {
784: MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
785: } else { /* MAT_REUSE_MATRIX */
786: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
787: MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
788: }
789: } else {
790: Mat AdB;
791: Vec diag;
793: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
794: MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
795: MatCreateVecs(A00,&diag,NULL);
796: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
797: MatGetRowSum(A00,diag);
798: } else {
799: MatGetDiagonal(A00,diag);
800: }
801: VecReciprocal(diag);
802: MatDiagonalScale(AdB,diag,NULL);
803: VecDestroy(&diag);
804: } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
805: Mat A00_inv;
806: MatType type;
807: MPI_Comm comm;
809: PetscObjectGetComm((PetscObject)A00,&comm);
810: MatGetType(A00,&type);
811: MatCreate(comm,&A00_inv);
812: MatSetType(A00_inv,type);
813: MatInvertBlockDiagonalMat(A00,A00_inv);
814: MatMatMult(A00_inv,A01,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AdB);
815: MatDestroy(&A00_inv);
816: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
817: /* Cannot really reuse Spmat in MatMatMult() because of MatAYPX() -->
818: MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult() */
819: MatDestroy(Spmat);
820: MatMatMult(A10,AdB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,Spmat);
821: if (!A11) {
822: MatScale(*Spmat,-1.0);
823: } else {
824: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
825: MatAYPX(*Spmat,-1,A11,DIFFERENT_NONZERO_PATTERN);
826: }
827: MatDestroy(&AdB);
828: }
829: return(0);
830: }
832: PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
833: {
834: Mat A,B,C,D;
835: Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
836: PetscErrorCode ierr;
839: if (preuse == MAT_IGNORE_MATRIX) return(0);
840: MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
841: if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
842: MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);
843: return(0);
844: }
846: /*@
847: MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
849: Collective on Mat
851: Input Parameters:
852: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
853: - 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
855: Output Parameter:
856: - Sp - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
858: Note:
859: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
860: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
861: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
862: before forming inv(diag(A00)).
864: Sometimes users would like to provide problem-specific data in the Schur complement, usually only
865: for special row and column index sets. In that case, the user should call PetscObjectComposeFunction() to set
866: "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
867: it should call MatSchurComplementGetPmat_Basic().
869: Developer Notes:
870: 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: }