Actual source code: precon.c
petsc-3.7.5 2017-01-01
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;
12: PetscInt PetscMGLevelId;
16: PetscErrorCode PCGetDefaultType_Private(PC pc,const char *type[])
17: {
19: PetscMPIInt size;
20: PetscBool flg1,flg2,set,flg3;
23: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
24: if (pc->pmat) {
25: PetscErrorCode (*f)(Mat,MatReuse,Mat*);
26: PetscObjectQueryFunction((PetscObject)pc->pmat,"MatGetDiagonalBlock_C",&f);
27: if (size == 1) {
28: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
29: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
30: MatIsSymmetricKnown(pc->pmat,&set,&flg3);
31: if (flg1 && (!flg2 || (set && flg3))) {
32: *type = PCICC;
33: } else if (flg2) {
34: *type = PCILU;
35: } else if (f) { /* likely is a parallel matrix run on one processor */
36: *type = PCBJACOBI;
37: } else {
38: *type = PCNONE;
39: }
40: } else {
41: if (f) {
42: *type = PCBJACOBI;
43: } else {
44: *type = PCNONE;
45: }
46: }
47: } else {
48: if (size == 1) {
49: *type = PCILU;
50: } else {
51: *type = PCBJACOBI;
52: }
53: }
54: return(0);
55: }
59: /*@
60: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
62: Collective on PC
64: Input Parameter:
65: . pc - the preconditioner context
67: Level: developer
69: 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
71: .keywords: PC, destroy
73: .seealso: PCCreate(), PCSetUp()
74: @*/
75: PetscErrorCode PCReset(PC pc)
76: {
81: if (pc->ops->reset) {
82: (*pc->ops->reset)(pc);
83: }
84: VecDestroy(&pc->diagonalscaleright);
85: VecDestroy(&pc->diagonalscaleleft);
86: MatDestroy(&pc->pmat);
87: MatDestroy(&pc->mat);
89: pc->setupcalled = 0;
90: return(0);
91: }
95: /*@
96: PCDestroy - Destroys PC context that was created with PCCreate().
98: Collective on PC
100: Input Parameter:
101: . pc - the preconditioner context
103: Level: developer
105: .keywords: PC, destroy
107: .seealso: PCCreate(), PCSetUp()
108: @*/
109: PetscErrorCode PCDestroy(PC *pc)
110: {
114: if (!*pc) return(0);
116: if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; return(0);}
118: PCReset(*pc);
120: /* if memory was published with SAWs then destroy it */
121: PetscObjectSAWsViewOff((PetscObject)*pc);
122: if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
123: DMDestroy(&(*pc)->dm);
124: PetscHeaderDestroy(pc);
125: return(0);
126: }
130: /*@C
131: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
132: scaling as needed by certain time-stepping codes.
134: Logically Collective on PC
136: Input Parameter:
137: . pc - the preconditioner context
139: Output Parameter:
140: . flag - PETSC_TRUE if it applies the scaling
142: Level: developer
144: Notes: If this returns PETSC_TRUE then the system solved via the Krylov method is
145: $ D M A D^{-1} y = D M b for left preconditioning or
146: $ D A M D^{-1} z = D b for right preconditioning
148: .keywords: PC
150: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
151: @*/
152: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
153: {
157: *flag = pc->diagonalscale;
158: return(0);
159: }
163: /*@
164: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
165: scaling as needed by certain time-stepping codes.
167: Logically Collective on PC
169: Input Parameters:
170: + pc - the preconditioner context
171: - s - scaling vector
173: Level: intermediate
175: Notes: The system solved via the Krylov method is
176: $ D M A D^{-1} y = D M b for left preconditioning or
177: $ D A M D^{-1} z = D b for right preconditioning
179: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
181: .keywords: PC
183: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
184: @*/
185: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
186: {
192: pc->diagonalscale = PETSC_TRUE;
194: PetscObjectReference((PetscObject)s);
195: VecDestroy(&pc->diagonalscaleleft);
197: pc->diagonalscaleleft = s;
199: VecDuplicate(s,&pc->diagonalscaleright);
200: VecCopy(s,pc->diagonalscaleright);
201: VecReciprocal(pc->diagonalscaleright);
202: return(0);
203: }
207: /*@
208: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
210: Logically Collective on PC
212: Input Parameters:
213: + pc - the preconditioner context
214: . in - input vector
215: + out - scaled vector (maybe the same as in)
217: Level: intermediate
219: Notes: The system solved via the Krylov method is
220: $ D M A D^{-1} y = D M b for left preconditioning or
221: $ D A M D^{-1} z = D b for right preconditioning
223: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
225: If diagonal scaling is turned off and in is not out then in is copied to out
227: .keywords: PC
229: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
230: @*/
231: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
232: {
239: if (pc->diagonalscale) {
240: VecPointwiseMult(out,pc->diagonalscaleleft,in);
241: } else if (in != out) {
242: VecCopy(in,out);
243: }
244: return(0);
245: }
249: /*@
250: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
252: Logically Collective on PC
254: Input Parameters:
255: + pc - the preconditioner context
256: . in - input vector
257: + out - scaled vector (maybe the same as in)
259: Level: intermediate
261: Notes: The system solved via the Krylov method is
262: $ D M A D^{-1} y = D M b for left preconditioning or
263: $ D A M D^{-1} z = D b for right preconditioning
265: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
267: If diagonal scaling is turned off and in is not out then in is copied to out
269: .keywords: PC
271: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
272: @*/
273: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
274: {
281: if (pc->diagonalscale) {
282: VecPointwiseMult(out,pc->diagonalscaleright,in);
283: } else if (in != out) {
284: VecCopy(in,out);
285: }
286: return(0);
287: }
291: /*@
292: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
293: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
294: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
296: Logically Collective on PC
298: Input Parameters:
299: + pc - the preconditioner context
300: - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
302: Options Database Key:
303: . -pc_use_amat <true,false>
305: Notes:
306: For the common case in which the linear system matrix and the matrix used to construct the
307: preconditioner are identical, this routine is does nothing.
309: Level: intermediate
311: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
312: @*/
313: PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
314: {
317: pc->useAmat = flg;
318: return(0);
319: }
323: /*@
324: PCSetErrorIfFailure - Causes PC to generate an error if a FPE, for example a zero pivot, is detected.
326: Logically Collective on PC
328: Input Parameters:
329: + pc - iterative context obtained from PCCreate()
330: - flg - PETSC_TRUE indicates you want the error generated
332: Level: advanced
334: Notes:
335: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
336: 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
337: detected.
339: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
341: .keywords: PC, set, initial guess, nonzero
343: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
344: @*/
345: PetscErrorCode PCSetErrorIfFailure(PC pc,PetscBool flg)
346: {
350: pc->erroriffailure = flg;
351: return(0);
352: }
356: /*@
357: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
358: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
359: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
361: Logically Collective on PC
363: Input Parameter:
364: . pc - the preconditioner context
366: Output Parameter:
367: . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
369: Notes:
370: For the common case in which the linear system matrix and the matrix used to construct the
371: preconditioner are identical, this routine is does nothing.
373: Level: intermediate
375: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
376: @*/
377: PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
378: {
381: *flg = pc->useAmat;
382: return(0);
383: }
387: /*@
388: PCCreate - Creates a preconditioner context.
390: Collective on MPI_Comm
392: Input Parameter:
393: . comm - MPI communicator
395: Output Parameter:
396: . pc - location to put the preconditioner context
398: Notes:
399: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or ICC
400: in parallel. For dense matrices it is always PCNONE.
402: Level: developer
404: .keywords: PC, create, context
406: .seealso: PCSetUp(), PCApply(), PCDestroy()
407: @*/
408: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
409: {
410: PC pc;
415: *newpc = 0;
416: PCInitializePackage();
418: PetscHeaderCreate(pc,PC_CLASSID,"PC","Preconditioner","PC",comm,PCDestroy,PCView);
420: pc->mat = 0;
421: pc->pmat = 0;
422: pc->setupcalled = 0;
423: pc->setfromoptionscalled = 0;
424: pc->data = 0;
425: pc->diagonalscale = PETSC_FALSE;
426: pc->diagonalscaleleft = 0;
427: pc->diagonalscaleright = 0;
429: pc->modifysubmatrices = 0;
430: pc->modifysubmatricesP = 0;
432: *newpc = pc;
433: return(0);
435: }
437: /* -------------------------------------------------------------------------------*/
441: /*@
442: PCApply - Applies the preconditioner to a vector.
444: Collective on PC and Vec
446: Input Parameters:
447: + pc - the preconditioner context
448: - x - input vector
450: Output Parameter:
451: . y - output vector
453: Level: developer
455: .keywords: PC, apply
457: .seealso: PCApplyTranspose(), PCApplyBAorAB()
458: @*/
459: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
460: {
462: PetscInt m,n,mv,nv;
468: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
469: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
470: /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
471: MatGetLocalSize(pc->pmat,&m,&n);
472: VecGetLocalSize(x,&nv);
473: VecGetLocalSize(y,&mv);
474: 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);
475: 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);
476: VecLocked(y,3);
478: PCSetUp(pc);
479: if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
480: VecLockPush(x);
481: PetscLogEventBegin(PC_Apply,pc,x,y,0);
482: (*pc->ops->apply)(pc,x,y);
483: PetscLogEventEnd(PC_Apply,pc,x,y,0);
484: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
485: VecLockPop(x);
486: return(0);
487: }
491: /*@
492: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
494: Collective on PC and Vec
496: Input Parameters:
497: + pc - the preconditioner context
498: - x - input vector
500: Output Parameter:
501: . y - output vector
503: Notes:
504: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
506: Level: developer
508: .keywords: PC, apply, symmetric, left
510: .seealso: PCApply(), PCApplySymmetricRight()
511: @*/
512: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
513: {
520: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
521: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
522: PCSetUp(pc);
523: if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
524: VecLockPush(x);
525: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
526: (*pc->ops->applysymmetricleft)(pc,x,y);
527: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
528: VecLockPop(x);
529: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
530: return(0);
531: }
535: /*@
536: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
538: Collective on PC and Vec
540: Input Parameters:
541: + pc - the preconditioner context
542: - x - input vector
544: Output Parameter:
545: . y - output vector
547: Level: developer
549: Notes:
550: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
552: .keywords: PC, apply, symmetric, right
554: .seealso: PCApply(), PCApplySymmetricLeft()
555: @*/
556: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
557: {
564: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
565: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
566: PCSetUp(pc);
567: if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
568: VecLockPush(x);
569: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
570: (*pc->ops->applysymmetricright)(pc,x,y);
571: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
572: VecLockPop(x);
573: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
574: return(0);
575: }
579: /*@
580: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
582: Collective on PC and Vec
584: Input Parameters:
585: + pc - the preconditioner context
586: - x - input vector
588: Output Parameter:
589: . y - output vector
591: Notes: For complex numbers this applies the non-Hermitian transpose.
593: Developer Notes: We need to implement a PCApplyHermitianTranspose()
595: Level: developer
597: .keywords: PC, apply, transpose
599: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
600: @*/
601: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
602: {
609: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
610: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
611: PCSetUp(pc);
612: if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
613: VecLockPush(x);
614: PetscLogEventBegin(PC_Apply,pc,x,y,0);
615: (*pc->ops->applytranspose)(pc,x,y);
616: PetscLogEventEnd(PC_Apply,pc,x,y,0);
617: VecLockPop(x);
618: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
619: return(0);
620: }
624: /*@
625: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
627: Collective on PC and Vec
629: Input Parameters:
630: . pc - the preconditioner context
632: Output Parameter:
633: . flg - PETSC_TRUE if a transpose operation is defined
635: Level: developer
637: .keywords: PC, apply, transpose
639: .seealso: PCApplyTranspose()
640: @*/
641: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
642: {
646: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
647: else *flg = PETSC_FALSE;
648: return(0);
649: }
653: /*@
654: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
656: Collective on PC and Vec
658: Input Parameters:
659: + pc - the preconditioner context
660: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
661: . x - input vector
662: - work - work vector
664: Output Parameter:
665: . y - output vector
667: Level: developer
669: 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
670: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
672: .keywords: PC, apply, operator
674: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
675: @*/
676: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
677: {
685: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
686: 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");
687: if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
688: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
690: PCSetUp(pc);
691: if (pc->diagonalscale) {
692: if (pc->ops->applyBA) {
693: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
694: VecDuplicate(x,&work2);
695: PCDiagonalScaleRight(pc,x,work2);
696: (*pc->ops->applyBA)(pc,side,work2,y,work);
697: PCDiagonalScaleLeft(pc,y,y);
698: VecDestroy(&work2);
699: } else if (side == PC_RIGHT) {
700: PCDiagonalScaleRight(pc,x,y);
701: PCApply(pc,y,work);
702: MatMult(pc->mat,work,y);
703: PCDiagonalScaleLeft(pc,y,y);
704: } else if (side == PC_LEFT) {
705: PCDiagonalScaleRight(pc,x,y);
706: MatMult(pc->mat,y,work);
707: PCApply(pc,work,y);
708: PCDiagonalScaleLeft(pc,y,y);
709: } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
710: } else {
711: if (pc->ops->applyBA) {
712: (*pc->ops->applyBA)(pc,side,x,y,work);
713: } else if (side == PC_RIGHT) {
714: PCApply(pc,x,work);
715: MatMult(pc->mat,work,y);
716: } else if (side == PC_LEFT) {
717: MatMult(pc->mat,x,work);
718: PCApply(pc,work,y);
719: } else if (side == PC_SYMMETRIC) {
720: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
721: PCApplySymmetricRight(pc,x,work);
722: MatMult(pc->mat,work,y);
723: VecCopy(y,work);
724: PCApplySymmetricLeft(pc,work,y);
725: }
726: }
727: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
728: return(0);
729: }
733: /*@
734: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
735: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
736: NOT tr(B*A) = tr(A)*tr(B).
738: Collective on PC and Vec
740: Input Parameters:
741: + pc - the preconditioner context
742: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
743: . x - input vector
744: - work - work vector
746: Output Parameter:
747: . y - output vector
750: 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
751: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
753: Level: developer
755: .keywords: PC, apply, operator, transpose
757: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
758: @*/
759: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
760: {
768: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
769: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
770: if (pc->ops->applyBAtranspose) {
771: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
772: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
773: return(0);
774: }
775: if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
777: PCSetUp(pc);
778: if (side == PC_RIGHT) {
779: PCApplyTranspose(pc,x,work);
780: MatMultTranspose(pc->mat,work,y);
781: } else if (side == PC_LEFT) {
782: MatMultTranspose(pc->mat,x,work);
783: PCApplyTranspose(pc,work,y);
784: }
785: /* add support for PC_SYMMETRIC */
786: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
787: return(0);
788: }
790: /* -------------------------------------------------------------------------------*/
794: /*@
795: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
796: built-in fast application of Richardson's method.
798: Not Collective
800: Input Parameter:
801: . pc - the preconditioner
803: Output Parameter:
804: . exists - PETSC_TRUE or PETSC_FALSE
806: Level: developer
808: .keywords: PC, apply, Richardson, exists
810: .seealso: PCApplyRichardson()
811: @*/
812: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
813: {
817: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
818: else *exists = PETSC_FALSE;
819: return(0);
820: }
824: /*@
825: PCApplyRichardson - Applies several steps of Richardson iteration with
826: the particular preconditioner. This routine is usually used by the
827: Krylov solvers and not the application code directly.
829: Collective on PC
831: Input Parameters:
832: + pc - the preconditioner context
833: . b - the right hand side
834: . w - one work vector
835: . rtol - relative decrease in residual norm convergence criteria
836: . abstol - absolute residual norm convergence criteria
837: . dtol - divergence residual norm increase criteria
838: . its - the number of iterations to apply.
839: - guesszero - if the input x contains nonzero initial guess
841: Output Parameter:
842: + outits - number of iterations actually used (for SOR this always equals its)
843: . reason - the reason the apply terminated
844: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
846: Notes:
847: Most preconditioners do not support this function. Use the command
848: PCApplyRichardsonExists() to determine if one does.
850: Except for the multigrid PC this routine ignores the convergence tolerances
851: and always runs for the number of iterations
853: Level: developer
855: .keywords: PC, apply, Richardson
857: .seealso: PCApplyRichardsonExists()
858: @*/
859: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
860: {
868: if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
869: PCSetUp(pc);
870: if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
871: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
872: return(0);
873: }
877: /*@
878: PCGetSetUpFailedReason - Gets the reason a PCSetUp() failed or 0 if it did not fail
880: Logically Collective on PC
882: Input Parameter:
883: . pc - the preconditioner context
885: Output Parameter:
886: . reason - the reason it failed, currently only -1
888: Level: advanced
890: .keywords: PC, setup
892: .seealso: PCCreate(), PCApply(), PCDestroy()
893: @*/
894: PetscErrorCode PCGetSetUpFailedReason(PC pc,PCFailedReason *reason)
895: {
897: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
898: else *reason = pc->failedreason;
899: return(0);
900: }
903: /*
904: a setupcall of 0 indicates never setup,
905: 1 indicates has been previously setup
906: -1 indicates a PCSetUp() was attempted and failed
907: */
910: /*@
911: PCSetUp - Prepares for the use of a preconditioner.
913: Collective on PC
915: Input Parameter:
916: . pc - the preconditioner context
918: Level: developer
920: .keywords: PC, setup
922: .seealso: PCCreate(), PCApply(), PCDestroy()
923: @*/
924: PetscErrorCode PCSetUp(PC pc)
925: {
926: PetscErrorCode ierr;
927: const char *def;
928: PetscObjectState matstate, matnonzerostate;
932: if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
934: if (pc->setupcalled && pc->reusepreconditioner) {
935: PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
936: return(0);
937: }
939: PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
940: MatGetNonzeroState(pc->pmat,&matnonzerostate);
941: if (!pc->setupcalled) {
942: PetscInfo(pc,"Setting up PC for first time\n");
943: pc->flag = DIFFERENT_NONZERO_PATTERN;
944: } else if (matstate == pc->matstate) {
945: PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
946: return(0);
947: } else {
948: if (matnonzerostate > pc->matnonzerostate) {
949: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
950: pc->flag = DIFFERENT_NONZERO_PATTERN;
951: } else {
952: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
953: pc->flag = SAME_NONZERO_PATTERN;
954: }
955: }
956: pc->matstate = matstate;
957: pc->matnonzerostate = matnonzerostate;
959: if (!((PetscObject)pc)->type_name) {
960: PCGetDefaultType_Private(pc,&def);
961: PCSetType(pc,def);
962: }
964: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
965: MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
966: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
967: if (pc->ops->setup) {
968: (*pc->ops->setup)(pc);
969: }
970: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
971: if (!pc->setupcalled) pc->setupcalled = 1;
972: return(0);
973: }
977: /*@
978: PCSetUpOnBlocks - Sets up the preconditioner for each block in
979: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
980: methods.
982: Collective on PC
984: Input Parameters:
985: . pc - the preconditioner context
987: Level: developer
989: .keywords: PC, setup, blocks
991: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
992: @*/
993: PetscErrorCode PCSetUpOnBlocks(PC pc)
994: {
999: if (!pc->ops->setuponblocks) return(0);
1000: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
1001: (*pc->ops->setuponblocks)(pc);
1002: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
1003: return(0);
1004: }
1008: /*@C
1009: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
1010: submatrices that arise within certain subdomain-based preconditioners.
1011: The basic submatrices are extracted from the preconditioner matrix as
1012: usual; the user can then alter these (for example, to set different boundary
1013: conditions for each submatrix) before they are used for the local solves.
1015: Logically Collective on PC
1017: Input Parameters:
1018: + pc - the preconditioner context
1019: . func - routine for modifying the submatrices
1020: - ctx - optional user-defined context (may be null)
1022: Calling sequence of func:
1023: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
1025: . row - an array of index sets that contain the global row numbers
1026: that comprise each local submatrix
1027: . col - an array of index sets that contain the global column numbers
1028: that comprise each local submatrix
1029: . submat - array of local submatrices
1030: - ctx - optional user-defined context for private data for the
1031: user-defined func routine (may be null)
1033: Notes:
1034: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
1035: KSPSolve().
1037: A routine set by PCSetModifySubMatrices() is currently called within
1038: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
1039: preconditioners. All other preconditioners ignore this routine.
1041: Level: advanced
1043: .keywords: PC, set, modify, submatrices
1045: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
1046: @*/
1047: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
1048: {
1051: pc->modifysubmatrices = func;
1052: pc->modifysubmatricesP = ctx;
1053: return(0);
1054: }
1058: /*@C
1059: PCModifySubMatrices - Calls an optional user-defined routine within
1060: certain preconditioners if one has been set with PCSetModifySubMarices().
1062: Collective on PC
1064: Input Parameters:
1065: + pc - the preconditioner context
1066: . nsub - the number of local submatrices
1067: . row - an array of index sets that contain the global row numbers
1068: that comprise each local submatrix
1069: . col - an array of index sets that contain the global column numbers
1070: that comprise each local submatrix
1071: . submat - array of local submatrices
1072: - ctx - optional user-defined context for private data for the
1073: user-defined routine (may be null)
1075: Output Parameter:
1076: . submat - array of local submatrices (the entries of which may
1077: have been modified)
1079: Notes:
1080: The user should NOT generally call this routine, as it will
1081: automatically be called within certain preconditioners (currently
1082: block Jacobi, additive Schwarz) if set.
1084: The basic submatrices are extracted from the preconditioner matrix
1085: as usual; the user can then alter these (for example, to set different
1086: boundary conditions for each submatrix) before they are used for the
1087: local solves.
1089: Level: developer
1091: .keywords: PC, modify, submatrices
1093: .seealso: PCSetModifySubMatrices()
1094: @*/
1095: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1096: {
1101: if (!pc->modifysubmatrices) return(0);
1102: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1103: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1104: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1105: return(0);
1106: }
1110: /*@
1111: PCSetOperators - Sets the matrix associated with the linear system and
1112: a (possibly) different one associated with the preconditioner.
1114: Logically Collective on PC and Mat
1116: Input Parameters:
1117: + pc - the preconditioner context
1118: . Amat - the matrix that defines the linear system
1119: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1121: Notes:
1122: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1124: If you wish to replace either Amat or Pmat but leave the other one untouched then
1125: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1126: on it and then pass it back in in your call to KSPSetOperators().
1128: More Notes about Repeated Solution of Linear Systems:
1129: PETSc does NOT reset the matrix entries of either Amat or Pmat
1130: to zero after a linear solve; the user is completely responsible for
1131: matrix assembly. See the routine MatZeroEntries() if desiring to
1132: zero all elements of a matrix.
1134: Level: intermediate
1136: .keywords: PC, set, operators, matrix, linear system
1138: .seealso: PCGetOperators(), MatZeroEntries()
1139: @*/
1140: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1141: {
1142: PetscErrorCode ierr;
1143: PetscInt m1,n1,m2,n2;
1151: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1152: MatGetLocalSize(Amat,&m1,&n1);
1153: MatGetLocalSize(pc->mat,&m2,&n2);
1154: 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);
1155: MatGetLocalSize(Pmat,&m1,&n1);
1156: MatGetLocalSize(pc->pmat,&m2,&n2);
1157: 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);
1158: }
1160: if (Pmat != pc->pmat) {
1161: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1162: pc->matnonzerostate = -1;
1163: pc->matstate = -1;
1164: }
1166: /* reference first in case the matrices are the same */
1167: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1168: MatDestroy(&pc->mat);
1169: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1170: MatDestroy(&pc->pmat);
1171: pc->mat = Amat;
1172: pc->pmat = Pmat;
1173: return(0);
1174: }
1178: /*@
1179: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1181: Logically Collective on PC
1183: Input Parameters:
1184: + pc - the preconditioner context
1185: - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1187: Level: intermediate
1189: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1190: @*/
1191: PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1192: {
1195: pc->reusepreconditioner = flag;
1196: return(0);
1197: }
1201: /*@
1202: PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1204: Not Collective
1206: Input Parameter:
1207: . pc - the preconditioner context
1209: Output Parameter:
1210: . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1212: Level: intermediate
1214: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1215: @*/
1216: PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1217: {
1220: *flag = pc->reusepreconditioner;
1221: return(0);
1222: }
1226: /*@C
1227: PCGetOperators - Gets the matrix associated with the linear system and
1228: possibly a different one associated with the preconditioner.
1230: Not collective, though parallel Mats are returned if the PC is parallel
1232: Input Parameter:
1233: . pc - the preconditioner context
1235: Output Parameters:
1236: + Amat - the matrix defining the linear system
1237: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1239: Level: intermediate
1241: Notes: Does not increase the reference count of the matrices, so you should not destroy them
1243: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1244: are created in PC and returned to the user. In this case, if both operators
1245: mat and pmat are requested, two DIFFERENT operators will be returned. If
1246: only one is requested both operators in the PC will be the same (i.e. as
1247: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1248: The user must set the sizes of the returned matrices and their type etc just
1249: as if the user created them with MatCreate(). For example,
1251: $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1252: $ set size, type, etc of Amat
1254: $ MatCreate(comm,&mat);
1255: $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1256: $ PetscObjectDereference((PetscObject)mat);
1257: $ set size, type, etc of Amat
1259: and
1261: $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1262: $ set size, type, etc of Amat and Pmat
1264: $ MatCreate(comm,&Amat);
1265: $ MatCreate(comm,&Pmat);
1266: $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1267: $ PetscObjectDereference((PetscObject)Amat);
1268: $ PetscObjectDereference((PetscObject)Pmat);
1269: $ set size, type, etc of Amat and Pmat
1271: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1272: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1273: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1274: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1275: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1276: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1277: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1278: it can be created for you?
1281: .keywords: PC, get, operators, matrix, linear system
1283: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1284: @*/
1285: PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1286: {
1291: if (Amat) {
1292: if (!pc->mat) {
1293: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1294: pc->mat = pc->pmat;
1295: PetscObjectReference((PetscObject)pc->mat);
1296: } else { /* both Amat and Pmat are empty */
1297: MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1298: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1299: pc->pmat = pc->mat;
1300: PetscObjectReference((PetscObject)pc->pmat);
1301: }
1302: }
1303: }
1304: *Amat = pc->mat;
1305: }
1306: if (Pmat) {
1307: if (!pc->pmat) {
1308: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1309: pc->pmat = pc->mat;
1310: PetscObjectReference((PetscObject)pc->pmat);
1311: } else {
1312: MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1313: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1314: pc->mat = pc->pmat;
1315: PetscObjectReference((PetscObject)pc->mat);
1316: }
1317: }
1318: }
1319: *Pmat = pc->pmat;
1320: }
1321: return(0);
1322: }
1326: /*@C
1327: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1328: possibly a different one associated with the preconditioner have been set in the PC.
1330: Not collective, though the results on all processes should be the same
1332: Input Parameter:
1333: . pc - the preconditioner context
1335: Output Parameters:
1336: + mat - the matrix associated with the linear system was set
1337: - pmat - matrix associated with the preconditioner was set, usually the same
1339: Level: intermediate
1341: .keywords: PC, get, operators, matrix, linear system
1343: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1344: @*/
1345: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1346: {
1349: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1350: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1351: return(0);
1352: }
1356: /*@
1357: PCFactorGetMatrix - Gets the factored matrix from the
1358: preconditioner context. This routine is valid only for the LU,
1359: incomplete LU, Cholesky, and incomplete Cholesky methods.
1361: Not Collective on PC though Mat is parallel if PC is parallel
1363: Input Parameters:
1364: . pc - the preconditioner context
1366: Output parameters:
1367: . mat - the factored matrix
1369: Level: advanced
1371: Notes: Does not increase the reference count for the matrix so DO NOT destroy it
1373: .keywords: PC, get, factored, matrix
1374: @*/
1375: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1376: {
1382: if (pc->ops->getfactoredmatrix) {
1383: (*pc->ops->getfactoredmatrix)(pc,mat);
1384: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1385: return(0);
1386: }
1390: /*@C
1391: PCSetOptionsPrefix - Sets the prefix used for searching for all
1392: PC options in the database.
1394: Logically Collective on PC
1396: Input Parameters:
1397: + pc - the preconditioner context
1398: - prefix - the prefix string to prepend to all PC option requests
1400: Notes:
1401: A hyphen (-) must NOT be given at the beginning of the prefix name.
1402: The first character of all runtime options is AUTOMATICALLY the
1403: hyphen.
1405: Level: advanced
1407: .keywords: PC, set, options, prefix, database
1409: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1410: @*/
1411: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1412: {
1417: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1418: return(0);
1419: }
1423: /*@C
1424: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1425: PC options in the database.
1427: Logically Collective on PC
1429: Input Parameters:
1430: + pc - the preconditioner context
1431: - prefix - the prefix string to prepend to all PC option requests
1433: Notes:
1434: A hyphen (-) must NOT be given at the beginning of the prefix name.
1435: The first character of all runtime options is AUTOMATICALLY the
1436: hyphen.
1438: Level: advanced
1440: .keywords: PC, append, options, prefix, database
1442: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1443: @*/
1444: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1445: {
1450: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1451: return(0);
1452: }
1456: /*@C
1457: PCGetOptionsPrefix - Gets the prefix used for searching for all
1458: PC options in the database.
1460: Not Collective
1462: Input Parameters:
1463: . pc - the preconditioner context
1465: Output Parameters:
1466: . prefix - pointer to the prefix string used, is returned
1468: Notes: On the fortran side, the user should pass in a string 'prifix' of
1469: sufficient length to hold the prefix.
1471: Level: advanced
1473: .keywords: PC, get, options, prefix, database
1475: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1476: @*/
1477: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1478: {
1484: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1485: return(0);
1486: }
1490: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc,PetscBool *change)
1491: {
1497: *change = PETSC_FALSE;
1498: PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1499: return(0);
1500: }
1504: /*@
1505: PCPreSolve - Optional pre-solve phase, intended for any
1506: preconditioner-specific actions that must be performed before
1507: the iterative solve itself.
1509: Collective on PC
1511: Input Parameters:
1512: + pc - the preconditioner context
1513: - ksp - the Krylov subspace context
1515: Level: developer
1517: Sample of Usage:
1518: .vb
1519: PCPreSolve(pc,ksp);
1520: KSPSolve(ksp,b,x);
1521: PCPostSolve(pc,ksp);
1522: .ve
1524: Notes:
1525: The pre-solve phase is distinct from the PCSetUp() phase.
1527: KSPSolve() calls this directly, so is rarely called by the user.
1529: .keywords: PC, pre-solve
1531: .seealso: PCPostSolve()
1532: @*/
1533: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1534: {
1536: Vec x,rhs;
1541: pc->presolvedone++;
1542: if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1543: KSPGetSolution(ksp,&x);
1544: KSPGetRhs(ksp,&rhs);
1546: if (pc->ops->presolve) {
1547: (*pc->ops->presolve)(pc,ksp,rhs,x);
1548: }
1549: return(0);
1550: }
1554: /*@
1555: PCPostSolve - Optional post-solve phase, intended for any
1556: preconditioner-specific actions that must be performed after
1557: the iterative solve itself.
1559: Collective on PC
1561: Input Parameters:
1562: + pc - the preconditioner context
1563: - ksp - the Krylov subspace context
1565: Sample of Usage:
1566: .vb
1567: PCPreSolve(pc,ksp);
1568: KSPSolve(ksp,b,x);
1569: PCPostSolve(pc,ksp);
1570: .ve
1572: Note:
1573: KSPSolve() calls this routine directly, so it is rarely called by the user.
1575: Level: developer
1577: .keywords: PC, post-solve
1579: .seealso: PCPreSolve(), KSPSolve()
1580: @*/
1581: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1582: {
1584: Vec x,rhs;
1589: pc->presolvedone--;
1590: KSPGetSolution(ksp,&x);
1591: KSPGetRhs(ksp,&rhs);
1592: if (pc->ops->postsolve) {
1593: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1594: }
1595: return(0);
1596: }
1600: /*@C
1601: PCLoad - Loads a PC that has been stored in binary with PCView().
1603: Collective on PetscViewer
1605: Input Parameters:
1606: + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1607: some related function before a call to PCLoad().
1608: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1610: Level: intermediate
1612: Notes:
1613: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1615: Notes for advanced users:
1616: Most users should not need to know the details of the binary storage
1617: format, since PCLoad() and PCView() completely hide these details.
1618: But for anyone who's interested, the standard binary matrix storage
1619: format is
1620: .vb
1621: has not yet been determined
1622: .ve
1624: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1625: @*/
1626: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1627: {
1629: PetscBool isbinary;
1630: PetscInt classid;
1631: char type[256];
1636: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1637: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1639: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1640: if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1641: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1642: PCSetType(newdm, type);
1643: if (newdm->ops->load) {
1644: (*newdm->ops->load)(newdm,viewer);
1645: }
1646: return(0);
1647: }
1649: #include <petscdraw.h>
1650: #if defined(PETSC_HAVE_SAWS)
1651: #include <petscviewersaws.h>
1652: #endif
1655: /*@C
1656: PCView - Prints the PC data structure.
1658: Collective on PC
1660: Input Parameters:
1661: + PC - the PC context
1662: - viewer - optional visualization context
1664: Note:
1665: The available visualization contexts include
1666: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1667: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1668: output where only the first processor opens
1669: the file. All other processors send their
1670: data to the first processor to print.
1672: The user can open an alternative visualization contexts with
1673: PetscViewerASCIIOpen() (output to a specified file).
1675: Level: developer
1677: .keywords: PC, view
1679: .seealso: KSPView(), PetscViewerASCIIOpen()
1680: @*/
1681: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1682: {
1683: PCType cstr;
1684: PetscErrorCode ierr;
1685: PetscBool iascii,isstring,isbinary,isdraw;
1686: PetscViewerFormat format;
1687: #if defined(PETSC_HAVE_SAWS)
1688: PetscBool issaws;
1689: #endif
1693: if (!viewer) {
1694: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1695: }
1699: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1700: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1701: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1702: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1703: #if defined(PETSC_HAVE_SAWS)
1704: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1705: #endif
1707: if (iascii) {
1708: PetscViewerGetFormat(viewer,&format);
1709: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1710: if (!pc->setupcalled) {
1711: PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");
1712: }
1713: if (pc->ops->view) {
1714: PetscViewerASCIIPushTab(viewer);
1715: (*pc->ops->view)(pc,viewer);
1716: PetscViewerASCIIPopTab(viewer);
1717: }
1718: if (pc->mat) {
1719: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1720: if (pc->pmat == pc->mat) {
1721: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1722: PetscViewerASCIIPushTab(viewer);
1723: MatView(pc->mat,viewer);
1724: PetscViewerASCIIPopTab(viewer);
1725: } else {
1726: if (pc->pmat) {
1727: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1728: } else {
1729: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1730: }
1731: PetscViewerASCIIPushTab(viewer);
1732: MatView(pc->mat,viewer);
1733: if (pc->pmat) {MatView(pc->pmat,viewer);}
1734: PetscViewerASCIIPopTab(viewer);
1735: }
1736: PetscViewerPopFormat(viewer);
1737: }
1738: } else if (isstring) {
1739: PCGetType(pc,&cstr);
1740: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1741: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1742: } else if (isbinary) {
1743: PetscInt classid = PC_FILE_CLASSID;
1744: MPI_Comm comm;
1745: PetscMPIInt rank;
1746: char type[256];
1748: PetscObjectGetComm((PetscObject)pc,&comm);
1749: MPI_Comm_rank(comm,&rank);
1750: if (!rank) {
1751: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1752: PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1753: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1754: }
1755: if (pc->ops->view) {
1756: (*pc->ops->view)(pc,viewer);
1757: }
1758: } else if (isdraw) {
1759: PetscDraw draw;
1760: char str[25];
1761: PetscReal x,y,bottom,h;
1762: PetscInt n;
1764: PetscViewerDrawGetDraw(viewer,0,&draw);
1765: PetscDrawGetCurrentPoint(draw,&x,&y);
1766: if (pc->mat) {
1767: MatGetSize(pc->mat,&n,NULL);
1768: PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1769: } else {
1770: PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1771: }
1772: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1773: bottom = y - h;
1774: PetscDrawPushCurrentPoint(draw,x,bottom);
1775: if (pc->ops->view) {
1776: (*pc->ops->view)(pc,viewer);
1777: }
1778: PetscDrawPopCurrentPoint(draw);
1779: #if defined(PETSC_HAVE_SAWS)
1780: } else if (issaws) {
1781: PetscMPIInt rank;
1783: PetscObjectName((PetscObject)pc);
1784: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1785: if (!((PetscObject)pc)->amsmem && !rank) {
1786: PetscObjectViewSAWs((PetscObject)pc,viewer);
1787: }
1788: if (pc->mat) {MatView(pc->mat,viewer);}
1789: if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1790: #endif
1791: }
1792: return(0);
1793: }
1797: /*@C
1798: PCRegister - Adds a method to the preconditioner package.
1800: Not collective
1802: Input Parameters:
1803: + name_solver - name of a new user-defined solver
1804: - routine_create - routine to create method context
1806: Notes:
1807: PCRegister() may be called multiple times to add several user-defined preconditioners.
1809: Sample usage:
1810: .vb
1811: PCRegister("my_solver", MySolverCreate);
1812: .ve
1814: Then, your solver can be chosen with the procedural interface via
1815: $ PCSetType(pc,"my_solver")
1816: or at runtime via the option
1817: $ -pc_type my_solver
1819: Level: advanced
1821: .keywords: PC, register
1823: .seealso: PCRegisterAll(), PCRegisterDestroy()
1824: @*/
1825: PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1826: {
1830: PetscFunctionListAdd(&PCList,sname,function);
1831: return(0);
1832: }
1836: /*@
1837: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1839: Collective on PC
1841: Input Parameter:
1842: . pc - the preconditioner object
1844: Output Parameter:
1845: . mat - the explict preconditioned operator
1847: Notes:
1848: This computation is done by applying the operators to columns of the
1849: identity matrix.
1851: Currently, this routine uses a dense matrix format when 1 processor
1852: is used and a sparse format otherwise. This routine is costly in general,
1853: and is recommended for use only with relatively small systems.
1855: Level: advanced
1857: .keywords: PC, compute, explicit, operator
1859: .seealso: KSPComputeExplicitOperator()
1861: @*/
1862: PetscErrorCode PCComputeExplicitOperator(PC pc,Mat *mat)
1863: {
1864: Vec in,out;
1866: PetscInt i,M,m,*rows,start,end;
1867: PetscMPIInt size;
1868: MPI_Comm comm;
1869: PetscScalar *array,one = 1.0;
1875: PetscObjectGetComm((PetscObject)pc,&comm);
1876: MPI_Comm_size(comm,&size);
1878: if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1879: MatCreateVecs(pc->pmat,&in,0);
1880: VecDuplicate(in,&out);
1881: VecGetOwnershipRange(in,&start,&end);
1882: VecGetSize(in,&M);
1883: VecGetLocalSize(in,&m);
1884: PetscMalloc1(m+1,&rows);
1885: for (i=0; i<m; i++) rows[i] = start + i;
1887: MatCreate(comm,mat);
1888: MatSetSizes(*mat,m,m,M,M);
1889: if (size == 1) {
1890: MatSetType(*mat,MATSEQDENSE);
1891: MatSeqDenseSetPreallocation(*mat,NULL);
1892: } else {
1893: MatSetType(*mat,MATMPIAIJ);
1894: MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1895: }
1896: MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
1898: for (i=0; i<M; i++) {
1900: VecSet(in,0.0);
1901: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1902: VecAssemblyBegin(in);
1903: VecAssemblyEnd(in);
1905: /* should fix, allowing user to choose side */
1906: PCApply(pc,in,out);
1908: VecGetArray(out,&array);
1909: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1910: VecRestoreArray(out,&array);
1912: }
1913: PetscFree(rows);
1914: VecDestroy(&out);
1915: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1916: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1917: return(0);
1918: }
1922: /*@
1923: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1925: Collective on PC
1927: Input Parameters:
1928: + pc - the solver context
1929: . dim - the dimension of the coordinates 1, 2, or 3
1930: - coords - the coordinates
1932: Level: intermediate
1934: Notes: coords is an array of the 3D coordinates for the nodes on
1935: the local processor. So if there are 108 equation on a processor
1936: for a displacement finite element discretization of elasticity (so
1937: that there are 36 = 108/3 nodes) then the array must have 108
1938: double precision values (ie, 3 * 36). These x y z coordinates
1939: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1940: ... , N-1.z ].
1942: .seealso: MatSetNearNullSpace
1943: @*/
1944: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1945: {
1949: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1950: return(0);
1951: }