Actual source code: schurm.c

petsc-master 2018-05-20
Report Typos and Errors
  1:  #include <petsc/private/matimpl.h>
  2:  #include <petscksp.h>
  3: const char *const MatSchurComplementAinvTypes[] = {"DIAG","LUMP","BLOCKDIAG","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;

 12: PetscErrorCode MatCreateVecs_SchurComplement(Mat N,Vec *right,Vec *left)
 13: {
 14:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 15:   PetscErrorCode      ierr;

 18:   if (Na->D) {
 19:     MatCreateVecs(Na->D,right,left);
 20:     return(0);
 21:   }
 22:   if (right) {
 23:     MatCreateVecs(Na->B,right,NULL);
 24:   }
 25:   if (left) {
 26:     MatCreateVecs(Na->C,NULL,left);
 27:   }
 28:   return(0);
 29: }

 31: PetscErrorCode MatView_SchurComplement(Mat N,PetscViewer viewer)
 32: {
 33:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 34:   PetscErrorCode      ierr;

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

 61: /*
 62:            A11^T - A01^T ksptrans(A00,Ap00) A10^T
 63: */
 64: PetscErrorCode MatMultTranspose_SchurComplement(Mat N,Vec x,Vec y)
 65: {
 66:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 67:   PetscErrorCode      ierr;

 70:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
 71:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
 72:   MatMultTranspose(Na->C,x,Na->work1);
 73:   KSPSolveTranspose(Na->ksp,Na->work1,Na->work2);
 74:   MatMultTranspose(Na->B,Na->work2,y);
 75:   VecScale(y,-1.0);
 76:   if (Na->D) {
 77:     MatMultTransposeAdd(Na->D,x,y,y);
 78:   }
 79:   return(0);
 80: }

 82: /*
 83:            A11 - A10 ksp(A00,Ap00) A01
 84: */
 85: PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
 86: {
 87:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 88:   PetscErrorCode      ierr;

 91:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
 92:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
 93:   MatMult(Na->B,x,Na->work1);
 94:   KSPSolve(Na->ksp,Na->work1,Na->work2);
 95:   MatMult(Na->C,Na->work2,y);
 96:   VecScale(y,-1.0);
 97:   if (Na->D) {
 98:     MatMultAdd(Na->D,x,y,y);
 99:   }
100:   return(0);
101: }

103: /*
104:            A11 - A10 ksp(A00,Ap00) A01
105: */
106: PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
107: {
108:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
109:   PetscErrorCode      ierr;

112:   if (!Na->work1) {MatCreateVecs(Na->A,&Na->work1,NULL);}
113:   if (!Na->work2) {MatCreateVecs(Na->A,&Na->work2,NULL);}
114:   MatMult(Na->B,x,Na->work1);
115:   KSPSolve(Na->ksp,Na->work1,Na->work2);
116:   if (y == z) {
117:     VecScale(Na->work2,-1.0);
118:     MatMultAdd(Na->C,Na->work2,z,z);
119:   } else {
120:     MatMult(Na->C,Na->work2,z);
121:     VecAYPX(z,-1.0,y);
122:   }
123:   if (Na->D) {
124:     MatMultAdd(Na->D,x,z,z);
125:   }
126:   return(0);
127: }

129: PetscErrorCode MatSetFromOptions_SchurComplement(PetscOptionItems *PetscOptionsObject,Mat N)
130: {
131:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
132:   PetscErrorCode      ierr;

135:   PetscOptionsHead(PetscOptionsObject,"MatSchurComplementOptions");
136:   Na->ainvtype = MAT_SCHUR_COMPLEMENT_AINV_DIAG;
137:   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);
138:   PetscOptionsTail();
139:   KSPSetFromOptions(Na->ksp);
140:   return(0);
141: }

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

149:   MatDestroy(&Na->A);
150:   MatDestroy(&Na->Ap);
151:   MatDestroy(&Na->B);
152:   MatDestroy(&Na->C);
153:   MatDestroy(&Na->D);
154:   VecDestroy(&Na->work1);
155:   VecDestroy(&Na->work2);
156:   KSPDestroy(&Na->ksp);
157:   PetscFree(N->data);
158:   return(0);
159: }

161: /*@C
162:       MatCreateSchurComplement - Creates a new matrix object that behaves like the Schur complement of a matrix

164:    Collective on Mat

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

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

173:    Level: intermediate

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

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

182:           A00 and  A11 must be square matrices.

184:           MatGetSchurComplement() takes as arguments the index sets for the submatrices and returns both the virtual Schur complement (what this returns) plus
185:           a sparse approximation to the true Schur complement (useful for building a preconditioner for the Schur complement).

187:           MatSchurComplementGetPmat() can be called on the output of this function to generate an explicit approximation to the Schur complement.

189:     Developer Notes:
190:     The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
191:     remove redundancy and be clearer and simplier.


194: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdateSubMatrices(), MatCreateTranspose(), MatGetSchurComplement(),
195:           MatSchurComplementGetPmat()

197: @*/
198: PetscErrorCode  MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *S)
199: {

203:   KSPInitializePackage();
204:   MatCreate(((PetscObject)A00)->comm,S);
205:   MatSetType(*S,MATSCHURCOMPLEMENT);
206:   MatSchurComplementSetSubMatrices(*S,A00,Ap00,A01,A10,A11);
207:   return(0);
208: }

210: /*@
211:       MatSchurComplementSetSubMatrices - Sets the matrices that define the Schur complement

213:    Collective on Mat

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

220:    Level: intermediate

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

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

229:           A00 and  A11 must be square matrices.

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

233: @*/
234: PetscErrorCode  MatSchurComplementSetSubMatrices(Mat S,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
235: {
236:   PetscErrorCode      ierr;
237:   PetscInt            m,n;
238:   Mat_SchurComplement *Na = (Mat_SchurComplement*)S->data;

241:   if (S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdateSubMatrices() for already used matrix");
249:   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);
250:   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);
251:   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);
252:   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);
253:   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);
254:   if (A11) {
257:     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);
258:   }

260:   MatGetLocalSize(A01,NULL,&n);
261:   MatGetLocalSize(A10,&m,NULL);
262:   MatSetSizes(S,m,n,PETSC_DECIDE,PETSC_DECIDE);
263:   PetscObjectReference((PetscObject)A00);
264:   PetscObjectReference((PetscObject)Ap00);
265:   PetscObjectReference((PetscObject)A01);
266:   PetscObjectReference((PetscObject)A10);
267:   Na->A  = A00;
268:   Na->Ap = Ap00;
269:   Na->B  = A01;
270:   Na->C  = A10;
271:   Na->D  = A11;
272:   if (A11) {
273:     PetscObjectReference((PetscObject)A11);
274:   }
275:   S->assembled    = PETSC_TRUE;
276:   S->preallocated = PETSC_TRUE;

278:   PetscLayoutSetUp((S)->rmap);
279:   PetscLayoutSetUp((S)->cmap);
280:   KSPSetOperators(Na->ksp,A00,Ap00);
281:   return(0);
282: }

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

287:   Not Collective

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

292:   Output Parameter:
293: . ksp - the linear solver object

295:   Options Database:
296: . -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.

298:   Level: intermediate

300: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
301: @*/
302: PetscErrorCode MatSchurComplementGetKSP(Mat S, KSP *ksp)
303: {
304:   Mat_SchurComplement *Na;

309:   Na   = (Mat_SchurComplement*) S->data;
310:   *ksp = Na->ksp;
311:   return(0);
312: }

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

317:   Not Collective

319:   Input Parameters:
320: + S   - matrix created with MatCreateSchurComplement()
321: - ksp - the linear solver object

323:   Level: developer

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

328: .seealso: MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate(), MATSCHURCOMPLEMENT
329: @*/
330: PetscErrorCode MatSchurComplementSetKSP(Mat S, KSP ksp)
331: {
332:   Mat_SchurComplement *Na;
333:   PetscErrorCode      ierr;

338:   Na      = (Mat_SchurComplement*) S->data;
339:   PetscObjectReference((PetscObject)ksp);
340:   KSPDestroy(&Na->ksp);
341:   Na->ksp = ksp;
342:   KSPSetOperators(Na->ksp, Na->A, Na->Ap);
343:   return(0);
344: }

346: /*@
347:       MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices

349:    Collective on Mat

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

356:    Level: intermediate

358:    Notes:
359:     All four matrices must have the same MPI communicator

361:           A00 and  A11 must be square matrices

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

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

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

375:   if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
382:   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);
383:   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);
384:   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);
385:   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);
386:   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);
387:   if (A11) {
390:     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);
391:   }

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

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

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

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


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

421:   Collective on Mat

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

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

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

432:   Level: intermediate

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

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

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

463: /*@
464:   MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

466:   Collective on Mat

468:   Input Parameter:
469: . M - the matrix obtained with MatCreateSchurComplement()

471:   Output Parameter:
472: . S - the Schur complement matrix

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

476:   Level: advanced

478: .seealso: MatCreateSchurComplement(), MatSchurComplementUpdate()
479: @*/
480: PetscErrorCode MatSchurComplementComputeExplicitOperator(Mat M, Mat *S)
481: {
482:   Mat            B, C, D;
483:   KSP            ksp;
484:   PC             pc;
485:   PetscBool      isLU, isILU;
486:   PetscReal      fill = 2.0;

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

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

514:     PCComputeExplicitOperator(pc, &Ainvd);
515:     MatConvert(Ainvd, MATAIJ, MAT_INITIAL_MATRIX, &Ainv);
516:     MatDestroy(&Ainvd);
517: #if 0
518:     /* Symmetric version */
519:     MatPtAP(Ainv, B, MAT_INITIAL_MATRIX, fill, S);
520: #else
521:     /* Nonsymmetric version */
522:     MatMatMatMult(C, Ainv, B, MAT_INITIAL_MATRIX, fill, S);
523: #endif
524:     MatDestroy(&Ainv);
525:   }
526:   if (D) {
527:     MatAXPY(*S, -1.0, D, DIFFERENT_NONZERO_PATTERN);
528:    }
529:   MatScale(*S,-1.0);
530:   return(0);
531: }

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

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

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

554:   reuse = MAT_INITIAL_MATRIX;
555:   if (mreuse == MAT_REUSE_MATRIX) {
556:     MatSchurComplementGetSubMatrices(*newmat,&A,&Ap,&B,&C,&D);
557:     if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
558:     if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
559:     MatDestroy(&Ap); /* get rid of extra reference */
560:     reuse = MAT_REUSE_MATRIX;
561:   }
562:   MatCreateSubMatrix(mat,isrow0,iscol0,reuse,&A);
563:   MatCreateSubMatrix(mat,isrow0,iscol1,reuse,&B);
564:   MatCreateSubMatrix(mat,isrow1,iscol0,reuse,&C);
565:   MatCreateSubMatrix(mat,isrow1,iscol1,reuse,&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:     if (mreuse != MAT_IGNORE_MATRIX) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse");
575:   }
576:   if (preuse != MAT_IGNORE_MATRIX) {
577:     MatCreateSchurComplementPmat(A,B,C,D,ainvtype,preuse,newpmat);
578:   }
579:   MatDestroy(&A);
580:   MatDestroy(&B);
581:   MatDestroy(&C);
582:   MatDestroy(&D);
583:   return(0);
584: }

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

589:     Collective on Mat

591:     Input Parameters:
592: +   A      - matrix in which the complement is to be taken
593: .   isrow0 - rows to eliminate
594: .   iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
595: .   isrow1 - rows in which the Schur complement is formed
596: .   iscol1 - columns in which the Schur complement is formed
597: .   mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in S
598: .   ainvtype - the type of approximation used for the inverse of the (0,0) block used in forming Sp:
599:                        MAT_SCHUR_COMPLEMENT_AINV_DIAG, MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG, or MAT_SCHUR_COMPLEMENT_AINV_LUMP
600: -   preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in Sp

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

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

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

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

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

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

624:     Developer Notes:
625:     The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
626:     remove redundancy and be clearer and simplier.

628:     Level: advanced

630:     Concepts: matrices^submatrices

632: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement(), MatSchurComplementAinvType
633: @*/
634: PetscErrorCode  MatGetSchurComplement(Mat A,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *S,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Sp)
635: {
636:   PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*) = NULL;

647:   if (A->factortype) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix");
648:   f = NULL;
649:   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. */
650:     PetscObjectQueryFunction((PetscObject)*S,"MatGetSchurComplement_C",&f);
651:   }
652:   if (f) {
653:       (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
654:   } else {
655:     MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
656:   }
657:   return(0);
658: }

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

663:     Not collective.

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

670:     Options database:
671:     -mat_schur_complement_ainv_type diag | lump | blockdiag

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

679:     Level: advanced

681:     Concepts: matrices^submatrices

683: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementGetAinvType()
684: @*/
685: PetscErrorCode  MatSchurComplementSetAinvType(Mat S,MatSchurComplementAinvType ainvtype)
686: {
687:   PetscErrorCode      ierr;
688:   const char*         t;
689:   PetscBool           isschur;
690:   Mat_SchurComplement *schur;

694:   PetscObjectGetType((PetscObject)S,&t);
695:   PetscStrcmp(t,MATSCHURCOMPLEMENT,&isschur);
696:   if (!isschur) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected Mat of type MATSCHURCOMPLEMENT, got %s instead",t);
697:   schur = (Mat_SchurComplement*)S->data;
698:   if (ainvtype != MAT_SCHUR_COMPLEMENT_AINV_DIAG && ainvtype != MAT_SCHUR_COMPLEMENT_AINV_LUMP && ainvtype != MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %d",(int)ainvtype);
699:   schur->ainvtype = ainvtype;
700:   return(0);
701: }

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

706:     Not collective.

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

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

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

721:     Level: advanced

723:     Concepts: matrices^submatrices

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

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

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

747:     Collective on Mat

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

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

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

763:     Level: advanced

765:     Concepts: matrices^submatrices

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

773:   PetscInt       N00;

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

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

782:   /* A zero size A00 or empty A01 or A10 imply S = A11. */
783:   MatGetSize(A00,&N00,NULL);
784:   if (!A01 || !A10 || !N00) {
785:     if (preuse == MAT_INITIAL_MATRIX) {
786:       MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
787:     } else { /* MAT_REUSE_MATRIX */
788:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
789:       MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
790:     }

792:   } else {
793:     Mat         AdB;
794:     Vec         diag;

796:     if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP || ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
797:       MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
798:       MatCreateVecs(A00,&diag,NULL);
799:       if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
800:         MatGetRowSum(A00,diag);
801:       } else {
802:         MatGetDiagonal(A00,diag);
803:       }
804:       VecReciprocal(diag);
805:       MatDiagonalScale(AdB,diag,NULL);
806:       VecDestroy(&diag);
807:     } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_BLOCK_DIAG) {
808:       Mat      A00_inv;
809:       MatType  type;
810:       MPI_Comm comm;

812:       PetscObjectGetComm((PetscObject)A00,&comm);
813:       MatGetType(A00,&type);
814:       MatCreate(comm,&A00_inv);
815:       MatSetType(A00_inv,type);
816:       MatInvertBlockDiagonalMat(A00,A00_inv);
817:       MatMatMult(A00_inv,A01,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&AdB);
818:       MatDestroy(&A00_inv);
819:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
820:     /* Cannot really reuse Spmat in MatMatMult() because of MatAYPX() -->
821:          MatAXPY() --> MatHeaderReplace() --> MatDestroy_XXX_MatMatMult()  */
822:     MatDestroy(Spmat);
823:     MatMatMult(A10,AdB,MAT_INITIAL_MATRIX,PETSC_DEFAULT,Spmat);
824:     if (!A11) {
825:       MatScale(*Spmat,-1.0);
826:     } else {
827:       /* TODO: when can we pass SAME_NONZERO_PATTERN? */
828:       MatAYPX(*Spmat,-1,A11,DIFFERENT_NONZERO_PATTERN);
829:     }
830:     MatDestroy(&AdB);
831:   }
832:   return(0);
833: }

835: PetscErrorCode  MatSchurComplementGetPmat_Basic(Mat S,MatReuse preuse,Mat *Spmat)
836: {
837:   Mat A,B,C,D;
838:   Mat_SchurComplement *schur = (Mat_SchurComplement *)S->data;
839:   PetscErrorCode      ierr;

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

844:   MatSchurComplementGetSubMatrices(S,&A,NULL,&B,&C,&D);
845:   if (!A) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Schur complement component matrices unset");
846:   MatCreateSchurComplementPmat(A,B,C,D,schur->ainvtype,preuse,Spmat);
847:   return(0);
848: }

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

853:     Collective on Mat

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

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

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

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

873:     Developer Notes:
874:     The API that includes MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementGetPmat() should be refactored to
875:     remove redundancy and be clearer and simplier.

877:     Level: advanced

879:     Concepts: matrices^submatrices

881: .seealso: MatCreateSubMatrix(), PCFIELDSPLIT, MatGetSchurComplement(), MatCreateSchurComplement(), MatSchurComplementSetAinvType()
882: @*/
883: PetscErrorCode  MatSchurComplementGetPmat(Mat S,MatReuse preuse,Mat *Sp)
884: {
885:   PetscErrorCode ierr,(*f)(Mat,MatReuse,Mat*);

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

893:   PetscObjectQueryFunction((PetscObject)S,"MatSchurComplementGetPmat_C",&f);
894:   if (f) {
895:     (*f)(S,preuse,Sp);
896:   } else {
897:     MatSchurComplementGetPmat_Basic(S,preuse,Sp);
898:   }
899:   return(0);
900: }

902: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
903: {
904:   PetscErrorCode      ierr;
905:   Mat_SchurComplement *Na;

908:   PetscNewLog(N,&Na);
909:   N->data = (void*) Na;

911:   N->ops->destroy        = MatDestroy_SchurComplement;
912:   N->ops->getvecs        = MatCreateVecs_SchurComplement;
913:   N->ops->view           = MatView_SchurComplement;
914:   N->ops->mult           = MatMult_SchurComplement;
915:   N->ops->multtranspose  = MatMultTranspose_SchurComplement;
916:   N->ops->multadd        = MatMultAdd_SchurComplement;
917:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
918:   N->assembled           = PETSC_FALSE;
919:   N->preallocated        = PETSC_FALSE;

921:   KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
922:   PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
923:   return(0);
924: }

926: static PetscBool KSPMatRegisterAllCalled;

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

931:   Not Collective

933:   Level: advanced

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

937: .seealso: MatRegisterAll(),  KSPInitializePackage()
938: @*/
939: PetscErrorCode KSPMatRegisterAll(void)
940: {

944:   if (KSPMatRegisterAllCalled) return(0);
945:   KSPMatRegisterAllCalled = PETSC_TRUE;
946:   MatRegister(MATSCHURCOMPLEMENT,MatCreate_SchurComplement);
947:   return(0);
948: }