Actual source code: schurm.c

petsc-3.5.3 2015-01-31
Report Typos and Errors
  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: }