Actual source code: schurm.c

petsc-3.5.2 2014-09-08
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 (A11->rmap->n != A11->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A11 %D do not equal local columns %D",A11->rmap->n,A11->cmap->n);
244:     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);
245:   }

247:   MatGetLocalSize(A01,NULL,&n);
248:   MatGetLocalSize(A10,&m,NULL);
249:   MatSetSizes(S,m,n,PETSC_DECIDE,PETSC_DECIDE);
250:   PetscObjectReference((PetscObject)A00);
251:   PetscObjectReference((PetscObject)Ap00);
252:   PetscObjectReference((PetscObject)A01);
253:   PetscObjectReference((PetscObject)A10);
254:   Na->A  = A00;
255:   Na->Ap = Ap00;
256:   Na->B  = A01;
257:   Na->C  = A10;
258:   Na->D  = A11;
259:   if (A11) {
260:     PetscObjectReference((PetscObject)A11);
261:   }
262:   S->assembled    = PETSC_TRUE;
263:   S->preallocated = PETSC_TRUE;

265:   PetscLayoutSetUp((S)->rmap);
266:   PetscLayoutSetUp((S)->cmap);
267:   KSPSetOperators(Na->ksp,A00,Ap00);
268:   return(0);
269: }

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;

298:   Na   = (Mat_SchurComplement*) S->data;
299:   *ksp = Na->ksp;
300:   return(0);
301: }

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

308:   Not Collective

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

314:   Level: developer

316:   Developer Notes:
317:     This is used in PCFieldSplit to reuse the 0-split KSP to implement ksp(A00,Ap00) in S.

319: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
320: @*/
321: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
322: {
323:   Mat_SchurComplement *Na;
324:   PetscErrorCode      ierr;

329:   Na      = (Mat_SchurComplement*) S->data;
330:   PetscObjectReference((PetscObject)ksp);
331:   KSPDestroy(&Na->ksp);
332:   Na->ksp = ksp;
333:   KSPSetOperators(Na->ksp, Na->A, Na->Ap);
334:   return(0);
335: }

339: /*@
340:       MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices

342:    Collective on Mat

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

349:    Level: intermediate

351:    Notes: All four matrices must have the same MPI communicator

353:           A00 and  A11 must be square matrices

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

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

360: @*/
361: PetscErrorCode  MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
362: {
363:   PetscErrorCode      ierr;
364:   Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;

367:   if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
374:   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);
375:   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);
376:   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);
377:   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);
378:   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);
379:   if (A11) {
382:     if (A11->rmap->n != A11->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A11 %D do not equal local columns %D",A11->rmap->n,A11->cmap->n);
383:     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);
384:   }

386:   PetscObjectReference((PetscObject)A00);
387:   PetscObjectReference((PetscObject)Ap00);
388:   PetscObjectReference((PetscObject)A01);
389:   PetscObjectReference((PetscObject)A10);
390:   if (A11) {
391:     PetscObjectReference((PetscObject)A11);
392:   }

394:   MatDestroy(&Na->A);
395:   MatDestroy(&Na->Ap);
396:   MatDestroy(&Na->B);
397:   MatDestroy(&Na->C);
398:   MatDestroy(&Na->D);

400:   Na->A  = A00;
401:   Na->Ap = Ap00;
402:   Na->B  = A01;
403:   Na->C  = A10;
404:   Na->D  = A11;

406:   KSPSetOperators(Na->ksp,A00,Ap00);
407:   return(0);
408: }


413: /*@C
414:   MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement

416:   Collective on Mat

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

421:   Output Paramters:
422: + A00,A01,A10,A11  - the four parts of the original matrix A = [A00 A01; A10 A11] (A11 is optional)
423: - Ap00             - preconditioning matrix for use in ksp(A00,Ap00) to approximate the action of A^{-1}.

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

427:   Level: intermediate

429: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
430: @*/
431: PetscErrorCode  MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
432: {
433:   Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
434:   PetscErrorCode      ierr;
435:   PetscBool           flg;

439:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);
440:   if (flg) {
441:     if (A00) *A00 = Na->A;
442:     if (Ap00) *Ap00 = Na->Ap;
443:     if (A01) *A01 = Na->B;
444:     if (A10) *A10 = Na->C;
445:     if (A11) *A11 = Na->D;
446:   } else {
447:     if (A00) *A00 = 0;
448:     if (Ap00) *Ap00 = 0;
449:     if (A01) *A01 = 0;
450:     if (A10) *A10 = 0;
451:     if (A11) *A11 = 0;
452:   }
453:   return(0);
454: }

456: #include <petsc-private/kspimpl.h>

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

463:   Collective on Mat

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 M, Mat *S)
478: {
479:   Mat            B, C, D;
480:   KSP            ksp;
481:   PC             pc;
482:   PetscBool      isLU, isILU;
483:   PetscReal      fill = 2.0;

487:   MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);
488:   MatSchurComplementGetKSP(M, &ksp);
489:   KSPGetPC(ksp, &pc);
490:   PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);
491:   PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);
492:   if (isLU || isILU) {
493:     Mat       fact, Bd, AinvB, AinvBd;
494:     PetscReal eps = 1.0e-10;

496:     /* This can be sped up for banded LU */
497:     KSPSetUp(ksp);
498:     PCFactorGetMatrix(pc, &fact);
499:     MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);
500:     MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);
501:     MatMatSolve(fact, Bd, AinvBd);
502:     MatDestroy(&Bd);
503:     MatChop(AinvBd, eps);
504:     MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);
505:     MatDestroy(&AinvBd);
506:     MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);
507:     MatDestroy(&AinvB);
508:   } else {
509:     Mat Ainvd, Ainv;

511:     PCComputeExplicitOperator(pc, &Ainvd);
512:     MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);
513:     MatDestroy(&Ainvd);
514: #if 0
515:     /* Symmetric version */
516:     MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);
517: #else
518:     /* Nonsymmetric version */
519:     MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);
520: #endif
521:     MatDestroy(&Ainv);
522:   }

524:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);
525:   MatView(*S, PETSC_VIEWER_STDOUT_WORLD);
526:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);

528:   if (D) {
529:     MatInfo info;

531:     MatGetInfo(D, MAT_GLOBAL_SUM, &info);
532:     if (info.nz_used) SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented");
533:   }
534:   return(0);
535: }

539: /* Developer Notes: This should be implemented with a MatCreate_SchurComplement() as that is the standard design for new Mat classes. */
540: PetscErrorCode MatGetSchurComplement_Basic(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatSchurComplementAinvType ainvtype, MatReuse preuse,Mat *newpmat)
541: {
543:   Mat            A=0,Ap=0,B=0,C=0,D=0;

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

556:   if (mreuse != MAT_IGNORE_MATRIX) {
557:     /* Use MatSchurComplement */
558:     if (mreuse == MAT_REUSE_MATRIX) {
559:       MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);
560:       if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
561:       if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
562:       MatDestroy(&Ap); /* get rid of extra reference */
563:     }
564:     MatGetSubMatrix(mat,isrow0,iscol0,mreuse,&A);
565:     MatGetSubMatrix(mat,isrow0,iscol1,mreuse,&B);
566:     MatGetSubMatrix(mat,isrow1,iscol0,mreuse,&C);
567:     MatGetSubMatrix(mat,isrow1,iscol1,mreuse,&D);
568:     switch (mreuse) {
569:     case MAT_INITIAL_MATRIX:
570:       MatCreateSchurComplement(A,A,B,C,D,newmat);
571:       break;
572:     case MAT_REUSE_MATRIX:
573:       MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D);
574:       break;
575:     default:
576:       SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse");
577:     }
578:   }
579:   if (preuse != MAT_IGNORE_MATRIX) {
580:     MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);
581:   }
582:   MatDestroy(&A);
583:   MatDestroy(&B);
584:   MatDestroy(&C);
585:   MatDestroy(&D);
586:   return(0);
587: }

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

594:     Collective on Mat

596:     Input Parameters:
597: +   A      - matrix in which the complement is to be taken
598: .   isrow0 - rows to eliminate
599: .   iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
600: .   isrow1 - rows in which the Schur complement is formed
601: .   iscol1 - columns in which the Schur complement is formed
602: .   mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in newmat
603: .   plump  - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
604: $                        MAT_SCHUR_COMPLEMENT_AINV_DIAG | MAT_SCHUR_COMPLEMENT_AINV_LUMP
605: -   preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in newpmat

607:     Output Parameters:
608: +   S      - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
609: -   Sp     - approximate Schur complement suitable for preconditioning

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

617:     Sometimes users would like to provide problem-specific data in the Schur complement, usually only for special row
618:     and column index sets.  In that case, the user should call PetscObjectComposeFunction() to set
619:     "MatNestGetSubMat_C" to their function.  If their function needs to fall back to the default implementation, it
620:     should call MatGetSchurComplement_Basic().

622:     Level: advanced

624:     Concepts: matrices^submatrices

626: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
627: @*/
628: PetscErrorCode  MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
629: {
630:   PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*);

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

643:   PetscObjectQueryFunction((PetscObject)S,"MatGetSchurComplement_C",&f);
644:   if (f) {
645:     (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
646:   } else {
647:     MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
648:   }
649:   return(0);
650: }

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

657:     Not collective.

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

664:     Options database:
665:     -mat_schur_complement_ainv_type diag | lump

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

673:     Level: advanced

675:     Concepts: matrices^submatrices

677: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
678: @*/
679: PetscErrorCode  MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
680: {
681:   PetscErrorCode      ierr;
682:   const char*         t;
683:   PetscBool           isschur;
684:   Mat_SchurComplement *schur;

688:   PetscObjectGetType((PetscObject)S,&t);
689:   PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
690:   if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
691:   schur = (Mat_SchurComplement*)S->data;
692:   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);
693:   schur->ainvtype = ainvtype;
694:   return(0);
695: }

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

702:     Not collective.

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

707:     Output Parameter:
708: .   ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
709: $                      MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP

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

717:     Level: advanced

719:     Concepts: matrices^submatrices

721: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
722: @*/
723: PetscErrorCode  MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
724: {
725:   PetscErrorCode      ierr;
726:   const char*         t;
727:   PetscBool           isschur;
728:   Mat_SchurComplement *schur;

732:   PetscObjectGetType((PetscObject)S,&t);
733:   PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
734:   if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
735:   schur = (Mat_SchurComplement*)S->data;
736:   if (ainvtype) *ainvtype = schur->ainvtype;
737:   return(0);
738: }

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

745:     Collective on Mat

747:     Input Parameters:
748: +   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)
749: .   ainvtype             - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
750: -   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

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

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

761:     Level: advanced

763:     Concepts: matrices^submatrices

765: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
766: @*/
767: PetscErrorCode  MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
768: {

771:   PetscInt       N00;

774:   /* Use an appropriate approximate inverse of A00 to form A11 - A10 inv(diag(A00)) A01; a NULL A01, A10 or A11 indicates a zero matrix. */
775:   /* TODO: Perhaps should create an appropriately-sized zero matrix of the same type as A00? */
776:   if ((!A01 || !A10) & !A11) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot assemble Spmat: A01, A10 and A11 are all NULL.");

778:   /* A zero size A00 or empty A01 or A10 imply S = A11. */
779:   MatGetSize(A00,&N00,NULL);
780:   if (!A01 || !A10 || !N00) {
781:     if (!preuse) {
782:       MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
783:     } else {
784:       MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
785:     }

787:   } else {
788:     Mat         AdB,Sp;
789:     Vec         diag;

791:     MatGetVecs(A00,&diag,NULL);
792:     if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
793:       MatGetRowSum(A00,diag);
794:     } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
795:       MatGetDiagonal(A00,diag);
796:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
797:     VecReciprocal(diag);
798:     MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
799:     MatDiagonalScale(AdB,diag,NULL);
800:     VecDestroy(&diag);
801:     Sp       = (preuse == MAT_REUSE_MATRIX) ? *Spmat : (Mat)0;
802:     MatMatMult(A10,AdB,preuse,PETSC_DEFAULT,&Sp);
803:     if (!A11) {
804:       MatScale(Sp,-1.0);
805:     } else {
806:       MatAYPX(Sp,-1,A11,DIFFERENT_NONZERO_PATTERN);
807:     }
808:     *Spmat    = Sp;
809:     MatDestroy(&AdB);
810:   }
811:   return(0);
812: }

816: PetscErrorCode  MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
817: {
818:   Mat A,B,C,D;
819:   Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
820:   PetscErrorCode      ierr;

823:   if (preuse == MAT_IGNORE_MATRIX) return(0);

825:   MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
826:   if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
827:   MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);
828:   return(0);
829: }

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

836:     Collective on Mat

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

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

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

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

856:     Level: advanced

858:     Concepts: matrices^submatrices

860: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
861: @*/
862: PetscErrorCode  MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
863: {
864:   PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);

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

872:   PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);
873:   if (f) {
874:     (*f)(S,preuse,Sp);
875:   } else {
876:     MatSchurComplementGetPmat_Basic(S,preuse,Sp);
877:   }
878:   return(0);
879: }

883: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
884: {
885:   PetscErrorCode      ierr;
886:   Mat_SchurComplement *Na;

889:   PetscNewLog(N,&Na);
890:   N->data = (void*) Na;

892:   N->ops->destroy        = MatDestroy_SchurComplement;
893:   N->ops->getvecs        = MatGetVecs_SchurComplement;
894:   N->ops->view           = MatView_SchurComplement;
895:   N->ops->mult           = MatMult_SchurComplement;
896:   N->ops->multadd        = MatMultAdd_SchurComplement;
897:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
898:   N->assembled           = PETSC_FALSE;
899:   N->preallocated        = PETSC_FALSE;

901:   KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
902:   PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
903:   return(0);
904: }

906: static PetscBool KSPMatRegisterAllCalled;

910: /*@C
911:   KSPMatRegisterAll - Registers all matrix implementations in the KSP package.

913:   Not Collective

915:   Level: advanced

917: .keywords: Mat, KSP, register, all

919: .seealso: MatRegisterAll(), MatRegisterDestroy(), KSPInitializePackage()
920: @*/
921: PetscErrorCode KSPMatRegisterAll()
922: {

926:   if (KSPMatRegisterAllCalled) return(0);
927:   KSPMatRegisterAllCalled = PETSC_TRUE;
928:   MatRegister(MATSCHURCOMPLEMENT,MatCreate_SchurComplement);
929:   return(0);
930: }