Actual source code: schurm.c

petsc-master 2016-05-24
Report Typos and Errors
  1: #include <petsc/private/matimpl.h>
  2: #include <petscksp.h>                 /*I "petscksp.h" I*/
  3: const char *const MatSchurComplementAinvTypes[] = {"DIAG","LUMP","MatSchurComplementAinvType","MAT_SCHUR_COMPLEMENT_AINV_",0};

  5: typedef struct {
  6:   Mat                        A,Ap,B,C,D;
  7:   KSP                        ksp;
  8:   Vec                        work1,work2;
  9:   MatSchurComplementAinvType ainvtype;
 10: } Mat_SchurComplement;



 16: PetscErrorCode MatCreateVecs_SchurComplement(Mat N,Vec *right,Vec *left)
 17: {
 18:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 19:   PetscErrorCode      ierr;

 22:   if (Na->D) {
 23:     MatCreateVecs(Na->D,right,left);
 24:     return(0);
 25:   }
 26:   if (right) {
 27:     MatCreateVecs(Na->B,right,NULL);
 28:   }
 29:   if (left) {
 30:     MatCreateVecs(Na->C,NULL,left);
 31:   }
 32:   return(0);
 33: }

 37: PetscErrorCode MatView_SchurComplement(Mat N,PetscViewer viewer)
 38: {
 39:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 40:   PetscErrorCode      ierr;

 43:   PetscViewerASCIIPrintf(viewer,"Schur complement A11 - A10 inv(A00) A01\n");
 44:   if (Na->D) {
 45:     PetscViewerASCIIPrintf(viewer,"A11\n");
 46:     PetscViewerASCIIPushTab(viewer);
 47:     MatView(Na->D,viewer);
 48:     PetscViewerASCIIPopTab(viewer);
 49:   } else {
 50:     PetscViewerASCIIPrintf(viewer,"A11 = 0\n");
 51:   }
 52:   PetscViewerASCIIPrintf(viewer,"A10\n");
 53:   PetscViewerASCIIPushTab(viewer);
 54:   MatView(Na->C,viewer);
 55:   PetscViewerASCIIPopTab(viewer);
 56:   PetscViewerASCIIPrintf(viewer,"KSP of A00\n");
 57:   PetscViewerASCIIPushTab(viewer);
 58:   KSPView(Na->ksp,viewer);
 59:   PetscViewerASCIIPopTab(viewer);
 60:   PetscViewerASCIIPrintf(viewer,"A01\n");
 61:   PetscViewerASCIIPushTab(viewer);
 62:   MatView(Na->B,viewer);
 63:   PetscViewerASCIIPopTab(viewer);
 64:   return(0);
 65: }

 67: /*
 68:            A11^T - A01^T ksptrans(A00,Ap00) A10^T
 69: */
 72: PetscErrorCode MatMultTranspose_SchurComplement(Mat N,Vec x,Vec y)
 73: {
 74:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 75:   PetscErrorCode      ierr;

 78:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
 79:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
 80:   MatMultTranspose(Na->C,x,Na->work1);
 81:   KSPSolveTranspose(Na->ksp,Na->work1,Na->work2);
 82:   MatMultTranspose(Na->B,Na->work2,y);
 83:   VecScale(y,-1.0);
 84:   if (Na->D) {
 85:     MatMultTransposeAdd(Na->D,x,y,y);
 86:   }
 87:   return(0);
 88: }

 90: /*
 91:            A11 - A10 ksp(A00,Ap00) A01
 92: */
 95: PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
 96: {
 97:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 98:   PetscErrorCode      ierr;

101:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
102:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
103:   MatMult(Na->B,x,Na->work1);
104:   KSPSolve(Na->ksp,Na->work1,Na->work2);
105:   MatMult(Na->C,Na->work2,y);
106:   VecScale(y,-1.0);
107:   if (Na->D) {
108:     MatMultAdd(Na->D,x,y,y);
109:   }
110:   return(0);
111: }

113: /*
114:            A11 - A10 ksp(A00,Ap00) A01
115: */
118: PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
119: {
120:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
121:   PetscErrorCode      ierr;

124:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
125:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
126:   MatMult(Na->B,x,Na->work1);
127:   KSPSolve(Na->ksp,Na->work1,Na->work2);
128:   if (y == z) {
129:     VecScale(Na->work2,-1.0);
130:     MatMultAdd(Na->C,Na->work2,z,z);
131:   } else {
132:     MatMult(Na->C,Na->work2,z);
133:     VecAYPX(z,-1.0,y);
134:   }
135:   if (Na->D) {
136:     MatMultAdd(Na->D,x,z,z);
137:   }
138:   return(0);
139: }

143: PetscErrorCode MatSetFromOptions_SchurComplement(PetscOptionItems *PetscOptionsObject,Mat N)
144: {
145:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
146:   PetscErrorCode      ierr;

149:   PetscOptionsHead(PetscOptionsObject,"MatSchurComplementOptions");
150:   Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
151:   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);
152:   PetscOptionsTail();
153:   KSPSetFromOptions(Na->ksp);
154:   return(0);
155: }

159: PetscErrorCode MatDestroy_SchurComplement(Mat N)
160: {
161:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
162:   PetscErrorCode      ierr;

165:   MatDestroy(&Na->A);
166:   MatDestroy(&Na->Ap);
167:   MatDestroy(&Na->B);
168:   MatDestroy(&Na->C);
169:   MatDestroy(&Na->D);
170:   VecDestroy(&Na->work1);
171:   VecDestroy(&Na->work2);
172:   KSPDestroy(&Na->ksp);
173:   PetscFree(N->data);
174:   return(0);
175: }

179: /*@
180:       MatCreateSchurComplement - Creates a new matrix object that behaves like the Schur complement of a matrix

182:    Collective on Mat

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

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

191:    Level: intermediate

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

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

199:           A00 and  A11 must be square matrices.

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

203: @*/
204: PetscErrorCode  MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
205: {

209:   KSPInitializePackage();
210:   MatCreate(((PetscObject)A00)->comm,S);
211:   MatSetType(*S,MATSCHURCOMPLEMENT);
212:   MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);
213:   return(0);
214: }

218: /*@
219:       MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement

221:    Collective on Mat

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

228:    Level: intermediate

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

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

236:           A00 and  A11 must be square matrices.

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

240: @*/
241: PetscErrorCode  MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
242: {
243:   PetscErrorCode      ierr;
244:   PetscInt            m,n;
245:   Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;

248:   if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
256:   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);
257:   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);
258:   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);
259:   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);
260:   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);
261:   if (A11) {
264:     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);
265:   }

267:   MatGetLocalSize(A01,NULL,&n);
268:   MatGetLocalSize(A10,&m,NULL);
269:   MatSetSizes(S,m,n,PETSC_DECIDE,PETSC_DECIDE);
270:   PetscObjectReference((PetscObject)A00);
271:   PetscObjectReference((PetscObject)Ap00);
272:   PetscObjectReference((PetscObject)A01);
273:   PetscObjectReference((PetscObject)A10);
274:   Na->A  = A00;
275:   Na->Ap = Ap00;
276:   Na->B  = A01;
277:   Na->C  = A10;
278:   Na->D  = A11;
279:   if (A11) {
280:     PetscObjectReference((PetscObject)A11);
281:   }
282:   S->assembled    = PETSC_TRUE;
283:   S->preallocated = PETSC_TRUE;

285:   PetscLayoutSetUp((S)->rmap);
286:   PetscLayoutSetUp((S)->cmap);
287:   KSPSetOperators(Na->ksp,A00,Ap00);
288:   return(0);
289: }

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

296:   Not Collective

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

301:   Output Parameter:
302: . ksp - the linear solver object

304:   Options Database:
305: . -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.

307:   Level: intermediate

309: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
310: @*/
311: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
312: {
313:   Mat_SchurComplement *Na;

318:   Na   = (Mat_SchurComplement*) S->data;
319:   *ksp = Na->ksp;
320:   return(0);
321: }

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

328:   Not Collective

330:   Input Parameters:
331: + S   - matrix created with MatCreateSchurComplement()
332: - ksp - the linear solver object

334:   Level: developer

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

339: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
340: @*/
341: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
342: {
343:   Mat_SchurComplement *Na;
344:   PetscErrorCode      ierr;

349:   Na      = (Mat_SchurComplement*) S->data;
350:   PetscObjectReference((PetscObject)ksp);
351:   KSPDestroy(&Na->ksp);
352:   Na->ksp = ksp;
353:   KSPSetOperators(Na->ksp, Na->A, Na->Ap);
354:   return(0);
355: }

359: /*@
360:       MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices

362:    Collective on Mat

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

369:    Level: intermediate

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

373:           A00 and  A11 must be square matrices

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

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

380: @*/
381: PetscErrorCode  MatSchurComplementUpdateSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
382: {
383:   PetscErrorCode      ierr;
384:   Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;

387:   if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
394:   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);
395:   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);
396:   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);
397:   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);
398:   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);
399:   if (A11) {
402:     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);
403:   }

405:   PetscObjectReference((PetscObject)A00);
406:   PetscObjectReference((PetscObject)Ap00);
407:   PetscObjectReference((PetscObject)A01);
408:   PetscObjectReference((PetscObject)A10);
409:   if (A11) {
410:     PetscObjectReference((PetscObject)A11);
411:   }

413:   MatDestroy(&Na->A);
414:   MatDestroy(&Na->Ap);
415:   MatDestroy(&Na->B);
416:   MatDestroy(&Na->C);
417:   MatDestroy(&Na->D);

419:   Na->A  = A00;
420:   Na->Ap = Ap00;
421:   Na->B  = A01;
422:   Na->C  = A10;
423:   Na->D  = A11;

425:   KSPSetOperators(Na->ksp,A00,Ap00);
426:   return(0);
427: }


432: /*@C
433:   MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement

435:   Collective on Mat

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

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

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

446:   Level: intermediate

448: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdateSubMatrices()
449: @*/
450: PetscErrorCode  MatSchurComplementGetSubMatrices(Mat S,Mat *A00,Mat *Ap00,Mat *A01,Mat *A10,Mat *A11)
451: {
452:   Mat_SchurComplement *Na = (Mat_SchurComplement*) S->data;
453:   PetscErrorCode      ierr;
454:   PetscBool           flg;

458:   PetscObjectTypeCompare((PetscObject)S,MATSCHURCOMPLEMENT,&flg);
459:   if (flg) {
460:     if (A00) *A00 = Na->A;
461:     if (Ap00) *Ap00 = Na->Ap;
462:     if (A01) *A01 = Na->B;
463:     if (A10) *A10 = Na->C;
464:     if (A11) *A11 = Na->D;
465:   } else {
466:     if (A00) *A00 = 0;
467:     if (Ap00) *Ap00 = 0;
468:     if (A01) *A01 = 0;
469:     if (A10) *A10 = 0;
470:     if (A11) *A11 = 0;
471:   }
472:   return(0);
473: }

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

479: /*@
480:   MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

482:   Collective on Mat

484:   Input Parameter:
485: . M - the matrix obtained with MatCreateSchurComplement()

487:   Output Parameter:
488: . S - the Schur complement matrix

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

492:   Level: advanced

494: .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
495: @*/
496: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
497: {
498:   Mat            B, C, D;
499:   KSP            ksp;
500:   PC             pc;
501:   PetscBool      isLU, isILU;
502:   PetscReal      fill = 2.0;

506:   MatSchurComplementGetSubMatrices(M, NULL, NULL, &B, &C, &D);
507:   MatSchurComplementGetKSP(M, &ksp);
508:   KSPGetPC(ksp, &pc);
509:   PetscObjectTypeCompare((PetscObject) pc, PCLU, &isLU);
510:   PetscObjectTypeCompare((PetscObject) pc, PCILU, &isILU);
511:   if (isLU || isILU) {
512:     Mat       fact, Bd, AinvB, AinvBd;
513:     PetscReal eps = 1.0e-10;

515:     /* This can be sped up for banded LU */
516:     KSPSetUp(ksp);
517:     PCFactorGetMatrix(pc, &fact);
518:     MatConvert(B, MATDENSE, MAT_INITIAL_MATRIX, &Bd);
519:     MatDuplicate(Bd, MAT_DO_NOT_COPY_VALUES, &AinvBd);
520:     MatMatSolve(fact, Bd, AinvBd);
521:     MatDestroy(&Bd);
522:     MatChop(AinvBd, eps);
523:     MatConvert(AinvBd, MATAIJ, MAT_INITIAL_MATRIX, &AinvB);
524:     MatDestroy(&AinvBd);
525:     MatMatMult(C, AinvB, MAT_INITIAL_MATRIX, fill, S);
526:     MatDestroy(&AinvB);
527:   } else {
528:     Mat Ainvd, Ainv;

530:     PCComputeExplicitOperator(pc, &Ainvd);
531:     MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);
532:     MatDestroy(&Ainvd);
533: #if 0
534:     /* Symmetric version */
535:     MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);
536: #else
537:     /* Nonsymmetric version */
538:     MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);
539: #endif
540:     MatDestroy(&Ainv);
541:   }

543:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);
544:   MatView(*S, PETSC_VIEWER_STDOUT_WORLD);
545:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);

547:   if (D) {
548:     MatInfo   info;
549:     PetscReal norm;

551:     MatGetInfo(D, MAT_GLOBAL_SUM, &info);
552:     if (info.nz_used) {
553:       MatNorm(D, NORM_INFINITY, &norm);
554:       if (norm > PETSC_MACHINE_EPSILON) SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented for Schur complements with non-vanishing D");
555:     }
556:   }
557:   return(0);
558: }

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

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

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

582:   reuse = MAT_INITIAL_MATRIX;
583:   if (mreuse == MAT_REUSE_MATRIX) {
584:     MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);
585:     if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
586:     if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
587:     MatDestroy(&Ap); /* get rid of extra reference */
588:     reuse = MAT_REUSE_MATRIX;
589:   }
590:   MatGetSubMatrix(mat,isrow0,iscol0,reuse,&A);
591:   MatGetSubMatrix(mat,isrow0,iscol1,reuse,&B);
592:   MatGetSubMatrix(mat,isrow1,iscol0,reuse,&C);
593:   MatGetSubMatrix(mat,isrow1,iscol1,reuse,&D);
594:   switch (mreuse) {
595:   case MAT_INITIAL_MATRIX:
596:     MatCreateSchurComplement(A,A,B,C,D,newmat);
597:     break;
598:   case MAT_REUSE_MATRIX:
599:     MatSchurComplementUpdateSubMatrices(*newmat,A,A,B,C,D);
600:     break;
601:   default:
602:     if (mreuse != MAT_IGNORE_MATRIX) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse");
603:   }
604:   if (preuse != MAT_IGNORE_MATRIX) {
605:     MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);
606:   }
607:   MatDestroy(&A);
608:   MatDestroy(&B);
609:   MatDestroy(&C);
610:   MatDestroy(&D);
611:   return(0);
612: }

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

619:     Collective on Mat

621:     Input Parameters:
622: +   A      - matrix in which the complement is to be taken
623: .   isrow0 - rows to eliminate
624: .   iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
625: .   isrow1 - rows in which the Schur complement is formed
626: .   iscol1 - columns in which the Schur complement is formed
627: .   mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in S
628: .   plump  - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
629:                        MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP
630: -   preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in Sp

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

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

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

647:     Level: advanced

649:     Concepts: matrices^submatrices

651: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
652: @*/
653: PetscErrorCode  MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
654: {
655:   PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*) = NULL;

666:   if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
667:   f = NULL;
668:   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. */
669:     PetscObjectQueryFunction((PetscObject)*S,"MatGetSchurComplement_C",&f);
670:   }
671:   if (f) {
672:       (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
673:   } else {
674:     MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
675:   }
676:   return(0);
677: }

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

684:     Not collective.

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

691:     Options database:
692:     -mat_schur_complement_ainv_type diag | lump

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

700:     Level: advanced

702:     Concepts: matrices^submatrices

704: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
705: @*/
706: PetscErrorCode  MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
707: {
708:   PetscErrorCode      ierr;
709:   const char*         t;
710:   PetscBool           isschur;
711:   Mat_SchurComplement *schur;

715:   PetscObjectGetType((PetscObject)S,&t);
716:   PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
717:   if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
718:   schur = (Mat_SchurComplement*)S->data;
719:   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);
720:   schur->ainvtype = ainvtype;
721:   return(0);
722: }

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

729:     Not collective.

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

734:     Output Parameter:
735: .   ainvtype - type of approximation used to form A00inv from A00 when assembling Sp = A11 - A10 A00inv A01:
736:                       MAT_SCHUR_COMPLEMENT_AINV_DIAG or MAT_SCHUR_COMPLEMENT_AINV_LUMP

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

744:     Level: advanced

746:     Concepts: matrices^submatrices

748: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
749: @*/
750: PetscErrorCode  MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
751: {
752:   PetscErrorCode      ierr;
753:   const char*         t;
754:   PetscBool           isschur;
755:   Mat_SchurComplement *schur;

759:   PetscObjectGetType((PetscObject)S,&t);
760:   PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
761:   if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
762:   schur = (Mat_SchurComplement*)S->data;
763:   if (ainvtype) *ainvtype = schur->ainvtype;
764:   return(0);
765: }

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

772:     Collective on Mat

774:     Input Parameters:
775: +   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)
776: .   ainvtype             - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
777: -   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

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

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

788:     Level: advanced

790:     Concepts: matrices^submatrices

792: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
793: @*/
794: PetscErrorCode  MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
795: {

798:   PetscInt       N00;

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

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

807:   /* A zero size A00 or empty A01 or A10 imply S = A11. */
808:   MatGetSize(A00,&N00,NULL);
809:   if (!A01 || !A10 || !N00) {
810:     if (preuse == MAT_INITIAL_MATRIX) {
811:       MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
812:     } else { /* MAT_REUSE_MATRIX */
813:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
814:       MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
815:     }

817:   } else {
818:     Mat         AdB;
819:     Vec         diag;

821:     MatCreateVecs(A00,&diag,NULL);
822:     if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
823:       MatGetRowSum(A00,diag);
824:     } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
825:       MatGetDiagonal(A00,diag);
826:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
827:     VecReciprocal(diag);
828:     MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
829:     MatDiagonalScale(AdB,diag,NULL);
830:     VecDestroy(&diag);
831:     /* Cannot really reuse Spmat in MatMatMult() because of MatAYPX() -->
832:          MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult()  */
833:     MatDestroy(Spmat);
834:     MatMatMult(A10,AdB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,Spmat);
835:     if (!A11) {
836:       MatScale(*Spmat,-1.0);
837:     } else {
838:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
839:       MatAYPX(*Spmat,-1,A11,DIFFERENT_NONZERO_PATTERN);
840:     }
841:     MatDestroy(&AdB);
842:   }
843:   return(0);
844: }

848: PetscErrorCode  MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
849: {
850:   Mat A,B,C,D;
851:   Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
852:   PetscErrorCode      ierr;

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

857:   MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
858:   if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
859:   MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);
860:   return(0);
861: }

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

868:     Collective on Mat

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

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

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

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

888:     Level: advanced

890:     Concepts: matrices^submatrices

892: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
893: @*/
894: PetscErrorCode  MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
895: {
896:   PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);

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

904:   PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);
905:   if (f) {
906:     (*f)(S,preuse,Sp);
907:   } else {
908:     MatSchurComplementGetPmat_Basic(S,preuse,Sp);
909:   }
910:   return(0);
911: }

915: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
916: {
917:   PetscErrorCode      ierr;
918:   Mat_SchurComplement *Na;

921:   PetscNewLog(N,&Na);
922:   N->data = (void*) Na;

924:   N->ops->destroy        = MatDestroy_SchurComplement;
925:   N->ops->getvecs        = MatCreateVecs_SchurComplement;
926:   N->ops->view           = MatView_SchurComplement;
927:   N->ops->mult           = MatMult_SchurComplement;
928:   N->ops->multtranspose  = MatMultTranspose_SchurComplement;
929:   N->ops->multadd        = MatMultAdd_SchurComplement;
930:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
931:   N->assembled           = PETSC_FALSE;
932:   N->preallocated        = PETSC_FALSE;

934:   KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
935:   PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
936:   return(0);
937: }

939: static PetscBool KSPMatRegisterAllCalled;

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

946:   Not Collective

948:   Level: advanced

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

952: .seealso: MatRegisterAll(),  KSPInitializePackage()
953: @*/
954: PetscErrorCode KSPMatRegisterAll()
955: {

959:   if (KSPMatRegisterAllCalled) return(0);
960:   KSPMatRegisterAllCalled = PETSC_TRUE;
961:   MatRegister(MATSCHURCOMPLEMENT,MatCreate_SchurComplement);
962:   return(0);
963: }