Actual source code: schurm.c

petsc-3.4.4 2014-03-13
  2: #include <petsc-private/matimpl.h>
  3: #include <petscksp.h>                 /*I "petscksp.h" I*/

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

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

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

 34: PetscErrorCode MatView_SchurComplement(Mat N,PetscViewer viewer)
 35: {
 36:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 37:   PetscErrorCode      ierr;

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


 65: /*
 66:            A11 - A10 ksp(A00,Ap00)  A01
 67: */
 70: PetscErrorCode MatMult_SchurComplement(Mat N,Vec x,Vec y)
 71: {
 72:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 73:   PetscErrorCode      ierr;

 76:   if (!Na->work1) {MatGetVecs(Na->A,&Na->work1,NULL);}
 77:   if (!Na->work2) {MatGetVecs(Na->A,&Na->work2,NULL);}
 78:   MatMult(Na->B,x,Na->work1);
 79:   KSPSolve(Na->ksp,Na->work1,Na->work2);
 80:   MatMult(Na->C,Na->work2,y);
 81:   VecScale(y,-1.0);
 82:   if (Na->D) {
 83:     MatMultAdd(Na->D,x,y,y);
 84:   }
 85:   return(0);
 86: }

 88: /*
 89:            A11 - A10 ksp(A00,Ap00)  A01
 90: */
 93: PetscErrorCode MatMultAdd_SchurComplement(Mat N,Vec x,Vec y,Vec z)
 94: {
 95:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
 96:   PetscErrorCode      ierr;

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

118: PetscErrorCode MatSetFromOptions_SchurComplement(Mat N)
119: {
120:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
121:   PetscErrorCode      ierr;

124:   KSPSetFromOptions(Na->ksp);
125:   return(0);
126: }

130: PetscErrorCode MatDestroy_SchurComplement(Mat N)
131: {
132:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;
133:   PetscErrorCode      ierr;

136:   MatDestroy(&Na->A);
137:   MatDestroy(&Na->Ap);
138:   MatDestroy(&Na->B);
139:   MatDestroy(&Na->C);
140:   MatDestroy(&Na->D);
141:   VecDestroy(&Na->work1);
142:   VecDestroy(&Na->work2);
143:   KSPDestroy(&Na->ksp);
144:   PetscFree(N->data);
145:   return(0);
146: }

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

153:    Collective on Mat

155:    Input Parameter:
156: .   A00,A01,A10,A11  - the four parts of the original matrix (A00 is optional)

158:    Output Parameter:
159: .   N - the matrix that the Schur complement A11 - A10 ksp(A00) A01

161:    Level: intermediate

163:    Notes: The Schur complement is NOT actually formed! Rather this
164:           object performs the matrix-vector product by using the the formula for
165:           the Schur complement and a KSP solver to approximate the action of inv(A)

167:           All four matrices must have the same MPI communicator

169:           A00 and  A11 must be square matrices

171: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdate(), MatCreateTranspose(), MatGetSchurComplement()

173: @*/
174: PetscErrorCode  MatCreateSchurComplement(Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11,Mat *N)
175: {

179: #if !defined(PETSC_USE_DYNAMIC_LIBRARIES)
180:   KSPInitializePackage();
181: #endif
182:   MatCreate(((PetscObject)A00)->comm,N);
183:   MatSetType(*N,MATSCHURCOMPLEMENT);
184:   MatSchurComplementSet(*N,A00,Ap00,A01,A10,A11);
185:   return(0);
186: }

190: /*@
191:       MatSchurComplementSet - Sets the matrices that define the Schur complement

193:    Collective on Mat

195:    Input Parameter:
196: +   N - matrix obtained with MatCreate() and MatSetType(MATSCHURCOMPLEMENT);
197: -   A00,A01,A10,A11  - the four parts of the original matrix (A00 is optional)

199:    Level: intermediate

201:    Notes: The Schur complement is NOT actually formed! Rather this
202:           object performs the matrix-vector product by using the the formula for
203:           the Schur complement and a KSP solver to approximate the action of inv(A)

205:           All four matrices must have the same MPI communicator

207:           A00 and  A11 must be square matrices

209: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatSchurComplementUpdate(), MatCreateTranspose(), MatGetSchurComplement()

211: @*/
212: PetscErrorCode  MatSchurComplementSet(Mat N,Mat A00,Mat Ap00,Mat A01,Mat A10,Mat A11)
213: {
214:   PetscErrorCode      ierr;
215:   PetscInt            m,n;
216:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;

219:   if (N->assembled) SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementUpdate() for already used matrix");
227:   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);
228:   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);
229:   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);
230:   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);
231:   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);
232:   if (A11) {
235:     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);
236:     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);
237:   }

239:   MatGetLocalSize(A01,NULL,&n);
240:   MatGetLocalSize(A10,&m,NULL);
241:   MatSetSizes(N,m,n,PETSC_DECIDE,PETSC_DECIDE);
242:   PetscObjectReference((PetscObject)A00);
243:   PetscObjectReference((PetscObject)Ap00);
244:   PetscObjectReference((PetscObject)A01);
245:   PetscObjectReference((PetscObject)A10);
246:   Na->A  = A00;
247:   Na->Ap = Ap00;
248:   Na->B  = A01;
249:   Na->C  = A10;
250:   Na->D  = A11;
251:   if (A11) {
252:     PetscObjectReference((PetscObject)A11);
253:   }
254:   N->assembled    = PETSC_TRUE;
255:   N->preallocated = PETSC_TRUE;

257:   PetscLayoutSetUp((N)->rmap);
258:   PetscLayoutSetUp((N)->cmap);
259:   KSPSetOperators(Na->ksp,A00,Ap00,SAME_NONZERO_PATTERN);
260:   return(0);
261: }

265: /*@
266:   MatSchurComplementGetKSP - Gets the KSP object that is used to invert A in the Schur complement matrix S = C A^{-1} B

268:   Not Collective

270:   Input Parameter:
271: . A - matrix created with MatCreateSchurComplement()

273:   Output Parameter:
274: . ksp - the linear solver object

276:   Options Database:
277: . -fieldsplit_0_XXX sets KSP and PC options for the A block solver inside the Schur complement

279:   Level: intermediate

281: .seealso: MatSchurComplementSetKSP(), MatCreateSchurComplement(), MatCreateNormal(), MatMult(), MatCreate()
282: @*/
283: PetscErrorCode MatSchurComplementGetKSP(Mat A, KSP *ksp)
284: {
285:   Mat_SchurComplement *Na;

290:   Na   = (Mat_SchurComplement*) A->data;
291:   *ksp = Na->ksp;
292:   return(0);
293: }

297: /*@
298:   MatSchurComplementSetKSP - Sets the KSP object that is used to invert A in the Schur complement matrix S = C A^{-1} B

300:   Not Collective

302:   Input Parameters:
303: + A - matrix created with MatCreateSchurComplement()
304: - ksp - the linear solver object

306:   Level: developer

308:   Developer Notes:
309:     There is likely no use case for this function.

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

321:   Na      = (Mat_SchurComplement*) A->data;
322:   PetscObjectReference((PetscObject)ksp);
323:   KSPDestroy(&Na->ksp);
324:   Na->ksp = ksp;
325:   KSPSetOperators(Na->ksp, Na->A, Na->Ap, SAME_NONZERO_PATTERN);
326:   return(0);
327: }

331: /*@
332:       MatSchurComplementUpdate - Updates the Schur complement matrix object with new submatrices

334:    Collective on Mat

336:    Input Parameters:
337: +   N - the matrix obtained with MatCreateSchurComplement()
338: .   A,B,C,D  - the four parts of the original matrix (D is optional)
339: -   str - either SAME_NONZERO_PATTERN,DIFFERENT_NONZERO_PATTERN,SAME_PRECONDITIONER


342:    Level: intermediate

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

346:           A and  D must be square matrices

348:           All of the matrices provided must have the same sizes as was used with MatCreateSchurComplement() or MatSchurComplementSet()
349:           though they need not be the same matrices

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

353: @*/
354: PetscErrorCode  MatSchurComplementUpdate(Mat N,Mat A,Mat Ap,Mat B,Mat C,Mat D,MatStructure str)
355: {
356:   PetscErrorCode      ierr;
357:   Mat_SchurComplement *Na = (Mat_SchurComplement*)N->data;

360:   if (!N->assembled) SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_ARG_WRONGSTATE,"Use MatSchurComplementSet() for new matrix");
367:   if (A->rmap->n != A->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A %D do not equal local columns %D",A->rmap->n,A->cmap->n);
368:   if (A->rmap->n != Ap->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of A %D do not equal local rows of Ap %D",A->rmap->n,Ap->rmap->n);
369:   if (Ap->rmap->n != Ap->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of Ap %D do not equal local columns %D",Ap->rmap->n,Ap->cmap->n);
370:   if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of A %D do not equal local rows of B %D",A->cmap->n,B->rmap->n);
371:   if (C->cmap->n != A->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local columns of C %D do not equal local rows of A %D",C->cmap->n,A->rmap->n);
372:   if (D) {
375:     if (D->rmap->n != D->cmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of D %D do not equal local columns %D",D->rmap->n,D->cmap->n);
376:     if (C->rmap->n != D->rmap->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local rows of C %D do not equal local rows D %D",C->rmap->n,D->rmap->n);
377:   }

379:   PetscObjectReference((PetscObject)A);
380:   PetscObjectReference((PetscObject)Ap);
381:   PetscObjectReference((PetscObject)B);
382:   PetscObjectReference((PetscObject)C);
383:   if (D) {
384:     PetscObjectReference((PetscObject)D);
385:   }

387:   MatDestroy(&Na->A);
388:   MatDestroy(&Na->Ap);
389:   MatDestroy(&Na->B);
390:   MatDestroy(&Na->C);
391:   MatDestroy(&Na->D);

393:   Na->A  = A;
394:   Na->Ap = Ap;
395:   Na->B  = B;
396:   Na->C  = C;
397:   Na->D  = D;

399:   KSPSetOperators(Na->ksp,A,Ap,str);
400:   return(0);
401: }


406: /*@C
407:   MatSchurComplementGetSubmatrices - Get the individual submatrices in the Schur complement

409:   Collective on Mat

411:   Input Parameters:
412: + N - the matrix obtained with MatCreateSchurComplement()
413: - A,B,C,D  - the four parts of the original matrix (D is optional)

415:   Note:
416:   D is optional, and thus can be NULL

418:   Level: intermediate

420: .seealso: MatCreateNormal(), MatMult(), MatCreate(), MatSchurComplementGetKSP(), MatCreateSchurComplement(), MatSchurComplementUpdate()
421: @*/
422: PetscErrorCode  MatSchurComplementGetSubmatrices(Mat N,Mat *A,Mat *Ap,Mat *B,Mat *C,Mat *D)
423: {
424:   Mat_SchurComplement *Na = (Mat_SchurComplement*) N->data;
425:   PetscErrorCode      ierr;
426:   PetscBool           flg;

430:   PetscObjectTypeCompare((PetscObject)N,MATSCHURCOMPLEMENT,&flg);
431:   if (flg) {
432:     if (A) *A = Na->A;
433:     if (Ap) *Ap = Na->Ap;
434:     if (B) *B = Na->B;
435:     if (C) *C = Na->C;
436:     if (D) *D = Na->D;
437:   } else {
438:     if (A) *A = 0;
439:     if (Ap) *Ap = 0;
440:     if (B) *B = 0;
441:     if (C) *C = 0;
442:     if (D) *D = 0;
443:   }
444:   return(0);
445: }

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

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

466:   if (mreuse != MAT_IGNORE_MATRIX) {
467:     /* Use MatSchurComplement */
468:     if (mreuse == MAT_REUSE_MATRIX) {
469:       MatSchurComplementGetSubmatrices(*newmat,&A,&Ap,&B,&C,&D);
470:       if (!A || !Ap || !B || !C) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Attempting to reuse matrix but Schur complement matrices unset");
471:       if (A != Ap) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"Preconditioning matrix does not match operator");
472:       MatDestroy(&Ap); /* get rid of extra reference */
473:     }
474:     MatGetSubMatrix(mat,isrow0,iscol0,mreuse,&A);
475:     MatGetSubMatrix(mat,isrow0,iscol1,mreuse,&B);
476:     MatGetSubMatrix(mat,isrow1,iscol0,mreuse,&C);
477:     MatGetSubMatrix(mat,isrow1,iscol1,mreuse,&D);
478:     switch (mreuse) {
479:     case MAT_INITIAL_MATRIX:
480:       MatCreateSchurComplement(A,A,B,C,D,newmat);
481:       break;
482:     case MAT_REUSE_MATRIX:
483:       MatSchurComplementUpdate(*newmat,A,A,B,C,D,DIFFERENT_NONZERO_PATTERN);
484:       break;
485:     default:
486:       SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Unrecognized value of mreuse");
487:     }
488:   }
489:   if (preuse != MAT_IGNORE_MATRIX) {
490:     /* Use the diagonal part of A to form D - C inv(diag(A)) B */
491:     Mat         Ad,AdB,S;
492:     Vec         diag;
493:     PetscInt    i,m,n,mstart,mend;
494:     PetscScalar *x;

496:     /* We could compose these with newpmat so that the matrices can be reused. */
497:     if (!A) {MatGetSubMatrix(mat,isrow0,iscol0,MAT_INITIAL_MATRIX,&A);}
498:     if (!B) {MatGetSubMatrix(mat,isrow0,iscol1,MAT_INITIAL_MATRIX,&B);}
499:     if (!C) {MatGetSubMatrix(mat,isrow1,iscol0,MAT_INITIAL_MATRIX,&C);}
500:     if (!D) {MatGetSubMatrix(mat,isrow1,iscol1,MAT_INITIAL_MATRIX,&D);}

502:     MatGetVecs(A,&diag,NULL);
503:     MatGetDiagonal(A,diag);
504:     VecReciprocal(diag);
505:     MatGetLocalSize(A,&m,&n);
506:     /* We need to compute S = D - C inv(diag(A)) B.  For row-oriented formats, it is easy to scale the rows of B and
507:      * for column-oriented formats the columns of C can be scaled.  Would skip creating a silly diagonal matrix. */
508:     MatCreate(PetscObjectComm((PetscObject)A),&Ad);
509:     MatSetSizes(Ad,m,n,PETSC_DETERMINE,PETSC_DETERMINE);
510:     MatSetOptionsPrefix(Ad,((PetscObject)mat)->prefix);
511:     MatAppendOptionsPrefix(Ad,"diag_");
512:     MatSetFromOptions(Ad);
513:     MatSeqAIJSetPreallocation(Ad,1,NULL);
514:     MatMPIAIJSetPreallocation(Ad,1,NULL,0,NULL);
515:     MatGetOwnershipRange(Ad,&mstart,&mend);
516:     VecGetArray(diag,&x);
517:     for (i=mstart; i<mend; i++) {
518:       MatSetValue(Ad,i,i,x[i-mstart],INSERT_VALUES);
519:     }
520:     VecRestoreArray(diag,&x);
521:     MatAssemblyBegin(Ad,MAT_FINAL_ASSEMBLY);
522:     MatAssemblyEnd(Ad,MAT_FINAL_ASSEMBLY);
523:     VecDestroy(&diag);

525:     MatMatMult(Ad,B,MAT_INITIAL_MATRIX,1,&AdB);
526:     S        = (preuse == MAT_REUSE_MATRIX) ? *newpmat : (Mat)0;
527:     MatMatMult(C,AdB,preuse,PETSC_DEFAULT,&S);
528:     MatAYPX(S,-1,D,DIFFERENT_NONZERO_PATTERN);
529:     *newpmat = S;
530:     MatDestroy(&Ad);
531:     MatDestroy(&AdB);
532:   }
533:   MatDestroy(&A);
534:   MatDestroy(&B);
535:   MatDestroy(&C);
536:   MatDestroy(&D);
537:   return(0);
538: }

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

545:     Collective on Mat

547:     Input Parameters:
548: +   mat - Matrix in which the complement is to be taken
549: .   isrow0 - rows to eliminate
550: .   iscol0 - columns to eliminate, (isrow0,iscol0) should be square and nonsingular
551: .   isrow1 - rows in which the Schur complement is formed
552: .   iscol1 - columns in which the Schur complement is formed
553: .   mreuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in newmat
554: -   preuse - MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX, use MAT_IGNORE_MATRIX to put nothing in newpmat

556:     Output Parameters:
557: +   newmat - exact Schur complement, often of type MATSCHURCOMPLEMENT which is difficult to use for preconditioning
558: -   newpmat - approximate Schur complement suitable for preconditioning

560:     Note:
561:     Since the real Schur complement is usually dense, providing a good approximation to newpmat usually requires
562:     application-specific information.  The default for assembled matrices is to use the diagonal of the (0,0) block
563:     which will rarely produce a scalable algorithm.

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

570:     Level: advanced

572:     Concepts: matrices^submatrices

574: .seealso: MatGetSubMatrix(), PCFIELDSPLIT, MatCreateSchurComplement()
575: @*/
576: PetscErrorCode  MatGetSchurComplement(Mat mat,IS isrow0,IS iscol0,IS isrow1,IS iscol1,MatReuse mreuse,Mat *newmat,MatReuse preuse,Mat *newpmat)
577: {
578:   PetscErrorCode ierr,(*f)(Mat,IS,IS,IS,IS,MatReuse,Mat*,MatReuse,Mat*);

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

591:   PetscObjectQueryFunction((PetscObject)mat,"MatGetSchurComplement_C",&f);
592:   if (f) {
593:     (*f)(mat,isrow0,iscol0,isrow1,iscol1,mreuse,newmat,preuse,newpmat);
594:   } else {
595:     MatGetSchurComplement_Basic(mat,isrow0,iscol0,isrow1,iscol1,mreuse,newmat,preuse,newpmat);
596:   }
597:   return(0);
598: }

602: PETSC_EXTERN PetscErrorCode MatCreate_SchurComplement(Mat N)
603: {
604:   PetscErrorCode      ierr;
605:   Mat_SchurComplement *Na;

608:   PetscNewLog(N,Mat_SchurComplement,&Na);
609:   N->data = (void*) Na;

611:   N->ops->destroy        = MatDestroy_SchurComplement;
612:   N->ops->getvecs        = MatGetVecs_SchurComplement;
613:   N->ops->view           = MatView_SchurComplement;
614:   N->ops->mult           = MatMult_SchurComplement;
615:   N->ops->multadd        = MatMultAdd_SchurComplement;
616:   N->ops->setfromoptions = MatSetFromOptions_SchurComplement;
617:   N->assembled           = PETSC_FALSE;
618:   N->preallocated        = PETSC_FALSE;

620:   KSPCreate(PetscObjectComm((PetscObject)N),&Na->ksp);
621:   PetscObjectChangeTypeName((PetscObject)N,MATSCHURCOMPLEMENT);
622:   return(0);
623: }

625: static PetscBool KSPMatRegisterAllCalled;

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

632:   Not Collective

634:   Level: advanced

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

638: .seealso: MatRegisterAll(), MatRegisterDestroy(), KSPInitializePackage()
639: @*/
640: PetscErrorCode KSPMatRegisterAll()
641: {

645:   if (KSPMatRegisterAllCalled) return(0);
646:   KSPMatRegisterAllCalled = PETSC_TRUE;
647:   MatRegister(MATSCHURCOMPLEMENT,MatCreate_SchurComplement);
648:   return(0);
649: }