Actual source code: schurm.c
petsc-3.4.5 2014-06-29
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: }