Actual source code: schurm.c

  1: #include <../src/ksp/ksp/utils/schurm/schurm.h>

  3: const char *const MatSchurComplementAinvTypes[] = {"DIAG","LUMP","BLOCKDIAG","MatSchurComplementAinvType","MAT_SCHUR_COMPLEMENT_AINV_",NULL};

  5: PetscErrorCode MatCreateVecs_SchurComplement(Mat N,Vec *right,Vec *left)
  6: {
  7:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
  8:   PetscErrorCode      ierr;

 11:   if (Na->D) {
 12:     MatCreateVecs(Na->D,right,left);
 13:     return(0);
 14:   }
 15:   if (right) {
 16:     MatCreateVecs(Na->B,right,NULL);
 17:   }
 18:   if (left) {
 19:     MatCreateVecs(Na->C,NULL,left);
 20:   }
 21:   return(0);
 22: }

 24: PetscErrorCode MatView_SchurComplement(Mat N,PetscViewer viewer)
 25: {
 26:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 27:   PetscErrorCode      ierr;

 30:   PetscViewerASCIIPrintf(viewer,"Schur complement A11 - A10 inv(A00) A01\n");
 31:   if (Na->D) {
 32:     PetscViewerASCIIPrintf(viewer,"A11\n");
 33:     PetscViewerASCIIPushTab(viewer);
 34:     MatView(Na->D,viewer);
 35:     PetscViewerASCIIPopTab(viewer);
 36:   } else {
 37:     PetscViewerASCIIPrintf(viewer,"A11 = 0\n");
 38:   }
 39:   PetscViewerASCIIPrintf(viewer,"A10\n");
 40:   PetscViewerASCIIPushTab(viewer);
 41:   MatView(Na->C,viewer);
 42:   PetscViewerASCIIPopTab(viewer);
 43:   PetscViewerASCIIPrintf(viewer,"KSP of A00\n");
 44:   PetscViewerASCIIPushTab(viewer);
 45:   KSPView(Na->ksp,viewer);
 46:   PetscViewerASCIIPopTab(viewer);
 47:   PetscViewerASCIIPrintf(viewer,"A01\n");
 48:   PetscViewerASCIIPushTab(viewer);
 49:   MatView(Na->B,viewer);
 50:   PetscViewerASCIIPopTab(viewer);
 51:   return(0);
 52: }

 54: /*
 55:            A11^T - A01^T ksptrans(A00,Ap00) A10^T
 56: */
 57: PetscErrorCode MatMultTranspose_SchurComplement(Mat N,Vec x,Vec y)
 58: {
 59:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 60:   PetscErrorCode      ierr;

 63:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
 64:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
 65:   MatMultTranspose(Na->C,x,Na->work1);
 66:   KSPSolveTranspose(Na->ksp,Na->work1,Na->work2);
 67:   MatMultTranspose(Na->B,Na->work2,y);
 68:   VecScale(y,-1.0);
 69:   if (Na->D) {
 70:     MatMultTransposeAdd(Na->D,x,y,y);
 71:   }
 72:   return(0);
 73: }

 75: /*
 76:            A11 - A10 ksp(A00,Ap00) A01
 77: */
 78: PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
 79: {
 80:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 81:   PetscErrorCode      ierr;

 84:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
 85:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
 86:   MatMult(Na->B,x,Na->work1);
 87:   KSPSolve(Na->ksp,Na->work1,Na->work2);
 88:   MatMult(Na->C,Na->work2,y);
 89:   VecScale(y,-1.0);
 90:   if (Na->D) {
 91:     MatMultAdd(Na->D,x,y,y);
 92:   }
 93:   return(0);
 94: }

 96: /*
 97:            A11 - A10 ksp(A00,Ap00) A01
 98: */
 99: PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
100: {
101:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
102:   PetscErrorCode      ierr;

105:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
106:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
107:   MatMult(Na->B,x,Na->work1);
108:   KSPSolve(Na->ksp,Na->work1,Na->work2);
109:   if (y == z) {
110:     VecScale(Na->work2,-1.0);
111:     MatMultAdd(Na->C,Na->work2,z,z);
112:   } else {
113:     MatMult(Na->C,Na->work2,z);
114:     VecAYPX(z,-1.0,y);
115:   }
116:   if (Na->D) {
117:     MatMultAdd(Na->D,x,z,z);
118:   }
119:   return(0);
120: }

122: PetscErrorCode MatSetFromOptions_SchurComplement(PetscOptionItems *PetscOptionsObject,Mat N)
123: {
124:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
125:   PetscErrorCode      ierr;

128:   PetscOptionsHead(PetscOptionsObject,"MatSchurComplementOptions");
129:   Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
130:   PetscOptionsEnum("-mat_schur_complement_ainv_type","Type of approximation for 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 A00

159:    Input Parameters:
160: +   A00,A01,A10,A11  - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
161: -   Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}

163:    Output Parameter:
164: .   S - the matrix that the Schur complement S = A11 - A10 ksp(A00,Ap00) A01

166:    Level: intermediate

168:    Notes:
169:     The Schur complement is NOT 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 simpler.

186: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement(),
187:           MatSchurComplementGetPmat()

189: @*/
190: PetscErrorCode  MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
191: {

195:   KSPInitializePackage();
196:   MatCreate(PetscObjectComm((PetscObject)A00),S);
197:   MatSetType(*S,MATSCHURCOMPLEMENT);
198:   MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);
199:   return(0);
200: }

202: /*@
203:       MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement

205:    Collective on S

207:    Input Parameters:
208: +   S                - matrix obtained with MatCreateSchurComplement (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
209: .   A00,A01,A10,A11  - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
210: -   Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.

212:    Level: intermediate

214:    Notes:
215:     The Schur complement is NOT actually formed! Rather, this
216:           object performs the matrix-vector product by using formula S = A11 - A10 A^{-1} A01
217:           for Schur complement S and a KSP solver to approximate the action of A^{-1}.

219:           All four matrices must have the same MPI communicator.

221:           A00 and  A11 must be square matrices.

223: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatCreateSchurComplement(), MatGetSchurComplement()

225: @*/
226: PetscErrorCode  MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
227: {
228:   PetscErrorCode      ierr;
229:   Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
230:   PetscBool           isschur;

233:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
234:   if (!isschur) return(0);
235:   if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
243:   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);
244:   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);
245:   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);
246:   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);
247:   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);
248:   if (A11) {
251:     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);
252:   }

254:   MatSetSizes(S,A10->rmap->n,A01->cmap->n,A10->rmap->N,A01->cmap->N);
255:   PetscObjectReference((PetscObject)A00);
256:   PetscObjectReference((PetscObject)Ap00);
257:   PetscObjectReference((PetscObject)A01);
258:   PetscObjectReference((PetscObject)A10);
259:   Na->A  = A00;
260:   Na->Ap = Ap00;
261:   Na->B  = A01;
262:   Na->C  = A10;
263:   Na->D  = A11;
264:   if (A11) {
265:     PetscObjectReference((PetscObject)A11);
266:   }
267:   MatSetUp(S);
268:   KSPSetOperators(Na->ksp,A00,Ap00);
269:   S->assembled = PETSC_TRUE;
270:   return(0);
271: }

273: /*@
274:   MatSchurComplementGetKSP - Gets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01

276:   Not Collective

278:   Input Parameter:
279: . S - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01

281:   Output Parameter:
282: . ksp - the linear solver object

284:   Options Database:
285: . -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.

287:   Level: intermediate

289: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
290: @*/
291: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
292: {
293:   Mat_SchurComplement *Na;
294:   PetscBool           isschur;
295:   PetscErrorCode      ierr;

299:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
300:   if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
302:   Na   = (Mat_SchurComplement*) S->data;
303:   *ksp = Na->ksp;
304:   return(0);
305: }

307: /*@
308:   MatSchurComplementSetKSP - Sets the KSP object that is used to invert A00 in the Schur complement matrix S = A11 - A10 ksp(A00,Ap00) A01

310:   Not Collective

312:   Input Parameters:
313: + S   - matrix created with MatCreateSchurComplement()
314: - ksp - the linear solver object

316:   Level: developer

318:   Developer Notes:
319:     This is used in PCFieldSplit to reuse the 0-split KSP to implement ksp(A00,Ap00) in S.
320:     The KSP operators are overwritten with A00 and Ap00 currently set in S.

322: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
323: @*/
324: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
325: {
326:   Mat_SchurComplement *Na;
327:   PetscErrorCode      ierr;
328:   PetscBool           isschur;

332:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
333:   if (!isschur) return(0);
335:   Na      = (Mat_SchurComplement*) S->data;
336:   PetscObjectReference((PetscObject)ksp);
337:   KSPDestroy(&Na->ksp);
338:   Na->ksp = ksp;
339:   KSPSetOperators(Na->ksp, Na->A, Na->Ap);
340:   return(0);
341: }

343: /*@
344:       MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices

346:    Collective on S

348:    Input Parameters:
349: +   S                - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
350: .   A00,A01,A10,A11  - the four parts of A = [A00 A01; A10 A11] (A11 is optional)
351: -   Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.

353:    Level: intermediate

355:    Notes:
356:     All four matrices must have the same MPI communicator

358:           A00 and  A11 must be square matrices

360:           All of the matrices provided must have the same sizes as was used with MatCreateSchurComplement() or MatSchurComplementSetSubMatrices()
361:           though they need not be the same matrices.

363: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement()

365: @*/
366: PetscErrorCode  MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
367: {
368:   PetscErrorCode      ierr;
369:   Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;
370:   PetscBool           isschur;

374:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
375:   if (!isschur) return(0);
376:   if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
384:   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);
385:   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);
386:   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);
387:   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);
388:   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);
389:   if (A11) {
392:     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);
393:   }

395:   PetscObjectReference((PetscObject)A00);
396:   PetscObjectReference((PetscObject)Ap00);
397:   PetscObjectReference((PetscObject)A01);
398:   PetscObjectReference((PetscObject)A10);
399:   if (A11) {
400:     PetscObjectReference((PetscObject)A11);
401:   }

403:   MatDestroy(&Na->A);
404:   MatDestroy(&Na->Ap);
405:   MatDestroy(&Na->B);
406:   MatDestroy(&Na->C);
407:   MatDestroy(&Na->D);

409:   Na->A  = A00;
410:   Na->Ap = Ap00;
411:   Na->B  = A01;
412:   Na->C  = A10;
413:   Na->D  = A11;

415:   KSPSetOperators(Na->ksp,A00,Ap00);
416:   return(0);
417: }

419: /*@C
420:   MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement

422:   Collective on S

424:   Input Parameter:
425: . S    - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01

427:   Output Parameters:
428: + A00  - the upper-left block of the original matrix A = [A00 A01; A10 A11]
429: . Ap00 - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}
430: . A01  - the upper-right block of the original matrix A = [A00 A01; A10 A11]
431: . A10  - the lower-left block of the original matrix A = [A00 A01; A10 A11]
432: - A11  - (optional) the lower-right block of the original matrix A = [A00 A01; A10 A11]

434:   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.

436:   Level: intermediate

438: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
439: @*/
440: PetscErrorCode  MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
441: {
442:   Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
443:   PetscErrorCode      ierr;
444:   PetscBool           flg;

448:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);
449:   if (!flg) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
450:   if (A00) *A00 = Na->A;
451:   if (Ap00) *Ap00 = Na->Ap;
452:   if (A01) *A01 = Na->B;
453:   if (A10) *A10 = Na->C;
454:   if (A11) *A11 = Na->D;
455:   return(0);
456: }

458: #include <petsc/private/kspimpl.h>

460: /*@
461:   MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

463:   Collective on M

465:   Input Parameter:
466: . M - the matrix obtained with MatCreateSchurComplement()

468:   Output Parameter:
469: . S - the Schur complement matrix

471:   Note: This can be expensive, so it is mainly for testing

473:   Level: advanced

475: .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
476: @*/
477: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat A, Mat *S)
478: {
479:   Mat            B, C, D, Bd, AinvBd;
480:   KSP            ksp;
481:   PetscInt       n,N,m,M;

485:   MatSchurComplementGetSubMatrices(A, NULL, NULL, &B, &C, &D);
486:   MatSchurComplementGetKSP(A, &ksp);
487:   KSPSetUp(ksp);
488:   MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);
489:   MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);
490:   KSPMatSolve(ksp, Bd, AinvBd);
491:   MatDestroy(&Bd);
492:   MatChop(AinvBd, PETSC_SMALL);
493:   if (D) {
494:     MatGetLocalSize(D, &m, &n);
495:     MatGetSize(D, &M, &N);
496:     MatCreateDense(PetscObjectComm((PetscObject)A), m, n, M, N, NULL, S);
497:   }
498:   MatMatMult(C, AinvBd, D ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, PETSC_DEFAULT, S);
499:   MatDestroy(&AinvBd);
500:   if (D) {
501:     MatAXPY(*S, -1.0, D, DIFFERENT_NONZERO_PATTERN);
502:   }
503:   MatConvert(*S, MATAIJ, MAT_INPLACE_MATRIX, S);
504:   MatScale(*S, -1.0);
505:   return(0);
506: }

508: /* Developer Notes:
509:     This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
510: PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *newpmat)
511: {
513:   Mat            A=NULL,Ap=NULL,B=NULL,C=NULL,D=NULL;
514:   MatReuse       reuse;

526:   if (mreuse == MAT_IGNORE_MATRIX && preuse == MAT_IGNORE_MATRIX) return(0);

530:   if (mat->factortype) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

532:   reuse = MAT_INITIAL_MATRIX;
533:   if (mreuse == MAT_REUSE_MATRIX) {
534:     MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);
535:     if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
536:     if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
537:     MatDestroy(&Ap); /* get rid of extra reference */
538:     reuse = MAT_REUSE_MATRIX;
539:   }
540:   MatCreateSubMatrix(mat,isrow0,iscol0,reuse,&A);
541:   MatCreateSubMatrix(mat,isrow0,iscol1,reuse,&B);
542:   MatCreateSubMatrix(mat,isrow1,iscol0,reuse,&C);
543:   MatCreateSubMatrix(mat,isrow1,iscol1,reuse,&D);
544:   switch (mreuse) {
545:   case MAT_INITIAL_MATRIX:
546:     MatCreateSchurComplement(A,A,B,C,D,newmat);
547:     break;
548:   case MAT_REUSE_MATRIX:
549:     MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D);
550:     break;
551:   default:
552:     if (mreuse != MAT_IGNORE_MATRIX) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse %d",(int)mreuse);
553:   }
554:   if (preuse != MAT_IGNORE_MATRIX) {
555:     MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);
556:   }
557:   MatDestroy(&A);
558:   MatDestroy(&B);
559:   MatDestroy(&C);
560:   MatDestroy(&D);
561:   return(0);
562: }

564: /*@
565:     MatGetSchurComplement - Obtain the Schur complement from eliminating part of the matrix in another part.

567:     Collective on A

569:     Input Parameters:
570: +   A      - matrix in which the complement is to be taken
571: .   isrow0 - rows to eliminate
572: .   iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
573: .   isrow1 - rows in which the Schur complement is formed
574: .   iscol1 - columns in which the Schur complement is formed
575: .   mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in S
576: .   ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
577:                        MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG, or MAT_SCHUR_COMPLEMENT_AINV_LUMP
578: -   preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in Sp

580:     Output Parameters:
581: +   S      - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
582: -   Sp     - approximate Schur complement from which a preconditioner can be built

584:     Note:
585:     Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
586:     application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
587:     the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
588:     before forming inv(diag(A00)).

590:     Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
591:     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
592:     "MatGetSchurComplement_C" to their function.  If their function needs to fall back to the default implementation, it
593:     should call MatGetSchurComplement_Basic().

595:     MatCreateSchurComplement() takes as arguments the four submatrices and returns the virtual Schur complement (what this returns in S).

597:     MatSchurComplementGetPmat() takes the virtual Schur complement and returns an explicit approximate Schur complement (what this returns in Sp).

599:     In other words calling MatCreateSchurComplement() followed by MatSchurComplementGetPmat() produces the same output as this function but with slightly different
600:     inputs. The actually submatrices of the original block matrix instead of index sets to the submatrices.

602:     Developer Notes:
603:     The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
604:     remove redundancy and be clearer and simpler.

606:     Level: advanced

608: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
609: @*/
610: PetscErrorCode  MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
611: {
612:   PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*) = NULL;

626:   if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
627:   f = NULL;
628:   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. */
629:     PetscObjectQueryFunction((PetscObject)*S,"MatGetSchurComplement_C",&f);
630:   }
631:   if (f) {
632:     (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
633:   } else {
634:     MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
635:   }
636:   return(0);
637: }

639: /*@
640:     MatSchurComplementSetAinvType - set the type of approximation used for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()

642:     Not collective.

644:     Input Parameters:
645: +   S        - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
646: -   ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
647:                       MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG

649:     Options database:
650:     -mat_schur_complement_ainv_type diag | lump | blockdiag

652:     Note:
653:     Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
654:     application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
655:     the (0,0) block A00 in place of A00^{-1}. This rarely produces a scalable algorithm. Optionally, A00 can be lumped
656:     before forming inv(diag(A00)).

658:     Level: advanced

660: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
661: @*/
662: PetscErrorCode  MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
663: {
664:   PetscErrorCode      ierr;
665:   PetscBool           isschur;
666:   Mat_SchurComplement *schur;

670:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
671:   if (!isschur) return(0);
673:   schur = (Mat_SchurComplement*)S->data;
674:   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);
675:   schur->ainvtype = ainvtype;
676:   return(0);
677: }

679: /*@
680:     MatSchurComplementGetAinvType - get the type of approximation for the inverse of the (0,0) block used in forming Sp in MatSchurComplementGetPmat()

682:     Not collective.

684:     Input Parameter:
685: .   S      - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01

687:     Output Parameter:
688: .   ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
689:                       MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_LUMP, or MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG

691:     Note:
692:     Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
693:     application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
694:     the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
695:     before forming inv(diag(A00)).

697:     Level: advanced

699: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
700: @*/
701: PetscErrorCode  MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
702: {
703:   PetscErrorCode      ierr;
704:   PetscBool           isschur;
705:   Mat_SchurComplement *schur;

709:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&isschur);
710:   if (!isschur) SETERRQ1(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONG,"Not for type %s",((PetscObject)S)->type_name);
711:   schur = (Mat_SchurComplement*)S->data;
712:   if (ainvtype) *ainvtype = schur->ainvtype;
713:   return(0);
714: }

716: /*@
717:     MatCreateSchurComplementPmat - create a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01

719:     Collective on A00

721:     Input Parameters:
722: +   A00      - the upper-left part of the original matrix A = [A00 A01; A10 A11]
723: .   A01      - (optional) the upper-right part of the original matrix A = [A00 A01; A10 A11]
724: .   A10      - (optional) the lower-left part of the original matrix A = [A00 A01; A10 A11]
725: .   A11      - (optional) the lower-right part of the original matrix A = [A00 A01; A10 A11]
726: .   ainvtype - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
727: -   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

729:     Output Parameter:
730: -   Spmat    - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01

732:     Note:
733:     Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
734:     application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
735:     the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
736:     before forming inv(diag(A00)).

738:     Level: advanced

740: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
741: @*/
742: PetscErrorCode  MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
743: {
745:   PetscInt       N00;

748:   /* 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. */
749:   /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
750:   if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Spmat: A01, A10 and A11 are all NULL.");
752:   if (preuse == MAT_IGNORE_MATRIX) return(0);

754:   /* A zero size A00 or empty A01 or A10 imply S = A11. */
755:   MatGetSize(A00,&N00,NULL);
756:   if (!A01 || !A10 || !N00) {
757:     if (preuse == MAT_INITIAL_MATRIX) {
758:       MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
759:     } else { /* MAT_REUSE_MATRIX */
760:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
761:       MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
762:     }
763:   } else {
764:     Mat AdB;
765:     Vec diag;

767:     if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
768:       MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
769:       MatCreateVecs(A00,&diag,NULL);
770:       if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
771:         MatGetRowSum(A00,diag);
772:       } else {
773:         MatGetDiagonal(A00,diag);
774:       }
775:       VecReciprocal(diag);
776:       MatDiagonalScale(AdB,diag,NULL);
777:       VecDestroy(&diag);
778:     } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
779:       Mat      A00_inv;
780:       MatType  type;
781:       MPI_Comm comm;

783:       PetscObjectGetComm((PetscObject)A00,&comm);
784:       MatGetType(A00,&type);
785:       MatCreate(comm,&A00_inv);
786:       MatSetType(A00_inv,type);
787:       MatInvertBlockDiagonalMat(A00,A00_inv);
788:       MatMatMult(A00_inv,A01,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AdB);
789:       MatDestroy(&A00_inv);
790:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
791:     /* Cannot really reuse Spmat in MatMatMult() because of MatAYPX() -->
792:          MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult()  */
793:     MatDestroy(Spmat);
794:     MatMatMult(A10,AdB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,Spmat);
795:     if (!A11) {
796:       MatScale(*Spmat,-1.0);
797:     } else {
798:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
799:       MatAYPX(*Spmat,-1,A11,DIFFERENT_NONZERO_PATTERN);
800:     }
801:     MatDestroy(&AdB);
802:   }
803:   return(0);
804: }

806: PetscErrorCode  MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
807: {
808:   Mat A,B,C,D;
809:   Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
810:   PetscErrorCode      ierr;

813:   if (preuse == MAT_IGNORE_MATRIX) return(0);
814:   MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
815:   if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
816:   MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);
817:   return(0);
818: }

820: /*@
821:     MatSchurComplementGetPmat - Obtain a preconditioning matrix for the Schur complement by assembling Sp = A11 - A10 inv(diag(A00)) A01

823:     Collective on S

825:     Input Parameters:
826: +   S      - matrix obtained with MatCreateSchurComplement() (or equivalent) and implementing the action of A11 - A10 ksp(A00,Ap00) A01
827: -   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

829:     Output Parameter:
830: -   Sp     - approximate Schur complement suitable for preconditioning S = A11 - A10 inv(diag(A00)) A01

832:     Note:
833:     Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
834:     application-specific information.  The default for assembled matrices is to use the inverse of the diagonal of
835:     the (0,0) block A00 in place of A00^{-1}. This rarely produce a scalable algorithm. Optionally, A00 can be lumped
836:     before forming inv(diag(A00)).

838:     Sometimes users would like to provide problem-specific data in the Schur complement, usually only
839:     for special row and column index sets.  In that case, the user should call PetscObjectComposeFunction() to set
840:     "MatSchurComplementGetPmat_C" to their function.  If their function needs to fall back to the default implementation,
841:     it should call MatSchurComplementGetPmat_Basic().

843:     Developer Notes:
844:     The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
845:     remove redundancy and be clearer and simpler.

847:     Level: advanced

849: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
850: @*/
851: PetscErrorCode  MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
852: {
853:   PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);

861:   if (S->factortype) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");

863:   PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);
864:   if (f) {
865:     (*f)(S,preuse,Sp);
866:   } else {
867:     MatSchurComplementGetPmat_Basic(S,preuse,Sp);
868:   }
869:   return(0);
870: }

872: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
873: {
874:   PetscErrorCode      ierr;
875:   Mat_SchurComplement *Na;

878:   PetscNewLog(N,&Na);
879:   N->data = (void*) Na;

881:   N->ops->destroy        = MatDestroy_SchurComplement;
882:   N->ops->getvecs        = MatCreateVecs_SchurComplement;
883:   N->ops->view           = MatView_SchurComplement;
884:   N->ops->mult           = MatMult_SchurComplement;
885:   N->ops->multtranspose  = MatMultTranspose_SchurComplement;
886:   N->ops->multadd        = MatMultAdd_SchurComplement;
887:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
888:   N->assembled           = PETSC_FALSE;
889:   N->preallocated        = PETSC_FALSE;

891:   KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
892:   PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
893:   return(0);
894: }