Actual source code: schurm.c

petsc-dev 2014-08-22
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: }

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

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

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

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

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

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

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

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

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

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

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

183:    Collective on Mat

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

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

192:    Level: intermediate

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

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

200:           A00 and  A11 must be square matrices.

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

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

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

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

222:    Collective on Mat

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

229:    Level: intermediate

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

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

237:           A00 and  A11 must be square matrices.

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

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

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

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

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

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

298:   Not Collective

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

303:   Output Parameter:
304: . ksp - the linear solver object

306:   Options Database:
307: . -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.

309:   Level: intermediate

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

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

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

330:   Not Collective

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

336:   Level: developer

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

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

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

361: /*@
362:       MatSchurComplementUpdateSubMatrices - Updates the Schur complement matrix object with new submatrices

364:    Collective on Mat

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

371:    Level: intermediate

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

375:           A00 and  A11 must be square matrices

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

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

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

389:   if (!S->assembled) SETERRQ(PetscObjectComm((PetscObject)S),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSetSubMatrices() for a new matrix");
396:   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);
397:   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);
398:   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);
399:   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);
400:   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);
401:   if (A11) {
404:     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);
405:     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);
406:   }

408:   PetscObjectReference((PetscObject)A00);
409:   PetscObjectReference((PetscObject)Ap00);
410:   PetscObjectReference((PetscObject)A01);
411:   PetscObjectReference((PetscObject)A10);
412:   if (A11) {
413:     PetscObjectReference((PetscObject)A11);
414:   }

416:   MatDestroy(&Na->A);
417:   MatDestroy(&Na->Ap);
418:   MatDestroy(&Na->B);
419:   MatDestroy(&Na->C);
420:   MatDestroy(&Na->D);

422:   Na->A  = A00;
423:   Na->Ap = Ap00;
424:   Na->B  = A01;
425:   Na->C  = A10;
426:   Na->D  = A11;

428:   KSPSetOperators(Na->ksp,A00,Ap00);
429:   return(0);
430: }


435: /*@C
436:   MatSchurComplementGetSubMatrices - Get the individual submatrices in the Schur complement

438:   Collective on Mat

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

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

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

449:   Level: intermediate

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

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

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

482: /*@
483:   MatSchurComplementComputeExplicitOperator - Compute the Schur complement matrix explicitly

485:   Collective on Mat

487:   Input Parameter:
488: . M - the matrix obtained with MatCreateSchurComplement()

490:   Output Parameter:
491: . S - the Schur complement matrix

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

495:   Level: advanced

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

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

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

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

546:   PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_INFO);
547:   MatView(*S, PETSC_VIEWER_STDOUT_WORLD);
548:   PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);

550:   if (D) {
551:     MatInfo info;

553:     MatGetInfo(D, MAT_GLOBAL_SUM, &info);
554:     if (info.nz_used) SETERRQ(PetscObjectComm((PetscObject) M), PETSC_ERR_SUP, "Not yet implemented");
555:   }
556:   return(0);
557: }

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

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

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

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

616:     Collective on Mat

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

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

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

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

644:     Level: advanced

646:     Concepts: matrices^submatrices

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

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

665:   PetscObjectQueryFunction((PetscObject)S,"MatGetSchurComplement_C",&f);
666:   if (f) {
667:     (*f)(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,preuse,Sp);
668:   } else {
669:     MatGetSchurComplement_Basic(A,isrow0,iscol0,isrow1,iscol1,mreuse,S,ainvtype,preuse,Sp);
670:   }
671:   return(0);
672: }

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

679:     Not collective.

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

686:     Options database:
687:     -mat_schur_complement_ainv_type diag | lump

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

695:     Level: advanced

697:     Concepts: matrices^submatrices

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

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

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

724:     Not collective.

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

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

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

739:     Level: advanced

741:     Concepts: matrices^submatrices

743: .seealso: MatSchurComplementAinvType, MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementSetAinvType()
744: @*/
745: PetscErrorCode  MatSchurComplementGetAinvType(Mat S,MatSchurComplementAinvType *ainvtype)
746: {
747:   PetscErrorCode      ierr;
748:   const char*         t;
749:   PetscBool           isschur;
750:   Mat_SchurComplement *schur;

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

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

767:     Collective on Mat

769:     Input Parameters:
770: +   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)
771: .   ainvtype             - type of approximation for inv(A00) used when forming Sp = A11 - A10 inv(A00) A01
772: -   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

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

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

783:     Level: advanced

785:     Concepts: matrices^submatrices

787: .seealso: MatCreateSchurComplement(), MatGetSchurComplement(), MatSchurComplementGetPmat(), MatSchurComplementAinvType
788: @*/
789: PetscErrorCode  MatCreateSchurComplementPmat(Mat A00,Mat A01,Mat A10,Mat A11,MatSchurComplementAinvType ainvtype,MatReuse preuse,Mat *Spmat)
790: {

793:   PetscInt       N00;

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

800:   /* A zero size A00 or empty A01 or A10 imply S = A11. */
801:   MatGetSize(A00,&N00,NULL);
802:   if (!A01 || !A10 || !N00) {
803:     if (!preuse) {
804:       MatDuplicate(A11,MAT_COPY_VALUES,Spmat);
805:     } else {
806:       MatCopy(A11,*Spmat,DIFFERENT_NONZERO_PATTERN);
807:     }

809:   } else {
810:     Mat         AdB,Sp;
811:     Vec         diag;

813:     MatGetVecs(A00,&diag,NULL);
814:     if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_LUMP) {
815:       MatGetRowSum(A00,diag);
816:     } else if (ainvtype == MAT_SCHUR_COMPLEMENT_AINV_DIAG) {
817:       MatGetDiagonal(A00,diag);
818:     } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown MatSchurComplementAinvType: %D", ainvtype);
819:     VecReciprocal(diag);
820:     MatDuplicate(A01,MAT_COPY_VALUES,&AdB);
821:     MatDiagonalScale(AdB,diag,NULL);
822:     VecDestroy(&diag);
823:     Sp       = (preuse == MAT_REUSE_MATRIX) ? *Spmat : (Mat)0;
824:     MatMatMult(A10,AdB,preuse,PETSC_DEFAULT,&Sp);
825:     if (!A11) {
826:       MatScale(Sp,-1.0);
827:     } else {
828:       MatAYPX(Sp,-1,A11,DIFFERENT_NONZERO_PATTERN);
829:     }
830:     *Spmat    = Sp;
831:     MatDestroy(&AdB);
832:   }
833:   return(0);
834: }

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

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

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

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

858:     Collective on Mat

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

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

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

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

878:     Level: advanced

880:     Concepts: matrices^submatrices

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

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

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

905: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
906: {
907:   PetscErrorCode      ierr;
908:   Mat_SchurComplement *Na;

911:   PetscNewLog(N,&Na);
912:   N->data = (void*) Na;

914:   N->ops->destroy        = MatDestroy_SchurComplement;
915:   N->ops->getvecs        = MatGetVecs_SchurComplement;
916:   N->ops->view           = MatView_SchurComplement;
917:   N->ops->mult           = MatMult_SchurComplement;
918:   N->ops->multtranspose  = MatMultTranspose_SchurComplement;
919:   N->ops->multadd        = MatMultAdd_SchurComplement;
920:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
921:   N->assembled           = PETSC_FALSE;
922:   N->preallocated        = PETSC_FALSE;

924:   KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
925:   PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
926:   return(0);
927: }

929: static PetscBool KSPMatRegisterAllCalled;

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

936:   Not Collective

938:   Level: advanced

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

942: .seealso: MatRegisterAll(), MatRegisterDestroy(), KSPInitializePackage()
943: @*/
944: PetscErrorCode KSPMatRegisterAll()
945: {

949:   if (KSPMatRegisterAllCalled) return(0);
950:   KSPMatRegisterAllCalled = PETSC_TRUE;
951:   MatRegister(MATSCHURCOMPLEMENT,MatCreate_SchurComplement);
952:   return(0);
953: }