Actual source code: schurm.c
petsc-3.5.4 2015-05-23
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 MatGetVecs_SchurComplement(Mat N,Vec *right,Vec *left)
18: {
19: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
20: PetscErrorCode ierr;
23: if (Na->D) {
24: MatGetVecs(Na->D,right,left);
25: return(0);
26: }
27: if (right) {
28: MatGetVecs(Na->B,right,NULL);
29: }
30: if (left) {
31: MatGetVecs(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: }
69: /*
70: A11 - A10 ksp(A00,Ap00) A01
71: */
74: PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
75: {
76: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
77: PetscErrorCode ierr;
80: if (!Na->work1) {MatGetVecs(Na->A,&Na->work1,NULL);}
81: if (!Na->work2) {MatGetVecs(Na->A,&Na->work2,NULL);}
82: MatMult(Na->B,x,Na->work1);
83: KSPSolve(Na->ksp,Na->work1,Na->work2);
84: MatMult(Na->C,Na->work2,y);
85: VecScale(y,-1.0);
86: if (Na->D) {
87: MatMultAdd(Na->D,x,y,y);
88: }
89: return(0);
90: }
92: /*
93: A11 - A10 ksp(A00,Ap00) A01
94: */
97: PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
98: {
99: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
100: PetscErrorCode ierr;
103: if (!Na->work1) {MatGetVecs(Na->A,&Na->work1,NULL);}
104: if (!Na->work2) {MatGetVecs(Na->A,&Na->work2,NULL);}
105: MatMult(Na->B,x,Na->work1);
106: KSPSolve(Na->ksp,Na->work1,Na->work2);
107: if (y == z) {
108: VecScale(Na->work2,-1.0);
109: MatMultAdd(Na->C,Na->work2,z,z);
110: } else {
111: MatMult(Na->C,Na->work2,z);
112: VecAYPX(z,-1.0,y);
113: }
114: if (Na->D) {
115: MatMultAdd(Na->D,x,z,z);
116: }
117: return(0);
118: }
122: PetscErrorCode MatSetFromOptions_SchurComplement(Mat N)
123: {
124: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
125: PetscErrorCode ierr;
128: PetscOptionsHead("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: }
138: PetscErrorCode MatDestroy_SchurComplement(Mat N)
139: {
140: Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
141: PetscErrorCode ierr;
144: MatDestroy(&Na->A);
145: MatDestroy(&Na->Ap);
146: MatDestroy(&Na->B);
147: MatDestroy(&Na->C);
148: MatDestroy(&Na->D);
149: VecDestroy(&Na->work1);
150: VecDestroy(&Na->work2);
151: KSPDestroy(&Na->ksp);
152: PetscFree(N->data);
153: return(0);
154: }
158: /*@
159: MatCreateSchurComplement - Creates a new matrix object that behaves like the Schur complement of a matrix
161: Collective on Mat
163: Input Parameters:
164: + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
165: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
167: Output Parameter:
168: . S - the matrix that the Schur complement S = A11 - A10 ksp(A00,Ap00) A01
170: Level: intermediate
172: Notes: The Schur complement is NOT actually formed! Rather, this
173: object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
174: for Schur complement S and a KSP solver to approximate the action of A^{-1}.
176: All four matrices must have the same MPI communicator.
178: A00 and A11 must be square matrices.
180: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement()
182: @*/
183: PetscErrorCode MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
184: {
188: KSPInitializePackage();
189: MatCreate(((PetscObject)A00)->comm,S);
190: MatSetType(*S,MATSCHURCOMPLEMENT);
191: MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);
192: return(0);
193: }
197: /*@
198: MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement
200: Collective on Mat
202: Input Parameter:
203: + S - matrix obtained with MatCreateSchurComplement (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
204: . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
205: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
207: Level: intermediate
209: Notes: The Schur complement is NOT actually formed! Rather, this
210: object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
211: for Schur complement S and a KSP solver to approximate the action of A^{-1}.
213: All four matrices must have the same MPI communicator.
215: A00 and A11 must be square matrices.
217: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatCreateSchurComplement(), MatGetSchurComplement()
219: @*/
220: PetscErrorCode MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
221: {
222: PetscErrorCode ierr;
223: PetscInt m,n;
224: Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
227: if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
235: 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);
236: 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);
237: 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);
238: 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);
239: 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);
240: if (A11) {
243: 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);
244: }
246: MatGetLocalSize(A01,NULL,&n);
247: MatGetLocalSize(A10,&m,NULL);
248: MatSetSizes(S,m,n,PETSC_DECIDE,PETSC_DECIDE);
249: PetscObjectReference((PetscObject)A00);
250: PetscObjectReference((PetscObject)Ap00);
251: PetscObjectReference((PetscObject)A01);
252: PetscObjectReference((PetscObject)A10);
253: Na->A = A00;
254: Na->Ap = Ap00;
255: Na->B = A01;
256: Na->C = A10;
257: Na->D = A11;
258: if (A11) {
259: PetscObjectReference((PetscObject)A11);
260: }
261: S->assembled = PETSC_TRUE;
262: S->preallocated = PETSC_TRUE;
264: PetscLayoutSetUp((S)->rmap);
265: PetscLayoutSetUp((S)->cmap);
266: KSPSetOperators(Na->ksp,A00,Ap00);
267: return(0);
268: }
272: /*@
273: MatSchurComplementGetKSP - Gets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
275: Not Collective
277: Input Parameter:
278: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
280: Output Parameter:
281: . ksp - the linear solver object
283: Options Database:
284: . -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.
286: Level: intermediate
288: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
289: @*/
290: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
291: {
292: Mat_SchurComplement *Na;
297: Na = (Mat_SchurComplement*) S->data;
298: *ksp = Na->ksp;
299: return(0);
300: }
304: /*@
305: MatSchurComplementSetKSP - Sets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01
307: Not Collective
309: Input Parameters:
310: + S - matrix created with MatCreateSchurComplement()
311: - ksp - the linear solver object
313: Level: developer
315: Developer Notes:
316: This is used in PCFieldSplit to reuse the 0-split KSP to implement ksp(A00,Ap00) in S.
318: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
319: @*/
320: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
321: {
322: Mat_SchurComplement *Na;
323: PetscErrorCode ierr;
328: Na = (Mat_SchurComplement*) S->data;
329: PetscObjectReference((PetscObject)ksp);
330: KSPDestroy(&Na->ksp);
331: Na->ksp = ksp;
332: KSPSetOperators(Na->ksp, Na->A, Na->Ap);
333: return(0);
334: }
338: /*@
339: MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices
341: Collective on Mat
343: Input Parameters:
344: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
345: . A00,A01,A10,A11 - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
346: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
348: Level: intermediate
350: Notes: All four matrices must have the same MPI communicator
352: A00 and A11 must be square matrices
354: All of the matrices provided must have the same sizes as was used with MatCreateSchurComplement() or MatSchurComplementSetSubMatrices()
355: though they need not be the same matrices.
357: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement()
359: @*/
360: PetscErrorCode MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
361: {
362: PetscErrorCode ierr;
363: Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
366: if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
373: 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);
374: 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);
375: 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);
376: 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);
377: 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);
378: if (A11) {
381: 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);
382: }
384: PetscObjectReference((PetscObject)A00);
385: PetscObjectReference((PetscObject)Ap00);
386: PetscObjectReference((PetscObject)A01);
387: PetscObjectReference((PetscObject)A10);
388: if (A11) {
389: PetscObjectReference((PetscObject)A11);
390: }
392: MatDestroy(&Na->A);
393: MatDestroy(&Na->Ap);
394: MatDestroy(&Na->B);
395: MatDestroy(&Na->C);
396: MatDestroy(&Na->D);
398: Na->A = A00;
399: Na->Ap = Ap00;
400: Na->B = A01;
401: Na->C = A10;
402: Na->D = A11;
404: KSPSetOperators(Na->ksp,A00,Ap00);
405: return(0);
406: }
411: /*@C
412: MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement
414: Collective on Mat
416: Input Parameter:
417: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
419: Output Paramters:
420: + A00,A01,A10,A11 - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
421: - Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.
423: 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.
425: Level: intermediate
427: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
428: @*/
429: PetscErrorCode MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
430: {
431: Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
432: PetscErrorCode ierr;
433: PetscBool flg;
437: PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);
438: if (flg) {
439: if (A00) *A00 = Na->A;
440: if (Ap00) *Ap00 = Na->Ap;
441: if (A01) *A01 = Na->B;
442: if (A10) *A10 = Na->C;
443: if (A11) *A11 = Na->D;
444: } else {
445: if (A00) *A00 = 0;
446: if (Ap00) *Ap00 = 0;
447: if (A01) *A01 = 0;
448: if (A10) *A10 = 0;
449: if (A11) *A11 = 0;
450: }
451: return(0);
452: }
454: #include <petsc-private/kspimpl.h>
458: /*@
459: MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly
461: Collective on Mat
463: Input Parameter:
464: . M - the matrix obtained with MatCreateSchurComplement()
466: Output Parameter:
467: . S - the Schur complement matrix
469: Note: This can be expensive, so it is mainly for testing
471: Level: advanced
473: .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
474: @*/
475: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
476: {
477: Mat B, C, D;
478: KSP ksp;
479: PC pc;
480: PetscBool isLU, isILU;
481: PetscReal fill = 2.0;
485: MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);
486: MatSchurComplementGetKSP(M, &ksp);
487: KSPGetPC(ksp, &pc);
488: PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);
489: PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);
490: if (isLU || isILU) {
491: Mat fact, Bd, AinvB, AinvBd;
492: PetscReal eps = 1.0e-10;
494: /* This can be sped up for banded LU */
495: KSPSetUp(ksp);
496: PCFactorGetMatrix(pc, &fact);
497: MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);
498: MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);
499: MatMatSolve(fact, Bd, AinvBd);
500: MatDestroy(&Bd);
501: MatChop(AinvBd, eps);
502: MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);
503: MatDestroy(&AinvBd);
504: MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);
505: MatDestroy(&AinvB);
506: } else {
507: Mat Ainvd, Ainv;
509: PCComputeExplicitOperator(pc, &Ainvd);
510: MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);
511: MatDestroy(&Ainvd);
512: #if 0
513: /* Symmetric version */
514: MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);
515: #else
516: /* Nonsymmetric version */
517: MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);
518: #endif
519: MatDestroy(&Ainv);
520: }
522: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);
523: MatView(*S, PETSC_VIEWER_STDOUT_WORLD);
524: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
526: if (D) {
527: MatInfo info;
529: MatGetInfo(D, MAT_GLOBAL_SUM, &info);
530: if (info.nz_used) SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented");
531: }
532: return(0);
533: }
537: /* Developer Notes: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
538: PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *newpmat)
539: {
541: Mat A=0,Ap=0,B=0,C=0,D=0;
552: if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
554: if (mreuse != MAT_IGNORE_MATRIX) {
555: /* Use MatSchurComplement */
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: }
562: MatGetSubMatrix(mat,isrow0,iscol0,mreuse,&A);
563: MatGetSubMatrix(mat,isrow0,iscol1,mreuse,&B);
564: MatGetSubMatrix(mat,isrow1,iscol0,mreuse,&C);
565: MatGetSubMatrix(mat,isrow1,iscol1,mreuse,&D);
566: switch (mreuse) {
567: case MAT_INITIAL_MATRIX:
568: MatCreateSchurComplement(A,A,B,C,D,newmat);
569: break;
570: case MAT_REUSE_MATRIX:
571: MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D);
572: break;
573: default:
574: SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse");
575: }
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: }
589: /*@
590: MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.
592: Collective on Mat
594: Input Parameters:
595: + A - matrix in which the complement is to be taken
596: . isrow0 - rows to eliminate
597: . iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
598: . isrow1 - rows in which the Schur complement is formed
599: . iscol1 - columns in which the Schur complement is formed
600: . mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in newmat
601: . plump - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
602: MAT_SCHUR_COMPLEMENT_AINV_DIAG | MAT_SCHUR_COMPLEMENT_AINV_LUMP
603: - preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in newpmat
605: Output Parameters:
606: + S - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
607: - Sp - approximate Schur complement suitable for preconditioning
609: Note:
610: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
611: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
612: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
613: before forming inv(diag(A00)).
615: Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
616: and column index sets. In that case, the user should call PetscObjectComposeFunction() to set
617: "MatNestGetSubMat_C" to their function. If their function needs to fall back to the default implementation, it
618: should call MatGetSchurComplement_Basic().
620: Level: advanced
622: Concepts: matrices^submatrices
624: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
625: @*/
626: PetscErrorCode MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
627: {
628: PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*);
639: if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
641: PetscObjectQueryFunction((PetscObject)S,"MatGetSchurComplement_C",&f);
642: if (f) {
643: (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
644: } else {
645: MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
646: }
647: return(0);
648: }
652: /*@
653: MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
655: Not collective.
657: Input Parameters:
658: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
659: - ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
660: MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP
662: Options database:
663: -mat_schur_complement_ainv_type diag | lump
665: Note:
666: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
667: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
668: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
669: before forming inv(diag(A00)).
671: Level: advanced
673: Concepts: matrices^submatrices
675: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
676: @*/
677: PetscErrorCode MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
678: {
679: PetscErrorCode ierr;
680: const char* t;
681: PetscBool isschur;
682: Mat_SchurComplement *schur;
686: PetscObjectGetType((PetscObject)S,&t);
687: PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
688: if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
689: schur = (Mat_SchurComplement*)S->data;
690: 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);
691: schur->ainvtype = ainvtype;
692: return(0);
693: }
697: /*@
698: MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()
700: Not collective.
702: Input Parameter:
703: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
705: Output Parameter:
706: . ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
707: MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP
709: Note:
710: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
711: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
712: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
713: before forming inv(diag(A00)).
715: Level: advanced
717: Concepts: matrices^submatrices
719: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
720: @*/
721: PetscErrorCode MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
722: {
723: PetscErrorCode ierr;
724: const char* t;
725: PetscBool isschur;
726: Mat_SchurComplement *schur;
730: PetscObjectGetType((PetscObject)S,&t);
731: PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
732: if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
733: schur = (Mat_SchurComplement*)S->data;
734: if (ainvtype) *ainvtype = schur->ainvtype;
735: return(0);
736: }
740: /*@
741: MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
743: Collective on Mat
745: Input Parameters:
746: + 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)
747: . ainvtype - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
748: - 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
750: Output Parameter:
751: - Spmat - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
753: Note:
754: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
755: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
756: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
757: before forming inv(diag(A00)).
759: Level: advanced
761: Concepts: matrices^submatrices
763: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
764: @*/
765: PetscErrorCode MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
766: {
769: PetscInt N00;
772: /* 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. */
773: /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
774: if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Spmat: A01, A10 and A11 are all NULL.");
776: /* A zero size A00 or empty A01 or A10 imply S = A11. */
777: MatGetSize(A00,&N00,NULL);
778: if (!A01 || !A10 || !N00) {
779: if (!preuse) {
780: MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
781: } else {
782: MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
783: }
785: } else {
786: Mat AdB,Sp;
787: Vec diag;
789: MatGetVecs(A00,&diag,NULL);
790: if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
791: MatGetRowSum(A00,diag);
792: } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
793: MatGetDiagonal(A00,diag);
794: } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
795: VecReciprocal(diag);
796: MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
797: MatDiagonalScale(AdB,diag,NULL);
798: VecDestroy(&diag);
799: Sp = (preuse == MAT_REUSE_MATRIX) ? *Spmat : (Mat)0;
800: MatMatMult(A10,AdB,preuse,PETSC_DEFAULT,&Sp);
801: if (!A11) {
802: MatScale(Sp,-1.0);
803: } else {
804: MatAYPX(Sp,-1,A11,DIFFERENT_NONZERO_PATTERN);
805: }
806: *Spmat = Sp;
807: MatDestroy(&AdB);
808: }
809: return(0);
810: }
814: PetscErrorCode MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
815: {
816: Mat A,B,C,D;
817: Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
818: PetscErrorCode ierr;
821: if (preuse == MAT_IGNORE_MATRIX) return(0);
823: MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
824: if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
825: MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);
826: return(0);
827: }
831: /*@
832: MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01
834: Collective on Mat
836: Input Parameters:
837: + S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
838: - 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
840: Output Parameter:
841: - Sp - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01
843: Note:
844: Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
845: application-specific information. The default for assembled matrices is to use the inverse of the diagonal of
846: the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
847: before forming inv(diag(A00)).
849: Sometimes users would like to provide problem-specific data in the Schur complement, usually only
850: for special row and column index sets. In that case, the user should call PetscObjectComposeFunction() to set
851: "MatSchurComplementGetPmat_C" to their function. If their function needs to fall back to the default implementation,
852: it should call MatSchurComplementGetPmat_Basic().
854: Level: advanced
856: Concepts: matrices^submatrices
858: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
859: @*/
860: PetscErrorCode MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
861: {
862: PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);
868: if (S->factortype) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
870: PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);
871: if (f) {
872: (*f)(S,preuse,Sp);
873: } else {
874: MatSchurComplementGetPmat_Basic(S,preuse,Sp);
875: }
876: return(0);
877: }
881: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
882: {
883: PetscErrorCode ierr;
884: Mat_SchurComplement *Na;
887: PetscNewLog(N,&Na);
888: N->data = (void*) Na;
890: N->ops->destroy = MatDestroy_SchurComplement;
891: N->ops->getvecs = MatGetVecs_SchurComplement;
892: N->ops->view = MatView_SchurComplement;
893: N->ops->mult = MatMult_SchurComplement;
894: N->ops->multadd = MatMultAdd_SchurComplement;
895: N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
896: N->assembled = PETSC_FALSE;
897: N->preallocated = PETSC_FALSE;
899: KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
900: PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
901: return(0);
902: }
904: static PetscBool KSPMatRegisterAllCalled;
908: /*@C
909: KSPMatRegisterAll - Registers all matrix implementations in the KSP package.
911: Not Collective
913: Level: advanced
915: .keywords: Mat, KSP, register, all
917: .seealso: MatRegisterAll(), MatRegisterDestroy(), KSPInitializePackage()
918: @*/
919: PetscErrorCode KSPMatRegisterAll()
920: {
924: if (KSPMatRegisterAllCalled) return(0);
925: KSPMatRegisterAllCalled = PETSC_TRUE;
926: MatRegister(MATSCHURCOMPLEMENT,MatCreate_SchurComplement);
927: return(0);
928: }