Actual source code: precon.c
petsc-3.10.5 2019-03-28
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_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: {
17: PetscMPIInt size;
18: PetscBool hasop,flg1,flg2,set,flg3;
21: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
22: if (pc->pmat) {
23: MatHasOperation(pc->pmat,MATOP_GET_DIAGONAL_BLOCK,&hasop);
24: if (size == 1) {
25: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
26: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
27: MatIsSymmetricKnown(pc->pmat,&set,&flg3);
28: if (flg1 && (!flg2 || (set && flg3))) {
29: *type = PCICC;
30: } else if (flg2) {
31: *type = PCILU;
32: } else if (hasop) { /* likely is a parallel matrix run on one processor */
33: *type = PCBJACOBI;
34: } else {
35: *type = PCNONE;
36: }
37: } else {
38: if (hasop) {
39: *type = PCBJACOBI;
40: } else {
41: *type = PCNONE;
42: }
43: }
44: } else {
45: if (size == 1) {
46: *type = PCILU;
47: } else {
48: *type = PCBJACOBI;
49: }
50: }
51: return(0);
52: }
54: /*@
55: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
57: Collective on PC
59: Input Parameter:
60: . pc - the preconditioner context
62: Level: developer
64: Notes:
65: 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
67: .keywords: PC, destroy
69: .seealso: PCCreate(), PCSetUp()
70: @*/
71: PetscErrorCode PCReset(PC pc)
72: {
77: if (pc->ops->reset) {
78: (*pc->ops->reset)(pc);
79: }
80: VecDestroy(&pc->diagonalscaleright);
81: VecDestroy(&pc->diagonalscaleleft);
82: MatDestroy(&pc->pmat);
83: MatDestroy(&pc->mat);
85: pc->setupcalled = 0;
86: return(0);
87: }
89: /*@
90: PCDestroy - Destroys PC context that was created with PCCreate().
92: Collective on PC
94: Input Parameter:
95: . pc - the preconditioner context
97: Level: developer
99: .keywords: PC, destroy
101: .seealso: PCCreate(), PCSetUp()
102: @*/
103: PetscErrorCode PCDestroy(PC *pc)
104: {
108: if (!*pc) return(0);
110: if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; return(0);}
112: PCReset(*pc);
114: /* if memory was published with SAWs then destroy it */
115: PetscObjectSAWsViewOff((PetscObject)*pc);
116: if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
117: DMDestroy(&(*pc)->dm);
118: PetscHeaderDestroy(pc);
119: return(0);
120: }
122: /*@C
123: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
124: scaling as needed by certain time-stepping codes.
126: Logically Collective on PC
128: Input Parameter:
129: . pc - the preconditioner context
131: Output Parameter:
132: . flag - PETSC_TRUE if it applies the scaling
134: Level: developer
136: Notes:
137: If this returns PETSC_TRUE then the system solved via the Krylov method is
138: $ D M A D^{-1} y = D M b for left preconditioning or
139: $ D A M D^{-1} z = D b for right preconditioning
141: .keywords: PC
143: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
144: @*/
145: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
146: {
150: *flag = pc->diagonalscale;
151: return(0);
152: }
154: /*@
155: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
156: scaling as needed by certain time-stepping codes.
158: Logically Collective on PC
160: Input Parameters:
161: + pc - the preconditioner context
162: - s - scaling vector
164: Level: intermediate
166: Notes:
167: The system solved via the Krylov method is
168: $ D M A D^{-1} y = D M b for left preconditioning or
169: $ D A M D^{-1} z = D b for right preconditioning
171: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
173: .keywords: PC
175: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
176: @*/
177: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
178: {
184: pc->diagonalscale = PETSC_TRUE;
186: PetscObjectReference((PetscObject)s);
187: VecDestroy(&pc->diagonalscaleleft);
189: pc->diagonalscaleleft = s;
191: VecDuplicate(s,&pc->diagonalscaleright);
192: VecCopy(s,pc->diagonalscaleright);
193: VecReciprocal(pc->diagonalscaleright);
194: return(0);
195: }
197: /*@
198: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
200: Logically Collective on PC
202: Input Parameters:
203: + pc - the preconditioner context
204: . in - input vector
205: + out - scaled vector (maybe the same as in)
207: Level: intermediate
209: Notes:
210: The system solved via the Krylov method is
211: $ D M A D^{-1} y = D M b for left preconditioning or
212: $ D A M D^{-1} z = D b for right preconditioning
214: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
216: If diagonal scaling is turned off and in is not out then in is copied to out
218: .keywords: PC
220: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
221: @*/
222: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
223: {
230: if (pc->diagonalscale) {
231: VecPointwiseMult(out,pc->diagonalscaleleft,in);
232: } else if (in != out) {
233: VecCopy(in,out);
234: }
235: return(0);
236: }
238: /*@
239: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
241: Logically Collective on PC
243: Input Parameters:
244: + pc - the preconditioner context
245: . in - input vector
246: + out - scaled vector (maybe the same as in)
248: Level: intermediate
250: Notes:
251: The system solved via the Krylov method is
252: $ D M A D^{-1} y = D M b for left preconditioning or
253: $ D A M D^{-1} z = D b for right preconditioning
255: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
257: If diagonal scaling is turned off and in is not out then in is copied to out
259: .keywords: PC
261: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
262: @*/
263: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
264: {
271: if (pc->diagonalscale) {
272: VecPointwiseMult(out,pc->diagonalscaleright,in);
273: } else if (in != out) {
274: VecCopy(in,out);
275: }
276: return(0);
277: }
279: /*@
280: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
281: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
282: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
284: Logically Collective on PC
286: Input Parameters:
287: + pc - the preconditioner context
288: - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
290: Options Database Key:
291: . -pc_use_amat <true,false>
293: Notes:
294: For the common case in which the linear system matrix and the matrix used to construct the
295: preconditioner are identical, this routine is does nothing.
297: Level: intermediate
299: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
300: @*/
301: PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
302: {
305: pc->useAmat = flg;
306: return(0);
307: }
309: /*@
310: PCSetErrorIfFailure - Causes PC to generate an error if a FPE, for example a zero pivot, is detected.
312: Logically Collective on PC
314: Input Parameters:
315: + pc - iterative context obtained from PCCreate()
316: - flg - PETSC_TRUE indicates you want the error generated
318: Level: advanced
320: Notes:
321: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
322: 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
323: detected.
325: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
327: .keywords: PC, set, initial guess, nonzero
329: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
330: @*/
331: PetscErrorCode PCSetErrorIfFailure(PC pc,PetscBool flg)
332: {
336: pc->erroriffailure = flg;
337: return(0);
338: }
340: /*@
341: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
342: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
343: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
345: Logically Collective on PC
347: Input Parameter:
348: . pc - the preconditioner context
350: Output Parameter:
351: . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
353: Notes:
354: For the common case in which the linear system matrix and the matrix used to construct the
355: preconditioner are identical, this routine is does nothing.
357: Level: intermediate
359: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
360: @*/
361: PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
362: {
365: *flg = pc->useAmat;
366: return(0);
367: }
369: /*@
370: PCCreate - Creates a preconditioner context.
372: Collective on MPI_Comm
374: Input Parameter:
375: . comm - MPI communicator
377: Output Parameter:
378: . pc - location to put the preconditioner context
380: Notes:
381: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or ICC
382: in parallel. For dense matrices it is always PCNONE.
384: Level: developer
386: .keywords: PC, create, context
388: .seealso: PCSetUp(), PCApply(), PCDestroy()
389: @*/
390: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
391: {
392: PC pc;
397: *newpc = 0;
398: PCInitializePackage();
400: PetscHeaderCreate(pc,PC_CLASSID,"PC","Preconditioner","PC",comm,PCDestroy,PCView);
402: pc->mat = 0;
403: pc->pmat = 0;
404: pc->setupcalled = 0;
405: pc->setfromoptionscalled = 0;
406: pc->data = 0;
407: pc->diagonalscale = PETSC_FALSE;
408: pc->diagonalscaleleft = 0;
409: pc->diagonalscaleright = 0;
411: pc->modifysubmatrices = 0;
412: pc->modifysubmatricesP = 0;
414: *newpc = pc;
415: return(0);
417: }
419: /* -------------------------------------------------------------------------------*/
421: /*@
422: PCApply - Applies the preconditioner to a vector.
424: Collective on PC and Vec
426: Input Parameters:
427: + pc - the preconditioner context
428: - x - input vector
430: Output Parameter:
431: . y - output vector
433: Level: developer
435: .keywords: PC, apply
437: .seealso: PCApplyTranspose(), PCApplyBAorAB()
438: @*/
439: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
440: {
442: PetscInt m,n,mv,nv;
448: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
449: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
450: /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
451: MatGetLocalSize(pc->pmat,&m,&n);
452: VecGetLocalSize(x,&nv);
453: VecGetLocalSize(y,&mv);
454: if (mv != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Preconditioner number of local rows %D does not equal resulting vector number of rows %D",m,mv);
455: if (nv != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Preconditioner number of local columns %D does not equal resulting vector number of rows %D",n,nv);
456: VecLocked(y,3);
458: PCSetUp(pc);
459: if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
460: VecLockPush(x);
461: PetscLogEventBegin(PC_Apply,pc,x,y,0);
462: (*pc->ops->apply)(pc,x,y);
463: PetscLogEventEnd(PC_Apply,pc,x,y,0);
464: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
465: VecLockPop(x);
466: return(0);
467: }
469: /*@
470: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
472: Collective on PC and Vec
474: Input Parameters:
475: + pc - the preconditioner context
476: - x - input vector
478: Output Parameter:
479: . y - output vector
481: Notes:
482: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
484: Level: developer
486: .keywords: PC, apply, symmetric, left
488: .seealso: PCApply(), PCApplySymmetricRight()
489: @*/
490: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
491: {
498: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
499: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
500: PCSetUp(pc);
501: if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
502: VecLockPush(x);
503: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
504: (*pc->ops->applysymmetricleft)(pc,x,y);
505: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
506: VecLockPop(x);
507: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
508: return(0);
509: }
511: /*@
512: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
514: Collective on PC and Vec
516: Input Parameters:
517: + pc - the preconditioner context
518: - x - input vector
520: Output Parameter:
521: . y - output vector
523: Level: developer
525: Notes:
526: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
528: .keywords: PC, apply, symmetric, right
530: .seealso: PCApply(), PCApplySymmetricLeft()
531: @*/
532: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
533: {
540: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
541: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
542: PCSetUp(pc);
543: if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
544: VecLockPush(x);
545: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
546: (*pc->ops->applysymmetricright)(pc,x,y);
547: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
548: VecLockPop(x);
549: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
550: return(0);
551: }
553: /*@
554: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
556: Collective on PC and Vec
558: Input Parameters:
559: + pc - the preconditioner context
560: - x - input vector
562: Output Parameter:
563: . y - output vector
565: Notes:
566: For complex numbers this applies the non-Hermitian transpose.
568: Developer Notes:
569: We need to implement a PCApplyHermitianTranspose()
571: Level: developer
573: .keywords: PC, apply, transpose
575: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
576: @*/
577: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
578: {
585: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
586: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
587: PCSetUp(pc);
588: if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
589: VecLockPush(x);
590: PetscLogEventBegin(PC_Apply,pc,x,y,0);
591: (*pc->ops->applytranspose)(pc,x,y);
592: PetscLogEventEnd(PC_Apply,pc,x,y,0);
593: VecLockPop(x);
594: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
595: return(0);
596: }
598: /*@
599: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
601: Collective on PC and Vec
603: Input Parameters:
604: . pc - the preconditioner context
606: Output Parameter:
607: . flg - PETSC_TRUE if a transpose operation is defined
609: Level: developer
611: .keywords: PC, apply, transpose
613: .seealso: PCApplyTranspose()
614: @*/
615: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
616: {
620: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
621: else *flg = PETSC_FALSE;
622: return(0);
623: }
625: /*@
626: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
628: Collective on PC and Vec
630: Input Parameters:
631: + pc - the preconditioner context
632: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
633: . x - input vector
634: - work - work vector
636: Output Parameter:
637: . y - output vector
639: Level: developer
641: Notes:
642: 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
643: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
645: .keywords: PC, apply, operator
647: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
648: @*/
649: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
650: {
658: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
659: if (side != PC_LEFT && side != PC_SYMMETRIC && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right, left, or symmetric");
660: if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
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);
682: } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
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 and Vec
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
721: Notes:
722: 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
723: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
725: Level: developer
727: .keywords: PC, apply, operator, transpose
729: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
730: @*/
731: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
732: {
740: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
741: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
742: if (pc->ops->applyBAtranspose) {
743: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
744: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
745: return(0);
746: }
747: if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
749: PCSetUp(pc);
750: if (side == PC_RIGHT) {
751: PCApplyTranspose(pc,x,work);
752: MatMultTranspose(pc->mat,work,y);
753: } else if (side == PC_LEFT) {
754: MatMultTranspose(pc->mat,x,work);
755: PCApplyTranspose(pc,work,y);
756: }
757: /* add support for PC_SYMMETRIC */
758: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
759: return(0);
760: }
762: /* -------------------------------------------------------------------------------*/
764: /*@
765: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
766: built-in fast application of Richardson's method.
768: Not Collective
770: Input Parameter:
771: . pc - the preconditioner
773: Output Parameter:
774: . exists - PETSC_TRUE or PETSC_FALSE
776: Level: developer
778: .keywords: PC, apply, Richardson, exists
780: .seealso: PCApplyRichardson()
781: @*/
782: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
783: {
787: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
788: else *exists = PETSC_FALSE;
789: return(0);
790: }
792: /*@
793: PCApplyRichardson - Applies several steps of Richardson iteration with
794: the particular preconditioner. This routine is usually used by the
795: Krylov solvers and not the application code directly.
797: Collective on PC
799: Input Parameters:
800: + pc - the preconditioner context
801: . b - the right hand side
802: . w - one work vector
803: . rtol - relative decrease in residual norm convergence criteria
804: . abstol - absolute residual norm convergence criteria
805: . dtol - divergence residual norm increase criteria
806: . its - the number of iterations to apply.
807: - guesszero - if the input x contains nonzero initial guess
809: Output Parameter:
810: + outits - number of iterations actually used (for SOR this always equals its)
811: . reason - the reason the apply terminated
812: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
814: Notes:
815: Most preconditioners do not support this function. Use the command
816: PCApplyRichardsonExists() to determine if one does.
818: Except for the multigrid PC this routine ignores the convergence tolerances
819: and always runs for the number of iterations
821: Level: developer
823: .keywords: PC, apply, Richardson
825: .seealso: PCApplyRichardsonExists()
826: @*/
827: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
828: {
836: if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
837: PCSetUp(pc);
838: if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
839: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
840: return(0);
841: }
843: /*@
844: PCGetSetUpFailedReason - Gets the reason a PCSetUp() failed or 0 if it did not fail
846: Logically Collective on PC
848: Input Parameter:
849: . pc - the preconditioner context
851: Output Parameter:
852: . reason - the reason it failed, currently only -1
854: Level: advanced
856: .keywords: PC, setup
858: .seealso: PCCreate(), PCApply(), PCDestroy()
859: @*/
860: PetscErrorCode PCGetSetUpFailedReason(PC pc,PCFailedReason *reason)
861: {
863: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
864: else *reason = pc->failedreason;
865: return(0);
866: }
869: /*
870: a setupcall of 0 indicates never setup,
871: 1 indicates has been previously setup
872: -1 indicates a PCSetUp() was attempted and failed
873: */
874: /*@
875: PCSetUp - Prepares for the use of a preconditioner.
877: Collective on PC
879: Input Parameter:
880: . pc - the preconditioner context
882: Level: developer
884: .keywords: PC, setup
886: .seealso: PCCreate(), PCApply(), PCDestroy()
887: @*/
888: PetscErrorCode PCSetUp(PC pc)
889: {
890: PetscErrorCode ierr;
891: const char *def;
892: PetscObjectState matstate, matnonzerostate;
896: if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
898: if (pc->setupcalled && pc->reusepreconditioner) {
899: PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
900: return(0);
901: }
903: PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
904: MatGetNonzeroState(pc->pmat,&matnonzerostate);
905: if (!pc->setupcalled) {
906: PetscInfo(pc,"Setting up PC for first time\n");
907: pc->flag = DIFFERENT_NONZERO_PATTERN;
908: } else if (matstate == pc->matstate) {
909: PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
910: return(0);
911: } else {
912: if (matnonzerostate > pc->matnonzerostate) {
913: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
914: pc->flag = DIFFERENT_NONZERO_PATTERN;
915: } else {
916: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
917: pc->flag = SAME_NONZERO_PATTERN;
918: }
919: }
920: pc->matstate = matstate;
921: pc->matnonzerostate = matnonzerostate;
923: if (!((PetscObject)pc)->type_name) {
924: PCGetDefaultType_Private(pc,&def);
925: PCSetType(pc,def);
926: }
928: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
929: MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
930: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
931: if (pc->ops->setup) {
932: (*pc->ops->setup)(pc);
933: }
934: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
935: if (!pc->setupcalled) pc->setupcalled = 1;
936: return(0);
937: }
939: /*@
940: PCSetUpOnBlocks - Sets up the preconditioner for each block in
941: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
942: methods.
944: Collective on PC
946: Input Parameters:
947: . pc - the preconditioner context
949: Level: developer
951: .keywords: PC, setup, blocks
953: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
954: @*/
955: PetscErrorCode PCSetUpOnBlocks(PC pc)
956: {
961: if (!pc->ops->setuponblocks) return(0);
962: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
963: (*pc->ops->setuponblocks)(pc);
964: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
965: return(0);
966: }
968: /*@C
969: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
970: submatrices that arise within certain subdomain-based preconditioners.
971: The basic submatrices are extracted from the preconditioner matrix as
972: usual; the user can then alter these (for example, to set different boundary
973: conditions for each submatrix) before they are used for the local solves.
975: Logically Collective on PC
977: Input Parameters:
978: + pc - the preconditioner context
979: . func - routine for modifying the submatrices
980: - ctx - optional user-defined context (may be null)
982: Calling sequence of func:
983: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
985: . row - an array of index sets that contain the global row numbers
986: that comprise each local submatrix
987: . col - an array of index sets that contain the global column numbers
988: that comprise each local submatrix
989: . submat - array of local submatrices
990: - ctx - optional user-defined context for private data for the
991: user-defined func routine (may be null)
993: Notes:
994: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
995: KSPSolve().
997: A routine set by PCSetModifySubMatrices() is currently called within
998: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
999: preconditioners. All other preconditioners ignore this routine.
1001: Level: advanced
1003: .keywords: PC, set, modify, submatrices
1005: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
1006: @*/
1007: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1008: {
1011: pc->modifysubmatrices = func;
1012: pc->modifysubmatricesP = ctx;
1013: return(0);
1014: }
1016: /*@C
1017: PCModifySubMatrices - Calls an optional user-defined routine within
1018: certain preconditioners if one has been set with PCSetModifySubMarices().
1020: Collective on PC
1022: Input Parameters:
1023: + pc - the preconditioner context
1024: . nsub - the number of local submatrices
1025: . row - an array of index sets that contain the global row numbers
1026: that comprise each local submatrix
1027: . col - an array of index sets that contain the global column numbers
1028: that comprise each local submatrix
1029: . submat - array of local submatrices
1030: - ctx - optional user-defined context for private data for the
1031: user-defined routine (may be null)
1033: Output Parameter:
1034: . submat - array of local submatrices (the entries of which may
1035: have been modified)
1037: Notes:
1038: The user should NOT generally call this routine, as it will
1039: automatically be called within certain preconditioners (currently
1040: block Jacobi, additive Schwarz) if set.
1042: The basic submatrices are extracted from the preconditioner matrix
1043: as usual; the user can then alter these (for example, to set different
1044: boundary conditions for each submatrix) before they are used for the
1045: local solves.
1047: Level: developer
1049: .keywords: PC, modify, submatrices
1051: .seealso: PCSetModifySubMatrices()
1052: @*/
1053: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1054: {
1059: if (!pc->modifysubmatrices) return(0);
1060: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1061: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1062: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1063: return(0);
1064: }
1066: /*@
1067: PCSetOperators - Sets the matrix associated with the linear system and
1068: a (possibly) different one associated with the preconditioner.
1070: Logically Collective on PC and Mat
1072: Input Parameters:
1073: + pc - the preconditioner context
1074: . Amat - the matrix that defines the linear system
1075: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1077: Notes:
1078: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1080: If you wish to replace either Amat or Pmat but leave the other one untouched then
1081: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1082: on it and then pass it back in in your call to KSPSetOperators().
1084: More Notes about Repeated Solution of Linear Systems:
1085: PETSc does NOT reset the matrix entries of either Amat or Pmat
1086: to zero after a linear solve; the user is completely responsible for
1087: matrix assembly. See the routine MatZeroEntries() if desiring to
1088: zero all elements of a matrix.
1090: Level: intermediate
1092: .keywords: PC, set, operators, matrix, linear system
1094: .seealso: PCGetOperators(), MatZeroEntries()
1095: @*/
1096: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1097: {
1098: PetscErrorCode ierr;
1099: PetscInt m1,n1,m2,n2;
1107: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1108: MatGetLocalSize(Amat,&m1,&n1);
1109: MatGetLocalSize(pc->mat,&m2,&n2);
1110: if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot change local size of Amat after use old sizes %D %D new sizes %D %D",m2,n2,m1,n1);
1111: MatGetLocalSize(Pmat,&m1,&n1);
1112: MatGetLocalSize(pc->pmat,&m2,&n2);
1113: if (m1 != m2 || n1 != n2) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot change local size of Pmat after use old sizes %D %D new sizes %D %D",m2,n2,m1,n1);
1114: }
1116: if (Pmat != pc->pmat) {
1117: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1118: pc->matnonzerostate = -1;
1119: pc->matstate = -1;
1120: }
1122: /* reference first in case the matrices are the same */
1123: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1124: MatDestroy(&pc->mat);
1125: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1126: MatDestroy(&pc->pmat);
1127: pc->mat = Amat;
1128: pc->pmat = Pmat;
1129: return(0);
1130: }
1132: /*@
1133: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1135: Logically Collective on PC
1137: Input Parameters:
1138: + pc - the preconditioner context
1139: - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1141: Level: intermediate
1143: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1144: @*/
1145: PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1146: {
1149: pc->reusepreconditioner = flag;
1150: return(0);
1151: }
1153: /*@
1154: PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1156: Not Collective
1158: Input Parameter:
1159: . pc - the preconditioner context
1161: Output Parameter:
1162: . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1164: Level: intermediate
1166: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1167: @*/
1168: PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1169: {
1172: *flag = pc->reusepreconditioner;
1173: return(0);
1174: }
1176: /*@
1177: PCGetOperators - Gets the matrix associated with the linear system and
1178: possibly a different one associated with the preconditioner.
1180: Not collective, though parallel Mats are returned if the PC is parallel
1182: Input Parameter:
1183: . pc - the preconditioner context
1185: Output Parameters:
1186: + Amat - the matrix defining the linear system
1187: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1189: Level: intermediate
1191: Notes:
1192: Does not increase the reference count of the matrices, so you should not destroy them
1194: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1195: are created in PC and returned to the user. In this case, if both operators
1196: mat and pmat are requested, two DIFFERENT operators will be returned. If
1197: only one is requested both operators in the PC will be the same (i.e. as
1198: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1199: The user must set the sizes of the returned matrices and their type etc just
1200: as if the user created them with MatCreate(). For example,
1202: $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1203: $ set size, type, etc of Amat
1205: $ MatCreate(comm,&mat);
1206: $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1207: $ PetscObjectDereference((PetscObject)mat);
1208: $ set size, type, etc of Amat
1210: and
1212: $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1213: $ set size, type, etc of Amat and Pmat
1215: $ MatCreate(comm,&Amat);
1216: $ MatCreate(comm,&Pmat);
1217: $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1218: $ PetscObjectDereference((PetscObject)Amat);
1219: $ PetscObjectDereference((PetscObject)Pmat);
1220: $ set size, type, etc of Amat and Pmat
1222: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1223: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1224: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1225: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1226: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1227: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1228: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1229: it can be created for you?
1232: .keywords: PC, get, operators, matrix, linear system
1234: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1235: @*/
1236: PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1237: {
1242: if (Amat) {
1243: if (!pc->mat) {
1244: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1245: pc->mat = pc->pmat;
1246: PetscObjectReference((PetscObject)pc->mat);
1247: } else { /* both Amat and Pmat are empty */
1248: MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1249: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1250: pc->pmat = pc->mat;
1251: PetscObjectReference((PetscObject)pc->pmat);
1252: }
1253: }
1254: }
1255: *Amat = pc->mat;
1256: }
1257: if (Pmat) {
1258: if (!pc->pmat) {
1259: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1260: pc->pmat = pc->mat;
1261: PetscObjectReference((PetscObject)pc->pmat);
1262: } else {
1263: MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1264: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1265: pc->mat = pc->pmat;
1266: PetscObjectReference((PetscObject)pc->mat);
1267: }
1268: }
1269: }
1270: *Pmat = pc->pmat;
1271: }
1272: return(0);
1273: }
1275: /*@C
1276: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1277: possibly a different one associated with the preconditioner have been set in the PC.
1279: Not collective, though the results on all processes should be the same
1281: Input Parameter:
1282: . pc - the preconditioner context
1284: Output Parameters:
1285: + mat - the matrix associated with the linear system was set
1286: - pmat - matrix associated with the preconditioner was set, usually the same
1288: Level: intermediate
1290: .keywords: PC, get, operators, matrix, linear system
1292: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1293: @*/
1294: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1295: {
1298: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1299: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1300: return(0);
1301: }
1303: /*@
1304: PCFactorGetMatrix - Gets the factored matrix from the
1305: preconditioner context. This routine is valid only for the LU,
1306: incomplete LU, Cholesky, and incomplete Cholesky methods.
1308: Not Collective on PC though Mat is parallel if PC is parallel
1310: Input Parameters:
1311: . pc - the preconditioner context
1313: Output parameters:
1314: . mat - the factored matrix
1316: Level: advanced
1318: Notes:
1319: Does not increase the reference count for the matrix so DO NOT destroy it
1321: .keywords: PC, get, factored, matrix
1322: @*/
1323: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1324: {
1330: if (pc->ops->getfactoredmatrix) {
1331: (*pc->ops->getfactoredmatrix)(pc,mat);
1332: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1333: return(0);
1334: }
1336: /*@C
1337: PCSetOptionsPrefix - Sets the prefix used for searching for all
1338: PC options in the database.
1340: Logically Collective on PC
1342: Input Parameters:
1343: + pc - the preconditioner context
1344: - prefix - the prefix string to prepend to all PC option requests
1346: Notes:
1347: A hyphen (-) must NOT be given at the beginning of the prefix name.
1348: The first character of all runtime options is AUTOMATICALLY the
1349: hyphen.
1351: Level: advanced
1353: .keywords: PC, set, options, prefix, database
1355: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1356: @*/
1357: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1358: {
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: .keywords: PC, append, options, prefix, database
1386: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1387: @*/
1388: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1389: {
1394: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1395: return(0);
1396: }
1398: /*@C
1399: PCGetOptionsPrefix - Gets the prefix used for searching for all
1400: PC options in the database.
1402: Not Collective
1404: Input Parameters:
1405: . pc - the preconditioner context
1407: Output Parameters:
1408: . prefix - pointer to the prefix string used, is returned
1410: Notes:
1411: On the fortran side, the user should pass in a string 'prifix' of
1412: sufficient length to hold the prefix.
1414: Level: advanced
1416: .keywords: PC, get, options, prefix, database
1418: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1419: @*/
1420: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1421: {
1427: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1428: return(0);
1429: }
1431: /*
1432: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1433: preconditioners including BDDC and Eisentat that transform the equations before applying
1434: the Krylov methods
1435: */
1436: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc,PetscBool *change)
1437: {
1443: *change = PETSC_FALSE;
1444: PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1445: return(0);
1446: }
1448: /*@
1449: PCPreSolve - Optional pre-solve phase, intended for any
1450: preconditioner-specific actions that must be performed before
1451: the iterative solve itself.
1453: Collective on PC
1455: Input Parameters:
1456: + pc - the preconditioner context
1457: - ksp - the Krylov subspace context
1459: Level: developer
1461: Sample of Usage:
1462: .vb
1463: PCPreSolve(pc,ksp);
1464: KSPSolve(ksp,b,x);
1465: PCPostSolve(pc,ksp);
1466: .ve
1468: Notes:
1469: The pre-solve phase is distinct from the PCSetUp() phase.
1471: KSPSolve() calls this directly, so is rarely called by the user.
1473: .keywords: PC, pre-solve
1475: .seealso: PCPostSolve()
1476: @*/
1477: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1478: {
1480: Vec x,rhs;
1485: pc->presolvedone++;
1486: if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1487: KSPGetSolution(ksp,&x);
1488: KSPGetRhs(ksp,&rhs);
1490: if (pc->ops->presolve) {
1491: (*pc->ops->presolve)(pc,ksp,rhs,x);
1492: }
1493: return(0);
1494: }
1496: /*@
1497: PCPostSolve - Optional post-solve phase, intended for any
1498: preconditioner-specific actions that must be performed after
1499: the iterative solve itself.
1501: Collective on PC
1503: Input Parameters:
1504: + pc - the preconditioner context
1505: - ksp - the Krylov subspace context
1507: Sample of Usage:
1508: .vb
1509: PCPreSolve(pc,ksp);
1510: KSPSolve(ksp,b,x);
1511: PCPostSolve(pc,ksp);
1512: .ve
1514: Note:
1515: KSPSolve() calls this routine directly, so it is rarely called by the user.
1517: Level: developer
1519: .keywords: PC, post-solve
1521: .seealso: PCPreSolve(), KSPSolve()
1522: @*/
1523: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1524: {
1526: Vec x,rhs;
1531: pc->presolvedone--;
1532: KSPGetSolution(ksp,&x);
1533: KSPGetRhs(ksp,&rhs);
1534: if (pc->ops->postsolve) {
1535: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1536: }
1537: return(0);
1538: }
1540: /*@C
1541: PCLoad - Loads a PC that has been stored in binary with PCView().
1543: Collective on PetscViewer
1545: Input Parameters:
1546: + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1547: some related function before a call to PCLoad().
1548: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1550: Level: intermediate
1552: Notes:
1553: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1555: Notes for advanced users:
1556: Most users should not need to know the details of the binary storage
1557: format, since PCLoad() and PCView() completely hide these details.
1558: But for anyone who's interested, the standard binary matrix storage
1559: format is
1560: .vb
1561: has not yet been determined
1562: .ve
1564: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1565: @*/
1566: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1567: {
1569: PetscBool isbinary;
1570: PetscInt classid;
1571: char type[256];
1576: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1577: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1579: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1580: if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1581: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1582: PCSetType(newdm, type);
1583: if (newdm->ops->load) {
1584: (*newdm->ops->load)(newdm,viewer);
1585: }
1586: return(0);
1587: }
1589: #include <petscdraw.h>
1590: #if defined(PETSC_HAVE_SAWS)
1591: #include <petscviewersaws.h>
1592: #endif
1593: /*@C
1594: PCView - Prints the PC data structure.
1596: Collective on PC
1598: Input Parameters:
1599: + PC - the PC context
1600: - viewer - optional visualization context
1602: Note:
1603: The available visualization contexts include
1604: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1605: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1606: output where only the first processor opens
1607: the file. All other processors send their
1608: data to the first processor to print.
1610: The user can open an alternative visualization contexts with
1611: PetscViewerASCIIOpen() (output to a specified file).
1613: Level: developer
1615: .keywords: PC, view
1617: .seealso: KSPView(), PetscViewerASCIIOpen()
1618: @*/
1619: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1620: {
1621: PCType cstr;
1623: PetscBool iascii,isstring,isbinary,isdraw;
1624: #if defined(PETSC_HAVE_SAWS)
1625: PetscBool issaws;
1626: #endif
1630: if (!viewer) {
1631: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1632: }
1636: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1637: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1638: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1639: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1640: #if defined(PETSC_HAVE_SAWS)
1641: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1642: #endif
1644: if (iascii) {
1645: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1646: if (!pc->setupcalled) {
1647: PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");
1648: }
1649: if (pc->ops->view) {
1650: PetscViewerASCIIPushTab(viewer);
1651: (*pc->ops->view)(pc,viewer);
1652: PetscViewerASCIIPopTab(viewer);
1653: }
1654: if (pc->mat) {
1655: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1656: if (pc->pmat == pc->mat) {
1657: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1658: PetscViewerASCIIPushTab(viewer);
1659: MatView(pc->mat,viewer);
1660: PetscViewerASCIIPopTab(viewer);
1661: } else {
1662: if (pc->pmat) {
1663: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1664: } else {
1665: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1666: }
1667: PetscViewerASCIIPushTab(viewer);
1668: MatView(pc->mat,viewer);
1669: if (pc->pmat) {MatView(pc->pmat,viewer);}
1670: PetscViewerASCIIPopTab(viewer);
1671: }
1672: PetscViewerPopFormat(viewer);
1673: }
1674: } else if (isstring) {
1675: PCGetType(pc,&cstr);
1676: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1677: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1678: } else if (isbinary) {
1679: PetscInt classid = PC_FILE_CLASSID;
1680: MPI_Comm comm;
1681: PetscMPIInt rank;
1682: char type[256];
1684: PetscObjectGetComm((PetscObject)pc,&comm);
1685: MPI_Comm_rank(comm,&rank);
1686: if (!rank) {
1687: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1688: PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1689: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1690: }
1691: if (pc->ops->view) {
1692: (*pc->ops->view)(pc,viewer);
1693: }
1694: } else if (isdraw) {
1695: PetscDraw draw;
1696: char str[25];
1697: PetscReal x,y,bottom,h;
1698: PetscInt n;
1700: PetscViewerDrawGetDraw(viewer,0,&draw);
1701: PetscDrawGetCurrentPoint(draw,&x,&y);
1702: if (pc->mat) {
1703: MatGetSize(pc->mat,&n,NULL);
1704: PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1705: } else {
1706: PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1707: }
1708: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1709: bottom = y - h;
1710: PetscDrawPushCurrentPoint(draw,x,bottom);
1711: if (pc->ops->view) {
1712: (*pc->ops->view)(pc,viewer);
1713: }
1714: PetscDrawPopCurrentPoint(draw);
1715: #if defined(PETSC_HAVE_SAWS)
1716: } else if (issaws) {
1717: PetscMPIInt rank;
1719: PetscObjectName((PetscObject)pc);
1720: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1721: if (!((PetscObject)pc)->amsmem && !rank) {
1722: PetscObjectViewSAWs((PetscObject)pc,viewer);
1723: }
1724: if (pc->mat) {MatView(pc->mat,viewer);}
1725: if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1726: #endif
1727: }
1728: return(0);
1729: }
1731: /*@C
1732: PCRegister - Adds a method to the preconditioner package.
1734: Not collective
1736: Input Parameters:
1737: + name_solver - name of a new user-defined solver
1738: - routine_create - routine to create method context
1740: Notes:
1741: PCRegister() may be called multiple times to add several user-defined preconditioners.
1743: Sample usage:
1744: .vb
1745: PCRegister("my_solver", MySolverCreate);
1746: .ve
1748: Then, your solver can be chosen with the procedural interface via
1749: $ PCSetType(pc,"my_solver")
1750: or at runtime via the option
1751: $ -pc_type my_solver
1753: Level: advanced
1755: .keywords: PC, register
1757: .seealso: PCRegisterAll(), PCRegisterDestroy()
1758: @*/
1759: PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1760: {
1764: PCInitializePackage();
1765: PetscFunctionListAdd(&PCList,sname,function);
1766: return(0);
1767: }
1769: /*@
1770: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1772: Collective on PC
1774: Input Parameter:
1775: . pc - the preconditioner object
1777: Output Parameter:
1778: . mat - the explict preconditioned operator
1780: Notes:
1781: This computation is done by applying the operators to columns of the
1782: identity matrix.
1784: Currently, this routine uses a dense matrix format when 1 processor
1785: is used and a sparse format otherwise. This routine is costly in general,
1786: and is recommended for use only with relatively small systems.
1788: Level: advanced
1790: .keywords: PC, compute, explicit, operator
1792: .seealso: KSPComputeExplicitOperator()
1794: @*/
1795: PetscErrorCode PCComputeExplicitOperator(PC pc,Mat *mat)
1796: {
1797: Vec in,out;
1799: PetscInt i,M,m,*rows,start,end;
1800: PetscMPIInt size;
1801: MPI_Comm comm;
1802: PetscScalar *array,one = 1.0;
1808: PetscObjectGetComm((PetscObject)pc,&comm);
1809: MPI_Comm_size(comm,&size);
1811: if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1812: MatCreateVecs(pc->pmat,&in,0);
1813: VecDuplicate(in,&out);
1814: VecGetOwnershipRange(in,&start,&end);
1815: VecGetSize(in,&M);
1816: VecGetLocalSize(in,&m);
1817: PetscMalloc1(m+1,&rows);
1818: for (i=0; i<m; i++) rows[i] = start + i;
1820: MatCreate(comm,mat);
1821: MatSetSizes(*mat,m,m,M,M);
1822: if (size == 1) {
1823: MatSetType(*mat,MATSEQDENSE);
1824: MatSeqDenseSetPreallocation(*mat,NULL);
1825: } else {
1826: MatSetType(*mat,MATMPIAIJ);
1827: MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1828: }
1829: MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
1831: for (i=0; i<M; i++) {
1833: VecSet(in,0.0);
1834: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1835: VecAssemblyBegin(in);
1836: VecAssemblyEnd(in);
1838: /* should fix, allowing user to choose side */
1839: PCApply(pc,in,out);
1841: VecGetArray(out,&array);
1842: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1843: VecRestoreArray(out,&array);
1845: }
1846: PetscFree(rows);
1847: VecDestroy(&out);
1848: VecDestroy(&in);
1849: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1850: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1851: return(0);
1852: }
1854: /*@
1855: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1857: Collective on PC
1859: Input Parameters:
1860: + pc - the solver context
1861: . dim - the dimension of the coordinates 1, 2, or 3
1862: - coords - the coordinates
1864: Level: intermediate
1866: Notes:
1867: coords is an array of the 3D coordinates for the nodes on
1868: the local processor. So if there are 108 equation on a processor
1869: for a displacement finite element discretization of elasticity (so
1870: that there are 36 = 108/3 nodes) then the array must have 108
1871: double precision values (ie, 3 * 36). These x y z coordinates
1872: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1873: ... , N-1.z ].
1875: .seealso: MatSetNearNullSpace()
1876: @*/
1877: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1878: {
1884: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1885: return(0);
1886: }