Actual source code: precon.c
petsc-3.6.4 2016-04-12
2: /*
3: The PC (preconditioner) interface routines, callable by users.
4: */
5: #include <petsc/private/pcimpl.h> /*I "petscksp.h" I*/
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, PC_ApplyOnMproc;
15: PetscErrorCode PCGetDefaultType_Private(PC pc,const char *type[])
16: {
18: PetscMPIInt size;
19: PetscBool flg1,flg2,set,flg3;
22: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
23: if (pc->pmat) {
24: PetscErrorCode (*f)(Mat,MatReuse,Mat*);
25: PetscObjectQueryFunction((PetscObject)pc->pmat,"MatGetDiagonalBlock_C",&f);
26: if (size == 1) {
27: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
28: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
29: MatIsSymmetricKnown(pc->pmat,&set,&flg3);
30: if (flg1 && (!flg2 || (set && flg3))) {
31: *type = PCICC;
32: } else if (flg2) {
33: *type = PCILU;
34: } else if (f) { /* likely is a parallel matrix run on one processor */
35: *type = PCBJACOBI;
36: } else {
37: *type = PCNONE;
38: }
39: } else {
40: if (f) {
41: *type = PCBJACOBI;
42: } else {
43: *type = PCNONE;
44: }
45: }
46: } else {
47: if (size == 1) {
48: *type = PCILU;
49: } else {
50: *type = PCBJACOBI;
51: }
52: }
53: return(0);
54: }
58: /*@
59: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
61: Collective on PC
63: Input Parameter:
64: . pc - the preconditioner context
66: Level: developer
68: Notes: 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
70: .keywords: PC, destroy
72: .seealso: PCCreate(), PCSetUp()
73: @*/
74: PetscErrorCode PCReset(PC pc)
75: {
80: if (pc->ops->reset) {
81: (*pc->ops->reset)(pc);
82: }
83: VecDestroy(&pc->diagonalscaleright);
84: VecDestroy(&pc->diagonalscaleleft);
85: MatDestroy(&pc->pmat);
86: MatDestroy(&pc->mat);
88: pc->setupcalled = 0;
89: return(0);
90: }
94: /*@
95: PCDestroy - Destroys PC context that was created with PCCreate().
97: Collective on PC
99: Input Parameter:
100: . pc - the preconditioner context
102: Level: developer
104: .keywords: PC, destroy
106: .seealso: PCCreate(), PCSetUp()
107: @*/
108: PetscErrorCode PCDestroy(PC *pc)
109: {
113: if (!*pc) return(0);
115: if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; return(0);}
117: PCReset(*pc);
119: /* if memory was published with SAWs then destroy it */
120: PetscObjectSAWsViewOff((PetscObject)*pc);
121: if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
122: DMDestroy(&(*pc)->dm);
123: PetscHeaderDestroy(pc);
124: return(0);
125: }
129: /*@C
130: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
131: scaling as needed by certain time-stepping codes.
133: Logically Collective on PC
135: Input Parameter:
136: . pc - the preconditioner context
138: Output Parameter:
139: . flag - PETSC_TRUE if it applies the scaling
141: Level: developer
143: Notes: If this returns PETSC_TRUE then the system solved via the Krylov method is
144: $ D M A D^{-1} y = D M b for left preconditioning or
145: $ D A M D^{-1} z = D b for right preconditioning
147: .keywords: PC
149: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
150: @*/
151: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
152: {
156: *flag = pc->diagonalscale;
157: return(0);
158: }
162: /*@
163: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
164: scaling as needed by certain time-stepping codes.
166: Logically Collective on PC
168: Input Parameters:
169: + pc - the preconditioner context
170: - s - scaling vector
172: Level: intermediate
174: Notes: The system solved via the Krylov method is
175: $ D M A D^{-1} y = D M b for left preconditioning or
176: $ D A M D^{-1} z = D b for right preconditioning
178: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
180: .keywords: PC
182: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
183: @*/
184: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
185: {
191: pc->diagonalscale = PETSC_TRUE;
193: PetscObjectReference((PetscObject)s);
194: VecDestroy(&pc->diagonalscaleleft);
196: pc->diagonalscaleleft = s;
198: VecDuplicate(s,&pc->diagonalscaleright);
199: VecCopy(s,pc->diagonalscaleright);
200: VecReciprocal(pc->diagonalscaleright);
201: return(0);
202: }
206: /*@
207: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
209: Logically Collective on PC
211: Input Parameters:
212: + pc - the preconditioner context
213: . in - input vector
214: + out - scaled vector (maybe the same as in)
216: Level: intermediate
218: Notes: The system solved via the Krylov method is
219: $ D M A D^{-1} y = D M b for left preconditioning or
220: $ D A M D^{-1} z = D b for right preconditioning
222: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
224: If diagonal scaling is turned off and in is not out then in is copied to out
226: .keywords: PC
228: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
229: @*/
230: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
231: {
238: if (pc->diagonalscale) {
239: VecPointwiseMult(out,pc->diagonalscaleleft,in);
240: } else if (in != out) {
241: VecCopy(in,out);
242: }
243: return(0);
244: }
248: /*@
249: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
251: Logically Collective on PC
253: Input Parameters:
254: + pc - the preconditioner context
255: . in - input vector
256: + out - scaled vector (maybe the same as in)
258: Level: intermediate
260: Notes: The system solved via the Krylov method is
261: $ D M A D^{-1} y = D M b for left preconditioning or
262: $ D A M D^{-1} z = D b for right preconditioning
264: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
266: If diagonal scaling is turned off and in is not out then in is copied to out
268: .keywords: PC
270: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
271: @*/
272: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
273: {
280: if (pc->diagonalscale) {
281: VecPointwiseMult(out,pc->diagonalscaleright,in);
282: } else if (in != out) {
283: VecCopy(in,out);
284: }
285: return(0);
286: }
290: /*@
291: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
292: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
293: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
295: Logically Collective on PC
297: Input Parameters:
298: + pc - the preconditioner context
299: - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
301: Options Database Key:
302: . -pc_use_amat <true,false>
304: Notes:
305: For the common case in which the linear system matrix and the matrix used to construct the
306: preconditioner are identical, this routine is does nothing.
308: Level: intermediate
310: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
311: @*/
312: PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
313: {
316: pc->useAmat = flg;
317: return(0);
318: }
322: /*@
323: PCSetErrorIfFailure - Causes PC to generate an error if a FPE, for example a zero pivot, is detected.
325: Logically Collective on PC
327: Input Parameters:
328: + pc - iterative context obtained from PCCreate()
329: - flg - PETSC_TRUE indicates you want the error generated
331: Level: advanced
333: Notes:
334: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
335: 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
336: detected.
338: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
340: .keywords: PC, set, initial guess, nonzero
342: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
343: @*/
344: PetscErrorCode PCSetErrorIfFailure(PC pc,PetscBool flg)
345: {
349: pc->erroriffailure = flg;
350: return(0);
351: }
355: /*@
356: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
357: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
358: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
360: Logically Collective on PC
362: Input Parameter:
363: . pc - the preconditioner context
365: Output Parameter:
366: . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
368: Notes:
369: For the common case in which the linear system matrix and the matrix used to construct the
370: preconditioner are identical, this routine is does nothing.
372: Level: intermediate
374: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
375: @*/
376: PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
377: {
380: *flg = pc->useAmat;
381: return(0);
382: }
386: /*@
387: PCCreate - Creates a preconditioner context.
389: Collective on MPI_Comm
391: Input Parameter:
392: . comm - MPI communicator
394: Output Parameter:
395: . pc - location to put the preconditioner context
397: Notes:
398: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or ICC
399: in parallel. For dense matrices it is always PCNONE.
401: Level: developer
403: .keywords: PC, create, context
405: .seealso: PCSetUp(), PCApply(), PCDestroy()
406: @*/
407: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
408: {
409: PC pc;
414: *newpc = 0;
415: PCInitializePackage();
417: PetscHeaderCreate(pc,PC_CLASSID,"PC","Preconditioner","PC",comm,PCDestroy,PCView);
419: pc->mat = 0;
420: pc->pmat = 0;
421: pc->setupcalled = 0;
422: pc->setfromoptionscalled = 0;
423: pc->data = 0;
424: pc->diagonalscale = PETSC_FALSE;
425: pc->diagonalscaleleft = 0;
426: pc->diagonalscaleright = 0;
428: pc->modifysubmatrices = 0;
429: pc->modifysubmatricesP = 0;
431: *newpc = pc;
432: return(0);
434: }
436: /* -------------------------------------------------------------------------------*/
440: /*@
441: PCApply - Applies the preconditioner to a vector.
443: Collective on PC and Vec
445: Input Parameters:
446: + pc - the preconditioner context
447: - x - input vector
449: Output Parameter:
450: . y - output vector
452: Level: developer
454: .keywords: PC, apply
456: .seealso: PCApplyTranspose(), PCApplyBAorAB()
457: @*/
458: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
459: {
461: PetscInt m,n,mv,nv;
467: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
468: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
469: /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
470: MatGetLocalSize(pc->pmat,&m,&n);
471: VecGetLocalSize(x,&nv);
472: VecGetLocalSize(y,&mv);
473: 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);
474: 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);
475: VecLocked(y,3);
477: if (pc->setupcalled < 2) {
478: PCSetUp(pc);
479: }
480: if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
481: VecLockPush(x);
482: PetscLogEventBegin(PC_Apply,pc,x,y,0);
483: (*pc->ops->apply)(pc,x,y);
484: PetscLogEventEnd(PC_Apply,pc,x,y,0);
485: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
486: VecLockPop(x);
487: return(0);
488: }
492: /*@
493: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
495: Collective on PC and Vec
497: Input Parameters:
498: + pc - the preconditioner context
499: - x - input vector
501: Output Parameter:
502: . y - output vector
504: Notes:
505: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
507: Level: developer
509: .keywords: PC, apply, symmetric, left
511: .seealso: PCApply(), PCApplySymmetricRight()
512: @*/
513: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
514: {
521: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
522: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
523: if (pc->setupcalled < 2) {
524: PCSetUp(pc);
525: }
526: if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
527: VecLockPush(x);
528: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
529: (*pc->ops->applysymmetricleft)(pc,x,y);
530: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
531: VecLockPop(x);
532: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
533: return(0);
534: }
538: /*@
539: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
541: Collective on PC and Vec
543: Input Parameters:
544: + pc - the preconditioner context
545: - x - input vector
547: Output Parameter:
548: . y - output vector
550: Level: developer
552: Notes:
553: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
555: .keywords: PC, apply, symmetric, right
557: .seealso: PCApply(), PCApplySymmetricLeft()
558: @*/
559: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
560: {
567: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
568: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
569: if (pc->setupcalled < 2) {
570: PCSetUp(pc);
571: }
572: if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
573: VecLockPush(x);
574: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
575: (*pc->ops->applysymmetricright)(pc,x,y);
576: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
577: VecLockPop(x);
578: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
579: return(0);
580: }
584: /*@
585: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
587: Collective on PC and Vec
589: Input Parameters:
590: + pc - the preconditioner context
591: - x - input vector
593: Output Parameter:
594: . y - output vector
596: Notes: For complex numbers this applies the non-Hermitian transpose.
598: Developer Notes: We need to implement a PCApplyHermitianTranspose()
600: Level: developer
602: .keywords: PC, apply, transpose
604: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
605: @*/
606: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
607: {
614: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
615: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
616: if (pc->setupcalled < 2) {
617: PCSetUp(pc);
618: }
619: if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
620: VecLockPush(x);
621: PetscLogEventBegin(PC_Apply,pc,x,y,0);
622: (*pc->ops->applytranspose)(pc,x,y);
623: PetscLogEventEnd(PC_Apply,pc,x,y,0);
624: VecLockPop(x);
625: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
626: return(0);
627: }
631: /*@
632: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
634: Collective on PC and Vec
636: Input Parameters:
637: . pc - the preconditioner context
639: Output Parameter:
640: . flg - PETSC_TRUE if a transpose operation is defined
642: Level: developer
644: .keywords: PC, apply, transpose
646: .seealso: PCApplyTranspose()
647: @*/
648: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
649: {
653: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
654: else *flg = PETSC_FALSE;
655: return(0);
656: }
660: /*@
661: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
663: Collective on PC and Vec
665: Input Parameters:
666: + pc - the preconditioner context
667: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
668: . x - input vector
669: - work - work vector
671: Output Parameter:
672: . y - output vector
674: Level: developer
676: Notes: 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
677: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
679: .keywords: PC, apply, operator
681: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
682: @*/
683: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
684: {
692: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
693: 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");
694: if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
695: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
697: if (pc->setupcalled < 2) {
698: PCSetUp(pc);
699: }
701: if (pc->diagonalscale) {
702: if (pc->ops->applyBA) {
703: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
704: VecDuplicate(x,&work2);
705: PCDiagonalScaleRight(pc,x,work2);
706: (*pc->ops->applyBA)(pc,side,work2,y,work);
707: PCDiagonalScaleLeft(pc,y,y);
708: VecDestroy(&work2);
709: } else if (side == PC_RIGHT) {
710: PCDiagonalScaleRight(pc,x,y);
711: PCApply(pc,y,work);
712: MatMult(pc->mat,work,y);
713: PCDiagonalScaleLeft(pc,y,y);
714: } else if (side == PC_LEFT) {
715: PCDiagonalScaleRight(pc,x,y);
716: MatMult(pc->mat,y,work);
717: PCApply(pc,work,y);
718: PCDiagonalScaleLeft(pc,y,y);
719: } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
720: } else {
721: if (pc->ops->applyBA) {
722: (*pc->ops->applyBA)(pc,side,x,y,work);
723: } else if (side == PC_RIGHT) {
724: PCApply(pc,x,work);
725: MatMult(pc->mat,work,y);
726: } else if (side == PC_LEFT) {
727: MatMult(pc->mat,x,work);
728: PCApply(pc,work,y);
729: } else if (side == PC_SYMMETRIC) {
730: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
731: PCApplySymmetricRight(pc,x,work);
732: MatMult(pc->mat,work,y);
733: VecCopy(y,work);
734: PCApplySymmetricLeft(pc,work,y);
735: }
736: }
737: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
738: return(0);
739: }
743: /*@
744: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
745: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
746: NOT tr(B*A) = tr(A)*tr(B).
748: Collective on PC and Vec
750: Input Parameters:
751: + pc - the preconditioner context
752: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
753: . x - input vector
754: - work - work vector
756: Output Parameter:
757: . y - output vector
760: Notes: 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
761: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
763: Level: developer
765: .keywords: PC, apply, operator, transpose
767: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
768: @*/
769: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
770: {
778: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
779: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
780: if (pc->ops->applyBAtranspose) {
781: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
782: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
783: return(0);
784: }
785: if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
787: if (pc->setupcalled < 2) {
788: PCSetUp(pc);
789: }
791: if (side == PC_RIGHT) {
792: PCApplyTranspose(pc,x,work);
793: MatMultTranspose(pc->mat,work,y);
794: } else if (side == PC_LEFT) {
795: MatMultTranspose(pc->mat,x,work);
796: PCApplyTranspose(pc,work,y);
797: }
798: /* add support for PC_SYMMETRIC */
799: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
800: return(0);
801: }
803: /* -------------------------------------------------------------------------------*/
807: /*@
808: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
809: built-in fast application of Richardson's method.
811: Not Collective
813: Input Parameter:
814: . pc - the preconditioner
816: Output Parameter:
817: . exists - PETSC_TRUE or PETSC_FALSE
819: Level: developer
821: .keywords: PC, apply, Richardson, exists
823: .seealso: PCApplyRichardson()
824: @*/
825: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
826: {
830: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
831: else *exists = PETSC_FALSE;
832: return(0);
833: }
837: /*@
838: PCApplyRichardson - Applies several steps of Richardson iteration with
839: the particular preconditioner. This routine is usually used by the
840: Krylov solvers and not the application code directly.
842: Collective on PC
844: Input Parameters:
845: + pc - the preconditioner context
846: . b - the right hand side
847: . w - one work vector
848: . rtol - relative decrease in residual norm convergence criteria
849: . abstol - absolute residual norm convergence criteria
850: . dtol - divergence residual norm increase criteria
851: . its - the number of iterations to apply.
852: - guesszero - if the input x contains nonzero initial guess
854: Output Parameter:
855: + outits - number of iterations actually used (for SOR this always equals its)
856: . reason - the reason the apply terminated
857: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
859: Notes:
860: Most preconditioners do not support this function. Use the command
861: PCApplyRichardsonExists() to determine if one does.
863: Except for the multigrid PC this routine ignores the convergence tolerances
864: and always runs for the number of iterations
866: Level: developer
868: .keywords: PC, apply, Richardson
870: .seealso: PCApplyRichardsonExists()
871: @*/
872: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
873: {
881: if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
882: if (pc->setupcalled < 2) {
883: PCSetUp(pc);
884: }
885: if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
886: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
887: return(0);
888: }
892: /*@
893: PCGetSetUpFailedReason - Gets the reason a PCSetUp() failed or 0 if it did not fail
895: Logically Collective on PC
897: Input Parameter:
898: . pc - the preconditioner context
900: Output Parameter:
901: . reason - the reason it failed, currently only -1
903: Level: advanced
905: .keywords: PC, setup
907: .seealso: PCCreate(), PCApply(), PCDestroy()
908: @*/
909: PetscErrorCode PCGetSetUpFailedReason(PC pc,PetscInt *reason)
910: {
912: if (pc->setupcalled < 0) *reason = pc->setupcalled;
913: else *reason = 0;
914: return(0);
915: }
918: /*
919: a setupcall of 0 indicates never setup,
920: 1 indicates has been previously setup
921: -1 indicates a PCSetUp() was attempted and failed
922: */
925: /*@
926: PCSetUp - Prepares for the use of a preconditioner.
928: Collective on PC
930: Input Parameter:
931: . pc - the preconditioner context
933: Level: developer
935: .keywords: PC, setup
937: .seealso: PCCreate(), PCApply(), PCDestroy()
938: @*/
939: PetscErrorCode PCSetUp(PC pc)
940: {
941: PetscErrorCode ierr;
942: const char *def;
943: PetscObjectState matstate, matnonzerostate;
947: if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
949: if (pc->setupcalled && pc->reusepreconditioner) {
950: PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
951: return(0);
952: }
954: PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
955: MatGetNonzeroState(pc->pmat,&matnonzerostate);
956: if (!pc->setupcalled) {
957: PetscInfo(pc,"Setting up PC for first time\n");
958: pc->flag = DIFFERENT_NONZERO_PATTERN;
959: } else if (matstate == pc->matstate) {
960: PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
961: return(0);
962: } else {
963: if (matnonzerostate > pc->matnonzerostate) {
964: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
965: pc->flag = DIFFERENT_NONZERO_PATTERN;
966: } else {
967: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
968: pc->flag = SAME_NONZERO_PATTERN;
969: }
970: }
971: pc->matstate = matstate;
972: pc->matnonzerostate = matnonzerostate;
974: if (!((PetscObject)pc)->type_name) {
975: PCGetDefaultType_Private(pc,&def);
976: PCSetType(pc,def);
977: }
979: MatSetErrorIfFPE(pc->pmat,pc->erroriffailure);
980: MatSetErrorIfFPE(pc->mat,pc->erroriffailure);
981: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
982: if (pc->ops->setup) {
983: (*pc->ops->setup)(pc);
984: }
985: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
986: if (!pc->setupcalled) pc->setupcalled = 1;
987: return(0);
988: }
992: /*@
993: PCSetUpOnBlocks - Sets up the preconditioner for each block in
994: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
995: methods.
997: Collective on PC
999: Input Parameters:
1000: . pc - the preconditioner context
1002: Level: developer
1004: .keywords: PC, setup, blocks
1006: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
1007: @*/
1008: PetscErrorCode PCSetUpOnBlocks(PC pc)
1009: {
1014: if (!pc->ops->setuponblocks) return(0);
1015: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
1016: (*pc->ops->setuponblocks)(pc);
1017: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
1018: return(0);
1019: }
1023: /*@C
1024: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1025: submatrices that arise within certain subdomain-based preconditioners.
1026: The basic submatrices are extracted from the preconditioner matrix as
1027: usual; the user can then alter these (for example, to set different boundary
1028: conditions for each submatrix) before they are used for the local solves.
1030: Logically Collective on PC
1032: Input Parameters:
1033: + pc - the preconditioner context
1034: . func - routine for modifying the submatrices
1035: - ctx - optional user-defined context (may be null)
1037: Calling sequence of func:
1038: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
1040: . row - an array of index sets that contain the global row numbers
1041: that comprise each local submatrix
1042: . col - an array of index sets that contain the global column numbers
1043: that comprise each local submatrix
1044: . submat - array of local submatrices
1045: - ctx - optional user-defined context for private data for the
1046: user-defined func routine (may be null)
1048: Notes:
1049: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
1050: KSPSolve().
1052: A routine set by PCSetModifySubMatrices() is currently called within
1053: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
1054: preconditioners. All other preconditioners ignore this routine.
1056: Level: advanced
1058: .keywords: PC, set, modify, submatrices
1060: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
1061: @*/
1062: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1063: {
1066: pc->modifysubmatrices = func;
1067: pc->modifysubmatricesP = ctx;
1068: return(0);
1069: }
1073: /*@C
1074: PCModifySubMatrices - Calls an optional user-defined routine within
1075: certain preconditioners if one has been set with PCSetModifySubMarices().
1077: Collective on PC
1079: Input Parameters:
1080: + pc - the preconditioner context
1081: . nsub - the number of local submatrices
1082: . row - an array of index sets that contain the global row numbers
1083: that comprise each local submatrix
1084: . col - an array of index sets that contain the global column numbers
1085: that comprise each local submatrix
1086: . submat - array of local submatrices
1087: - ctx - optional user-defined context for private data for the
1088: user-defined routine (may be null)
1090: Output Parameter:
1091: . submat - array of local submatrices (the entries of which may
1092: have been modified)
1094: Notes:
1095: The user should NOT generally call this routine, as it will
1096: automatically be called within certain preconditioners (currently
1097: block Jacobi, additive Schwarz) if set.
1099: The basic submatrices are extracted from the preconditioner matrix
1100: as usual; the user can then alter these (for example, to set different
1101: boundary conditions for each submatrix) before they are used for the
1102: local solves.
1104: Level: developer
1106: .keywords: PC, modify, submatrices
1108: .seealso: PCSetModifySubMatrices()
1109: @*/
1110: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1111: {
1116: if (!pc->modifysubmatrices) return(0);
1117: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1118: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1119: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1120: return(0);
1121: }
1125: /*@
1126: PCSetOperators - Sets the matrix associated with the linear system and
1127: a (possibly) different one associated with the preconditioner.
1129: Logically Collective on PC and Mat
1131: Input Parameters:
1132: + pc - the preconditioner context
1133: . Amat - the matrix that defines the linear system
1134: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1136: Notes:
1137: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1139: If you wish to replace either Amat or Pmat but leave the other one untouched then
1140: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1141: on it and then pass it back in in your call to KSPSetOperators().
1143: More Notes about Repeated Solution of Linear Systems:
1144: PETSc does NOT reset the matrix entries of either Amat or Pmat
1145: to zero after a linear solve; the user is completely responsible for
1146: matrix assembly. See the routine MatZeroEntries() if desiring to
1147: zero all elements of a matrix.
1149: Level: intermediate
1151: .keywords: PC, set, operators, matrix, linear system
1153: .seealso: PCGetOperators(), MatZeroEntries()
1154: @*/
1155: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1156: {
1157: PetscErrorCode ierr;
1158: PetscInt m1,n1,m2,n2;
1166: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1167: MatGetLocalSize(Amat,&m1,&n1);
1168: MatGetLocalSize(pc->mat,&m2,&n2);
1169: 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);
1170: MatGetLocalSize(Pmat,&m1,&n1);
1171: MatGetLocalSize(pc->pmat,&m2,&n2);
1172: 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);
1173: }
1175: if (Pmat != pc->pmat) {
1176: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1177: pc->matnonzerostate = -1;
1178: pc->matstate = -1;
1179: }
1181: /* reference first in case the matrices are the same */
1182: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1183: MatDestroy(&pc->mat);
1184: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1185: MatDestroy(&pc->pmat);
1186: pc->mat = Amat;
1187: pc->pmat = Pmat;
1188: return(0);
1189: }
1193: /*@
1194: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1196: Logically Collective on PC
1198: Input Parameters:
1199: + pc - the preconditioner context
1200: - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1202: Level: intermediate
1204: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1205: @*/
1206: PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1207: {
1210: pc->reusepreconditioner = flag;
1211: return(0);
1212: }
1216: /*@
1217: PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1219: Not Collective
1221: Input Parameter:
1222: . pc - the preconditioner context
1224: Output Parameter:
1225: . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1227: Level: intermediate
1229: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1230: @*/
1231: PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1232: {
1235: *flag = pc->reusepreconditioner;
1236: return(0);
1237: }
1241: /*@C
1242: PCGetOperators - Gets the matrix associated with the linear system and
1243: possibly a different one associated with the preconditioner.
1245: Not collective, though parallel Mats are returned if the PC is parallel
1247: Input Parameter:
1248: . pc - the preconditioner context
1250: Output Parameters:
1251: + Amat - the matrix defining the linear system
1252: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1254: Level: intermediate
1256: Notes: Does not increase the reference count of the matrices, so you should not destroy them
1258: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1259: are created in PC and returned to the user. In this case, if both operators
1260: mat and pmat are requested, two DIFFERENT operators will be returned. If
1261: only one is requested both operators in the PC will be the same (i.e. as
1262: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1263: The user must set the sizes of the returned matrices and their type etc just
1264: as if the user created them with MatCreate(). For example,
1266: $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1267: $ set size, type, etc of Amat
1269: $ MatCreate(comm,&mat);
1270: $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1271: $ PetscObjectDereference((PetscObject)mat);
1272: $ set size, type, etc of Amat
1274: and
1276: $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1277: $ set size, type, etc of Amat and Pmat
1279: $ MatCreate(comm,&Amat);
1280: $ MatCreate(comm,&Pmat);
1281: $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1282: $ PetscObjectDereference((PetscObject)Amat);
1283: $ PetscObjectDereference((PetscObject)Pmat);
1284: $ set size, type, etc of Amat and Pmat
1286: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1287: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1288: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1289: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1290: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1291: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1292: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1293: it can be created for you?
1296: .keywords: PC, get, operators, matrix, linear system
1298: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1299: @*/
1300: PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1301: {
1306: if (Amat) {
1307: if (!pc->mat) {
1308: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1309: pc->mat = pc->pmat;
1310: PetscObjectReference((PetscObject)pc->mat);
1311: } else { /* both Amat and Pmat are empty */
1312: MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1313: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1314: pc->pmat = pc->mat;
1315: PetscObjectReference((PetscObject)pc->pmat);
1316: }
1317: }
1318: }
1319: *Amat = pc->mat;
1320: }
1321: if (Pmat) {
1322: if (!pc->pmat) {
1323: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1324: pc->pmat = pc->mat;
1325: PetscObjectReference((PetscObject)pc->pmat);
1326: } else {
1327: MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1328: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1329: pc->mat = pc->pmat;
1330: PetscObjectReference((PetscObject)pc->mat);
1331: }
1332: }
1333: }
1334: *Pmat = pc->pmat;
1335: }
1336: return(0);
1337: }
1341: /*@C
1342: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1343: possibly a different one associated with the preconditioner have been set in the PC.
1345: Not collective, though the results on all processes should be the same
1347: Input Parameter:
1348: . pc - the preconditioner context
1350: Output Parameters:
1351: + mat - the matrix associated with the linear system was set
1352: - pmat - matrix associated with the preconditioner was set, usually the same
1354: Level: intermediate
1356: .keywords: PC, get, operators, matrix, linear system
1358: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1359: @*/
1360: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1361: {
1364: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1365: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1366: return(0);
1367: }
1371: /*@
1372: PCFactorGetMatrix - Gets the factored matrix from the
1373: preconditioner context. This routine is valid only for the LU,
1374: incomplete LU, Cholesky, and incomplete Cholesky methods.
1376: Not Collective on PC though Mat is parallel if PC is parallel
1378: Input Parameters:
1379: . pc - the preconditioner context
1381: Output parameters:
1382: . mat - the factored matrix
1384: Level: advanced
1386: Notes: Does not increase the reference count for the matrix so DO NOT destroy it
1388: .keywords: PC, get, factored, matrix
1389: @*/
1390: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1391: {
1397: if (pc->ops->getfactoredmatrix) {
1398: (*pc->ops->getfactoredmatrix)(pc,mat);
1399: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1400: return(0);
1401: }
1405: /*@C
1406: PCSetOptionsPrefix - Sets the prefix used for searching for all
1407: PC options in the database.
1409: Logically Collective on PC
1411: Input Parameters:
1412: + pc - the preconditioner context
1413: - prefix - the prefix string to prepend to all PC option requests
1415: Notes:
1416: A hyphen (-) must NOT be given at the beginning of the prefix name.
1417: The first character of all runtime options is AUTOMATICALLY the
1418: hyphen.
1420: Level: advanced
1422: .keywords: PC, set, options, prefix, database
1424: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1425: @*/
1426: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1427: {
1432: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1433: return(0);
1434: }
1438: /*@C
1439: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1440: PC options in the database.
1442: Logically Collective on PC
1444: Input Parameters:
1445: + pc - the preconditioner context
1446: - prefix - the prefix string to prepend to all PC option requests
1448: Notes:
1449: A hyphen (-) must NOT be given at the beginning of the prefix name.
1450: The first character of all runtime options is AUTOMATICALLY the
1451: hyphen.
1453: Level: advanced
1455: .keywords: PC, append, options, prefix, database
1457: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1458: @*/
1459: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1460: {
1465: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1466: return(0);
1467: }
1471: /*@C
1472: PCGetOptionsPrefix - Gets the prefix used for searching for all
1473: PC options in the database.
1475: Not Collective
1477: Input Parameters:
1478: . pc - the preconditioner context
1480: Output Parameters:
1481: . prefix - pointer to the prefix string used, is returned
1483: Notes: On the fortran side, the user should pass in a string 'prifix' of
1484: sufficient length to hold the prefix.
1486: Level: advanced
1488: .keywords: PC, get, options, prefix, database
1490: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1491: @*/
1492: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1493: {
1499: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1500: return(0);
1501: }
1505: /*@
1506: PCPreSolve - Optional pre-solve phase, intended for any
1507: preconditioner-specific actions that must be performed before
1508: the iterative solve itself.
1510: Collective on PC
1512: Input Parameters:
1513: + pc - the preconditioner context
1514: - ksp - the Krylov subspace context
1516: Level: developer
1518: Sample of Usage:
1519: .vb
1520: PCPreSolve(pc,ksp);
1521: KSPSolve(ksp,b,x);
1522: PCPostSolve(pc,ksp);
1523: .ve
1525: Notes:
1526: The pre-solve phase is distinct from the PCSetUp() phase.
1528: KSPSolve() calls this directly, so is rarely called by the user.
1530: .keywords: PC, pre-solve
1532: .seealso: PCPostSolve()
1533: @*/
1534: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1535: {
1537: Vec x,rhs;
1542: pc->presolvedone++;
1543: if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1544: KSPGetSolution(ksp,&x);
1545: KSPGetRhs(ksp,&rhs);
1547: if (pc->ops->presolve) {
1548: (*pc->ops->presolve)(pc,ksp,rhs,x);
1549: }
1550: return(0);
1551: }
1555: /*@
1556: PCPostSolve - Optional post-solve phase, intended for any
1557: preconditioner-specific actions that must be performed after
1558: the iterative solve itself.
1560: Collective on PC
1562: Input Parameters:
1563: + pc - the preconditioner context
1564: - ksp - the Krylov subspace context
1566: Sample of Usage:
1567: .vb
1568: PCPreSolve(pc,ksp);
1569: KSPSolve(ksp,b,x);
1570: PCPostSolve(pc,ksp);
1571: .ve
1573: Note:
1574: KSPSolve() calls this routine directly, so it is rarely called by the user.
1576: Level: developer
1578: .keywords: PC, post-solve
1580: .seealso: PCPreSolve(), KSPSolve()
1581: @*/
1582: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1583: {
1585: Vec x,rhs;
1590: pc->presolvedone--;
1591: KSPGetSolution(ksp,&x);
1592: KSPGetRhs(ksp,&rhs);
1593: if (pc->ops->postsolve) {
1594: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1595: }
1596: return(0);
1597: }
1601: /*@C
1602: PCLoad - Loads a PC that has been stored in binary with PCView().
1604: Collective on PetscViewer
1606: Input Parameters:
1607: + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1608: some related function before a call to PCLoad().
1609: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1611: Level: intermediate
1613: Notes:
1614: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1616: Notes for advanced users:
1617: Most users should not need to know the details of the binary storage
1618: format, since PCLoad() and PCView() completely hide these details.
1619: But for anyone who's interested, the standard binary matrix storage
1620: format is
1621: .vb
1622: has not yet been determined
1623: .ve
1625: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1626: @*/
1627: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1628: {
1630: PetscBool isbinary;
1631: PetscInt classid;
1632: char type[256];
1637: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1638: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1640: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1641: if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1642: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1643: PCSetType(newdm, type);
1644: if (newdm->ops->load) {
1645: (*newdm->ops->load)(newdm,viewer);
1646: }
1647: return(0);
1648: }
1650: #include <petscdraw.h>
1651: #if defined(PETSC_HAVE_SAWS)
1652: #include <petscviewersaws.h>
1653: #endif
1656: /*@C
1657: PCView - Prints the PC data structure.
1659: Collective on PC
1661: Input Parameters:
1662: + PC - the PC context
1663: - viewer - optional visualization context
1665: Note:
1666: The available visualization contexts include
1667: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1668: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1669: output where only the first processor opens
1670: the file. All other processors send their
1671: data to the first processor to print.
1673: The user can open an alternative visualization contexts with
1674: PetscViewerASCIIOpen() (output to a specified file).
1676: Level: developer
1678: .keywords: PC, view
1680: .seealso: KSPView(), PetscViewerASCIIOpen()
1681: @*/
1682: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1683: {
1684: PCType cstr;
1685: PetscErrorCode ierr;
1686: PetscBool iascii,isstring,isbinary,isdraw;
1687: PetscViewerFormat format;
1688: #if defined(PETSC_HAVE_SAWS)
1689: PetscBool issaws;
1690: #endif
1694: if (!viewer) {
1695: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1696: }
1700: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1701: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1702: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1703: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1704: #if defined(PETSC_HAVE_SAWS)
1705: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1706: #endif
1708: if (iascii) {
1709: PetscViewerGetFormat(viewer,&format);
1710: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1711: if (!pc->setupcalled) {
1712: PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");
1713: }
1714: if (pc->ops->view) {
1715: PetscViewerASCIIPushTab(viewer);
1716: (*pc->ops->view)(pc,viewer);
1717: PetscViewerASCIIPopTab(viewer);
1718: }
1719: if (pc->mat) {
1720: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1721: if (pc->pmat == pc->mat) {
1722: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1723: PetscViewerASCIIPushTab(viewer);
1724: MatView(pc->mat,viewer);
1725: PetscViewerASCIIPopTab(viewer);
1726: } else {
1727: if (pc->pmat) {
1728: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1729: } else {
1730: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1731: }
1732: PetscViewerASCIIPushTab(viewer);
1733: MatView(pc->mat,viewer);
1734: if (pc->pmat) {MatView(pc->pmat,viewer);}
1735: PetscViewerASCIIPopTab(viewer);
1736: }
1737: PetscViewerPopFormat(viewer);
1738: }
1739: } else if (isstring) {
1740: PCGetType(pc,&cstr);
1741: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1742: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1743: } else if (isbinary) {
1744: PetscInt classid = PC_FILE_CLASSID;
1745: MPI_Comm comm;
1746: PetscMPIInt rank;
1747: char type[256];
1749: PetscObjectGetComm((PetscObject)pc,&comm);
1750: MPI_Comm_rank(comm,&rank);
1751: if (!rank) {
1752: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1753: PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1754: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1755: }
1756: if (pc->ops->view) {
1757: (*pc->ops->view)(pc,viewer);
1758: }
1759: } else if (isdraw) {
1760: PetscDraw draw;
1761: char str[25];
1762: PetscReal x,y,bottom,h;
1763: PetscInt n;
1765: PetscViewerDrawGetDraw(viewer,0,&draw);
1766: PetscDrawGetCurrentPoint(draw,&x,&y);
1767: if (pc->mat) {
1768: MatGetSize(pc->mat,&n,NULL);
1769: PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1770: } else {
1771: PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1772: }
1773: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1774: bottom = y - h;
1775: PetscDrawPushCurrentPoint(draw,x,bottom);
1776: if (pc->ops->view) {
1777: (*pc->ops->view)(pc,viewer);
1778: }
1779: PetscDrawPopCurrentPoint(draw);
1780: #if defined(PETSC_HAVE_SAWS)
1781: } else if (issaws) {
1782: PetscMPIInt rank;
1784: PetscObjectName((PetscObject)pc);
1785: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1786: if (!((PetscObject)pc)->amsmem && !rank) {
1787: PetscObjectViewSAWs((PetscObject)pc,viewer);
1788: }
1789: if (pc->mat) {MatView(pc->mat,viewer);}
1790: if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1791: #endif
1792: }
1793: return(0);
1794: }
1799: /*@
1800: PCSetInitialGuessNonzero - Tells the iterative solver that the
1801: initial guess is nonzero; otherwise PC assumes the initial guess
1802: is to be zero (and thus zeros it out before solving).
1804: Logically Collective on PC
1806: Input Parameters:
1807: + pc - iterative context obtained from PCCreate()
1808: - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
1810: Level: Developer
1812: Notes:
1813: This is a weird function. Since PC's are linear operators on the right hand side they
1814: CANNOT use an initial guess. This function is for the "pass-through" preconditioners
1815: PCKSP and PCREDUNDANT and causes the inner KSP object to use the nonzero
1816: initial guess. Not currently working for PCREDUNDANT, that has to be rewritten to use KSP.
1819: .keywords: PC, set, initial guess, nonzero
1821: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
1822: @*/
1823: PetscErrorCode PCSetInitialGuessNonzero(PC pc,PetscBool flg)
1824: {
1827: pc->nonzero_guess = flg;
1828: return(0);
1829: }
1833: /*@
1834: PCGetInitialGuessNonzero - Determines if the iterative solver assumes that the
1835: initial guess is nonzero; otherwise PC assumes the initial guess
1836: is to be zero (and thus zeros it out before solving).
1838: Logically Collective on PC
1840: Input Parameter:
1841: . pc - iterative context obtained from PCCreate()
1843: Output Parameter:
1844: . flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
1846: Level: Developer
1848: .keywords: PC, set, initial guess, nonzero
1850: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll(), PCSetInitialGuessNonzero()
1851: @*/
1852: PetscErrorCode PCGetInitialGuessNonzero(PC pc,PetscBool *flg)
1853: {
1855: *flg = pc->nonzero_guess;
1856: return(0);
1857: }
1861: /*@C
1862: PCRegister - Adds a method to the preconditioner package.
1864: Not collective
1866: Input Parameters:
1867: + name_solver - name of a new user-defined solver
1868: - routine_create - routine to create method context
1870: Notes:
1871: PCRegister() may be called multiple times to add several user-defined preconditioners.
1873: Sample usage:
1874: .vb
1875: PCRegister("my_solver", MySolverCreate);
1876: .ve
1878: Then, your solver can be chosen with the procedural interface via
1879: $ PCSetType(pc,"my_solver")
1880: or at runtime via the option
1881: $ -pc_type my_solver
1883: Level: advanced
1885: .keywords: PC, register
1887: .seealso: PCRegisterAll(), PCRegisterDestroy()
1888: @*/
1889: PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1890: {
1894: PetscFunctionListAdd(&PCList,sname,function);
1895: return(0);
1896: }
1900: /*@
1901: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1903: Collective on PC
1905: Input Parameter:
1906: . pc - the preconditioner object
1908: Output Parameter:
1909: . mat - the explict preconditioned operator
1911: Notes:
1912: This computation is done by applying the operators to columns of the
1913: identity matrix.
1915: Currently, this routine uses a dense matrix format when 1 processor
1916: is used and a sparse format otherwise. This routine is costly in general,
1917: and is recommended for use only with relatively small systems.
1919: Level: advanced
1921: .keywords: PC, compute, explicit, operator
1923: .seealso: KSPComputeExplicitOperator()
1925: @*/
1926: PetscErrorCode PCComputeExplicitOperator(PC pc,Mat *mat)
1927: {
1928: Vec in,out;
1930: PetscInt i,M,m,*rows,start,end;
1931: PetscMPIInt size;
1932: MPI_Comm comm;
1933: PetscScalar *array,one = 1.0;
1939: PetscObjectGetComm((PetscObject)pc,&comm);
1940: MPI_Comm_size(comm,&size);
1942: if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1943: MatCreateVecs(pc->pmat,&in,0);
1944: VecDuplicate(in,&out);
1945: VecGetOwnershipRange(in,&start,&end);
1946: VecGetSize(in,&M);
1947: VecGetLocalSize(in,&m);
1948: PetscMalloc1(m+1,&rows);
1949: for (i=0; i<m; i++) rows[i] = start + i;
1951: MatCreate(comm,mat);
1952: MatSetSizes(*mat,m,m,M,M);
1953: if (size == 1) {
1954: MatSetType(*mat,MATSEQDENSE);
1955: MatSeqDenseSetPreallocation(*mat,NULL);
1956: } else {
1957: MatSetType(*mat,MATMPIAIJ);
1958: MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1959: }
1960: MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
1962: for (i=0; i<M; i++) {
1964: VecSet(in,0.0);
1965: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1966: VecAssemblyBegin(in);
1967: VecAssemblyEnd(in);
1969: /* should fix, allowing user to choose side */
1970: PCApply(pc,in,out);
1972: VecGetArray(out,&array);
1973: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1974: VecRestoreArray(out,&array);
1976: }
1977: PetscFree(rows);
1978: VecDestroy(&out);
1979: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1980: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1981: return(0);
1982: }
1986: /*@
1987: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1989: Collective on PC
1991: Input Parameters:
1992: + pc - the solver context
1993: . dim - the dimension of the coordinates 1, 2, or 3
1994: - coords - the coordinates
1996: Level: intermediate
1998: Notes: coords is an array of the 3D coordinates for the nodes on
1999: the local processor. So if there are 108 equation on a processor
2000: for a displacement finite element discretization of elasticity (so
2001: that there are 36 = 108/3 nodes) then the array must have 108
2002: double precision values (ie, 3 * 36). These x y z coordinates
2003: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
2004: ... , N-1.z ].
2006: .seealso: MatSetNearNullSpace
2007: @*/
2008: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2009: {
2013: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
2014: return(0);
2015: }