Actual source code: schurm.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/matimpl.h>
3: #include <petscksp.h> /*I "petscksp.h" I*/
4: const char *const MatSchurComplementAinvTypes[] = {"DIAG","LUMP","MatSchurComplementAinvType","MAT_SCHUR_COMPLEMENT_AINV_",0};
6: typedef struct {
7: Mat A,Ap,B,C,D;
8: KSP ksp;
9: Vec work1,work2;
10: MatSchurComplementAinvType ainvtype;
11: } Mat_SchurComplement;
17: PetscErrorCode MatCreateVecs_SchurComplement(Mat N,Vec *right,Vec *left)
18: {
19: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
20: PetscErrorCode ierr;
23: if (Na->D) {
24: MatCreateVecs(Na->D,right,left);
25: return(0);
26: }
27: if (right) {
28: MatCreateVecs(Na->B,right,NULL);
29: }
30: if (left) {
31: MatCreateVecs(Na->C,NULL,left);
32: }
33: return(0);
34: }
38: PetscErrorCode MatView_SchurComplement(Mat N,PetscViewer viewer)
39: {
40: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
41: PetscErrorCode ierr;
44: PetscViewerASCIIPrintf(viewer,"Schur complement A11 - A10 inv(A00) A01\n");
45: if (Na->D) {
46: PetscViewerASCIIPrintf(viewer,"A11\n");
47: PetscViewerASCIIPushTab(viewer);
48: MatView(Na->D,viewer);
49: PetscViewerASCIIPopTab(viewer);
50: } else {
51: PetscViewerASCIIPrintf(viewer,"A11 = 0\n");
52: }
53: PetscViewerASCIIPrintf(viewer,"A10\n");
54: PetscViewerASCIIPushTab(viewer);
55: MatView(Na->C,viewer);
56: PetscViewerASCIIPopTab(viewer);
57: PetscViewerASCIIPrintf(viewer,"KSP of A00\n");
58: PetscViewerASCIIPushTab(viewer);
59: KSPView(Na->ksp,viewer);
60: PetscViewerASCIIPopTab(viewer);
61: PetscViewerASCIIPrintf(viewer,"A01\n");
62: PetscViewerASCIIPushTab(viewer);
63: MatView(Na->B,viewer);
64: PetscViewerASCIIPopTab(viewer);
65: return(0);
66: }
68: /*
69: A11^T - A01^T ksptrans(A00,Ap00) A10^T
70: */
73: PetscErrorCode MatMultTranspose_SchurComplement(Mat N,Vec x,Vec y)
74: {
75: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
76: PetscErrorCode ierr;
79: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
80: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
81: MatMultTranspose(Na->C,x,Na->work1);
82: KSPSolveTranspose(Na->ksp,Na->work1,Na->work2);
83: MatMultTranspose(Na->B,Na->work2,y);
84: VecScale(y,-1.0);
85: if (Na->D) {
86: MatMultTransposeAdd(Na->D,x,y,y);
87: }
88: return(0);
89: }
91: /*
92: A11 - A10 ksp(A00,Ap00) A01
93: */
96: PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
97: {
98: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
99: PetscErrorCode ierr;
102: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
103: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
104: MatMult(Na->B,x,Na->work1);
105: KSPSolve(Na->ksp,Na->work1,Na->work2);
106: MatMult(Na->C,Na->work2,y);
107: VecScale(y,-1.0);
108: if (Na->D) {
109: MatMultAdd(Na->D,x,y,y);
110: }
111: return(0);
112: }
114: /*
115: A11 - A10 ksp(A00,Ap00) A01
116: */
119: PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
120: {
121: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
122: PetscErrorCode ierr;
125: if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
126: if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
127: MatMult(Na->B,x,Na->work1);
128: KSPSolve(Na->ksp,Na->work1,Na->work2);
129: if (y == z) {
130: VecScale(Na->work2,-1.0);
131: MatMultAdd(Na->C,Na->work2,z,z);
132: } else {
133: MatMult(Na->C,Na->work2,z);
134: VecAYPX(z,-1.0,y);
135: }
136: if (Na->D) {
137: MatMultAdd(Na->D,x,z,z);
138: }
139: return(0);
140: }
144: PetscErrorCode MatSetFromOptions_SchurComplement(PetscOptionItems *PetscOptionsObject,Mat N)
145: {
146: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
147: PetscErrorCode ierr;
150: PetscOptionsHead(PetscOptionsObject,"MatSchurComplementOptions");
151: Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
152: 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);
153: PetscOptionsTail();
154: KSPSetFromOptions(Na->ksp);
155: return(0);
156: }
160: PetscErrorCode MatDestroy_SchurComplement(Mat N)
161: {
162: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
163: PetscErrorCode ierr;
166: MatDestroy(&Na->A);
167: MatDestroy(&Na->Ap);
168: MatDestroy(&Na->B);
169: MatDestroy(&Na->C);
170: MatDestroy(&Na->D);
171: VecDestroy(&Na->work1);
172: VecDestroy(&Na->work2);
173: KSPDestroy(&Na->ksp);
174: PetscFree(N->data);
175: return(0);
176: }
180: /*@
181: MatCreateSchurComplement - Creates a new matrix object that behaves like the Schur complement of a matrix
183: Collective on Mat
185: Input Parameters:
186: + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
187: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
189: Output Parameter:
190: . S - the matrix that the Schur complement S = A11 - A10 ksp(A00,Ap00) A01
192: Level: intermediate
194: Notes: The Schur complement is NOT actually formed! Rather, this
195: object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
196: for Schur complement S and a KSP solver to approximate the action of A^{-1}.
198: All four matrices must have the same MPI communicator.
200: A00 and A11 must be square matrices.
202: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement()
204: @*/
205: PetscErrorCode MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
206: {
210: KSPInitializePackage();
211: MatCreate(((PetscObject)A00)->comm,S);
212: MatSetType(*S,MATSCHURCOMPLEMENT);
213: MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);
214: return(0);
215: }
219: /*@
220: MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
222: Collective on Mat
224: Input Parameter:
225: + S - matrix obtained with MatCreateSchurComplement (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
226: . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
227: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
229: Level: intermediate
231: Notes: The Schur complement is NOT actually formed! Rather, this
232: object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
233: for Schur complement S and a KSP solver to approximate the action of A^{-1}.
235: All four matrices must have the same MPI communicator.
237: A00 and A11 must be square matrices.
239: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatCreateSchurComplement(), MatGetSchurComplement()
241: @*/
242: PetscErrorCode MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
243: {
244: PetscErrorCode ierr;
245: PetscInt m,n;
246: Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
249: if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
257: 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);
258: 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);
259: 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);
260: 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);
261: 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);
262: if (A11) {
265: 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);
266: }
268: MatGetLocalSize(A01,NULL,&n);
269: MatGetLocalSize(A10,&m,NULL);
270: MatSetSizes(S,m,n,PETSC_DECIDE,PETSC_DECIDE);
271: PetscObjectReference((PetscObject)A00);
272: PetscObjectReference((PetscObject)Ap00);
273: PetscObjectReference((PetscObject)A01);
274: PetscObjectReference((PetscObject)A10);
275: Na->A = A00;
276: Na->Ap = Ap00;
277: Na->B = A01;
278: Na->C = A10;
279: Na->D = A11;
280: if (A11) {
281: PetscObjectReference((PetscObject)A11);
282: }
283: S->assembled = PETSC_TRUE;
284: S->preallocated = PETSC_TRUE;
286: PetscLayoutSetUp((S)->rmap);
287: PetscLayoutSetUp((S)->cmap);
288: KSPSetOperators(Na->ksp,A00,Ap00);
289: return(0);
290: }
294: /*@
295: MatSchurComplementGetKSP - Gets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
297: Not Collective
299: Input Parameter:
300: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
302: Output Parameter:
303: . ksp - the linear solver object
305: Options Database:
306: . -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.
308: Level: intermediate
310: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
311: @*/
312: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
313: {
314: Mat_SchurComplement *Na;
319: Na = (Mat_SchurComplement*) S->data;
320: *ksp = Na->ksp;
321: return(0);
322: }
326: /*@
327: MatSchurComplementSetKSP - Sets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
329: Not Collective
331: Input Parameters:
332: + S - matrix created with MatCreateSchurComplement()
333: - ksp - the linear solver object
335: Level: developer
337: Developer Notes:
338: This is used in PCFieldSplit to reuse the 0-split KSP to implement ksp(A00,Ap00) in S.
340: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
341: @*/
342: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
343: {
344: Mat_SchurComplement *Na;
345: PetscErrorCode ierr;
350: Na = (Mat_SchurComplement*) S->data;
351: PetscObjectReference((PetscObject)ksp);
352: KSPDestroy(&Na->ksp);
353: Na->ksp = ksp;
354: KSPSetOperators(Na->ksp, Na->A, Na->Ap);
355: return(0);
356: }
360: /*@
361: MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
363: Collective on Mat
365: Input Parameters:
366: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
367: . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
368: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
370: Level: intermediate
372: Notes: All four matrices must have the same MPI communicator
374: A00 and A11 must be square matrices
376: All of the matrices provided must have the same sizes as was used with MatCreateSchurComplement() or MatSchurComplementSetSubMatrices()
377: though they need not be the same matrices.
379: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement()
381: @*/
382: PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
383: {
384: PetscErrorCode ierr;
385: Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
388: if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
395: 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);
396: 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);
397: 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);
398: 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);
399: 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);
400: if (A11) {
403: 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);
404: }
406: PetscObjectReference((PetscObject)A00);
407: PetscObjectReference((PetscObject)Ap00);
408: PetscObjectReference((PetscObject)A01);
409: PetscObjectReference((PetscObject)A10);
410: if (A11) {
411: PetscObjectReference((PetscObject)A11);
412: }
414: MatDestroy(&Na->A);
415: MatDestroy(&Na->Ap);
416: MatDestroy(&Na->B);
417: MatDestroy(&Na->C);
418: MatDestroy(&Na->D);
420: Na->A = A00;
421: Na->Ap = Ap00;
422: Na->B = A01;
423: Na->C = A10;
424: Na->D = A11;
426: KSPSetOperators(Na->ksp,A00,Ap00);
427: return(0);
428: }
433: /*@C
434: MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
436: Collective on Mat
438: Input Parameter:
439: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
441: Output Paramters:
442: + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
443: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
445: 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.
447: Level: intermediate
449: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
450: @*/
451: PetscErrorCode MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
452: {
453: Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
454: PetscErrorCode ierr;
455: PetscBool flg;
459: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);
460: if (flg) {
461: if (A00) *A00 = Na->A;
462: if (Ap00) *Ap00 = Na->Ap;
463: if (A01) *A01 = Na->B;
464: if (A10) *A10 = Na->C;
465: if (A11) *A11 = Na->D;
466: } else {
467: if (A00) *A00 = 0;
468: if (Ap00) *Ap00 = 0;
469: if (A01) *A01 = 0;
470: if (A10) *A10 = 0;
471: if (A11) *A11 = 0;
472: }
473: return(0);
474: }
476: #include <petsc/private/kspimpl.h>
480: /*@
481: MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
483: Collective on Mat
485: Input Parameter:
486: . M - the matrix obtained with MatCreateSchurComplement()
488: Output Parameter:
489: . S - the Schur complement matrix
491: Note: This can be expensive, so it is mainly for testing
493: Level: advanced
495: .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
496: @*/
497: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
498: {
499: Mat B, C, D;
500: KSP ksp;
501: PC pc;
502: PetscBool isLU, isILU;
503: PetscReal fill = 2.0;
507: MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);
508: MatSchurComplementGetKSP(M, &ksp);
509: KSPGetPC(ksp, &pc);
510: PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);
511: PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);
512: if (isLU || isILU) {
513: Mat fact, Bd, AinvB, AinvBd;
514: PetscReal eps = 1.0e-10;
516: /* This can be sped up for banded LU */
517: KSPSetUp(ksp);
518: PCFactorGetMatrix(pc, &fact);
519: MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);
520: MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);
521: MatMatSolve(fact, Bd, AinvBd);
522: MatDestroy(&Bd);
523: MatChop(AinvBd, eps);
524: MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);
525: MatDestroy(&AinvBd);
526: MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);
527: MatDestroy(&AinvB);
528: } else {
529: Mat Ainvd, Ainv;
531: PCComputeExplicitOperator(pc, &Ainvd);
532: MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);
533: MatDestroy(&Ainvd);
534: #if 0
535: /* Symmetric version */
536: MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);
537: #else
538: /* Nonsymmetric version */
539: MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);
540: #endif
541: MatDestroy(&Ainv);
542: }
544: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);
545: MatView(*S, PETSC_VIEWER_STDOUT_WORLD);
546: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
548: if (D) {
549: MatInfo info;
550: PetscReal norm;
552: MatGetInfo(D, MAT_GLOBAL_SUM, &info);
553: if (info.nz_used) {
554: MatNorm(D, NORM_INFINITY, &norm);
555: if (norm > PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented for Schur complements with non-vanishing D");
556: }
557: }
558: return(0);
559: }
563: /* Developer Notes: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
564: PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *newpmat)
565: {
567: Mat A=0,Ap=0,B=0,C=0,D=0;
568: MatReuse reuse;
577: if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) return(0);
581: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
583: reuse = MAT_INITIAL_MATRIX;
584: if (mreuse == MAT_REUSE_MATRIX) {
585: MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);
586: if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
587: if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
588: MatDestroy(&Ap); /* get rid of extra reference */
589: reuse = MAT_REUSE_MATRIX;
590: }
591: MatGetSubMatrix(mat,isrow0,iscol0,reuse,&A);
592: MatGetSubMatrix(mat,isrow0,iscol1,reuse,&B);
593: MatGetSubMatrix(mat,isrow1,iscol0,reuse,&C);
594: MatGetSubMatrix(mat,isrow1,iscol1,reuse,&D);
595: switch (mreuse) {
596: case MAT_INITIAL_MATRIX:
597: MatCreateSchurComplement(A,A,B,C,D,newmat);
598: break;
599: case MAT_REUSE_MATRIX:
600: MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D);
601: break;
602: default:
603: if (mreuse != MAT_IGNORE_MATRIX) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse");
604: }
605: if (preuse != MAT_IGNORE_MATRIX) {
606: MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);
607: }
608: MatDestroy(&A);
609: MatDestroy(&B);
610: MatDestroy(&C);
611: MatDestroy(&D);
612: return(0);
613: }
617: /*@
618: MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
620: Collective on Mat
622: Input Parameters:
623: + A - matrix in which the complement is to be taken
624: . isrow0 - rows to eliminate
625: . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
626: . isrow1 - rows in which the Schur complement is formed
627: . iscol1 - columns in which the Schur complement is formed
628: . mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in S
629: . plump - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
630: MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP
631: - preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in Sp
633: Output Parameters:
634: + S - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
635: - Sp - approximate Schur complement suitable for preconditioning
637: Note:
638: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
639: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
640: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
641: before forming inv(diag(A00)).
643: Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
644: 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
645: "MatGetSchurComplement_C" to their function. If their function needs to fall back to the default implementation, it
646: should call MatGetSchurComplement_Basic().
648: Level: advanced
650: Concepts: matrices^submatrices
652: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
653: @*/
654: PetscErrorCode MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
655: {
656: PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*) = NULL;
667: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
668: f = NULL;
669: 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. */
670: PetscObjectQueryFunction((PetscObject)*S,"MatGetSchurComplement_C",&f);
671: }
672: if (f) {
673: (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
674: } else {
675: MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
676: }
677: return(0);
678: }
682: /*@
683: MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
685: Not collective.
687: Input Parameters:
688: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
689: - ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
690: MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP
692: Options database:
693: -mat_schur_complement_ainv_type diag | lump
695: Note:
696: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
697: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
698: the (0,0) block A00 in place of A00^{-1}. This rarely produces a scalable algorithm. Optionally, A00 can be lumped
699: before forming inv(diag(A00)).
701: Level: advanced
703: Concepts: matrices^submatrices
705: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
706: @*/
707: PetscErrorCode MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
708: {
709: PetscErrorCode ierr;
710: const char* t;
711: PetscBool isschur;
712: Mat_SchurComplement *schur;
716: PetscObjectGetType((PetscObject)S,&t);
717: PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
718: if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
719: schur = (Mat_SchurComplement*)S->data;
720: if (ainvtype != MAT_SCHUR_COMPLEMENT_AINV_DIAG && ainvtype != MAT_SCHUR_COMPLEMENT_AINV_LUMP) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D",ainvtype);
721: schur->ainvtype = ainvtype;
722: return(0);
723: }
727: /*@
728: MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
730: Not collective.
732: Input Parameter:
733: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
735: Output Parameter:
736: . ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
737: MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP
739: Note:
740: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
741: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
742: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
743: before forming inv(diag(A00)).
745: Level: advanced
747: Concepts: matrices^submatrices
749: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
750: @*/
751: PetscErrorCode MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
752: {
753: PetscErrorCode ierr;
754: const char* t;
755: PetscBool isschur;
756: Mat_SchurComplement *schur;
760: PetscObjectGetType((PetscObject)S,&t);
761: PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
762: if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
763: schur = (Mat_SchurComplement*)S->data;
764: if (ainvtype) *ainvtype = schur->ainvtype;
765: return(0);
766: }
770: /*@
771: MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
773: Collective on Mat
775: Input Parameters:
776: + 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)
777: . ainvtype - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
778: - 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
780: Output Parameter:
781: - Spmat - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
783: Note:
784: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
785: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
786: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
787: before forming inv(diag(A00)).
789: Level: advanced
791: Concepts: matrices^submatrices
793: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
794: @*/
795: PetscErrorCode MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
796: {
799: PetscInt N00;
802: /* 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. */
803: /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
804: if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Spmat: A01, A10 and A11 are all NULL.");
806: if (preuse == MAT_IGNORE_MATRIX) return(0);
808: /* A zero size A00 or empty A01 or A10 imply S = A11. */
809: MatGetSize(A00,&N00,NULL);
810: if (!A01 || !A10 || !N00) {
811: if (preuse == MAT_INITIAL_MATRIX) {
812: MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
813: } else { /* MAT_REUSE_MATRIX */
814: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
815: MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
816: }
818: } else {
819: Mat AdB;
820: Vec diag;
822: MatCreateVecs(A00,&diag,NULL);
823: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
824: MatGetRowSum(A00,diag);
825: } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
826: MatGetDiagonal(A00,diag);
827: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
828: VecReciprocal(diag);
829: MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
830: MatDiagonalScale(AdB,diag,NULL);
831: VecDestroy(&diag);
832: /* Cannot really reuse Spmat in MatMatMult() because of MatAYPX() -->
833: MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult() */
834: MatDestroy(Spmat);
835: MatMatMult(A10,AdB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,Spmat);
836: if (!A11) {
837: MatScale(*Spmat,-1.0);
838: } else {
839: /* TODO: when can we pass SAME_NONZERO_PATTERN? */
840: MatAYPX(*Spmat,-1,A11,DIFFERENT_NONZERO_PATTERN);
841: }
842: MatDestroy(&AdB);
843: }
844: return(0);
845: }
849: PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
850: {
851: Mat A,B,C,D;
852: Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
853: PetscErrorCode ierr;
856: if (preuse == MAT_IGNORE_MATRIX) return(0);
858: MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
859: if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
860: MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);
861: return(0);
862: }
866: /*@
867: MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
869: Collective on Mat
871: Input Parameters:
872: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
873: - 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
875: Output Parameter:
876: - Sp - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
878: Note:
879: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
880: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
881: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
882: before forming inv(diag(A00)).
884: Sometimes users would like to provide problem-specific data in the Schur complement, usually only
885: for special row and column index sets. In that case, the user should call PetscObjectComposeFunction() to set
886: "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
887: it should call MatSchurComplementGetPmat_Basic().
889: Level: advanced
891: Concepts: matrices^submatrices
893: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
894: @*/
895: PetscErrorCode MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
896: {
897: PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);
903: if (S->factortype) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
905: PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);
906: if (f) {
907: (*f)(S,preuse,Sp);
908: } else {
909: MatSchurComplementGetPmat_Basic(S,preuse,Sp);
910: }
911: return(0);
912: }
916: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
917: {
918: PetscErrorCode ierr;
919: Mat_SchurComplement *Na;
922: PetscNewLog(N,&Na);
923: N->data = (void*) Na;
925: N->ops->destroy = MatDestroy_SchurComplement;
926: N->ops->getvecs = MatCreateVecs_SchurComplement;
927: N->ops->view = MatView_SchurComplement;
928: N->ops->mult = MatMult_SchurComplement;
929: N->ops->multtranspose = MatMultTranspose_SchurComplement;
930: N->ops->multadd = MatMultAdd_SchurComplement;
931: N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
932: N->assembled = PETSC_FALSE;
933: N->preallocated = PETSC_FALSE;
935: KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
936: PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
937: return(0);
938: }
940: static PetscBool KSPMatRegisterAllCalled;
944: /*@C
945: KSPMatRegisterAll - Registers all matrix implementations in the KSP package.
947: Not Collective
949: Level: advanced
951: .keywords: Mat, KSP, register, all
953: .seealso: MatRegisterAll(), KSPInitializePackage()
954: @*/
955: PetscErrorCode KSPMatRegisterAll()
956: {
960: if (KSPMatRegisterAllCalled) return(0);
961: KSPMatRegisterAllCalled = PETSC_TRUE;
962: MatRegister(MATSCHURCOMPLEMENT,MatCreate_SchurComplement);
963: return(0);
964: }