Actual source code: precon.c
petsc-3.13.6 2020-09-29
2: /*
3: The PC (preconditioner) interface routines, callable by users.
4: */
5: #include <petsc/private/pcimpl.h>
6: #include <petscdm.h>
8: /* Logging support */
9: PetscClassId PC_CLASSID;
10: PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
11: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks;
12: PetscInt PetscMGLevelId;
14: PetscErrorCode PCGetDefaultType_Private(PC pc,const char *type[])
15: {
17: PetscMPIInt size;
18: PetscBool hasop,flg1,flg2,set,flg3;
21: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
22: if (pc->pmat) {
23: MatHasOperation(pc->pmat,MATOP_GET_DIAGONAL_BLOCK,&hasop);
24: if (size == 1) {
25: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
26: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
27: MatIsSymmetricKnown(pc->pmat,&set,&flg3);
28: if (flg1 && (!flg2 || (set && flg3))) {
29: *type = PCICC;
30: } else if (flg2) {
31: *type = PCILU;
32: } else if (hasop) { /* likely is a parallel matrix run on one processor */
33: *type = PCBJACOBI;
34: } else {
35: *type = PCNONE;
36: }
37: } else {
38: if (hasop) {
39: *type = PCBJACOBI;
40: } else {
41: *type = PCNONE;
42: }
43: }
44: } else {
45: if (size == 1) {
46: *type = PCILU;
47: } else {
48: *type = PCBJACOBI;
49: }
50: }
51: return(0);
52: }
54: /*@
55: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
57: Collective on PC
59: Input Parameter:
60: . pc - the preconditioner context
62: Level: developer
64: Notes:
65: This allows a PC to be reused for a different sized linear system but using the same options that have been previously set in the PC
67: .seealso: PCCreate(), PCSetUp()
68: @*/
69: PetscErrorCode PCReset(PC pc)
70: {
75: if (pc->ops->reset) {
76: (*pc->ops->reset)(pc);
77: }
78: VecDestroy(&pc->diagonalscaleright);
79: VecDestroy(&pc->diagonalscaleleft);
80: MatDestroy(&pc->pmat);
81: MatDestroy(&pc->mat);
83: pc->setupcalled = 0;
84: return(0);
85: }
87: /*@
88: PCDestroy - Destroys PC context that was created with PCCreate().
90: Collective on PC
92: Input Parameter:
93: . pc - the preconditioner context
95: Level: developer
97: .seealso: PCCreate(), PCSetUp()
98: @*/
99: PetscErrorCode PCDestroy(PC *pc)
100: {
104: if (!*pc) return(0);
106: if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; return(0);}
108: PCReset(*pc);
110: /* if memory was published with SAWs then destroy it */
111: PetscObjectSAWsViewOff((PetscObject)*pc);
112: if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
113: DMDestroy(&(*pc)->dm);
114: PetscHeaderDestroy(pc);
115: return(0);
116: }
118: /*@C
119: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
120: scaling as needed by certain time-stepping codes.
122: Logically Collective on PC
124: Input Parameter:
125: . pc - the preconditioner context
127: Output Parameter:
128: . flag - PETSC_TRUE if it applies the scaling
130: Level: developer
132: Notes:
133: If this returns PETSC_TRUE then the system solved via the Krylov method is
134: $ D M A D^{-1} y = D M b for left preconditioning or
135: $ D A M D^{-1} z = D b for right preconditioning
137: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
138: @*/
139: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
140: {
144: *flag = pc->diagonalscale;
145: return(0);
146: }
148: /*@
149: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
150: scaling as needed by certain time-stepping codes.
152: Logically Collective on PC
154: Input Parameters:
155: + pc - the preconditioner context
156: - s - scaling vector
158: Level: intermediate
160: Notes:
161: The system solved via the Krylov method is
162: $ D M A D^{-1} y = D M b for left preconditioning or
163: $ D A M D^{-1} z = D b for right preconditioning
165: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
167: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
168: @*/
169: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
170: {
176: pc->diagonalscale = PETSC_TRUE;
178: PetscObjectReference((PetscObject)s);
179: VecDestroy(&pc->diagonalscaleleft);
181: pc->diagonalscaleleft = s;
183: VecDuplicate(s,&pc->diagonalscaleright);
184: VecCopy(s,pc->diagonalscaleright);
185: VecReciprocal(pc->diagonalscaleright);
186: return(0);
187: }
189: /*@
190: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
192: Logically Collective on PC
194: Input Parameters:
195: + pc - the preconditioner context
196: . in - input vector
197: - out - scaled vector (maybe the same as in)
199: Level: intermediate
201: Notes:
202: The system solved via the Krylov method is
203: $ D M A D^{-1} y = D M b for left preconditioning or
204: $ D A M D^{-1} z = D b for right preconditioning
206: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
208: If diagonal scaling is turned off and in is not out then in is copied to out
210: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
211: @*/
212: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
213: {
220: if (pc->diagonalscale) {
221: VecPointwiseMult(out,pc->diagonalscaleleft,in);
222: } else if (in != out) {
223: VecCopy(in,out);
224: }
225: return(0);
226: }
228: /*@
229: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
231: Logically Collective on PC
233: Input Parameters:
234: + pc - the preconditioner context
235: . in - input vector
236: - out - scaled vector (maybe the same as in)
238: Level: intermediate
240: Notes:
241: The system solved via the Krylov method is
242: $ D M A D^{-1} y = D M b for left preconditioning or
243: $ D A M D^{-1} z = D b for right preconditioning
245: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
247: If diagonal scaling is turned off and in is not out then in is copied to out
249: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
250: @*/
251: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
252: {
259: if (pc->diagonalscale) {
260: VecPointwiseMult(out,pc->diagonalscaleright,in);
261: } else if (in != out) {
262: VecCopy(in,out);
263: }
264: return(0);
265: }
267: /*@
268: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
269: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
270: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
272: Logically Collective on PC
274: Input Parameters:
275: + pc - the preconditioner context
276: - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
278: Options Database Key:
279: . -pc_use_amat <true,false>
281: Notes:
282: For the common case in which the linear system matrix and the matrix used to construct the
283: preconditioner are identical, this routine is does nothing.
285: Level: intermediate
287: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
288: @*/
289: PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
290: {
293: pc->useAmat = flg;
294: return(0);
295: }
297: /*@
298: PCSetErrorIfFailure - Causes PC to generate an error if a FPE, for example a zero pivot, is detected.
300: Logically Collective on PC
302: Input Parameters:
303: + pc - iterative context obtained from PCCreate()
304: - flg - PETSC_TRUE indicates you want the error generated
306: Level: advanced
308: Notes:
309: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
310: 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
311: detected.
313: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
315: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
316: @*/
317: PetscErrorCode PCSetErrorIfFailure(PC pc,PetscBool flg)
318: {
322: pc->erroriffailure = flg;
323: return(0);
324: }
326: /*@
327: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
328: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
329: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
331: Logically Collective on PC
333: Input Parameter:
334: . pc - the preconditioner context
336: Output Parameter:
337: . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
339: Notes:
340: For the common case in which the linear system matrix and the matrix used to construct the
341: preconditioner are identical, this routine is does nothing.
343: Level: intermediate
345: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
346: @*/
347: PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
348: {
351: *flg = pc->useAmat;
352: return(0);
353: }
355: /*@
356: PCCreate - Creates a preconditioner context.
358: Collective
360: Input Parameter:
361: . comm - MPI communicator
363: Output Parameter:
364: . pc - location to put the preconditioner context
366: Notes:
367: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or PCICC
368: in parallel. For dense matrices it is always PCNONE.
370: Level: developer
372: .seealso: PCSetUp(), PCApply(), PCDestroy()
373: @*/
374: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
375: {
376: PC pc;
381: *newpc = 0;
382: PCInitializePackage();
384: PetscHeaderCreate(pc,PC_CLASSID,"PC","Preconditioner","PC",comm,PCDestroy,PCView);
386: pc->mat = 0;
387: pc->pmat = 0;
388: pc->setupcalled = 0;
389: pc->setfromoptionscalled = 0;
390: pc->data = 0;
391: pc->diagonalscale = PETSC_FALSE;
392: pc->diagonalscaleleft = 0;
393: pc->diagonalscaleright = 0;
395: pc->modifysubmatrices = 0;
396: pc->modifysubmatricesP = 0;
398: *newpc = pc;
399: return(0);
401: }
403: /* -------------------------------------------------------------------------------*/
405: /*@
406: PCApply - Applies the preconditioner to a vector.
408: Collective on PC
410: Input Parameters:
411: + pc - the preconditioner context
412: - x - input vector
414: Output Parameter:
415: . y - output vector
417: Level: developer
419: .seealso: PCApplyTranspose(), PCApplyBAorAB()
420: @*/
421: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
422: {
424: PetscInt m,n,mv,nv;
430: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
431: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
432: /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
433: MatGetLocalSize(pc->pmat,&m,&n);
434: VecGetLocalSize(x,&nv);
435: VecGetLocalSize(y,&mv);
436: 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);
437: 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);
438: VecSetErrorIfLocked(y,3);
440: PCSetUp(pc);
441: if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
442: VecLockReadPush(x);
443: PetscLogEventBegin(PC_Apply,pc,x,y,0);
444: (*pc->ops->apply)(pc,x,y);
445: PetscLogEventEnd(PC_Apply,pc,x,y,0);
446: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
447: VecLockReadPop(x);
448: return(0);
449: }
451: /*@
452: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
454: Collective on PC
456: Input Parameters:
457: + pc - the preconditioner context
458: - x - input vector
460: Output Parameter:
461: . y - output vector
463: Notes:
464: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
466: Level: developer
468: .seealso: PCApply(), PCApplySymmetricRight()
469: @*/
470: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
471: {
478: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
479: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
480: PCSetUp(pc);
481: if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
482: VecLockReadPush(x);
483: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
484: (*pc->ops->applysymmetricleft)(pc,x,y);
485: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
486: VecLockReadPop(x);
487: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
488: return(0);
489: }
491: /*@
492: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
494: Collective on PC
496: Input Parameters:
497: + pc - the preconditioner context
498: - x - input vector
500: Output Parameter:
501: . y - output vector
503: Level: developer
505: Notes:
506: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
508: .seealso: PCApply(), PCApplySymmetricLeft()
509: @*/
510: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
511: {
518: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
519: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
520: PCSetUp(pc);
521: if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
522: VecLockReadPush(x);
523: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
524: (*pc->ops->applysymmetricright)(pc,x,y);
525: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
526: VecLockReadPop(x);
527: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
528: return(0);
529: }
531: /*@
532: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
534: Collective on PC
536: Input Parameters:
537: + pc - the preconditioner context
538: - x - input vector
540: Output Parameter:
541: . y - output vector
543: Notes:
544: For complex numbers this applies the non-Hermitian transpose.
546: Developer Notes:
547: We need to implement a PCApplyHermitianTranspose()
549: Level: developer
551: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
552: @*/
553: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
554: {
561: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
562: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
563: PCSetUp(pc);
564: if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
565: VecLockReadPush(x);
566: PetscLogEventBegin(PC_Apply,pc,x,y,0);
567: (*pc->ops->applytranspose)(pc,x,y);
568: PetscLogEventEnd(PC_Apply,pc,x,y,0);
569: VecLockReadPop(x);
570: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
571: return(0);
572: }
574: /*@
575: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
577: Collective on PC
579: Input Parameters:
580: . pc - the preconditioner context
582: Output Parameter:
583: . flg - PETSC_TRUE if a transpose operation is defined
585: Level: developer
587: .seealso: PCApplyTranspose()
588: @*/
589: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
590: {
594: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
595: else *flg = PETSC_FALSE;
596: return(0);
597: }
599: /*@
600: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
602: Collective on PC
604: Input Parameters:
605: + pc - the preconditioner context
606: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
607: . x - input vector
608: - work - work vector
610: Output Parameter:
611: . y - output vector
613: Level: developer
615: Notes:
616: 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
617: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
619: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
620: @*/
621: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
622: {
634: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
635: 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");
636: if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner Section 1.5 Writing Application Codes with PETSc");
637: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
639: PCSetUp(pc);
640: if (pc->diagonalscale) {
641: if (pc->ops->applyBA) {
642: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
643: VecDuplicate(x,&work2);
644: PCDiagonalScaleRight(pc,x,work2);
645: (*pc->ops->applyBA)(pc,side,work2,y,work);
646: PCDiagonalScaleLeft(pc,y,y);
647: VecDestroy(&work2);
648: } else if (side == PC_RIGHT) {
649: PCDiagonalScaleRight(pc,x,y);
650: PCApply(pc,y,work);
651: MatMult(pc->mat,work,y);
652: PCDiagonalScaleLeft(pc,y,y);
653: } else if (side == PC_LEFT) {
654: PCDiagonalScaleRight(pc,x,y);
655: MatMult(pc->mat,y,work);
656: PCApply(pc,work,y);
657: PCDiagonalScaleLeft(pc,y,y);
658: } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric Section 1.5 Writing Application Codes with PETSc of preconditioner");
659: } else {
660: if (pc->ops->applyBA) {
661: (*pc->ops->applyBA)(pc,side,x,y,work);
662: } else if (side == PC_RIGHT) {
663: PCApply(pc,x,work);
664: MatMult(pc->mat,work,y);
665: } else if (side == PC_LEFT) {
666: MatMult(pc->mat,x,work);
667: PCApply(pc,work,y);
668: } else if (side == PC_SYMMETRIC) {
669: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
670: PCApplySymmetricRight(pc,x,work);
671: MatMult(pc->mat,work,y);
672: VecCopy(y,work);
673: PCApplySymmetricLeft(pc,work,y);
674: }
675: }
676: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
677: return(0);
678: }
680: /*@
681: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
682: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
683: NOT tr(B*A) = tr(A)*tr(B).
685: Collective on PC
687: Input Parameters:
688: + pc - the preconditioner context
689: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
690: . x - input vector
691: - work - work vector
693: Output Parameter:
694: . y - output vector
697: Notes:
698: 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
699: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
701: Level: developer
703: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
704: @*/
705: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
706: {
714: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
715: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
716: if (pc->ops->applyBAtranspose) {
717: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
718: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
719: return(0);
720: }
721: if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
723: PCSetUp(pc);
724: if (side == PC_RIGHT) {
725: PCApplyTranspose(pc,x,work);
726: MatMultTranspose(pc->mat,work,y);
727: } else if (side == PC_LEFT) {
728: MatMultTranspose(pc->mat,x,work);
729: PCApplyTranspose(pc,work,y);
730: }
731: /* add support for PC_SYMMETRIC */
732: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
733: return(0);
734: }
736: /* -------------------------------------------------------------------------------*/
738: /*@
739: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
740: built-in fast Section 1.5 Writing Application Codes with PETSc of Richardson's method.
742: Not Collective
744: Input Parameter:
745: . pc - the preconditioner
747: Output Parameter:
748: . exists - PETSC_TRUE or PETSC_FALSE
750: Level: developer
752: .seealso: PCApplyRichardson()
753: @*/
754: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
755: {
759: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
760: else *exists = PETSC_FALSE;
761: return(0);
762: }
764: /*@
765: PCApplyRichardson - Applies several steps of Richardson iteration with
766: the particular preconditioner. This routine is usually used by the
767: Krylov solvers and not the Section 1.5 Writing Application Codes with PETSc code directly.
769: Collective on PC
771: Input Parameters:
772: + pc - the preconditioner context
773: . b - the right hand side
774: . w - one work vector
775: . rtol - relative decrease in residual norm convergence criteria
776: . abstol - absolute residual norm convergence criteria
777: . dtol - divergence residual norm increase criteria
778: . its - the number of iterations to apply.
779: - guesszero - if the input x contains nonzero initial guess
781: Output Parameter:
782: + outits - number of iterations actually used (for SOR this always equals its)
783: . reason - the reason the apply terminated
784: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
786: Notes:
787: Most preconditioners do not support this function. Use the command
788: PCApplyRichardsonExists() to determine if one does.
790: Except for the multigrid PC this routine ignores the convergence tolerances
791: and always runs for the number of iterations
793: Level: developer
795: .seealso: PCApplyRichardsonExists()
796: @*/
797: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
798: {
806: if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
807: PCSetUp(pc);
808: if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
809: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
810: return(0);
811: }
813: /*@
814: PCGetFailedReason - Gets the reason a PCSetUp() failed or 0 if it did not fail
816: Logically Collective on PC
818: Input Parameter:
819: . pc - the preconditioner context
821: Output Parameter:
822: . reason - the reason it failed, currently only -1
824: Level: advanced
826: .seealso: PCCreate(), PCApply(), PCDestroy()
827: @*/
828: PetscErrorCode PCGetFailedReason(PC pc,PCFailedReason *reason)
829: {
831: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
832: else *reason = pc->failedreason;
833: return(0);
834: }
837: /*
838: a setupcall of 0 indicates never setup,
839: 1 indicates has been previously setup
840: -1 indicates a PCSetUp() was attempted and failed
841: */
842: /*@
843: PCSetUp - Prepares for the use of a preconditioner.
845: Collective on PC
847: Input Parameter:
848: . pc - the preconditioner context
850: Level: developer
852: .seealso: PCCreate(), PCApply(), PCDestroy()
853: @*/
854: PetscErrorCode PCSetUp(PC pc)
855: {
856: PetscErrorCode ierr;
857: const char *def;
858: PetscObjectState matstate, matnonzerostate;
862: if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
864: if (pc->setupcalled && pc->reusepreconditioner) {
865: PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
866: return(0);
867: }
869: PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
870: MatGetNonzeroState(pc->pmat,&matnonzerostate);
871: if (!pc->setupcalled) {
872: PetscInfo(pc,"Setting up PC for first time\n");
873: pc->flag = DIFFERENT_NONZERO_PATTERN;
874: } else if (matstate == pc->matstate) {
875: PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
876: return(0);
877: } else {
878: if (matnonzerostate > pc->matnonzerostate) {
879: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
880: pc->flag = DIFFERENT_NONZERO_PATTERN;
881: } else {
882: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
883: pc->flag = SAME_NONZERO_PATTERN;
884: }
885: }
886: pc->matstate = matstate;
887: pc->matnonzerostate = matnonzerostate;
889: if (!((PetscObject)pc)->type_name) {
890: PCGetDefaultType_Private(pc,&def);
891: PCSetType(pc,def);
892: }
894: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
895: MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
896: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
897: if (pc->ops->setup) {
898: (*pc->ops->setup)(pc);
899: }
900: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
901: if (!pc->setupcalled) pc->setupcalled = 1;
902: return(0);
903: }
905: /*@
906: PCSetUpOnBlocks - Sets up the preconditioner for each block in
907: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
908: methods.
910: Collective on PC
912: Input Parameters:
913: . pc - the preconditioner context
915: Level: developer
917: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
918: @*/
919: PetscErrorCode PCSetUpOnBlocks(PC pc)
920: {
925: if (!pc->ops->setuponblocks) return(0);
926: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
927: (*pc->ops->setuponblocks)(pc);
928: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
929: return(0);
930: }
932: /*@C
933: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
934: submatrices that arise within certain subdomain-based preconditioners.
935: The basic submatrices are extracted from the preconditioner matrix as
936: usual; the user can then alter these (for example, to set different boundary
937: conditions for each submatrix) before they are used for the local solves.
939: Logically Collective on PC
941: Input Parameters:
942: + pc - the preconditioner context
943: . func - routine for modifying the submatrices
944: - ctx - optional user-defined context (may be null)
946: Calling sequence of func:
947: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
949: + row - an array of index sets that contain the global row numbers
950: that comprise each local submatrix
951: . col - an array of index sets that contain the global column numbers
952: that comprise each local submatrix
953: . submat - array of local submatrices
954: - ctx - optional user-defined context for private data for the
955: user-defined func routine (may be null)
957: Notes:
958: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
959: KSPSolve().
961: A routine set by PCSetModifySubMatrices() is currently called within
962: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
963: preconditioners. All other preconditioners ignore this routine.
965: Level: advanced
967: .seealso: PCModifySubMatrices()
968: @*/
969: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
970: {
973: pc->modifysubmatrices = func;
974: pc->modifysubmatricesP = ctx;
975: return(0);
976: }
978: /*@C
979: PCModifySubMatrices - Calls an optional user-defined routine within
980: certain preconditioners if one has been set with PCSetModifySubMatrices().
982: Collective on PC
984: Input Parameters:
985: + pc - the preconditioner context
986: . nsub - the number of local submatrices
987: . row - an array of index sets that contain the global row numbers
988: that comprise each local submatrix
989: . col - an array of index sets that contain the global column numbers
990: that comprise each local submatrix
991: . submat - array of local submatrices
992: - ctx - optional user-defined context for private data for the
993: user-defined routine (may be null)
995: Output Parameter:
996: . submat - array of local submatrices (the entries of which may
997: have been modified)
999: Notes:
1000: The user should NOT generally call this routine, as it will
1001: automatically be called within certain preconditioners (currently
1002: block Jacobi, additive Schwarz) if set.
1004: The basic submatrices are extracted from the preconditioner matrix
1005: as usual; the user can then alter these (for example, to set different
1006: boundary conditions for each submatrix) before they are used for the
1007: local solves.
1009: Level: developer
1011: .seealso: PCSetModifySubMatrices()
1012: @*/
1013: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1014: {
1019: if (!pc->modifysubmatrices) return(0);
1020: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1021: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1022: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1023: return(0);
1024: }
1026: /*@
1027: PCSetOperators - Sets the matrix associated with the linear system and
1028: a (possibly) different one associated with the preconditioner.
1030: Logically Collective on PC
1032: Input Parameters:
1033: + pc - the preconditioner context
1034: . Amat - the matrix that defines the linear system
1035: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1037: Notes:
1038: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1040: If you wish to replace either Amat or Pmat but leave the other one untouched then
1041: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1042: on it and then pass it back in in your call to KSPSetOperators().
1044: More Notes about Repeated Solution of Linear Systems:
1045: PETSc does NOT reset the matrix entries of either Amat or Pmat
1046: to zero after a linear solve; the user is completely responsible for
1047: matrix assembly. See the routine MatZeroEntries() if desiring to
1048: zero all elements of a matrix.
1050: Level: intermediate
1052: .seealso: PCGetOperators(), MatZeroEntries()
1053: @*/
1054: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1055: {
1056: PetscErrorCode ierr;
1057: PetscInt m1,n1,m2,n2;
1065: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1066: MatGetLocalSize(Amat,&m1,&n1);
1067: MatGetLocalSize(pc->mat,&m2,&n2);
1068: 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);
1069: MatGetLocalSize(Pmat,&m1,&n1);
1070: MatGetLocalSize(pc->pmat,&m2,&n2);
1071: 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);
1072: }
1074: if (Pmat != pc->pmat) {
1075: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1076: pc->matnonzerostate = -1;
1077: pc->matstate = -1;
1078: }
1080: /* reference first in case the matrices are the same */
1081: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1082: MatDestroy(&pc->mat);
1083: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1084: MatDestroy(&pc->pmat);
1085: pc->mat = Amat;
1086: pc->pmat = Pmat;
1087: return(0);
1088: }
1090: /*@
1091: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1093: Logically Collective on PC
1095: Input Parameters:
1096: + pc - the preconditioner context
1097: - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1099: Level: intermediate
1101: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1102: @*/
1103: PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1104: {
1107: pc->reusepreconditioner = flag;
1108: return(0);
1109: }
1111: /*@
1112: PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1114: Not Collective
1116: Input Parameter:
1117: . pc - the preconditioner context
1119: Output Parameter:
1120: . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1122: Level: intermediate
1124: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1125: @*/
1126: PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1127: {
1130: *flag = pc->reusepreconditioner;
1131: return(0);
1132: }
1134: /*@
1135: PCGetOperators - Gets the matrix associated with the linear system and
1136: possibly a different one associated with the preconditioner.
1138: Not collective, though parallel Mats are returned if the PC is parallel
1140: Input Parameter:
1141: . pc - the preconditioner context
1143: Output Parameters:
1144: + Amat - the matrix defining the linear system
1145: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1147: Level: intermediate
1149: Notes:
1150: Does not increase the reference count of the matrices, so you should not destroy them
1152: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1153: are created in PC and returned to the user. In this case, if both operators
1154: mat and pmat are requested, two DIFFERENT operators will be returned. If
1155: only one is requested both operators in the PC will be the same (i.e. as
1156: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1157: The user must set the sizes of the returned matrices and their type etc just
1158: as if the user created them with MatCreate(). For example,
1160: $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1161: $ set size, type, etc of Amat
1163: $ MatCreate(comm,&mat);
1164: $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1165: $ PetscObjectDereference((PetscObject)mat);
1166: $ set size, type, etc of Amat
1168: and
1170: $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1171: $ set size, type, etc of Amat and Pmat
1173: $ MatCreate(comm,&Amat);
1174: $ MatCreate(comm,&Pmat);
1175: $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1176: $ PetscObjectDereference((PetscObject)Amat);
1177: $ PetscObjectDereference((PetscObject)Pmat);
1178: $ set size, type, etc of Amat and Pmat
1180: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1181: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1182: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1183: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1184: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1185: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1186: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1187: it can be created for you?
1190: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1191: @*/
1192: PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1193: {
1198: if (Amat) {
1199: if (!pc->mat) {
1200: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1201: pc->mat = pc->pmat;
1202: PetscObjectReference((PetscObject)pc->mat);
1203: } else { /* both Amat and Pmat are empty */
1204: MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1205: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1206: pc->pmat = pc->mat;
1207: PetscObjectReference((PetscObject)pc->pmat);
1208: }
1209: }
1210: }
1211: *Amat = pc->mat;
1212: }
1213: if (Pmat) {
1214: if (!pc->pmat) {
1215: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1216: pc->pmat = pc->mat;
1217: PetscObjectReference((PetscObject)pc->pmat);
1218: } else {
1219: MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1220: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1221: pc->mat = pc->pmat;
1222: PetscObjectReference((PetscObject)pc->mat);
1223: }
1224: }
1225: }
1226: *Pmat = pc->pmat;
1227: }
1228: return(0);
1229: }
1231: /*@C
1232: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1233: possibly a different one associated with the preconditioner have been set in the PC.
1235: Not collective, though the results on all processes should be the same
1237: Input Parameter:
1238: . pc - the preconditioner context
1240: Output Parameters:
1241: + mat - the matrix associated with the linear system was set
1242: - pmat - matrix associated with the preconditioner was set, usually the same
1244: Level: intermediate
1246: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1247: @*/
1248: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1249: {
1252: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1253: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1254: return(0);
1255: }
1257: /*@
1258: PCFactorGetMatrix - Gets the factored matrix from the
1259: preconditioner context. This routine is valid only for the LU,
1260: incomplete LU, Cholesky, and incomplete Cholesky methods.
1262: Not Collective on PC though Mat is parallel if PC is parallel
1264: Input Parameters:
1265: . pc - the preconditioner context
1267: Output parameters:
1268: . mat - the factored matrix
1270: Level: advanced
1272: Notes:
1273: Does not increase the reference count for the matrix so DO NOT destroy it
1275: @*/
1276: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1277: {
1283: if (pc->ops->getfactoredmatrix) {
1284: (*pc->ops->getfactoredmatrix)(pc,mat);
1285: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1286: return(0);
1287: }
1289: /*@C
1290: PCSetOptionsPrefix - Sets the prefix used for searching for all
1291: PC options in the database.
1293: Logically Collective on PC
1295: Input Parameters:
1296: + pc - the preconditioner context
1297: - prefix - the prefix string to prepend to all PC option requests
1299: Notes:
1300: A hyphen (-) must NOT be given at the beginning of the prefix name.
1301: The first character of all runtime options is AUTOMATICALLY the
1302: hyphen.
1304: Level: advanced
1306: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1307: @*/
1308: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1309: {
1314: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1315: return(0);
1316: }
1318: /*@C
1319: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1320: PC options in the database.
1322: Logically Collective on PC
1324: Input Parameters:
1325: + pc - the preconditioner context
1326: - prefix - the prefix string to prepend to all PC option requests
1328: Notes:
1329: A hyphen (-) must NOT be given at the beginning of the prefix name.
1330: The first character of all runtime options is AUTOMATICALLY the
1331: hyphen.
1333: Level: advanced
1335: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1336: @*/
1337: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1338: {
1343: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1344: return(0);
1345: }
1347: /*@C
1348: PCGetOptionsPrefix - Gets the prefix used for searching for all
1349: PC options in the database.
1351: Not Collective
1353: Input Parameters:
1354: . pc - the preconditioner context
1356: Output Parameters:
1357: . prefix - pointer to the prefix string used, is returned
1359: Notes:
1360: On the fortran side, the user should pass in a string 'prifix' of
1361: sufficient length to hold the prefix.
1363: Level: advanced
1365: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1366: @*/
1367: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1368: {
1374: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1375: return(0);
1376: }
1378: /*
1379: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1380: preconditioners including BDDC and Eisentat that transform the equations before applying
1381: the Krylov methods
1382: */
1383: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc,PetscBool *change)
1384: {
1390: *change = PETSC_FALSE;
1391: PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1392: return(0);
1393: }
1395: /*@
1396: PCPreSolve - Optional pre-solve phase, intended for any
1397: preconditioner-specific actions that must be performed before
1398: the iterative solve itself.
1400: Collective on PC
1402: Input Parameters:
1403: + pc - the preconditioner context
1404: - ksp - the Krylov subspace context
1406: Level: developer
1408: Sample of Usage:
1409: .vb
1410: PCPreSolve(pc,ksp);
1411: KSPSolve(ksp,b,x);
1412: PCPostSolve(pc,ksp);
1413: .ve
1415: Notes:
1416: The pre-solve phase is distinct from the PCSetUp() phase.
1418: KSPSolve() calls this directly, so is rarely called by the user.
1420: .seealso: PCPostSolve()
1421: @*/
1422: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1423: {
1425: Vec x,rhs;
1430: pc->presolvedone++;
1431: if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1432: KSPGetSolution(ksp,&x);
1433: KSPGetRhs(ksp,&rhs);
1435: if (pc->ops->presolve) {
1436: (*pc->ops->presolve)(pc,ksp,rhs,x);
1437: }
1438: return(0);
1439: }
1441: /*@
1442: PCPostSolve - Optional post-solve phase, intended for any
1443: preconditioner-specific actions that must be performed after
1444: the iterative solve itself.
1446: Collective on PC
1448: Input Parameters:
1449: + pc - the preconditioner context
1450: - ksp - the Krylov subspace context
1452: Sample of Usage:
1453: .vb
1454: PCPreSolve(pc,ksp);
1455: KSPSolve(ksp,b,x);
1456: PCPostSolve(pc,ksp);
1457: .ve
1459: Note:
1460: KSPSolve() calls this routine directly, so it is rarely called by the user.
1462: Level: developer
1464: .seealso: PCPreSolve(), KSPSolve()
1465: @*/
1466: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1467: {
1469: Vec x,rhs;
1474: pc->presolvedone--;
1475: KSPGetSolution(ksp,&x);
1476: KSPGetRhs(ksp,&rhs);
1477: if (pc->ops->postsolve) {
1478: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1479: }
1480: return(0);
1481: }
1483: /*@C
1484: PCLoad - Loads a PC that has been stored in binary with PCView().
1486: Collective on PetscViewer
1488: Input Parameters:
1489: + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1490: some related function before a call to PCLoad().
1491: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1493: Level: intermediate
1495: Notes:
1496: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1498: Notes for advanced users:
1499: Most users should not need to know the details of the binary storage
1500: format, since PCLoad() and PCView() completely hide these details.
1501: But for anyone who's interested, the standard binary matrix storage
1502: format is
1503: .vb
1504: has not yet been determined
1505: .ve
1507: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1508: @*/
1509: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1510: {
1512: PetscBool isbinary;
1513: PetscInt classid;
1514: char type[256];
1519: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1520: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1522: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1523: if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1524: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1525: PCSetType(newdm, type);
1526: if (newdm->ops->load) {
1527: (*newdm->ops->load)(newdm,viewer);
1528: }
1529: return(0);
1530: }
1532: #include <petscdraw.h>
1533: #if defined(PETSC_HAVE_SAWS)
1534: #include <petscviewersaws.h>
1535: #endif
1537: /*@C
1538: PCViewFromOptions - View from Options
1540: Collective on PC
1542: Input Parameters:
1543: + A - the PC context
1544: . obj - Optional object
1545: - name - command line option
1547: Level: intermediate
1548: .seealso: PC, PCView, PetscObjectViewFromOptions(), PCCreate()
1549: @*/
1550: PetscErrorCode PCViewFromOptions(PC A,PetscObject obj,const char name[])
1551: {
1556: PetscObjectViewFromOptions((PetscObject)A,obj,name);
1557: return(0);
1558: }
1560: /*@C
1561: PCView - Prints the PC data structure.
1563: Collective on PC
1565: Input Parameters:
1566: + PC - the PC context
1567: - viewer - optional visualization context
1569: Note:
1570: The available visualization contexts include
1571: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1572: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1573: output where only the first processor opens
1574: the file. All other processors send their
1575: data to the first processor to print.
1577: The user can open an alternative visualization contexts with
1578: PetscViewerASCIIOpen() (output to a specified file).
1580: Level: developer
1582: .seealso: KSPView(), PetscViewerASCIIOpen()
1583: @*/
1584: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1585: {
1586: PCType cstr;
1588: PetscBool iascii,isstring,isbinary,isdraw;
1589: #if defined(PETSC_HAVE_SAWS)
1590: PetscBool issaws;
1591: #endif
1595: if (!viewer) {
1596: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1597: }
1601: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1602: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1603: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1604: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1605: #if defined(PETSC_HAVE_SAWS)
1606: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1607: #endif
1609: if (iascii) {
1610: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1611: if (!pc->setupcalled) {
1612: PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");
1613: }
1614: if (pc->ops->view) {
1615: PetscViewerASCIIPushTab(viewer);
1616: (*pc->ops->view)(pc,viewer);
1617: PetscViewerASCIIPopTab(viewer);
1618: }
1619: if (pc->mat) {
1620: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1621: if (pc->pmat == pc->mat) {
1622: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1623: PetscViewerASCIIPushTab(viewer);
1624: MatView(pc->mat,viewer);
1625: PetscViewerASCIIPopTab(viewer);
1626: } else {
1627: if (pc->pmat) {
1628: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1629: } else {
1630: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1631: }
1632: PetscViewerASCIIPushTab(viewer);
1633: MatView(pc->mat,viewer);
1634: if (pc->pmat) {MatView(pc->pmat,viewer);}
1635: PetscViewerASCIIPopTab(viewer);
1636: }
1637: PetscViewerPopFormat(viewer);
1638: }
1639: } else if (isstring) {
1640: PCGetType(pc,&cstr);
1641: PetscViewerStringSPrintf(viewer," PCType: %-7.7s",cstr);
1642: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1643: if (pc->mat) {MatView(pc->mat,viewer);}
1644: if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1645: } else if (isbinary) {
1646: PetscInt classid = PC_FILE_CLASSID;
1647: MPI_Comm comm;
1648: PetscMPIInt rank;
1649: char type[256];
1651: PetscObjectGetComm((PetscObject)pc,&comm);
1652: MPI_Comm_rank(comm,&rank);
1653: if (!rank) {
1654: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT);
1655: PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1656: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR);
1657: }
1658: if (pc->ops->view) {
1659: (*pc->ops->view)(pc,viewer);
1660: }
1661: } else if (isdraw) {
1662: PetscDraw draw;
1663: char str[25];
1664: PetscReal x,y,bottom,h;
1665: PetscInt n;
1667: PetscViewerDrawGetDraw(viewer,0,&draw);
1668: PetscDrawGetCurrentPoint(draw,&x,&y);
1669: if (pc->mat) {
1670: MatGetSize(pc->mat,&n,NULL);
1671: PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1672: } else {
1673: PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1674: }
1675: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1676: bottom = y - h;
1677: PetscDrawPushCurrentPoint(draw,x,bottom);
1678: if (pc->ops->view) {
1679: (*pc->ops->view)(pc,viewer);
1680: }
1681: PetscDrawPopCurrentPoint(draw);
1682: #if defined(PETSC_HAVE_SAWS)
1683: } else if (issaws) {
1684: PetscMPIInt rank;
1686: PetscObjectName((PetscObject)pc);
1687: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1688: if (!((PetscObject)pc)->amsmem && !rank) {
1689: PetscObjectViewSAWs((PetscObject)pc,viewer);
1690: }
1691: if (pc->mat) {MatView(pc->mat,viewer);}
1692: if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1693: #endif
1694: }
1695: return(0);
1696: }
1698: /*@C
1699: PCRegister - Adds a method to the preconditioner package.
1701: Not collective
1703: Input Parameters:
1704: + name_solver - name of a new user-defined solver
1705: - routine_create - routine to create method context
1707: Notes:
1708: PCRegister() may be called multiple times to add several user-defined preconditioners.
1710: Sample usage:
1711: .vb
1712: PCRegister("my_solver", MySolverCreate);
1713: .ve
1715: Then, your solver can be chosen with the procedural interface via
1716: $ PCSetType(pc,"my_solver")
1717: or at runtime via the option
1718: $ -pc_type my_solver
1720: Level: advanced
1722: .seealso: PCRegisterAll()
1723: @*/
1724: PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1725: {
1729: PCInitializePackage();
1730: PetscFunctionListAdd(&PCList,sname,function);
1731: return(0);
1732: }
1734: static PetscErrorCode MatMult_PC(Mat A,Vec X,Vec Y)
1735: {
1736: PC pc;
1740: MatShellGetContext(A,&pc);
1741: PCApply(pc,X,Y);
1742: return(0);
1743: }
1745: /*@
1746: PCComputeOperator - Computes the explicit preconditioned operator.
1748: Collective on PC
1750: Input Parameter:
1751: + pc - the preconditioner object
1752: - mattype - the matrix type to be used for the operator
1754: Output Parameter:
1755: . mat - the explict preconditioned operator
1757: Notes:
1758: This computation is done by applying the operators to columns of the identity matrix.
1759: This routine is costly in general, and is recommended for use only with relatively small systems.
1760: Currently, this routine uses a dense matrix format when mattype == NULL
1762: Level: advanced
1764: .seealso: KSPComputeOperator(), MatType
1766: @*/
1767: PetscErrorCode PCComputeOperator(PC pc,MatType mattype,Mat *mat)
1768: {
1770: PetscInt N,M,m,n;
1771: Mat A,Apc;
1776: PCGetOperators(pc,&A,NULL);
1777: MatGetLocalSize(A,&m,&n);
1778: MatGetSize(A,&M,&N);
1779: MatCreateShell(PetscObjectComm((PetscObject)pc),m,n,M,N,pc,&Apc);
1780: MatShellSetOperation(Apc,MATOP_MULT,(void (*)(void))MatMult_PC);
1781: MatComputeOperator(Apc,mattype,mat);
1782: MatDestroy(&Apc);
1783: return(0);
1784: }
1786: /*@
1787: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1789: Collective on PC
1791: Input Parameters:
1792: + pc - the solver context
1793: . dim - the dimension of the coordinates 1, 2, or 3
1794: . nloc - the blocked size of the coordinates array
1795: - coords - the coordinates array
1797: Level: intermediate
1799: Notes:
1800: coords is an array of the dim coordinates for the nodes on
1801: the local processor, of size dim*nloc.
1802: If there are 108 equation on a processor
1803: for a displacement finite element discretization of elasticity (so
1804: that there are nloc = 36 = 108/3 nodes) then the array must have 108
1805: double precision values (ie, 3 * 36). These x y z coordinates
1806: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1807: ... , N-1.z ].
1809: .seealso: MatSetNearNullSpace()
1810: @*/
1811: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])
1812: {
1818: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1819: return(0);
1820: }
1822: /*@
1823: PCGetInterpolations - Gets interpolation matrices for all levels (except level 0)
1825: Logically Collective on PC
1827: Input Parameters:
1828: + pc - the precondition context
1830: Output Parameter:
1831: - num_levels - the number of levels
1832: . interpolations - the interpolation matrices (size of num_levels-1)
1834: Level: advanced
1836: .keywords: MG, GAMG, BoomerAMG, multigrid, interpolation, level
1838: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetInterpolation(), PCGetCoarseOperators()
1839: @*/
1840: PetscErrorCode PCGetInterpolations(PC pc,PetscInt *num_levels,Mat *interpolations[])
1841: {
1848: PetscUseMethod(pc,"PCGetInterpolations_C",(PC,PetscInt*,Mat*[]),(pc,num_levels,interpolations));
1849: return(0);
1850: }
1852: /*@
1853: PCGetCoarseOperators - Gets coarse operator matrices for all levels (except the finest level)
1855: Logically Collective on PC
1857: Input Parameters:
1858: + pc - the precondition context
1860: Output Parameter:
1861: - num_levels - the number of levels
1862: . coarseOperators - the coarse operator matrices (size of num_levels-1)
1864: Level: advanced
1866: .keywords: MG, GAMG, BoomerAMG, get, multigrid, interpolation, level
1868: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale(), PCMGGetInterpolation(), PCGetInterpolations()
1869: @*/
1870: PetscErrorCode PCGetCoarseOperators(PC pc,PetscInt *num_levels,Mat *coarseOperators[])
1871: {
1878: PetscUseMethod(pc,"PCGetCoarseOperators_C",(PC,PetscInt*,Mat*[]),(pc,num_levels,coarseOperators));
1879: return(0);
1880: }