Actual source code: precon.c
petsc-3.8.4 2018-03-24
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, PC_ApplyOnMproc;
12: PetscInt PetscMGLevelId;
14: PetscErrorCode PCGetDefaultType_Private(PC pc,const char *type[])
15: {
17: PetscMPIInt size;
18: PetscBool flg1,flg2,set,flg3;
21: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
22: if (pc->pmat) {
23: void (*f)(void);
24: MatShellGetOperation(pc->pmat,MATOP_GET_DIAGONAL_BLOCK,&f);
25: if (size == 1) {
26: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
27: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
28: MatIsSymmetricKnown(pc->pmat,&set,&flg3);
29: if (flg1 && (!flg2 || (set && flg3))) {
30: *type = PCICC;
31: } else if (flg2) {
32: *type = PCILU;
33: } else if (f) { /* likely is a parallel matrix run on one processor */
34: *type = PCBJACOBI;
35: } else {
36: *type = PCNONE;
37: }
38: } else {
39: if (f) {
40: *type = PCBJACOBI;
41: } else {
42: *type = PCNONE;
43: }
44: }
45: } else {
46: if (size == 1) {
47: *type = PCILU;
48: } else {
49: *type = PCBJACOBI;
50: }
51: }
52: return(0);
53: }
55: /*@
56: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
58: Collective on PC
60: Input Parameter:
61: . pc - the preconditioner context
63: Level: developer
65: 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
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: If this returns PETSC_TRUE then the system solved via the Krylov method is
137: $ D M A D^{-1} y = D M b for left preconditioning or
138: $ D A M D^{-1} z = D b for right preconditioning
140: .keywords: PC
142: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
143: @*/
144: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
145: {
149: *flag = pc->diagonalscale;
150: return(0);
151: }
153: /*@
154: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
155: scaling as needed by certain time-stepping codes.
157: Logically Collective on PC
159: Input Parameters:
160: + pc - the preconditioner context
161: - s - scaling vector
163: Level: intermediate
165: Notes: The system solved via the Krylov method is
166: $ D M A D^{-1} y = D M b for left preconditioning or
167: $ D A M D^{-1} z = D b for right preconditioning
169: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
171: .keywords: PC
173: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
174: @*/
175: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
176: {
182: pc->diagonalscale = PETSC_TRUE;
184: PetscObjectReference((PetscObject)s);
185: VecDestroy(&pc->diagonalscaleleft);
187: pc->diagonalscaleleft = s;
189: VecDuplicate(s,&pc->diagonalscaleright);
190: VecCopy(s,pc->diagonalscaleright);
191: VecReciprocal(pc->diagonalscaleright);
192: return(0);
193: }
195: /*@
196: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
198: Logically Collective on PC
200: Input Parameters:
201: + pc - the preconditioner context
202: . in - input vector
203: + out - scaled vector (maybe the same as in)
205: Level: intermediate
207: Notes: The system solved via the Krylov method is
208: $ D M A D^{-1} y = D M b for left preconditioning or
209: $ D A M D^{-1} z = D b for right preconditioning
211: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
213: If diagonal scaling is turned off and in is not out then in is copied to out
215: .keywords: PC
217: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
218: @*/
219: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
220: {
227: if (pc->diagonalscale) {
228: VecPointwiseMult(out,pc->diagonalscaleleft,in);
229: } else if (in != out) {
230: VecCopy(in,out);
231: }
232: return(0);
233: }
235: /*@
236: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
238: Logically Collective on PC
240: Input Parameters:
241: + pc - the preconditioner context
242: . in - input vector
243: + out - scaled vector (maybe the same as in)
245: Level: intermediate
247: Notes: The system solved via the Krylov method is
248: $ D M A D^{-1} y = D M b for left preconditioning or
249: $ D A M D^{-1} z = D b for right preconditioning
251: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
253: If diagonal scaling is turned off and in is not out then in is copied to out
255: .keywords: PC
257: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
258: @*/
259: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
260: {
267: if (pc->diagonalscale) {
268: VecPointwiseMult(out,pc->diagonalscaleright,in);
269: } else if (in != out) {
270: VecCopy(in,out);
271: }
272: return(0);
273: }
275: /*@
276: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
277: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
278: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
280: Logically Collective on PC
282: Input Parameters:
283: + pc - the preconditioner context
284: - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
286: Options Database Key:
287: . -pc_use_amat <true,false>
289: Notes:
290: For the common case in which the linear system matrix and the matrix used to construct the
291: preconditioner are identical, this routine is does nothing.
293: Level: intermediate
295: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
296: @*/
297: PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
298: {
301: pc->useAmat = flg;
302: return(0);
303: }
305: /*@
306: PCSetErrorIfFailure - Causes PC to generate an error if a FPE, for example a zero pivot, is detected.
308: Logically Collective on PC
310: Input Parameters:
311: + pc - iterative context obtained from PCCreate()
312: - flg - PETSC_TRUE indicates you want the error generated
314: Level: advanced
316: Notes:
317: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
318: 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
319: detected.
321: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
323: .keywords: PC, set, initial guess, nonzero
325: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
326: @*/
327: PetscErrorCode PCSetErrorIfFailure(PC pc,PetscBool flg)
328: {
332: pc->erroriffailure = flg;
333: return(0);
334: }
336: /*@
337: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
338: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
339: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
341: Logically Collective on PC
343: Input Parameter:
344: . pc - the preconditioner context
346: Output Parameter:
347: . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
349: Notes:
350: For the common case in which the linear system matrix and the matrix used to construct the
351: preconditioner are identical, this routine is does nothing.
353: Level: intermediate
355: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
356: @*/
357: PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
358: {
361: *flg = pc->useAmat;
362: return(0);
363: }
365: /*@
366: PCCreate - Creates a preconditioner context.
368: Collective on MPI_Comm
370: Input Parameter:
371: . comm - MPI communicator
373: Output Parameter:
374: . pc - location to put the preconditioner context
376: Notes:
377: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or ICC
378: in parallel. For dense matrices it is always PCNONE.
380: Level: developer
382: .keywords: PC, create, context
384: .seealso: PCSetUp(), PCApply(), PCDestroy()
385: @*/
386: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
387: {
388: PC pc;
393: *newpc = 0;
394: PCInitializePackage();
396: PetscHeaderCreate(pc,PC_CLASSID,"PC","Preconditioner","PC",comm,PCDestroy,PCView);
398: pc->mat = 0;
399: pc->pmat = 0;
400: pc->setupcalled = 0;
401: pc->setfromoptionscalled = 0;
402: pc->data = 0;
403: pc->diagonalscale = PETSC_FALSE;
404: pc->diagonalscaleleft = 0;
405: pc->diagonalscaleright = 0;
407: pc->modifysubmatrices = 0;
408: pc->modifysubmatricesP = 0;
410: *newpc = pc;
411: return(0);
413: }
415: /* -------------------------------------------------------------------------------*/
417: /*@
418: PCApply - Applies the preconditioner to a vector.
420: Collective on PC and Vec
422: Input Parameters:
423: + pc - the preconditioner context
424: - x - input vector
426: Output Parameter:
427: . y - output vector
429: Level: developer
431: .keywords: PC, apply
433: .seealso: PCApplyTranspose(), PCApplyBAorAB()
434: @*/
435: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
436: {
438: PetscInt m,n,mv,nv;
444: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
445: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
446: /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
447: MatGetLocalSize(pc->pmat,&m,&n);
448: VecGetLocalSize(x,&nv);
449: VecGetLocalSize(y,&mv);
450: 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);
451: 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);
452: VecLocked(y,3);
454: PCSetUp(pc);
455: if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
456: VecLockPush(x);
457: PetscLogEventBegin(PC_Apply,pc,x,y,0);
458: (*pc->ops->apply)(pc,x,y);
459: PetscLogEventEnd(PC_Apply,pc,x,y,0);
460: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
461: VecLockPop(x);
462: return(0);
463: }
465: /*@
466: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
468: Collective on PC and Vec
470: Input Parameters:
471: + pc - the preconditioner context
472: - x - input vector
474: Output Parameter:
475: . y - output vector
477: Notes:
478: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
480: Level: developer
482: .keywords: PC, apply, symmetric, left
484: .seealso: PCApply(), PCApplySymmetricRight()
485: @*/
486: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
487: {
494: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
495: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
496: PCSetUp(pc);
497: if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
498: VecLockPush(x);
499: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
500: (*pc->ops->applysymmetricleft)(pc,x,y);
501: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
502: VecLockPop(x);
503: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
504: return(0);
505: }
507: /*@
508: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
510: Collective on PC and Vec
512: Input Parameters:
513: + pc - the preconditioner context
514: - x - input vector
516: Output Parameter:
517: . y - output vector
519: Level: developer
521: Notes:
522: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
524: .keywords: PC, apply, symmetric, right
526: .seealso: PCApply(), PCApplySymmetricLeft()
527: @*/
528: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
529: {
536: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
537: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
538: PCSetUp(pc);
539: if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
540: VecLockPush(x);
541: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
542: (*pc->ops->applysymmetricright)(pc,x,y);
543: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
544: VecLockPop(x);
545: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
546: return(0);
547: }
549: /*@
550: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
552: Collective on PC and Vec
554: Input Parameters:
555: + pc - the preconditioner context
556: - x - input vector
558: Output Parameter:
559: . y - output vector
561: Notes: For complex numbers this applies the non-Hermitian transpose.
563: Developer Notes: We need to implement a PCApplyHermitianTranspose()
565: Level: developer
567: .keywords: PC, apply, transpose
569: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
570: @*/
571: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
572: {
579: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
580: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
581: PCSetUp(pc);
582: if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
583: VecLockPush(x);
584: PetscLogEventBegin(PC_Apply,pc,x,y,0);
585: (*pc->ops->applytranspose)(pc,x,y);
586: PetscLogEventEnd(PC_Apply,pc,x,y,0);
587: VecLockPop(x);
588: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
589: return(0);
590: }
592: /*@
593: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
595: Collective on PC and Vec
597: Input Parameters:
598: . pc - the preconditioner context
600: Output Parameter:
601: . flg - PETSC_TRUE if a transpose operation is defined
603: Level: developer
605: .keywords: PC, apply, transpose
607: .seealso: PCApplyTranspose()
608: @*/
609: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
610: {
614: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
615: else *flg = PETSC_FALSE;
616: return(0);
617: }
619: /*@
620: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
622: Collective on PC and Vec
624: Input Parameters:
625: + pc - the preconditioner context
626: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
627: . x - input vector
628: - work - work vector
630: Output Parameter:
631: . y - output vector
633: Level: developer
635: 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
636: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
638: .keywords: PC, apply, operator
640: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
641: @*/
642: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
643: {
651: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
652: 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");
653: if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
654: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
656: PCSetUp(pc);
657: if (pc->diagonalscale) {
658: if (pc->ops->applyBA) {
659: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
660: VecDuplicate(x,&work2);
661: PCDiagonalScaleRight(pc,x,work2);
662: (*pc->ops->applyBA)(pc,side,work2,y,work);
663: PCDiagonalScaleLeft(pc,y,y);
664: VecDestroy(&work2);
665: } else if (side == PC_RIGHT) {
666: PCDiagonalScaleRight(pc,x,y);
667: PCApply(pc,y,work);
668: MatMult(pc->mat,work,y);
669: PCDiagonalScaleLeft(pc,y,y);
670: } else if (side == PC_LEFT) {
671: PCDiagonalScaleRight(pc,x,y);
672: MatMult(pc->mat,y,work);
673: PCApply(pc,work,y);
674: PCDiagonalScaleLeft(pc,y,y);
675: } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
676: } else {
677: if (pc->ops->applyBA) {
678: (*pc->ops->applyBA)(pc,side,x,y,work);
679: } else if (side == PC_RIGHT) {
680: PCApply(pc,x,work);
681: MatMult(pc->mat,work,y);
682: } else if (side == PC_LEFT) {
683: MatMult(pc->mat,x,work);
684: PCApply(pc,work,y);
685: } else if (side == PC_SYMMETRIC) {
686: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
687: PCApplySymmetricRight(pc,x,work);
688: MatMult(pc->mat,work,y);
689: VecCopy(y,work);
690: PCApplySymmetricLeft(pc,work,y);
691: }
692: }
693: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
694: return(0);
695: }
697: /*@
698: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
699: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
700: NOT tr(B*A) = tr(A)*tr(B).
702: Collective on PC and Vec
704: Input Parameters:
705: + pc - the preconditioner context
706: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
707: . x - input vector
708: - work - work vector
710: Output Parameter:
711: . y - output vector
714: 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
715: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
717: Level: developer
719: .keywords: PC, apply, operator, transpose
721: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
722: @*/
723: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
724: {
732: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
733: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
734: if (pc->ops->applyBAtranspose) {
735: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
736: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
737: return(0);
738: }
739: if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
741: PCSetUp(pc);
742: if (side == PC_RIGHT) {
743: PCApplyTranspose(pc,x,work);
744: MatMultTranspose(pc->mat,work,y);
745: } else if (side == PC_LEFT) {
746: MatMultTranspose(pc->mat,x,work);
747: PCApplyTranspose(pc,work,y);
748: }
749: /* add support for PC_SYMMETRIC */
750: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
751: return(0);
752: }
754: /* -------------------------------------------------------------------------------*/
756: /*@
757: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
758: built-in fast application of Richardson's method.
760: Not Collective
762: Input Parameter:
763: . pc - the preconditioner
765: Output Parameter:
766: . exists - PETSC_TRUE or PETSC_FALSE
768: Level: developer
770: .keywords: PC, apply, Richardson, exists
772: .seealso: PCApplyRichardson()
773: @*/
774: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
775: {
779: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
780: else *exists = PETSC_FALSE;
781: return(0);
782: }
784: /*@
785: PCApplyRichardson - Applies several steps of Richardson iteration with
786: the particular preconditioner. This routine is usually used by the
787: Krylov solvers and not the application code directly.
789: Collective on PC
791: Input Parameters:
792: + pc - the preconditioner context
793: . b - the right hand side
794: . w - one work vector
795: . rtol - relative decrease in residual norm convergence criteria
796: . abstol - absolute residual norm convergence criteria
797: . dtol - divergence residual norm increase criteria
798: . its - the number of iterations to apply.
799: - guesszero - if the input x contains nonzero initial guess
801: Output Parameter:
802: + outits - number of iterations actually used (for SOR this always equals its)
803: . reason - the reason the apply terminated
804: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
806: Notes:
807: Most preconditioners do not support this function. Use the command
808: PCApplyRichardsonExists() to determine if one does.
810: Except for the multigrid PC this routine ignores the convergence tolerances
811: and always runs for the number of iterations
813: Level: developer
815: .keywords: PC, apply, Richardson
817: .seealso: PCApplyRichardsonExists()
818: @*/
819: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
820: {
828: if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
829: PCSetUp(pc);
830: if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
831: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
832: return(0);
833: }
835: /*@
836: PCGetSetUpFailedReason - Gets the reason a PCSetUp() failed or 0 if it did not fail
838: Logically Collective on PC
840: Input Parameter:
841: . pc - the preconditioner context
843: Output Parameter:
844: . reason - the reason it failed, currently only -1
846: Level: advanced
848: .keywords: PC, setup
850: .seealso: PCCreate(), PCApply(), PCDestroy()
851: @*/
852: PetscErrorCode PCGetSetUpFailedReason(PC pc,PCFailedReason *reason)
853: {
855: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
856: else *reason = pc->failedreason;
857: return(0);
858: }
861: /*
862: a setupcall of 0 indicates never setup,
863: 1 indicates has been previously setup
864: -1 indicates a PCSetUp() was attempted and failed
865: */
866: /*@
867: PCSetUp - Prepares for the use of a preconditioner.
869: Collective on PC
871: Input Parameter:
872: . pc - the preconditioner context
874: Level: developer
876: .keywords: PC, setup
878: .seealso: PCCreate(), PCApply(), PCDestroy()
879: @*/
880: PetscErrorCode PCSetUp(PC pc)
881: {
882: PetscErrorCode ierr;
883: const char *def;
884: PetscObjectState matstate, matnonzerostate;
888: if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
890: if (pc->setupcalled && pc->reusepreconditioner) {
891: PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
892: return(0);
893: }
895: PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
896: MatGetNonzeroState(pc->pmat,&matnonzerostate);
897: if (!pc->setupcalled) {
898: PetscInfo(pc,"Setting up PC for first time\n");
899: pc->flag = DIFFERENT_NONZERO_PATTERN;
900: } else if (matstate == pc->matstate) {
901: PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
902: return(0);
903: } else {
904: if (matnonzerostate > pc->matnonzerostate) {
905: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
906: pc->flag = DIFFERENT_NONZERO_PATTERN;
907: } else {
908: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
909: pc->flag = SAME_NONZERO_PATTERN;
910: }
911: }
912: pc->matstate = matstate;
913: pc->matnonzerostate = matnonzerostate;
915: if (!((PetscObject)pc)->type_name) {
916: PCGetDefaultType_Private(pc,&def);
917: PCSetType(pc,def);
918: }
920: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
921: MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
922: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
923: if (pc->ops->setup) {
924: (*pc->ops->setup)(pc);
925: }
926: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
927: if (!pc->setupcalled) pc->setupcalled = 1;
928: return(0);
929: }
931: /*@
932: PCSetUpOnBlocks - Sets up the preconditioner for each block in
933: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
934: methods.
936: Collective on PC
938: Input Parameters:
939: . pc - the preconditioner context
941: Level: developer
943: .keywords: PC, setup, blocks
945: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
946: @*/
947: PetscErrorCode PCSetUpOnBlocks(PC pc)
948: {
953: if (!pc->ops->setuponblocks) return(0);
954: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
955: (*pc->ops->setuponblocks)(pc);
956: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
957: return(0);
958: }
960: /*@C
961: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
962: submatrices that arise within certain subdomain-based preconditioners.
963: The basic submatrices are extracted from the preconditioner matrix as
964: usual; the user can then alter these (for example, to set different boundary
965: conditions for each submatrix) before they are used for the local solves.
967: Logically Collective on PC
969: Input Parameters:
970: + pc - the preconditioner context
971: . func - routine for modifying the submatrices
972: - ctx - optional user-defined context (may be null)
974: Calling sequence of func:
975: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
977: . row - an array of index sets that contain the global row numbers
978: that comprise each local submatrix
979: . col - an array of index sets that contain the global column numbers
980: that comprise each local submatrix
981: . submat - array of local submatrices
982: - ctx - optional user-defined context for private data for the
983: user-defined func routine (may be null)
985: Notes:
986: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
987: KSPSolve().
989: A routine set by PCSetModifySubMatrices() is currently called within
990: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
991: preconditioners. All other preconditioners ignore this routine.
993: Level: advanced
995: .keywords: PC, set, modify, submatrices
997: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
998: @*/
999: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1000: {
1003: pc->modifysubmatrices = func;
1004: pc->modifysubmatricesP = ctx;
1005: return(0);
1006: }
1008: /*@C
1009: PCModifySubMatrices - Calls an optional user-defined routine within
1010: certain preconditioners if one has been set with PCSetModifySubMarices().
1012: Collective on PC
1014: Input Parameters:
1015: + pc - the preconditioner context
1016: . nsub - the number of local submatrices
1017: . row - an array of index sets that contain the global row numbers
1018: that comprise each local submatrix
1019: . col - an array of index sets that contain the global column numbers
1020: that comprise each local submatrix
1021: . submat - array of local submatrices
1022: - ctx - optional user-defined context for private data for the
1023: user-defined routine (may be null)
1025: Output Parameter:
1026: . submat - array of local submatrices (the entries of which may
1027: have been modified)
1029: Notes:
1030: The user should NOT generally call this routine, as it will
1031: automatically be called within certain preconditioners (currently
1032: block Jacobi, additive Schwarz) if set.
1034: The basic submatrices are extracted from the preconditioner matrix
1035: as usual; the user can then alter these (for example, to set different
1036: boundary conditions for each submatrix) before they are used for the
1037: local solves.
1039: Level: developer
1041: .keywords: PC, modify, submatrices
1043: .seealso: PCSetModifySubMatrices()
1044: @*/
1045: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1046: {
1051: if (!pc->modifysubmatrices) return(0);
1052: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1053: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1054: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1055: return(0);
1056: }
1058: /*@
1059: PCSetOperators - Sets the matrix associated with the linear system and
1060: a (possibly) different one associated with the preconditioner.
1062: Logically Collective on PC and Mat
1064: Input Parameters:
1065: + pc - the preconditioner context
1066: . Amat - the matrix that defines the linear system
1067: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1069: Notes:
1070: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1072: If you wish to replace either Amat or Pmat but leave the other one untouched then
1073: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1074: on it and then pass it back in in your call to KSPSetOperators().
1076: More Notes about Repeated Solution of Linear Systems:
1077: PETSc does NOT reset the matrix entries of either Amat or Pmat
1078: to zero after a linear solve; the user is completely responsible for
1079: matrix assembly. See the routine MatZeroEntries() if desiring to
1080: zero all elements of a matrix.
1082: Level: intermediate
1084: .keywords: PC, set, operators, matrix, linear system
1086: .seealso: PCGetOperators(), MatZeroEntries()
1087: @*/
1088: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1089: {
1090: PetscErrorCode ierr;
1091: PetscInt m1,n1,m2,n2;
1099: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1100: MatGetLocalSize(Amat,&m1,&n1);
1101: MatGetLocalSize(pc->mat,&m2,&n2);
1102: 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);
1103: MatGetLocalSize(Pmat,&m1,&n1);
1104: MatGetLocalSize(pc->pmat,&m2,&n2);
1105: 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);
1106: }
1108: if (Pmat != pc->pmat) {
1109: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1110: pc->matnonzerostate = -1;
1111: pc->matstate = -1;
1112: }
1114: /* reference first in case the matrices are the same */
1115: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1116: MatDestroy(&pc->mat);
1117: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1118: MatDestroy(&pc->pmat);
1119: pc->mat = Amat;
1120: pc->pmat = Pmat;
1121: return(0);
1122: }
1124: /*@
1125: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1127: Logically Collective on PC
1129: Input Parameters:
1130: + pc - the preconditioner context
1131: - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1133: Level: intermediate
1135: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1136: @*/
1137: PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1138: {
1141: pc->reusepreconditioner = flag;
1142: return(0);
1143: }
1145: /*@
1146: PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1148: Not Collective
1150: Input Parameter:
1151: . pc - the preconditioner context
1153: Output Parameter:
1154: . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1156: Level: intermediate
1158: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1159: @*/
1160: PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1161: {
1164: *flag = pc->reusepreconditioner;
1165: return(0);
1166: }
1168: /*@C
1169: PCGetOperators - Gets the matrix associated with the linear system and
1170: possibly a different one associated with the preconditioner.
1172: Not collective, though parallel Mats are returned if the PC is parallel
1174: Input Parameter:
1175: . pc - the preconditioner context
1177: Output Parameters:
1178: + Amat - the matrix defining the linear system
1179: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1181: Level: intermediate
1183: Notes: Does not increase the reference count of the matrices, so you should not destroy them
1185: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1186: are created in PC and returned to the user. In this case, if both operators
1187: mat and pmat are requested, two DIFFERENT operators will be returned. If
1188: only one is requested both operators in the PC will be the same (i.e. as
1189: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1190: The user must set the sizes of the returned matrices and their type etc just
1191: as if the user created them with MatCreate(). For example,
1193: $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1194: $ set size, type, etc of Amat
1196: $ MatCreate(comm,&mat);
1197: $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1198: $ PetscObjectDereference((PetscObject)mat);
1199: $ set size, type, etc of Amat
1201: and
1203: $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1204: $ set size, type, etc of Amat and Pmat
1206: $ MatCreate(comm,&Amat);
1207: $ MatCreate(comm,&Pmat);
1208: $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1209: $ PetscObjectDereference((PetscObject)Amat);
1210: $ PetscObjectDereference((PetscObject)Pmat);
1211: $ set size, type, etc of Amat and Pmat
1213: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1214: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1215: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1216: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1217: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1218: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1219: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1220: it can be created for you?
1223: .keywords: PC, get, operators, matrix, linear system
1225: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1226: @*/
1227: PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1228: {
1233: if (Amat) {
1234: if (!pc->mat) {
1235: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1236: pc->mat = pc->pmat;
1237: PetscObjectReference((PetscObject)pc->mat);
1238: } else { /* both Amat and Pmat are empty */
1239: MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1240: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1241: pc->pmat = pc->mat;
1242: PetscObjectReference((PetscObject)pc->pmat);
1243: }
1244: }
1245: }
1246: *Amat = pc->mat;
1247: }
1248: if (Pmat) {
1249: if (!pc->pmat) {
1250: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1251: pc->pmat = pc->mat;
1252: PetscObjectReference((PetscObject)pc->pmat);
1253: } else {
1254: MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1255: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1256: pc->mat = pc->pmat;
1257: PetscObjectReference((PetscObject)pc->mat);
1258: }
1259: }
1260: }
1261: *Pmat = pc->pmat;
1262: }
1263: return(0);
1264: }
1266: /*@C
1267: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1268: possibly a different one associated with the preconditioner have been set in the PC.
1270: Not collective, though the results on all processes should be the same
1272: Input Parameter:
1273: . pc - the preconditioner context
1275: Output Parameters:
1276: + mat - the matrix associated with the linear system was set
1277: - pmat - matrix associated with the preconditioner was set, usually the same
1279: Level: intermediate
1281: .keywords: PC, get, operators, matrix, linear system
1283: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1284: @*/
1285: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1286: {
1289: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1290: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1291: return(0);
1292: }
1294: /*@
1295: PCFactorGetMatrix - Gets the factored matrix from the
1296: preconditioner context. This routine is valid only for the LU,
1297: incomplete LU, Cholesky, and incomplete Cholesky methods.
1299: Not Collective on PC though Mat is parallel if PC is parallel
1301: Input Parameters:
1302: . pc - the preconditioner context
1304: Output parameters:
1305: . mat - the factored matrix
1307: Level: advanced
1309: Notes: Does not increase the reference count for the matrix so DO NOT destroy it
1311: .keywords: PC, get, factored, matrix
1312: @*/
1313: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1314: {
1320: if (pc->ops->getfactoredmatrix) {
1321: (*pc->ops->getfactoredmatrix)(pc,mat);
1322: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1323: return(0);
1324: }
1326: /*@C
1327: PCSetOptionsPrefix - Sets the prefix used for searching for all
1328: PC options in the database.
1330: Logically Collective on PC
1332: Input Parameters:
1333: + pc - the preconditioner context
1334: - prefix - the prefix string to prepend to all PC option requests
1336: Notes:
1337: A hyphen (-) must NOT be given at the beginning of the prefix name.
1338: The first character of all runtime options is AUTOMATICALLY the
1339: hyphen.
1341: Level: advanced
1343: .keywords: PC, set, options, prefix, database
1345: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1346: @*/
1347: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1348: {
1353: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1354: return(0);
1355: }
1357: /*@C
1358: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1359: PC options in the database.
1361: Logically Collective on PC
1363: Input Parameters:
1364: + pc - the preconditioner context
1365: - prefix - the prefix string to prepend to all PC option requests
1367: Notes:
1368: A hyphen (-) must NOT be given at the beginning of the prefix name.
1369: The first character of all runtime options is AUTOMATICALLY the
1370: hyphen.
1372: Level: advanced
1374: .keywords: PC, append, options, prefix, database
1376: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1377: @*/
1378: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1379: {
1384: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1385: return(0);
1386: }
1388: /*@C
1389: PCGetOptionsPrefix - Gets the prefix used for searching for all
1390: PC options in the database.
1392: Not Collective
1394: Input Parameters:
1395: . pc - the preconditioner context
1397: Output Parameters:
1398: . prefix - pointer to the prefix string used, is returned
1400: Notes: On the fortran side, the user should pass in a string 'prifix' of
1401: sufficient length to hold the prefix.
1403: Level: advanced
1405: .keywords: PC, get, options, prefix, database
1407: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1408: @*/
1409: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1410: {
1416: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1417: return(0);
1418: }
1420: /*
1421: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1422: preconditioners including BDDC and Eisentat that transform the equations before applying
1423: the Krylov methods
1424: */
1425: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc,PetscBool *change)
1426: {
1432: *change = PETSC_FALSE;
1433: PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1434: return(0);
1435: }
1437: /*@
1438: PCPreSolve - Optional pre-solve phase, intended for any
1439: preconditioner-specific actions that must be performed before
1440: the iterative solve itself.
1442: Collective on PC
1444: Input Parameters:
1445: + pc - the preconditioner context
1446: - ksp - the Krylov subspace context
1448: Level: developer
1450: Sample of Usage:
1451: .vb
1452: PCPreSolve(pc,ksp);
1453: KSPSolve(ksp,b,x);
1454: PCPostSolve(pc,ksp);
1455: .ve
1457: Notes:
1458: The pre-solve phase is distinct from the PCSetUp() phase.
1460: KSPSolve() calls this directly, so is rarely called by the user.
1462: .keywords: PC, pre-solve
1464: .seealso: PCPostSolve()
1465: @*/
1466: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1467: {
1469: Vec x,rhs;
1474: pc->presolvedone++;
1475: if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1476: KSPGetSolution(ksp,&x);
1477: KSPGetRhs(ksp,&rhs);
1479: if (pc->ops->presolve) {
1480: (*pc->ops->presolve)(pc,ksp,rhs,x);
1481: }
1482: return(0);
1483: }
1485: /*@
1486: PCPostSolve - Optional post-solve phase, intended for any
1487: preconditioner-specific actions that must be performed after
1488: the iterative solve itself.
1490: Collective on PC
1492: Input Parameters:
1493: + pc - the preconditioner context
1494: - ksp - the Krylov subspace context
1496: Sample of Usage:
1497: .vb
1498: PCPreSolve(pc,ksp);
1499: KSPSolve(ksp,b,x);
1500: PCPostSolve(pc,ksp);
1501: .ve
1503: Note:
1504: KSPSolve() calls this routine directly, so it is rarely called by the user.
1506: Level: developer
1508: .keywords: PC, post-solve
1510: .seealso: PCPreSolve(), KSPSolve()
1511: @*/
1512: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1513: {
1515: Vec x,rhs;
1520: pc->presolvedone--;
1521: KSPGetSolution(ksp,&x);
1522: KSPGetRhs(ksp,&rhs);
1523: if (pc->ops->postsolve) {
1524: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1525: }
1526: return(0);
1527: }
1529: /*@C
1530: PCLoad - Loads a PC that has been stored in binary with PCView().
1532: Collective on PetscViewer
1534: Input Parameters:
1535: + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1536: some related function before a call to PCLoad().
1537: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1539: Level: intermediate
1541: Notes:
1542: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1544: Notes for advanced users:
1545: Most users should not need to know the details of the binary storage
1546: format, since PCLoad() and PCView() completely hide these details.
1547: But for anyone who's interested, the standard binary matrix storage
1548: format is
1549: .vb
1550: has not yet been determined
1551: .ve
1553: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1554: @*/
1555: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1556: {
1558: PetscBool isbinary;
1559: PetscInt classid;
1560: char type[256];
1565: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1566: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1568: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1569: if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1570: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1571: PCSetType(newdm, type);
1572: if (newdm->ops->load) {
1573: (*newdm->ops->load)(newdm,viewer);
1574: }
1575: return(0);
1576: }
1578: #include <petscdraw.h>
1579: #if defined(PETSC_HAVE_SAWS)
1580: #include <petscviewersaws.h>
1581: #endif
1582: /*@C
1583: PCView - Prints the PC data structure.
1585: Collective on PC
1587: Input Parameters:
1588: + PC - the PC context
1589: - viewer - optional visualization context
1591: Note:
1592: The available visualization contexts include
1593: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1594: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1595: output where only the first processor opens
1596: the file. All other processors send their
1597: data to the first processor to print.
1599: The user can open an alternative visualization contexts with
1600: PetscViewerASCIIOpen() (output to a specified file).
1602: Level: developer
1604: .keywords: PC, view
1606: .seealso: KSPView(), PetscViewerASCIIOpen()
1607: @*/
1608: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1609: {
1610: PCType cstr;
1611: PetscErrorCode ierr;
1612: PetscBool iascii,isstring,isbinary,isdraw;
1613: PetscViewerFormat format;
1614: #if defined(PETSC_HAVE_SAWS)
1615: PetscBool issaws;
1616: #endif
1620: if (!viewer) {
1621: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1622: }
1626: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1627: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1628: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1629: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1630: #if defined(PETSC_HAVE_SAWS)
1631: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1632: #endif
1634: if (iascii) {
1635: PetscViewerGetFormat(viewer,&format);
1636: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1637: if (!pc->setupcalled) {
1638: PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");
1639: }
1640: if (pc->ops->view) {
1641: PetscViewerASCIIPushTab(viewer);
1642: (*pc->ops->view)(pc,viewer);
1643: PetscViewerASCIIPopTab(viewer);
1644: }
1645: if (pc->mat) {
1646: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1647: if (pc->pmat == pc->mat) {
1648: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1649: PetscViewerASCIIPushTab(viewer);
1650: MatView(pc->mat,viewer);
1651: PetscViewerASCIIPopTab(viewer);
1652: } else {
1653: if (pc->pmat) {
1654: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1655: } else {
1656: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1657: }
1658: PetscViewerASCIIPushTab(viewer);
1659: MatView(pc->mat,viewer);
1660: if (pc->pmat) {MatView(pc->pmat,viewer);}
1661: PetscViewerASCIIPopTab(viewer);
1662: }
1663: PetscViewerPopFormat(viewer);
1664: }
1665: } else if (isstring) {
1666: PCGetType(pc,&cstr);
1667: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1668: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1669: } else if (isbinary) {
1670: PetscInt classid = PC_FILE_CLASSID;
1671: MPI_Comm comm;
1672: PetscMPIInt rank;
1673: char type[256];
1675: PetscObjectGetComm((PetscObject)pc,&comm);
1676: MPI_Comm_rank(comm,&rank);
1677: if (!rank) {
1678: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1679: PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1680: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1681: }
1682: if (pc->ops->view) {
1683: (*pc->ops->view)(pc,viewer);
1684: }
1685: } else if (isdraw) {
1686: PetscDraw draw;
1687: char str[25];
1688: PetscReal x,y,bottom,h;
1689: PetscInt n;
1691: PetscViewerDrawGetDraw(viewer,0,&draw);
1692: PetscDrawGetCurrentPoint(draw,&x,&y);
1693: if (pc->mat) {
1694: MatGetSize(pc->mat,&n,NULL);
1695: PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1696: } else {
1697: PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1698: }
1699: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1700: bottom = y - h;
1701: PetscDrawPushCurrentPoint(draw,x,bottom);
1702: if (pc->ops->view) {
1703: (*pc->ops->view)(pc,viewer);
1704: }
1705: PetscDrawPopCurrentPoint(draw);
1706: #if defined(PETSC_HAVE_SAWS)
1707: } else if (issaws) {
1708: PetscMPIInt rank;
1710: PetscObjectName((PetscObject)pc);
1711: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1712: if (!((PetscObject)pc)->amsmem && !rank) {
1713: PetscObjectViewSAWs((PetscObject)pc,viewer);
1714: }
1715: if (pc->mat) {MatView(pc->mat,viewer);}
1716: if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1717: #endif
1718: }
1719: return(0);
1720: }
1722: /*@C
1723: PCRegister - Adds a method to the preconditioner package.
1725: Not collective
1727: Input Parameters:
1728: + name_solver - name of a new user-defined solver
1729: - routine_create - routine to create method context
1731: Notes:
1732: PCRegister() may be called multiple times to add several user-defined preconditioners.
1734: Sample usage:
1735: .vb
1736: PCRegister("my_solver", MySolverCreate);
1737: .ve
1739: Then, your solver can be chosen with the procedural interface via
1740: $ PCSetType(pc,"my_solver")
1741: or at runtime via the option
1742: $ -pc_type my_solver
1744: Level: advanced
1746: .keywords: PC, register
1748: .seealso: PCRegisterAll(), PCRegisterDestroy()
1749: @*/
1750: PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1751: {
1755: PetscFunctionListAdd(&PCList,sname,function);
1756: return(0);
1757: }
1759: /*@
1760: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1762: Collective on PC
1764: Input Parameter:
1765: . pc - the preconditioner object
1767: Output Parameter:
1768: . mat - the explict preconditioned operator
1770: Notes:
1771: This computation is done by applying the operators to columns of the
1772: identity matrix.
1774: Currently, this routine uses a dense matrix format when 1 processor
1775: is used and a sparse format otherwise. This routine is costly in general,
1776: and is recommended for use only with relatively small systems.
1778: Level: advanced
1780: .keywords: PC, compute, explicit, operator
1782: .seealso: KSPComputeExplicitOperator()
1784: @*/
1785: PetscErrorCode PCComputeExplicitOperator(PC pc,Mat *mat)
1786: {
1787: Vec in,out;
1789: PetscInt i,M,m,*rows,start,end;
1790: PetscMPIInt size;
1791: MPI_Comm comm;
1792: PetscScalar *array,one = 1.0;
1798: PetscObjectGetComm((PetscObject)pc,&comm);
1799: MPI_Comm_size(comm,&size);
1801: if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1802: MatCreateVecs(pc->pmat,&in,0);
1803: VecDuplicate(in,&out);
1804: VecGetOwnershipRange(in,&start,&end);
1805: VecGetSize(in,&M);
1806: VecGetLocalSize(in,&m);
1807: PetscMalloc1(m+1,&rows);
1808: for (i=0; i<m; i++) rows[i] = start + i;
1810: MatCreate(comm,mat);
1811: MatSetSizes(*mat,m,m,M,M);
1812: if (size == 1) {
1813: MatSetType(*mat,MATSEQDENSE);
1814: MatSeqDenseSetPreallocation(*mat,NULL);
1815: } else {
1816: MatSetType(*mat,MATMPIAIJ);
1817: MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1818: }
1819: MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
1821: for (i=0; i<M; i++) {
1823: VecSet(in,0.0);
1824: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1825: VecAssemblyBegin(in);
1826: VecAssemblyEnd(in);
1828: /* should fix, allowing user to choose side */
1829: PCApply(pc,in,out);
1831: VecGetArray(out,&array);
1832: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1833: VecRestoreArray(out,&array);
1835: }
1836: PetscFree(rows);
1837: VecDestroy(&out);
1838: VecDestroy(&in);
1839: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1840: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1841: return(0);
1842: }
1844: /*@
1845: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1847: Collective on PC
1849: Input Parameters:
1850: + pc - the solver context
1851: . dim - the dimension of the coordinates 1, 2, or 3
1852: - coords - the coordinates
1854: Level: intermediate
1856: Notes: coords is an array of the 3D coordinates for the nodes on
1857: the local processor. So if there are 108 equation on a processor
1858: for a displacement finite element discretization of elasticity (so
1859: that there are 36 = 108/3 nodes) then the array must have 108
1860: double precision values (ie, 3 * 36). These x y z coordinates
1861: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1862: ... , N-1.z ].
1864: .seealso: MatSetNearNullSpace
1865: @*/
1866: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1867: {
1871: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1872: return(0);
1873: }