Actual source code: precon.c
2: /*
3: The PC (preconditioner) interface routines, callable by users.
4: */
5: #include <petsc/private/pcimpl.h>
6: #include <petscdm.h>
8: /* Logging support */
9: PetscClassId PC_CLASSID;
10: PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_MatApply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
11: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks;
12: PetscInt PetscMGLevelId;
14: PetscErrorCode PCGetDefaultType_Private(PC pc,const char *type[])
15: {
16: PetscMPIInt size;
17: PetscBool hasop,flg1,flg2,set,flg3;
19: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
20: if (pc->pmat) {
21: MatHasOperation(pc->pmat,MATOP_GET_DIAGONAL_BLOCK,&hasop);
22: if (size == 1) {
23: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
24: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
25: MatIsSymmetricKnown(pc->pmat,&set,&flg3);
26: if (flg1 && (!flg2 || (set && flg3))) {
27: *type = PCICC;
28: } else if (flg2) {
29: *type = PCILU;
30: } else if (hasop) { /* likely is a parallel matrix run on one processor */
31: *type = PCBJACOBI;
32: } else {
33: *type = PCNONE;
34: }
35: } else {
36: if (hasop) {
37: *type = PCBJACOBI;
38: } else {
39: *type = PCNONE;
40: }
41: }
42: } else {
43: if (size == 1) {
44: *type = PCILU;
45: } else {
46: *type = PCBJACOBI;
47: }
48: }
49: return 0;
50: }
52: /*@
53: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
55: Collective on PC
57: Input Parameter:
58: . pc - the preconditioner context
60: Level: developer
62: Notes:
63: This allows a PC to be reused for a different sized linear system but using the same options that have been previously set in the PC
65: .seealso: PCCreate(), PCSetUp()
66: @*/
67: PetscErrorCode PCReset(PC pc)
68: {
70: if (pc->ops->reset) {
71: (*pc->ops->reset)(pc);
72: }
73: VecDestroy(&pc->diagonalscaleright);
74: VecDestroy(&pc->diagonalscaleleft);
75: MatDestroy(&pc->pmat);
76: MatDestroy(&pc->mat);
78: pc->setupcalled = 0;
79: return 0;
80: }
82: /*@C
83: PCDestroy - Destroys PC context that was created with PCCreate().
85: Collective on PC
87: Input Parameter:
88: . pc - the preconditioner context
90: Level: developer
92: .seealso: PCCreate(), PCSetUp()
93: @*/
94: PetscErrorCode PCDestroy(PC *pc)
95: {
96: if (!*pc) return 0;
98: if (--((PetscObject)(*pc))->refct > 0) {*pc = NULL; return 0;}
100: PCReset(*pc);
102: /* if memory was published with SAWs then destroy it */
103: PetscObjectSAWsViewOff((PetscObject)*pc);
104: if ((*pc)->ops->destroy) (*(*pc)->ops->destroy)((*pc));
105: DMDestroy(&(*pc)->dm);
106: PetscHeaderDestroy(pc);
107: return 0;
108: }
110: /*@C
111: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
112: scaling as needed by certain time-stepping codes.
114: Logically Collective on PC
116: Input Parameter:
117: . pc - the preconditioner context
119: Output Parameter:
120: . flag - PETSC_TRUE if it applies the scaling
122: Level: developer
124: Notes:
125: If this returns PETSC_TRUE then the system solved via the Krylov method is
126: $ D M A D^{-1} y = D M b for left preconditioning or
127: $ D A M D^{-1} z = D b for right preconditioning
129: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
130: @*/
131: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
132: {
135: *flag = pc->diagonalscale;
136: return 0;
137: }
139: /*@
140: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
141: scaling as needed by certain time-stepping codes.
143: Logically Collective on PC
145: Input Parameters:
146: + pc - the preconditioner context
147: - s - scaling vector
149: Level: intermediate
151: Notes:
152: The system solved via the Krylov method is
153: $ D M A D^{-1} y = D M b for left preconditioning or
154: $ D A M D^{-1} z = D b for right preconditioning
156: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
158: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
159: @*/
160: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
161: {
164: pc->diagonalscale = PETSC_TRUE;
166: PetscObjectReference((PetscObject)s);
167: VecDestroy(&pc->diagonalscaleleft);
169: pc->diagonalscaleleft = s;
171: VecDuplicate(s,&pc->diagonalscaleright);
172: VecCopy(s,pc->diagonalscaleright);
173: VecReciprocal(pc->diagonalscaleright);
174: return 0;
175: }
177: /*@
178: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
180: Logically Collective on PC
182: Input Parameters:
183: + pc - the preconditioner context
184: . in - input vector
185: - out - scaled vector (maybe the same as in)
187: Level: intermediate
189: Notes:
190: The system solved via the Krylov method is
191: $ D M A D^{-1} y = D M b for left preconditioning or
192: $ D A M D^{-1} z = D b for right preconditioning
194: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
196: If diagonal scaling is turned off and in is not out then in is copied to out
198: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
199: @*/
200: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
201: {
205: if (pc->diagonalscale) {
206: VecPointwiseMult(out,pc->diagonalscaleleft,in);
207: } else if (in != out) {
208: VecCopy(in,out);
209: }
210: return 0;
211: }
213: /*@
214: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
216: Logically Collective on PC
218: Input Parameters:
219: + pc - the preconditioner context
220: . in - input vector
221: - out - scaled vector (maybe the same as in)
223: Level: intermediate
225: Notes:
226: The system solved via the Krylov method is
227: $ D M A D^{-1} y = D M b for left preconditioning or
228: $ D A M D^{-1} z = D b for right preconditioning
230: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
232: If diagonal scaling is turned off and in is not out then in is copied to out
234: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
235: @*/
236: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
237: {
241: if (pc->diagonalscale) {
242: VecPointwiseMult(out,pc->diagonalscaleright,in);
243: } else if (in != out) {
244: VecCopy(in,out);
245: }
246: return 0;
247: }
249: /*@
250: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
251: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
252: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperators() or PCSetOperators() not the Pmat.
254: Logically Collective on PC
256: Input Parameters:
257: + pc - the preconditioner context
258: - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
260: Options Database Key:
261: . -pc_use_amat <true,false> - use the amat to apply the operator
263: Notes:
264: For the common case in which the linear system matrix and the matrix used to construct the
265: preconditioner are identical, this routine is does nothing.
267: Level: intermediate
269: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
270: @*/
271: PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
272: {
274: pc->useAmat = flg;
275: return 0;
276: }
278: /*@
279: PCSetErrorIfFailure - Causes PC to generate an error if a FPE, for example a zero pivot, is detected.
281: Logically Collective on PC
283: Input Parameters:
284: + pc - iterative context obtained from PCCreate()
285: - flg - PETSC_TRUE indicates you want the error generated
287: Level: advanced
289: Notes:
290: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
291: to determine if it has converged or failed. Or use -ksp_error_if_not_converged to cause the program to terminate as soon as lack of convergence is
292: detected.
294: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
296: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
297: @*/
298: PetscErrorCode PCSetErrorIfFailure(PC pc,PetscBool flg)
299: {
302: pc->erroriffailure = flg;
303: return 0;
304: }
306: /*@
307: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
308: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
309: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperators() or PCSetOperators() not the Pmat.
311: Logically Collective on PC
313: Input Parameter:
314: . pc - the preconditioner context
316: Output Parameter:
317: . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
319: Notes:
320: For the common case in which the linear system matrix and the matrix used to construct the
321: preconditioner are identical, this routine is does nothing.
323: Level: intermediate
325: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
326: @*/
327: PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
328: {
330: *flg = pc->useAmat;
331: return 0;
332: }
334: /*@
335: PCCreate - Creates a preconditioner context.
337: Collective
339: Input Parameter:
340: . comm - MPI communicator
342: Output Parameter:
343: . pc - location to put the preconditioner context
345: Notes:
346: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or PCICC
347: in parallel. For dense matrices it is always PCNONE.
349: Level: developer
351: .seealso: PCSetUp(), PCApply(), PCDestroy()
352: @*/
353: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
354: {
355: PC pc;
358: *newpc = NULL;
359: PCInitializePackage();
361: PetscHeaderCreate(pc,PC_CLASSID,"PC","Preconditioner","PC",comm,PCDestroy,PCView);
363: pc->mat = NULL;
364: pc->pmat = NULL;
365: pc->setupcalled = 0;
366: pc->setfromoptionscalled = 0;
367: pc->data = NULL;
368: pc->diagonalscale = PETSC_FALSE;
369: pc->diagonalscaleleft = NULL;
370: pc->diagonalscaleright = NULL;
372: pc->modifysubmatrices = NULL;
373: pc->modifysubmatricesP = NULL;
375: *newpc = pc;
376: return 0;
378: }
380: /* -------------------------------------------------------------------------------*/
382: /*@
383: PCApply - Applies the preconditioner to a vector.
385: Collective on PC
387: Input Parameters:
388: + pc - the preconditioner context
389: - x - input vector
391: Output Parameter:
392: . y - output vector
394: Level: developer
396: .seealso: PCApplyTranspose(), PCApplyBAorAB()
397: @*/
398: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
399: {
400: PetscInt m,n,mv,nv;
406: if (pc->erroriffailure) VecValidValues(x,2,PETSC_TRUE);
407: /* use pmat to check vector sizes since for KSPLSQR the pmat may be of a different size than mat */
408: MatGetLocalSize(pc->pmat,&m,&n);
409: VecGetLocalSize(x,&mv);
410: VecGetLocalSize(y,&nv);
411: /* check pmat * y = x is feasible */
414: VecSetErrorIfLocked(y,3);
416: PCSetUp(pc);
418: VecLockReadPush(x);
419: PetscLogEventBegin(PC_Apply,pc,x,y,0);
420: (*pc->ops->apply)(pc,x,y);
421: PetscLogEventEnd(PC_Apply,pc,x,y,0);
422: if (pc->erroriffailure) VecValidValues(y,3,PETSC_FALSE);
423: VecLockReadPop(x);
424: return 0;
425: }
427: /*@
428: PCMatApply - Applies the preconditioner to multiple vectors stored as a MATDENSE. Like PCApply(), Y and X must be different matrices.
430: Collective on PC
432: Input Parameters:
433: + pc - the preconditioner context
434: - X - block of input vectors
436: Output Parameter:
437: . Y - block of output vectors
439: Level: developer
441: .seealso: PCApply(), KSPMatSolve()
442: @*/
443: PetscErrorCode PCMatApply(PC pc,Mat X,Mat Y)
444: {
445: Mat A;
446: Vec cy, cx;
447: PetscInt m1, M1, m2, M2, n1, N1, n2, N2, m3, M3, n3, N3;
448: PetscBool match;
456: PCGetOperators(pc, NULL, &A);
457: MatGetLocalSize(A, &m3, &n3);
458: MatGetLocalSize(X, &m2, &n2);
459: MatGetLocalSize(Y, &m1, &n1);
460: MatGetSize(A, &M3, &N3);
461: MatGetSize(X, &M2, &N2);
462: MatGetSize(Y, &M1, &N1);
466: PetscObjectBaseTypeCompareAny((PetscObject)Y, &match, MATSEQDENSE, MATMPIDENSE, "");
468: PetscObjectBaseTypeCompareAny((PetscObject)X, &match, MATSEQDENSE, MATMPIDENSE, "");
470: PCSetUp(pc);
471: if (pc->ops->matapply) {
472: PetscLogEventBegin(PC_MatApply, pc, X, Y, 0);
473: (*pc->ops->matapply)(pc, X, Y);
474: PetscLogEventEnd(PC_MatApply, pc, X, Y, 0);
475: } else {
476: PetscInfo(pc, "PC type %s applying column by column\n", ((PetscObject)pc)->type_name);
477: for (n1 = 0; n1 < N1; ++n1) {
478: MatDenseGetColumnVecRead(X, n1, &cx);
479: MatDenseGetColumnVecWrite(Y, n1, &cy);
480: PCApply(pc, cx, cy);
481: MatDenseRestoreColumnVecWrite(Y, n1, &cy);
482: MatDenseRestoreColumnVecRead(X, n1, &cx);
483: }
484: }
485: return 0;
486: }
488: /*@
489: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
491: Collective on PC
493: Input Parameters:
494: + pc - the preconditioner context
495: - x - input vector
497: Output Parameter:
498: . y - output vector
500: Notes:
501: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
503: Level: developer
505: .seealso: PCApply(), PCApplySymmetricRight()
506: @*/
507: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
508: {
513: if (pc->erroriffailure) VecValidValues(x,2,PETSC_TRUE);
514: PCSetUp(pc);
516: VecLockReadPush(x);
517: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
518: (*pc->ops->applysymmetricleft)(pc,x,y);
519: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
520: VecLockReadPop(x);
521: if (pc->erroriffailure) VecValidValues(y,3,PETSC_FALSE);
522: return 0;
523: }
525: /*@
526: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
528: Collective on PC
530: Input Parameters:
531: + pc - the preconditioner context
532: - x - input vector
534: Output Parameter:
535: . y - output vector
537: Level: developer
539: Notes:
540: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
542: .seealso: PCApply(), PCApplySymmetricLeft()
543: @*/
544: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
545: {
550: if (pc->erroriffailure) VecValidValues(x,2,PETSC_TRUE);
551: PCSetUp(pc);
553: VecLockReadPush(x);
554: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
555: (*pc->ops->applysymmetricright)(pc,x,y);
556: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
557: VecLockReadPop(x);
558: if (pc->erroriffailure) VecValidValues(y,3,PETSC_FALSE);
559: return 0;
560: }
562: /*@
563: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
565: Collective on PC
567: Input Parameters:
568: + pc - the preconditioner context
569: - x - input vector
571: Output Parameter:
572: . y - output vector
574: Notes:
575: For complex numbers this applies the non-Hermitian transpose.
577: Developer Notes:
578: We need to implement a PCApplyHermitianTranspose()
580: Level: developer
582: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
583: @*/
584: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
585: {
590: if (pc->erroriffailure) VecValidValues(x,2,PETSC_TRUE);
591: PCSetUp(pc);
593: VecLockReadPush(x);
594: PetscLogEventBegin(PC_Apply,pc,x,y,0);
595: (*pc->ops->applytranspose)(pc,x,y);
596: PetscLogEventEnd(PC_Apply,pc,x,y,0);
597: VecLockReadPop(x);
598: if (pc->erroriffailure) VecValidValues(y,3,PETSC_FALSE);
599: return 0;
600: }
602: /*@
603: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
605: Collective on PC
607: Input Parameters:
608: . pc - the preconditioner context
610: Output Parameter:
611: . flg - PETSC_TRUE if a transpose operation is defined
613: Level: developer
615: .seealso: PCApplyTranspose()
616: @*/
617: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
618: {
621: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
622: else *flg = PETSC_FALSE;
623: return 0;
624: }
626: /*@
627: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
629: Collective on PC
631: Input Parameters:
632: + pc - the preconditioner context
633: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
634: . x - input vector
635: - work - work vector
637: Output Parameter:
638: . y - output vector
640: Level: developer
642: Notes:
643: If the PC has had PCSetDiagonalScale() set then D M A D^{-1} for left preconditioning or D A M D^{-1} is actually applied. Note that the
644: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
646: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
647: @*/
648: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
649: {
661: if (pc->erroriffailure) VecValidValues(x,3,PETSC_TRUE);
663: PCSetUp(pc);
664: if (pc->diagonalscale) {
665: if (pc->ops->applyBA) {
666: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
667: VecDuplicate(x,&work2);
668: PCDiagonalScaleRight(pc,x,work2);
669: (*pc->ops->applyBA)(pc,side,work2,y,work);
670: PCDiagonalScaleLeft(pc,y,y);
671: VecDestroy(&work2);
672: } else if (side == PC_RIGHT) {
673: PCDiagonalScaleRight(pc,x,y);
674: PCApply(pc,y,work);
675: MatMult(pc->mat,work,y);
676: PCDiagonalScaleLeft(pc,y,y);
677: } else if (side == PC_LEFT) {
678: PCDiagonalScaleRight(pc,x,y);
679: MatMult(pc->mat,y,work);
680: PCApply(pc,work,y);
681: PCDiagonalScaleLeft(pc,y,y);
683: } else {
684: if (pc->ops->applyBA) {
685: (*pc->ops->applyBA)(pc,side,x,y,work);
686: } else if (side == PC_RIGHT) {
687: PCApply(pc,x,work);
688: MatMult(pc->mat,work,y);
689: } else if (side == PC_LEFT) {
690: MatMult(pc->mat,x,work);
691: PCApply(pc,work,y);
692: } else if (side == PC_SYMMETRIC) {
693: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
694: PCApplySymmetricRight(pc,x,work);
695: MatMult(pc->mat,work,y);
696: VecCopy(y,work);
697: PCApplySymmetricLeft(pc,work,y);
698: }
699: }
700: if (pc->erroriffailure) VecValidValues(y,4,PETSC_FALSE);
701: return 0;
702: }
704: /*@
705: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
706: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
707: NOT tr(B*A) = tr(A)*tr(B).
709: Collective on PC
711: Input Parameters:
712: + pc - the preconditioner context
713: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
714: . x - input vector
715: - work - work vector
717: Output Parameter:
718: . y - output vector
720: Notes:
721: this routine is used internally so that the same Krylov code can be used to solve A x = b and A' x = b, with a preconditioner
722: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
724: Level: developer
726: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
727: @*/
728: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
729: {
735: if (pc->erroriffailure) VecValidValues(x,3,PETSC_TRUE);
736: if (pc->ops->applyBAtranspose) {
737: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
738: if (pc->erroriffailure) VecValidValues(y,4,PETSC_FALSE);
739: return 0;
740: }
743: PCSetUp(pc);
744: if (side == PC_RIGHT) {
745: PCApplyTranspose(pc,x,work);
746: MatMultTranspose(pc->mat,work,y);
747: } else if (side == PC_LEFT) {
748: MatMultTranspose(pc->mat,x,work);
749: PCApplyTranspose(pc,work,y);
750: }
751: /* add support for PC_SYMMETRIC */
752: if (pc->erroriffailure) VecValidValues(y,4,PETSC_FALSE);
753: return 0;
754: }
756: /* -------------------------------------------------------------------------------*/
758: /*@
759: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
760: built-in fast application of Richardson's method.
762: Not Collective
764: Input Parameter:
765: . pc - the preconditioner
767: Output Parameter:
768: . exists - PETSC_TRUE or PETSC_FALSE
770: Level: developer
772: .seealso: PCApplyRichardson()
773: @*/
774: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
775: {
778: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
779: else *exists = PETSC_FALSE;
780: return 0;
781: }
783: /*@
784: PCApplyRichardson - Applies several steps of Richardson iteration with
785: the particular preconditioner. This routine is usually used by the
786: Krylov solvers and not the application code directly.
788: Collective on PC
790: Input Parameters:
791: + pc - the preconditioner context
792: . b - the right hand side
793: . w - one work vector
794: . rtol - relative decrease in residual norm convergence criteria
795: . abstol - absolute residual norm convergence criteria
796: . dtol - divergence residual norm increase criteria
797: . its - the number of iterations to apply.
798: - guesszero - if the input x contains nonzero initial guess
800: Output Parameters:
801: + outits - number of iterations actually used (for SOR this always equals its)
802: . reason - the reason the apply terminated
803: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
805: Notes:
806: Most preconditioners do not support this function. Use the command
807: PCApplyRichardsonExists() to determine if one does.
809: Except for the multigrid PC this routine ignores the convergence tolerances
810: and always runs for the number of iterations
812: Level: developer
814: .seealso: PCApplyRichardsonExists()
815: @*/
816: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
817: {
823: PCSetUp(pc);
825: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
826: return 0;
827: }
829: /*@
830: PCSetFailedReason - Sets the reason a PCSetUp() failed or PC_NOERROR if it did not fail
832: Logically Collective on PC
834: Input Parameters:
835: + pc - the preconditioner context
836: - reason - the reason it failedx
838: Level: advanced
840: .seealso: PCCreate(), PCApply(), PCDestroy(), PCFailedReason
841: @*/
842: PetscErrorCode PCSetFailedReason(PC pc,PCFailedReason reason)
843: {
844: pc->failedreason = reason;
845: return 0;
846: }
848: /*@
849: PCGetFailedReason - Gets the reason a PCSetUp() failed or PC_NOERROR if it did not fail
851: Logically Collective on PC
853: Input Parameter:
854: . pc - the preconditioner context
856: Output Parameter:
857: . reason - the reason it failed
859: Level: advanced
861: Notes: This is the maximum over reason over all ranks in the PC communicator. It is only valid after
862: a call KSPCheckDot() or KSPCheckNorm() inside a KSPSolve(). It is not valid immediately after a PCSetUp()
863: or PCApply(), then use PCGetFailedReasonRank()
865: .seealso: PCCreate(), PCApply(), PCDestroy(), PCGetFailedReasonRank(), PCSetFailedReason()
866: @*/
867: PetscErrorCode PCGetFailedReason(PC pc,PCFailedReason *reason)
868: {
869: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
870: else *reason = pc->failedreason;
871: return 0;
872: }
874: /*@
875: PCGetFailedReasonRank - Gets the reason a PCSetUp() failed or PC_NOERROR if it did not fail on this MPI rank
877: Not Collective on PC
879: Input Parameter:
880: . pc - the preconditioner context
882: Output Parameter:
883: . reason - the reason it failed
885: Notes:
886: Different ranks may have different reasons or no reason, see PCGetFailedReason()
888: Level: advanced
890: .seealso: PCCreate(), PCApply(), PCDestroy(), PCGetFailedReason(), PCSetFailedReason()
891: @*/
892: PetscErrorCode PCGetFailedReasonRank(PC pc,PCFailedReason *reason)
893: {
894: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
895: else *reason = pc->failedreason;
896: return 0;
897: }
899: /* Next line needed to deactivate KSP_Solve logging */
900: #include <petsc/private/kspimpl.h>
902: /*
903: a setupcall of 0 indicates never setup,
904: 1 indicates has been previously setup
905: -1 indicates a PCSetUp() was attempted and failed
906: */
907: /*@
908: PCSetUp - Prepares for the use of a preconditioner.
910: Collective on PC
912: Input Parameter:
913: . pc - the preconditioner context
915: Level: developer
917: .seealso: PCCreate(), PCApply(), PCDestroy()
918: @*/
919: PetscErrorCode PCSetUp(PC pc)
920: {
921: const char *def;
922: PetscObjectState matstate, matnonzerostate;
927: if (pc->setupcalled && pc->reusepreconditioner) {
928: PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
929: return 0;
930: }
932: PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
933: MatGetNonzeroState(pc->pmat,&matnonzerostate);
934: if (!pc->setupcalled) {
935: PetscInfo(pc,"Setting up PC for first time\n");
936: pc->flag = DIFFERENT_NONZERO_PATTERN;
937: } else if (matstate == pc->matstate) {
938: PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
939: return 0;
940: } else {
941: if (matnonzerostate > pc->matnonzerostate) {
942: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
943: pc->flag = DIFFERENT_NONZERO_PATTERN;
944: } else {
945: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
946: pc->flag = SAME_NONZERO_PATTERN;
947: }
948: }
949: pc->matstate = matstate;
950: pc->matnonzerostate = matnonzerostate;
952: if (!((PetscObject)pc)->type_name) {
953: PCGetDefaultType_Private(pc,&def);
954: PCSetType(pc,def);
955: }
957: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
958: MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
959: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
960: if (pc->ops->setup) {
961: /* do not log solves and applications of preconditioners while constructing preconditioners; perhaps they should be logged separately from the regular solves */
962: KSPInitializePackage();
963: PetscLogEventDeactivatePush(KSP_Solve);
964: PetscLogEventDeactivatePush(PC_Apply);
965: (*pc->ops->setup)(pc);
966: PetscLogEventDeactivatePop(KSP_Solve);
967: PetscLogEventDeactivatePop(PC_Apply);
968: }
969: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
970: if (!pc->setupcalled) pc->setupcalled = 1;
971: return 0;
972: }
974: /*@
975: PCSetUpOnBlocks - Sets up the preconditioner for each block in
976: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
977: methods.
979: Collective on PC
981: Input Parameters:
982: . pc - the preconditioner context
984: Level: developer
986: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
987: @*/
988: PetscErrorCode PCSetUpOnBlocks(PC pc)
989: {
991: if (!pc->ops->setuponblocks) return 0;
992: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
993: (*pc->ops->setuponblocks)(pc);
994: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
995: return 0;
996: }
998: /*@C
999: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1000: submatrices that arise within certain subdomain-based preconditioners.
1001: The basic submatrices are extracted from the preconditioner matrix as
1002: usual; the user can then alter these (for example, to set different boundary
1003: conditions for each submatrix) before they are used for the local solves.
1005: Logically Collective on PC
1007: Input Parameters:
1008: + pc - the preconditioner context
1009: . func - routine for modifying the submatrices
1010: - ctx - optional user-defined context (may be null)
1012: Calling sequence of func:
1013: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
1015: + row - an array of index sets that contain the global row numbers
1016: that comprise each local submatrix
1017: . col - an array of index sets that contain the global column numbers
1018: that comprise each local submatrix
1019: . submat - array of local submatrices
1020: - ctx - optional user-defined context for private data for the
1021: user-defined func routine (may be null)
1023: Notes:
1024: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
1025: KSPSolve().
1027: A routine set by PCSetModifySubMatrices() is currently called within
1028: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
1029: preconditioners. All other preconditioners ignore this routine.
1031: Level: advanced
1033: .seealso: PCModifySubMatrices()
1034: @*/
1035: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1036: {
1038: pc->modifysubmatrices = func;
1039: pc->modifysubmatricesP = ctx;
1040: return 0;
1041: }
1043: /*@C
1044: PCModifySubMatrices - Calls an optional user-defined routine within
1045: certain preconditioners if one has been set with PCSetModifySubMatrices().
1047: Collective on PC
1049: Input Parameters:
1050: + pc - the preconditioner context
1051: . nsub - the number of local submatrices
1052: . row - an array of index sets that contain the global row numbers
1053: that comprise each local submatrix
1054: . col - an array of index sets that contain the global column numbers
1055: that comprise each local submatrix
1056: . submat - array of local submatrices
1057: - ctx - optional user-defined context for private data for the
1058: user-defined routine (may be null)
1060: Output Parameter:
1061: . submat - array of local submatrices (the entries of which may
1062: have been modified)
1064: Notes:
1065: The user should NOT generally call this routine, as it will
1066: automatically be called within certain preconditioners (currently
1067: block Jacobi, additive Schwarz) if set.
1069: The basic submatrices are extracted from the preconditioner matrix
1070: as usual; the user can then alter these (for example, to set different
1071: boundary conditions for each submatrix) before they are used for the
1072: local solves.
1074: Level: developer
1076: .seealso: PCSetModifySubMatrices()
1077: @*/
1078: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1079: {
1081: if (!pc->modifysubmatrices) return 0;
1082: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1083: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1084: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1085: return 0;
1086: }
1088: /*@
1089: PCSetOperators - Sets the matrix associated with the linear system and
1090: a (possibly) different one associated with the preconditioner.
1092: Logically Collective on PC
1094: Input Parameters:
1095: + pc - the preconditioner context
1096: . Amat - the matrix that defines the linear system
1097: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1099: Notes:
1100: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1102: If you wish to replace either Amat or Pmat but leave the other one untouched then
1103: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1104: on it and then pass it back in in your call to KSPSetOperators().
1106: More Notes about Repeated Solution of Linear Systems:
1107: PETSc does NOT reset the matrix entries of either Amat or Pmat
1108: to zero after a linear solve; the user is completely responsible for
1109: matrix assembly. See the routine MatZeroEntries() if desiring to
1110: zero all elements of a matrix.
1112: Level: intermediate
1114: .seealso: PCGetOperators(), MatZeroEntries()
1115: @*/
1116: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1117: {
1118: PetscInt m1,n1,m2,n2;
1125: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1126: MatGetLocalSize(Amat,&m1,&n1);
1127: MatGetLocalSize(pc->mat,&m2,&n2);
1129: MatGetLocalSize(Pmat,&m1,&n1);
1130: MatGetLocalSize(pc->pmat,&m2,&n2);
1132: }
1134: if (Pmat != pc->pmat) {
1135: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1136: pc->matnonzerostate = -1;
1137: pc->matstate = -1;
1138: }
1140: /* reference first in case the matrices are the same */
1141: if (Amat) PetscObjectReference((PetscObject)Amat);
1142: MatDestroy(&pc->mat);
1143: if (Pmat) PetscObjectReference((PetscObject)Pmat);
1144: MatDestroy(&pc->pmat);
1145: pc->mat = Amat;
1146: pc->pmat = Pmat;
1147: return 0;
1148: }
1150: /*@
1151: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1153: Logically Collective on PC
1155: Input Parameters:
1156: + pc - the preconditioner context
1157: - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1159: Level: intermediate
1161: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1162: @*/
1163: PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1164: {
1167: pc->reusepreconditioner = flag;
1168: return 0;
1169: }
1171: /*@
1172: PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1174: Not Collective
1176: Input Parameter:
1177: . pc - the preconditioner context
1179: Output Parameter:
1180: . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1182: Level: intermediate
1184: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1185: @*/
1186: PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1187: {
1190: *flag = pc->reusepreconditioner;
1191: return 0;
1192: }
1194: /*@
1195: PCGetOperators - Gets the matrix associated with the linear system and
1196: possibly a different one associated with the preconditioner.
1198: Not collective, though parallel Mats are returned if the PC is parallel
1200: Input Parameter:
1201: . pc - the preconditioner context
1203: Output Parameters:
1204: + Amat - the matrix defining the linear system
1205: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1207: Level: intermediate
1209: Notes:
1210: Does not increase the reference count of the matrices, so you should not destroy them
1212: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1213: are created in PC and returned to the user. In this case, if both operators
1214: mat and pmat are requested, two DIFFERENT operators will be returned. If
1215: only one is requested both operators in the PC will be the same (i.e. as
1216: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1217: The user must set the sizes of the returned matrices and their type etc just
1218: as if the user created them with MatCreate(). For example,
1220: $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1221: $ set size, type, etc of Amat
1223: $ MatCreate(comm,&mat);
1224: $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1225: $ PetscObjectDereference((PetscObject)mat);
1226: $ set size, type, etc of Amat
1228: and
1230: $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1231: $ set size, type, etc of Amat and Pmat
1233: $ MatCreate(comm,&Amat);
1234: $ MatCreate(comm,&Pmat);
1235: $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1236: $ PetscObjectDereference((PetscObject)Amat);
1237: $ PetscObjectDereference((PetscObject)Pmat);
1238: $ set size, type, etc of Amat and Pmat
1240: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1241: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1242: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1243: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1244: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1245: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1246: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1247: it can be created for you?
1249: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1250: @*/
1251: PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1252: {
1254: if (Amat) {
1255: if (!pc->mat) {
1256: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1257: pc->mat = pc->pmat;
1258: PetscObjectReference((PetscObject)pc->mat);
1259: } else { /* both Amat and Pmat are empty */
1260: MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1261: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1262: pc->pmat = pc->mat;
1263: PetscObjectReference((PetscObject)pc->pmat);
1264: }
1265: }
1266: }
1267: *Amat = pc->mat;
1268: }
1269: if (Pmat) {
1270: if (!pc->pmat) {
1271: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1272: pc->pmat = pc->mat;
1273: PetscObjectReference((PetscObject)pc->pmat);
1274: } else {
1275: MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1276: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1277: pc->mat = pc->pmat;
1278: PetscObjectReference((PetscObject)pc->mat);
1279: }
1280: }
1281: }
1282: *Pmat = pc->pmat;
1283: }
1284: return 0;
1285: }
1287: /*@C
1288: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1289: possibly a different one associated with the preconditioner have been set in the PC.
1291: Not collective, though the results on all processes should be the same
1293: Input Parameter:
1294: . pc - the preconditioner context
1296: Output Parameters:
1297: + mat - the matrix associated with the linear system was set
1298: - pmat - matrix associated with the preconditioner was set, usually the same
1300: Level: intermediate
1302: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1303: @*/
1304: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1305: {
1307: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1308: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1309: return 0;
1310: }
1312: /*@
1313: PCFactorGetMatrix - Gets the factored matrix from the
1314: preconditioner context. This routine is valid only for the LU,
1315: incomplete LU, Cholesky, and incomplete Cholesky methods.
1317: Not Collective on PC though Mat is parallel if PC is parallel
1319: Input Parameters:
1320: . pc - the preconditioner context
1322: Output parameters:
1323: . mat - the factored matrix
1325: Level: advanced
1327: Notes:
1328: Does not increase the reference count for the matrix so DO NOT destroy it
1330: @*/
1331: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1332: {
1335: if (pc->ops->getfactoredmatrix) {
1336: (*pc->ops->getfactoredmatrix)(pc,mat);
1337: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1338: return 0;
1339: }
1341: /*@C
1342: PCSetOptionsPrefix - Sets the prefix used for searching for all
1343: PC options in the database.
1345: Logically Collective on PC
1347: Input Parameters:
1348: + pc - the preconditioner context
1349: - prefix - the prefix string to prepend to all PC option requests
1351: Notes:
1352: A hyphen (-) must NOT be given at the beginning of the prefix name.
1353: The first character of all runtime options is AUTOMATICALLY the
1354: hyphen.
1356: Level: advanced
1358: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1359: @*/
1360: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1361: {
1363: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1364: return 0;
1365: }
1367: /*@C
1368: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1369: PC options in the database.
1371: Logically Collective on PC
1373: Input Parameters:
1374: + pc - the preconditioner context
1375: - prefix - the prefix string to prepend to all PC option requests
1377: Notes:
1378: A hyphen (-) must NOT be given at the beginning of the prefix name.
1379: The first character of all runtime options is AUTOMATICALLY the
1380: hyphen.
1382: Level: advanced
1384: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1385: @*/
1386: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1387: {
1389: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1390: return 0;
1391: }
1393: /*@C
1394: PCGetOptionsPrefix - Gets the prefix used for searching for all
1395: PC options in the database.
1397: Not Collective
1399: Input Parameters:
1400: . pc - the preconditioner context
1402: Output Parameters:
1403: . prefix - pointer to the prefix string used, is returned
1405: Notes:
1406: On the fortran side, the user should pass in a string 'prifix' of
1407: sufficient length to hold the prefix.
1409: Level: advanced
1411: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1412: @*/
1413: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1414: {
1417: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1418: return 0;
1419: }
1421: /*
1422: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1423: preconditioners including BDDC and Eisentat that transform the equations before applying
1424: the Krylov methods
1425: */
1426: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc,PetscBool *change)
1427: {
1430: *change = PETSC_FALSE;
1431: PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1432: return 0;
1433: }
1435: /*@
1436: PCPreSolve - Optional pre-solve phase, intended for any
1437: preconditioner-specific actions that must be performed before
1438: the iterative solve itself.
1440: Collective on PC
1442: Input Parameters:
1443: + pc - the preconditioner context
1444: - ksp - the Krylov subspace context
1446: Level: developer
1448: Sample of Usage:
1449: .vb
1450: PCPreSolve(pc,ksp);
1451: KSPSolve(ksp,b,x);
1452: PCPostSolve(pc,ksp);
1453: .ve
1455: Notes:
1456: The pre-solve phase is distinct from the PCSetUp() phase.
1458: KSPSolve() calls this directly, so is rarely called by the user.
1460: .seealso: PCPostSolve()
1461: @*/
1462: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1463: {
1464: Vec x,rhs;
1468: pc->presolvedone++;
1470: KSPGetSolution(ksp,&x);
1471: KSPGetRhs(ksp,&rhs);
1473: if (pc->ops->presolve) {
1474: (*pc->ops->presolve)(pc,ksp,rhs,x);
1475: } else if (pc->presolve) {
1476: (pc->presolve)(pc,ksp);
1477: }
1478: return 0;
1479: }
1481: /*@C
1482: PCSetPreSolve - Sets function PCPreSolve() which is intended for any
1483: preconditioner-specific actions that must be performed before
1484: the iterative solve itself.
1486: Logically Collective on pc
1488: Input Parameters:
1489: + pc - the preconditioner object
1490: - presolve - the function to call before the solve
1492: Calling sequence of presolve:
1493: $ func(PC pc,KSP ksp)
1495: + pc - the PC context
1496: - ksp - the KSP context
1498: Level: developer
1500: .seealso: PC, PCSetUp(), PCPreSolve()
1501: @*/
1502: PetscErrorCode PCSetPreSolve(PC pc,PetscErrorCode (*presolve)(PC,KSP))
1503: {
1505: pc->presolve = presolve;
1506: return 0;
1507: }
1509: /*@
1510: PCPostSolve - Optional post-solve phase, intended for any
1511: preconditioner-specific actions that must be performed after
1512: the iterative solve itself.
1514: Collective on PC
1516: Input Parameters:
1517: + pc - the preconditioner context
1518: - ksp - the Krylov subspace context
1520: Sample of Usage:
1521: .vb
1522: PCPreSolve(pc,ksp);
1523: KSPSolve(ksp,b,x);
1524: PCPostSolve(pc,ksp);
1525: .ve
1527: Note:
1528: KSPSolve() calls this routine directly, so it is rarely called by the user.
1530: Level: developer
1532: .seealso: PCPreSolve(), KSPSolve()
1533: @*/
1534: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1535: {
1536: Vec x,rhs;
1540: pc->presolvedone--;
1541: KSPGetSolution(ksp,&x);
1542: KSPGetRhs(ksp,&rhs);
1543: if (pc->ops->postsolve) {
1544: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1545: }
1546: return 0;
1547: }
1549: /*@C
1550: PCLoad - Loads a PC that has been stored in binary with PCView().
1552: Collective on PetscViewer
1554: Input Parameters:
1555: + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1556: some related function before a call to PCLoad().
1557: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1559: Level: intermediate
1561: Notes:
1562: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1564: Notes for advanced users:
1565: Most users should not need to know the details of the binary storage
1566: format, since PCLoad() and PCView() completely hide these details.
1567: But for anyone who's interested, the standard binary matrix storage
1568: format is
1569: .vb
1570: has not yet been determined
1571: .ve
1573: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1574: @*/
1575: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1576: {
1577: PetscBool isbinary;
1578: PetscInt classid;
1579: char type[256];
1583: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1586: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1588: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1589: PCSetType(newdm, type);
1590: if (newdm->ops->load) {
1591: (*newdm->ops->load)(newdm,viewer);
1592: }
1593: return 0;
1594: }
1596: #include <petscdraw.h>
1597: #if defined(PETSC_HAVE_SAWS)
1598: #include <petscviewersaws.h>
1599: #endif
1601: /*@C
1602: PCViewFromOptions - View from Options
1604: Collective on PC
1606: Input Parameters:
1607: + A - the PC context
1608: . obj - Optional object
1609: - name - command line option
1611: Level: intermediate
1612: .seealso: PC, PCView, PetscObjectViewFromOptions(), PCCreate()
1613: @*/
1614: PetscErrorCode PCViewFromOptions(PC A,PetscObject obj,const char name[])
1615: {
1617: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1618: return 0;
1619: }
1621: /*@C
1622: PCView - Prints the PC data structure.
1624: Collective on PC
1626: Input Parameters:
1627: + PC - the PC context
1628: - viewer - optional visualization context
1630: Note:
1631: The available visualization contexts include
1632: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1633: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1634: output where only the first processor opens
1635: the file. All other processors send their
1636: data to the first processor to print.
1638: The user can open an alternative visualization contexts with
1639: PetscViewerASCIIOpen() (output to a specified file).
1641: Level: developer
1643: .seealso: KSPView(), PetscViewerASCIIOpen()
1644: @*/
1645: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1646: {
1647: PCType cstr;
1648: PetscBool iascii,isstring,isbinary,isdraw;
1649: #if defined(PETSC_HAVE_SAWS)
1650: PetscBool issaws;
1651: #endif
1654: if (!viewer) {
1655: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1656: }
1660: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1661: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1662: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1663: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1664: #if defined(PETSC_HAVE_SAWS)
1665: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1666: #endif
1668: if (iascii) {
1669: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1670: if (!pc->setupcalled) {
1671: PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");
1672: }
1673: if (pc->ops->view) {
1674: PetscViewerASCIIPushTab(viewer);
1675: (*pc->ops->view)(pc,viewer);
1676: PetscViewerASCIIPopTab(viewer);
1677: }
1678: if (pc->mat) {
1679: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1680: if (pc->pmat == pc->mat) {
1681: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1682: PetscViewerASCIIPushTab(viewer);
1683: MatView(pc->mat,viewer);
1684: PetscViewerASCIIPopTab(viewer);
1685: } else {
1686: if (pc->pmat) {
1687: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1688: } else {
1689: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1690: }
1691: PetscViewerASCIIPushTab(viewer);
1692: MatView(pc->mat,viewer);
1693: if (pc->pmat) MatView(pc->pmat,viewer);
1694: PetscViewerASCIIPopTab(viewer);
1695: }
1696: PetscViewerPopFormat(viewer);
1697: }
1698: } else if (isstring) {
1699: PCGetType(pc,&cstr);
1700: PetscViewerStringSPrintf(viewer," PCType: %-7.7s",cstr);
1701: if (pc->ops->view) (*pc->ops->view)(pc,viewer);
1702: if (pc->mat) MatView(pc->mat,viewer);
1703: if (pc->pmat && pc->pmat != pc->mat) MatView(pc->pmat,viewer);
1704: } else if (isbinary) {
1705: PetscInt classid = PC_FILE_CLASSID;
1706: MPI_Comm comm;
1707: PetscMPIInt rank;
1708: char type[256];
1710: PetscObjectGetComm((PetscObject)pc,&comm);
1711: MPI_Comm_rank(comm,&rank);
1712: if (rank == 0) {
1713: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);
1714: PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1715: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);
1716: }
1717: if (pc->ops->view) {
1718: (*pc->ops->view)(pc,viewer);
1719: }
1720: } else if (isdraw) {
1721: PetscDraw draw;
1722: char str[25];
1723: PetscReal x,y,bottom,h;
1724: PetscInt n;
1726: PetscViewerDrawGetDraw(viewer,0,&draw);
1727: PetscDrawGetCurrentPoint(draw,&x,&y);
1728: if (pc->mat) {
1729: MatGetSize(pc->mat,&n,NULL);
1730: PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1731: } else {
1732: PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1733: }
1734: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1735: bottom = y - h;
1736: PetscDrawPushCurrentPoint(draw,x,bottom);
1737: if (pc->ops->view) {
1738: (*pc->ops->view)(pc,viewer);
1739: }
1740: PetscDrawPopCurrentPoint(draw);
1741: #if defined(PETSC_HAVE_SAWS)
1742: } else if (issaws) {
1743: PetscMPIInt rank;
1745: PetscObjectName((PetscObject)pc);
1746: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1747: if (!((PetscObject)pc)->amsmem && rank == 0) {
1748: PetscObjectViewSAWs((PetscObject)pc,viewer);
1749: }
1750: if (pc->mat) MatView(pc->mat,viewer);
1751: if (pc->pmat && pc->pmat != pc->mat) MatView(pc->pmat,viewer);
1752: #endif
1753: }
1754: return 0;
1755: }
1757: /*@C
1758: PCRegister - Adds a method to the preconditioner package.
1760: Not collective
1762: Input Parameters:
1763: + name_solver - name of a new user-defined solver
1764: - routine_create - routine to create method context
1766: Notes:
1767: PCRegister() may be called multiple times to add several user-defined preconditioners.
1769: Sample usage:
1770: .vb
1771: PCRegister("my_solver", MySolverCreate);
1772: .ve
1774: Then, your solver can be chosen with the procedural interface via
1775: $ PCSetType(pc,"my_solver")
1776: or at runtime via the option
1777: $ -pc_type my_solver
1779: Level: advanced
1781: .seealso: PCRegisterAll()
1782: @*/
1783: PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1784: {
1785: PCInitializePackage();
1786: PetscFunctionListAdd(&PCList,sname,function);
1787: return 0;
1788: }
1790: static PetscErrorCode MatMult_PC(Mat A,Vec X,Vec Y)
1791: {
1792: PC pc;
1794: MatShellGetContext(A,&pc);
1795: PCApply(pc,X,Y);
1796: return 0;
1797: }
1799: /*@
1800: PCComputeOperator - Computes the explicit preconditioned operator.
1802: Collective on PC
1804: Input Parameters:
1805: + pc - the preconditioner object
1806: - mattype - the matrix type to be used for the operator
1808: Output Parameter:
1809: . mat - the explicit preconditioned operator
1811: Notes:
1812: This computation is done by applying the operators to columns of the identity matrix.
1813: This routine is costly in general, and is recommended for use only with relatively small systems.
1814: Currently, this routine uses a dense matrix format when mattype == NULL
1816: Level: advanced
1818: .seealso: KSPComputeOperator(), MatType
1820: @*/
1821: PetscErrorCode PCComputeOperator(PC pc,MatType mattype,Mat *mat)
1822: {
1823: PetscInt N,M,m,n;
1824: Mat A,Apc;
1828: PCGetOperators(pc,&A,NULL);
1829: MatGetLocalSize(A,&m,&n);
1830: MatGetSize(A,&M,&N);
1831: MatCreateShell(PetscObjectComm((PetscObject)pc),m,n,M,N,pc,&Apc);
1832: MatShellSetOperation(Apc,MATOP_MULT,(void (*)(void))MatMult_PC);
1833: MatComputeOperator(Apc,mattype,mat);
1834: MatDestroy(&Apc);
1835: return 0;
1836: }
1838: /*@
1839: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1841: Collective on PC
1843: Input Parameters:
1844: + pc - the solver context
1845: . dim - the dimension of the coordinates 1, 2, or 3
1846: . nloc - the blocked size of the coordinates array
1847: - coords - the coordinates array
1849: Level: intermediate
1851: Notes:
1852: coords is an array of the dim coordinates for the nodes on
1853: the local processor, of size dim*nloc.
1854: If there are 108 equation on a processor
1855: for a displacement finite element discretization of elasticity (so
1856: that there are nloc = 36 = 108/3 nodes) then the array must have 108
1857: double precision values (ie, 3 * 36). These x y z coordinates
1858: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1859: ... , N-1.z ].
1861: .seealso: MatSetNearNullSpace()
1862: @*/
1863: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1864: {
1867: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1868: return 0;
1869: }
1871: /*@
1872: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
1874: Logically Collective on PC
1876: Input Parameter:
1877: . pc - the precondition context
1879: Output Parameters:
1880: + num_levels - the number of levels
1881: - interpolations - the interpolation matrices (size of num_levels-1)
1883: Level: advanced
1885: .keywords: MG, GAMG, BoomerAMG, multigrid, interpolation, level
1887: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetInterpolation(), PCGetCoarseOperators()
1888: @*/
1889: PetscErrorCode PCGetInterpolations(PC pc,PetscInt *num_levels,Mat *interpolations[])
1890: {
1894: PetscUseMethod(pc,"PCGetInterpolations_C",(PC,PetscInt*,Mat*[]),(pc,num_levels,interpolations));
1895: return 0;
1896: }
1898: /*@
1899: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
1901: Logically Collective on PC
1903: Input Parameter:
1904: . pc - the precondition context
1906: Output Parameters:
1907: + num_levels - the number of levels
1908: - coarseOperators - the coarse operator matrices (size of num_levels-1)
1910: Level: advanced
1912: .keywords: MG, GAMG, BoomerAMG, get, multigrid, interpolation, level
1914: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale(), PCMGGetInterpolation(), PCGetInterpolations()
1915: @*/
1916: PetscErrorCode PCGetCoarseOperators(PC pc,PetscInt *num_levels,Mat *coarseOperators[])
1917: {
1921: PetscUseMethod(pc,"PCGetCoarseOperators_C",(PC,PetscInt*,Mat*[]),(pc,num_levels,coarseOperators));
1922: return 0;
1923: }