Actual source code: precon.c
petsc-3.9.4 2018-09-11
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: 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
66: .keywords: PC, destroy
68: .seealso: PCCreate(), PCSetUp()
69: @*/
70: PetscErrorCode PCReset(PC pc)
71: {
76: if (pc->ops->reset) {
77: (*pc->ops->reset)(pc);
78: }
79: VecDestroy(&pc->diagonalscaleright);
80: VecDestroy(&pc->diagonalscaleleft);
81: MatDestroy(&pc->pmat);
82: MatDestroy(&pc->mat);
84: pc->setupcalled = 0;
85: return(0);
86: }
88: /*@
89: PCDestroy - Destroys PC context that was created with PCCreate().
91: Collective on PC
93: Input Parameter:
94: . pc - the preconditioner context
96: Level: developer
98: .keywords: PC, destroy
100: .seealso: PCCreate(), PCSetUp()
101: @*/
102: PetscErrorCode PCDestroy(PC *pc)
103: {
107: if (!*pc) return(0);
109: if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; return(0);}
111: PCReset(*pc);
113: /* if memory was published with SAWs then destroy it */
114: PetscObjectSAWsViewOff((PetscObject)*pc);
115: if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
116: DMDestroy(&(*pc)->dm);
117: PetscHeaderDestroy(pc);
118: return(0);
119: }
121: /*@C
122: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
123: scaling as needed by certain time-stepping codes.
125: Logically Collective on PC
127: Input Parameter:
128: . pc - the preconditioner context
130: Output Parameter:
131: . flag - PETSC_TRUE if it applies the scaling
133: Level: developer
135: Notes: If this returns PETSC_TRUE then the system solved via the Krylov method is
136: $ D M A D^{-1} y = D M b for left preconditioning or
137: $ D A M D^{-1} z = D b for right preconditioning
139: .keywords: PC
141: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
142: @*/
143: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
144: {
148: *flag = pc->diagonalscale;
149: return(0);
150: }
152: /*@
153: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
154: scaling as needed by certain time-stepping codes.
156: Logically Collective on PC
158: Input Parameters:
159: + pc - the preconditioner context
160: - s - scaling vector
162: Level: intermediate
164: Notes: The system solved via the Krylov method is
165: $ D M A D^{-1} y = D M b for left preconditioning or
166: $ D A M D^{-1} z = D b for right preconditioning
168: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
170: .keywords: PC
172: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
173: @*/
174: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
175: {
181: pc->diagonalscale = PETSC_TRUE;
183: PetscObjectReference((PetscObject)s);
184: VecDestroy(&pc->diagonalscaleleft);
186: pc->diagonalscaleleft = s;
188: VecDuplicate(s,&pc->diagonalscaleright);
189: VecCopy(s,pc->diagonalscaleright);
190: VecReciprocal(pc->diagonalscaleright);
191: return(0);
192: }
194: /*@
195: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
197: Logically Collective on PC
199: Input Parameters:
200: + pc - the preconditioner context
201: . in - input vector
202: + out - scaled vector (maybe the same as in)
204: Level: intermediate
206: Notes: The system solved via the Krylov method is
207: $ D M A D^{-1} y = D M b for left preconditioning or
208: $ D A M D^{-1} z = D b for right preconditioning
210: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
212: If diagonal scaling is turned off and in is not out then in is copied to out
214: .keywords: PC
216: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
217: @*/
218: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
219: {
226: if (pc->diagonalscale) {
227: VecPointwiseMult(out,pc->diagonalscaleleft,in);
228: } else if (in != out) {
229: VecCopy(in,out);
230: }
231: return(0);
232: }
234: /*@
235: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
237: Logically Collective on PC
239: Input Parameters:
240: + pc - the preconditioner context
241: . in - input vector
242: + out - scaled vector (maybe the same as in)
244: Level: intermediate
246: Notes: The system solved via the Krylov method is
247: $ D M A D^{-1} y = D M b for left preconditioning or
248: $ D A M D^{-1} z = D b for right preconditioning
250: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
252: If diagonal scaling is turned off and in is not out then in is copied to out
254: .keywords: PC
256: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
257: @*/
258: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
259: {
266: if (pc->diagonalscale) {
267: VecPointwiseMult(out,pc->diagonalscaleright,in);
268: } else if (in != out) {
269: VecCopy(in,out);
270: }
271: return(0);
272: }
274: /*@
275: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
276: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
277: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
279: Logically Collective on PC
281: Input Parameters:
282: + pc - the preconditioner context
283: - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
285: Options Database Key:
286: . -pc_use_amat <true,false>
288: Notes:
289: For the common case in which the linear system matrix and the matrix used to construct the
290: preconditioner are identical, this routine is does nothing.
292: Level: intermediate
294: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
295: @*/
296: PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
297: {
300: pc->useAmat = flg;
301: return(0);
302: }
304: /*@
305: PCSetErrorIfFailure - Causes PC to generate an error if a FPE, for example a zero pivot, is detected.
307: Logically Collective on PC
309: Input Parameters:
310: + pc - iterative context obtained from PCCreate()
311: - flg - PETSC_TRUE indicates you want the error generated
313: Level: advanced
315: Notes:
316: Normally PETSc continues if a linear solver fails due to a failed setup of a preconditioner, you can call KSPGetConvergedReason() after a KSPSolve()
317: 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
318: detected.
320: This is propagated into KSPs used by this PC, which then propagate it into PCs used by those KSPs
322: .keywords: PC, set, initial guess, nonzero
324: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
325: @*/
326: PetscErrorCode PCSetErrorIfFailure(PC pc,PetscBool flg)
327: {
331: pc->erroriffailure = flg;
332: return(0);
333: }
335: /*@
336: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
337: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
338: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
340: Logically Collective on PC
342: Input Parameter:
343: . pc - the preconditioner context
345: Output Parameter:
346: . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
348: Notes:
349: For the common case in which the linear system matrix and the matrix used to construct the
350: preconditioner are identical, this routine is does nothing.
352: Level: intermediate
354: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
355: @*/
356: PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
357: {
360: *flg = pc->useAmat;
361: return(0);
362: }
364: /*@
365: PCCreate - Creates a preconditioner context.
367: Collective on MPI_Comm
369: Input Parameter:
370: . comm - MPI communicator
372: Output Parameter:
373: . pc - location to put the preconditioner context
375: Notes:
376: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or ICC
377: in parallel. For dense matrices it is always PCNONE.
379: Level: developer
381: .keywords: PC, create, context
383: .seealso: PCSetUp(), PCApply(), PCDestroy()
384: @*/
385: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
386: {
387: PC pc;
392: *newpc = 0;
393: PCInitializePackage();
395: PetscHeaderCreate(pc,PC_CLASSID,"PC","Preconditioner","PC",comm,PCDestroy,PCView);
397: pc->mat = 0;
398: pc->pmat = 0;
399: pc->setupcalled = 0;
400: pc->setfromoptionscalled = 0;
401: pc->data = 0;
402: pc->diagonalscale = PETSC_FALSE;
403: pc->diagonalscaleleft = 0;
404: pc->diagonalscaleright = 0;
406: pc->modifysubmatrices = 0;
407: pc->modifysubmatricesP = 0;
409: *newpc = pc;
410: return(0);
412: }
414: /* -------------------------------------------------------------------------------*/
416: /*@
417: PCApply - Applies the preconditioner to a vector.
419: Collective on PC and Vec
421: Input Parameters:
422: + pc - the preconditioner context
423: - x - input vector
425: Output Parameter:
426: . y - output vector
428: Level: developer
430: .keywords: PC, apply
432: .seealso: PCApplyTranspose(), PCApplyBAorAB()
433: @*/
434: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
435: {
437: PetscInt m,n,mv,nv;
443: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
444: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
445: /* use pmat to check vector sizes since for KSPLQR the pmat may be of a different size than mat */
446: MatGetLocalSize(pc->pmat,&m,&n);
447: VecGetLocalSize(x,&nv);
448: VecGetLocalSize(y,&mv);
449: 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);
450: 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);
451: VecLocked(y,3);
453: PCSetUp(pc);
454: if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
455: VecLockPush(x);
456: PetscLogEventBegin(PC_Apply,pc,x,y,0);
457: (*pc->ops->apply)(pc,x,y);
458: PetscLogEventEnd(PC_Apply,pc,x,y,0);
459: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
460: VecLockPop(x);
461: return(0);
462: }
464: /*@
465: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
467: Collective on PC and Vec
469: Input Parameters:
470: + pc - the preconditioner context
471: - x - input vector
473: Output Parameter:
474: . y - output vector
476: Notes:
477: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
479: Level: developer
481: .keywords: PC, apply, symmetric, left
483: .seealso: PCApply(), PCApplySymmetricRight()
484: @*/
485: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
486: {
493: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
494: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
495: PCSetUp(pc);
496: if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
497: VecLockPush(x);
498: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
499: (*pc->ops->applysymmetricleft)(pc,x,y);
500: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
501: VecLockPop(x);
502: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
503: return(0);
504: }
506: /*@
507: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
509: Collective on PC and Vec
511: Input Parameters:
512: + pc - the preconditioner context
513: - x - input vector
515: Output Parameter:
516: . y - output vector
518: Level: developer
520: Notes:
521: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
523: .keywords: PC, apply, symmetric, right
525: .seealso: PCApply(), PCApplySymmetricLeft()
526: @*/
527: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
528: {
535: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
536: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
537: PCSetUp(pc);
538: if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
539: VecLockPush(x);
540: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
541: (*pc->ops->applysymmetricright)(pc,x,y);
542: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
543: VecLockPop(x);
544: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
545: return(0);
546: }
548: /*@
549: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
551: Collective on PC and Vec
553: Input Parameters:
554: + pc - the preconditioner context
555: - x - input vector
557: Output Parameter:
558: . y - output vector
560: Notes: For complex numbers this applies the non-Hermitian transpose.
562: Developer Notes: We need to implement a PCApplyHermitianTranspose()
564: Level: developer
566: .keywords: PC, apply, transpose
568: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
569: @*/
570: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
571: {
578: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
579: if (pc->erroriffailure) {VecValidValues(x,2,PETSC_TRUE);}
580: PCSetUp(pc);
581: if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
582: VecLockPush(x);
583: PetscLogEventBegin(PC_Apply,pc,x,y,0);
584: (*pc->ops->applytranspose)(pc,x,y);
585: PetscLogEventEnd(PC_Apply,pc,x,y,0);
586: VecLockPop(x);
587: if (pc->erroriffailure) {VecValidValues(y,3,PETSC_FALSE);}
588: return(0);
589: }
591: /*@
592: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
594: Collective on PC and Vec
596: Input Parameters:
597: . pc - the preconditioner context
599: Output Parameter:
600: . flg - PETSC_TRUE if a transpose operation is defined
602: Level: developer
604: .keywords: PC, apply, transpose
606: .seealso: PCApplyTranspose()
607: @*/
608: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
609: {
613: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
614: else *flg = PETSC_FALSE;
615: return(0);
616: }
618: /*@
619: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
621: Collective on PC and Vec
623: Input Parameters:
624: + pc - the preconditioner context
625: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
626: . x - input vector
627: - work - work vector
629: Output Parameter:
630: . y - output vector
632: Level: developer
634: 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
635: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
637: .keywords: PC, apply, operator
639: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
640: @*/
641: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
642: {
650: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
651: 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");
652: if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
653: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
655: PCSetUp(pc);
656: if (pc->diagonalscale) {
657: if (pc->ops->applyBA) {
658: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
659: VecDuplicate(x,&work2);
660: PCDiagonalScaleRight(pc,x,work2);
661: (*pc->ops->applyBA)(pc,side,work2,y,work);
662: PCDiagonalScaleLeft(pc,y,y);
663: VecDestroy(&work2);
664: } else if (side == PC_RIGHT) {
665: PCDiagonalScaleRight(pc,x,y);
666: PCApply(pc,y,work);
667: MatMult(pc->mat,work,y);
668: PCDiagonalScaleLeft(pc,y,y);
669: } else if (side == PC_LEFT) {
670: PCDiagonalScaleRight(pc,x,y);
671: MatMult(pc->mat,y,work);
672: PCApply(pc,work,y);
673: PCDiagonalScaleLeft(pc,y,y);
674: } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
675: } else {
676: if (pc->ops->applyBA) {
677: (*pc->ops->applyBA)(pc,side,x,y,work);
678: } else if (side == PC_RIGHT) {
679: PCApply(pc,x,work);
680: MatMult(pc->mat,work,y);
681: } else if (side == PC_LEFT) {
682: MatMult(pc->mat,x,work);
683: PCApply(pc,work,y);
684: } else if (side == PC_SYMMETRIC) {
685: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
686: PCApplySymmetricRight(pc,x,work);
687: MatMult(pc->mat,work,y);
688: VecCopy(y,work);
689: PCApplySymmetricLeft(pc,work,y);
690: }
691: }
692: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
693: return(0);
694: }
696: /*@
697: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
698: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
699: NOT tr(B*A) = tr(A)*tr(B).
701: Collective on PC and Vec
703: Input Parameters:
704: + pc - the preconditioner context
705: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
706: . x - input vector
707: - work - work vector
709: Output Parameter:
710: . y - output vector
713: 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
714: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
716: Level: developer
718: .keywords: PC, apply, operator, transpose
720: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
721: @*/
722: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
723: {
731: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
732: if (pc->erroriffailure) {VecValidValues(x,3,PETSC_TRUE);}
733: if (pc->ops->applyBAtranspose) {
734: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
735: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
736: return(0);
737: }
738: if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
740: PCSetUp(pc);
741: if (side == PC_RIGHT) {
742: PCApplyTranspose(pc,x,work);
743: MatMultTranspose(pc->mat,work,y);
744: } else if (side == PC_LEFT) {
745: MatMultTranspose(pc->mat,x,work);
746: PCApplyTranspose(pc,work,y);
747: }
748: /* add support for PC_SYMMETRIC */
749: if (pc->erroriffailure) {VecValidValues(y,4,PETSC_FALSE);}
750: return(0);
751: }
753: /* -------------------------------------------------------------------------------*/
755: /*@
756: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
757: built-in fast application of Richardson's method.
759: Not Collective
761: Input Parameter:
762: . pc - the preconditioner
764: Output Parameter:
765: . exists - PETSC_TRUE or PETSC_FALSE
767: Level: developer
769: .keywords: PC, apply, Richardson, exists
771: .seealso: PCApplyRichardson()
772: @*/
773: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
774: {
778: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
779: else *exists = PETSC_FALSE;
780: return(0);
781: }
783: /*@
784: PCApplyRichardson - Applies several steps of Richardson iteration with
785: the particular preconditioner. This routine is usually used by the
786: Krylov solvers and not the application code directly.
788: Collective on PC
790: Input Parameters:
791: + pc - the preconditioner context
792: . b - the right hand side
793: . w - one work vector
794: . rtol - relative decrease in residual norm convergence criteria
795: . abstol - absolute residual norm convergence criteria
796: . dtol - divergence residual norm increase criteria
797: . its - the number of iterations to apply.
798: - guesszero - if the input x contains nonzero initial guess
800: Output Parameter:
801: + outits - number of iterations actually used (for SOR this always equals its)
802: . reason - the reason the apply terminated
803: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
805: Notes:
806: Most preconditioners do not support this function. Use the command
807: PCApplyRichardsonExists() to determine if one does.
809: Except for the multigrid PC this routine ignores the convergence tolerances
810: and always runs for the number of iterations
812: Level: developer
814: .keywords: PC, apply, Richardson
816: .seealso: PCApplyRichardsonExists()
817: @*/
818: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
819: {
827: if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
828: PCSetUp(pc);
829: if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
830: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
831: return(0);
832: }
834: /*@
835: PCGetSetUpFailedReason - Gets the reason a PCSetUp() failed or 0 if it did not fail
837: Logically Collective on PC
839: Input Parameter:
840: . pc - the preconditioner context
842: Output Parameter:
843: . reason - the reason it failed, currently only -1
845: Level: advanced
847: .keywords: PC, setup
849: .seealso: PCCreate(), PCApply(), PCDestroy()
850: @*/
851: PetscErrorCode PCGetSetUpFailedReason(PC pc,PCFailedReason *reason)
852: {
854: if (pc->setupcalled < 0) *reason = (PCFailedReason)pc->setupcalled;
855: else *reason = pc->failedreason;
856: return(0);
857: }
860: /*
861: a setupcall of 0 indicates never setup,
862: 1 indicates has been previously setup
863: -1 indicates a PCSetUp() was attempted and failed
864: */
865: /*@
866: PCSetUp - Prepares for the use of a preconditioner.
868: Collective on PC
870: Input Parameter:
871: . pc - the preconditioner context
873: Level: developer
875: .keywords: PC, setup
877: .seealso: PCCreate(), PCApply(), PCDestroy()
878: @*/
879: PetscErrorCode PCSetUp(PC pc)
880: {
881: PetscErrorCode ierr;
882: const char *def;
883: PetscObjectState matstate, matnonzerostate;
887: if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
889: if (pc->setupcalled && pc->reusepreconditioner) {
890: PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
891: return(0);
892: }
894: PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
895: MatGetNonzeroState(pc->pmat,&matnonzerostate);
896: if (!pc->setupcalled) {
897: PetscInfo(pc,"Setting up PC for first time\n");
898: pc->flag = DIFFERENT_NONZERO_PATTERN;
899: } else if (matstate == pc->matstate) {
900: PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
901: return(0);
902: } else {
903: if (matnonzerostate > pc->matnonzerostate) {
904: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
905: pc->flag = DIFFERENT_NONZERO_PATTERN;
906: } else {
907: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
908: pc->flag = SAME_NONZERO_PATTERN;
909: }
910: }
911: pc->matstate = matstate;
912: pc->matnonzerostate = matnonzerostate;
914: if (!((PetscObject)pc)->type_name) {
915: PCGetDefaultType_Private(pc,&def);
916: PCSetType(pc,def);
917: }
919: MatSetErrorIfFailure(pc->pmat,pc->erroriffailure);
920: MatSetErrorIfFailure(pc->mat,pc->erroriffailure);
921: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
922: if (pc->ops->setup) {
923: (*pc->ops->setup)(pc);
924: }
925: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
926: if (!pc->setupcalled) pc->setupcalled = 1;
927: return(0);
928: }
930: /*@
931: PCSetUpOnBlocks - Sets up the preconditioner for each block in
932: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
933: methods.
935: Collective on PC
937: Input Parameters:
938: . pc - the preconditioner context
940: Level: developer
942: .keywords: PC, setup, blocks
944: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
945: @*/
946: PetscErrorCode PCSetUpOnBlocks(PC pc)
947: {
952: if (!pc->ops->setuponblocks) return(0);
953: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
954: (*pc->ops->setuponblocks)(pc);
955: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
956: return(0);
957: }
959: /*@C
960: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
961: submatrices that arise within certain subdomain-based preconditioners.
962: The basic submatrices are extracted from the preconditioner matrix as
963: usual; the user can then alter these (for example, to set different boundary
964: conditions for each submatrix) before they are used for the local solves.
966: Logically Collective on PC
968: Input Parameters:
969: + pc - the preconditioner context
970: . func - routine for modifying the submatrices
971: - ctx - optional user-defined context (may be null)
973: Calling sequence of func:
974: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
976: . row - an array of index sets that contain the global row numbers
977: that comprise each local submatrix
978: . col - an array of index sets that contain the global column numbers
979: that comprise each local submatrix
980: . submat - array of local submatrices
981: - ctx - optional user-defined context for private data for the
982: user-defined func routine (may be null)
984: Notes:
985: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
986: KSPSolve().
988: A routine set by PCSetModifySubMatrices() is currently called within
989: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
990: preconditioners. All other preconditioners ignore this routine.
992: Level: advanced
994: .keywords: PC, set, modify, submatrices
996: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
997: @*/
998: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
999: {
1002: pc->modifysubmatrices = func;
1003: pc->modifysubmatricesP = ctx;
1004: return(0);
1005: }
1007: /*@C
1008: PCModifySubMatrices - Calls an optional user-defined routine within
1009: certain preconditioners if one has been set with PCSetModifySubMarices().
1011: Collective on PC
1013: Input Parameters:
1014: + pc - the preconditioner context
1015: . nsub - the number of local submatrices
1016: . row - an array of index sets that contain the global row numbers
1017: that comprise each local submatrix
1018: . col - an array of index sets that contain the global column numbers
1019: that comprise each local submatrix
1020: . submat - array of local submatrices
1021: - ctx - optional user-defined context for private data for the
1022: user-defined routine (may be null)
1024: Output Parameter:
1025: . submat - array of local submatrices (the entries of which may
1026: have been modified)
1028: Notes:
1029: The user should NOT generally call this routine, as it will
1030: automatically be called within certain preconditioners (currently
1031: block Jacobi, additive Schwarz) if set.
1033: The basic submatrices are extracted from the preconditioner matrix
1034: as usual; the user can then alter these (for example, to set different
1035: boundary conditions for each submatrix) before they are used for the
1036: local solves.
1038: Level: developer
1040: .keywords: PC, modify, submatrices
1042: .seealso: PCSetModifySubMatrices()
1043: @*/
1044: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1045: {
1050: if (!pc->modifysubmatrices) return(0);
1051: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1052: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1053: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1054: return(0);
1055: }
1057: /*@
1058: PCSetOperators - Sets the matrix associated with the linear system and
1059: a (possibly) different one associated with the preconditioner.
1061: Logically Collective on PC and Mat
1063: Input Parameters:
1064: + pc - the preconditioner context
1065: . Amat - the matrix that defines the linear system
1066: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1068: Notes:
1069: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1071: If you wish to replace either Amat or Pmat but leave the other one untouched then
1072: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1073: on it and then pass it back in in your call to KSPSetOperators().
1075: More Notes about Repeated Solution of Linear Systems:
1076: PETSc does NOT reset the matrix entries of either Amat or Pmat
1077: to zero after a linear solve; the user is completely responsible for
1078: matrix assembly. See the routine MatZeroEntries() if desiring to
1079: zero all elements of a matrix.
1081: Level: intermediate
1083: .keywords: PC, set, operators, matrix, linear system
1085: .seealso: PCGetOperators(), MatZeroEntries()
1086: @*/
1087: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1088: {
1089: PetscErrorCode ierr;
1090: PetscInt m1,n1,m2,n2;
1098: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1099: MatGetLocalSize(Amat,&m1,&n1);
1100: MatGetLocalSize(pc->mat,&m2,&n2);
1101: 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);
1102: MatGetLocalSize(Pmat,&m1,&n1);
1103: MatGetLocalSize(pc->pmat,&m2,&n2);
1104: 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);
1105: }
1107: if (Pmat != pc->pmat) {
1108: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1109: pc->matnonzerostate = -1;
1110: pc->matstate = -1;
1111: }
1113: /* reference first in case the matrices are the same */
1114: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1115: MatDestroy(&pc->mat);
1116: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1117: MatDestroy(&pc->pmat);
1118: pc->mat = Amat;
1119: pc->pmat = Pmat;
1120: return(0);
1121: }
1123: /*@
1124: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1126: Logically Collective on PC
1128: Input Parameters:
1129: + pc - the preconditioner context
1130: - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1132: Level: intermediate
1134: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1135: @*/
1136: PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1137: {
1140: pc->reusepreconditioner = flag;
1141: return(0);
1142: }
1144: /*@
1145: PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1147: Not Collective
1149: Input Parameter:
1150: . pc - the preconditioner context
1152: Output Parameter:
1153: . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1155: Level: intermediate
1157: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1158: @*/
1159: PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1160: {
1163: *flag = pc->reusepreconditioner;
1164: return(0);
1165: }
1167: /*@C
1168: PCGetOperators - Gets the matrix associated with the linear system and
1169: possibly a different one associated with the preconditioner.
1171: Not collective, though parallel Mats are returned if the PC is parallel
1173: Input Parameter:
1174: . pc - the preconditioner context
1176: Output Parameters:
1177: + Amat - the matrix defining the linear system
1178: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1180: Level: intermediate
1182: Notes: Does not increase the reference count of the matrices, so you should not destroy them
1184: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1185: are created in PC and returned to the user. In this case, if both operators
1186: mat and pmat are requested, two DIFFERENT operators will be returned. If
1187: only one is requested both operators in the PC will be the same (i.e. as
1188: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1189: The user must set the sizes of the returned matrices and their type etc just
1190: as if the user created them with MatCreate(). For example,
1192: $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1193: $ set size, type, etc of Amat
1195: $ MatCreate(comm,&mat);
1196: $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1197: $ PetscObjectDereference((PetscObject)mat);
1198: $ set size, type, etc of Amat
1200: and
1202: $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1203: $ set size, type, etc of Amat and Pmat
1205: $ MatCreate(comm,&Amat);
1206: $ MatCreate(comm,&Pmat);
1207: $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1208: $ PetscObjectDereference((PetscObject)Amat);
1209: $ PetscObjectDereference((PetscObject)Pmat);
1210: $ set size, type, etc of Amat and Pmat
1212: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1213: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1214: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1215: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1216: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1217: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1218: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1219: it can be created for you?
1222: .keywords: PC, get, operators, matrix, linear system
1224: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1225: @*/
1226: PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1227: {
1232: if (Amat) {
1233: if (!pc->mat) {
1234: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1235: pc->mat = pc->pmat;
1236: PetscObjectReference((PetscObject)pc->mat);
1237: } else { /* both Amat and Pmat are empty */
1238: MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1239: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1240: pc->pmat = pc->mat;
1241: PetscObjectReference((PetscObject)pc->pmat);
1242: }
1243: }
1244: }
1245: *Amat = pc->mat;
1246: }
1247: if (Pmat) {
1248: if (!pc->pmat) {
1249: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1250: pc->pmat = pc->mat;
1251: PetscObjectReference((PetscObject)pc->pmat);
1252: } else {
1253: MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1254: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1255: pc->mat = pc->pmat;
1256: PetscObjectReference((PetscObject)pc->mat);
1257: }
1258: }
1259: }
1260: *Pmat = pc->pmat;
1261: }
1262: return(0);
1263: }
1265: /*@C
1266: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1267: possibly a different one associated with the preconditioner have been set in the PC.
1269: Not collective, though the results on all processes should be the same
1271: Input Parameter:
1272: . pc - the preconditioner context
1274: Output Parameters:
1275: + mat - the matrix associated with the linear system was set
1276: - pmat - matrix associated with the preconditioner was set, usually the same
1278: Level: intermediate
1280: .keywords: PC, get, operators, matrix, linear system
1282: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1283: @*/
1284: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1285: {
1288: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1289: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1290: return(0);
1291: }
1293: /*@
1294: PCFactorGetMatrix - Gets the factored matrix from the
1295: preconditioner context. This routine is valid only for the LU,
1296: incomplete LU, Cholesky, and incomplete Cholesky methods.
1298: Not Collective on PC though Mat is parallel if PC is parallel
1300: Input Parameters:
1301: . pc - the preconditioner context
1303: Output parameters:
1304: . mat - the factored matrix
1306: Level: advanced
1308: Notes: Does not increase the reference count for the matrix so DO NOT destroy it
1310: .keywords: PC, get, factored, matrix
1311: @*/
1312: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1313: {
1319: if (pc->ops->getfactoredmatrix) {
1320: (*pc->ops->getfactoredmatrix)(pc,mat);
1321: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1322: return(0);
1323: }
1325: /*@C
1326: PCSetOptionsPrefix - Sets the prefix used for searching for all
1327: PC options in the database.
1329: Logically Collective on PC
1331: Input Parameters:
1332: + pc - the preconditioner context
1333: - prefix - the prefix string to prepend to all PC option requests
1335: Notes:
1336: A hyphen (-) must NOT be given at the beginning of the prefix name.
1337: The first character of all runtime options is AUTOMATICALLY the
1338: hyphen.
1340: Level: advanced
1342: .keywords: PC, set, options, prefix, database
1344: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1345: @*/
1346: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1347: {
1352: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1353: return(0);
1354: }
1356: /*@C
1357: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1358: PC options in the database.
1360: Logically Collective on PC
1362: Input Parameters:
1363: + pc - the preconditioner context
1364: - prefix - the prefix string to prepend to all PC option requests
1366: Notes:
1367: A hyphen (-) must NOT be given at the beginning of the prefix name.
1368: The first character of all runtime options is AUTOMATICALLY the
1369: hyphen.
1371: Level: advanced
1373: .keywords: PC, append, options, prefix, database
1375: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1376: @*/
1377: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1378: {
1383: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1384: return(0);
1385: }
1387: /*@C
1388: PCGetOptionsPrefix - Gets the prefix used for searching for all
1389: PC options in the database.
1391: Not Collective
1393: Input Parameters:
1394: . pc - the preconditioner context
1396: Output Parameters:
1397: . prefix - pointer to the prefix string used, is returned
1399: Notes: On the fortran side, the user should pass in a string 'prifix' of
1400: sufficient length to hold the prefix.
1402: Level: advanced
1404: .keywords: PC, get, options, prefix, database
1406: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1407: @*/
1408: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1409: {
1415: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1416: return(0);
1417: }
1419: /*
1420: Indicates the right hand side will be changed by KSPSolve(), this occurs for a few
1421: preconditioners including BDDC and Eisentat that transform the equations before applying
1422: the Krylov methods
1423: */
1424: PETSC_INTERN PetscErrorCode PCPreSolveChangeRHS(PC pc,PetscBool *change)
1425: {
1431: *change = PETSC_FALSE;
1432: PetscTryMethod(pc,"PCPreSolveChangeRHS_C",(PC,PetscBool*),(pc,change));
1433: return(0);
1434: }
1436: /*@
1437: PCPreSolve - Optional pre-solve phase, intended for any
1438: preconditioner-specific actions that must be performed before
1439: the iterative solve itself.
1441: Collective on PC
1443: Input Parameters:
1444: + pc - the preconditioner context
1445: - ksp - the Krylov subspace context
1447: Level: developer
1449: Sample of Usage:
1450: .vb
1451: PCPreSolve(pc,ksp);
1452: KSPSolve(ksp,b,x);
1453: PCPostSolve(pc,ksp);
1454: .ve
1456: Notes:
1457: The pre-solve phase is distinct from the PCSetUp() phase.
1459: KSPSolve() calls this directly, so is rarely called by the user.
1461: .keywords: PC, pre-solve
1463: .seealso: PCPostSolve()
1464: @*/
1465: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1466: {
1468: Vec x,rhs;
1473: pc->presolvedone++;
1474: if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1475: KSPGetSolution(ksp,&x);
1476: KSPGetRhs(ksp,&rhs);
1478: if (pc->ops->presolve) {
1479: (*pc->ops->presolve)(pc,ksp,rhs,x);
1480: }
1481: return(0);
1482: }
1484: /*@
1485: PCPostSolve - Optional post-solve phase, intended for any
1486: preconditioner-specific actions that must be performed after
1487: the iterative solve itself.
1489: Collective on PC
1491: Input Parameters:
1492: + pc - the preconditioner context
1493: - ksp - the Krylov subspace context
1495: Sample of Usage:
1496: .vb
1497: PCPreSolve(pc,ksp);
1498: KSPSolve(ksp,b,x);
1499: PCPostSolve(pc,ksp);
1500: .ve
1502: Note:
1503: KSPSolve() calls this routine directly, so it is rarely called by the user.
1505: Level: developer
1507: .keywords: PC, post-solve
1509: .seealso: PCPreSolve(), KSPSolve()
1510: @*/
1511: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1512: {
1514: Vec x,rhs;
1519: pc->presolvedone--;
1520: KSPGetSolution(ksp,&x);
1521: KSPGetRhs(ksp,&rhs);
1522: if (pc->ops->postsolve) {
1523: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1524: }
1525: return(0);
1526: }
1528: /*@C
1529: PCLoad - Loads a PC that has been stored in binary with PCView().
1531: Collective on PetscViewer
1533: Input Parameters:
1534: + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1535: some related function before a call to PCLoad().
1536: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1538: Level: intermediate
1540: Notes:
1541: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1543: Notes for advanced users:
1544: Most users should not need to know the details of the binary storage
1545: format, since PCLoad() and PCView() completely hide these details.
1546: But for anyone who's interested, the standard binary matrix storage
1547: format is
1548: .vb
1549: has not yet been determined
1550: .ve
1552: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1553: @*/
1554: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1555: {
1557: PetscBool isbinary;
1558: PetscInt classid;
1559: char type[256];
1564: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1565: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1567: PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
1568: if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1569: PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
1570: PCSetType(newdm, type);
1571: if (newdm->ops->load) {
1572: (*newdm->ops->load)(newdm,viewer);
1573: }
1574: return(0);
1575: }
1577: #include <petscdraw.h>
1578: #if defined(PETSC_HAVE_SAWS)
1579: #include <petscviewersaws.h>
1580: #endif
1581: /*@C
1582: PCView - Prints the PC data structure.
1584: Collective on PC
1586: Input Parameters:
1587: + PC - the PC context
1588: - viewer - optional visualization context
1590: Note:
1591: The available visualization contexts include
1592: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1593: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1594: output where only the first processor opens
1595: the file. All other processors send their
1596: data to the first processor to print.
1598: The user can open an alternative visualization contexts with
1599: PetscViewerASCIIOpen() (output to a specified file).
1601: Level: developer
1603: .keywords: PC, view
1605: .seealso: KSPView(), PetscViewerASCIIOpen()
1606: @*/
1607: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1608: {
1609: PCType cstr;
1611: PetscBool iascii,isstring,isbinary,isdraw;
1612: #if defined(PETSC_HAVE_SAWS)
1613: PetscBool issaws;
1614: #endif
1618: if (!viewer) {
1619: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1620: }
1624: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1625: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1626: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1627: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1628: #if defined(PETSC_HAVE_SAWS)
1629: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&issaws);
1630: #endif
1632: if (iascii) {
1633: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1634: if (!pc->setupcalled) {
1635: PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");
1636: }
1637: if (pc->ops->view) {
1638: PetscViewerASCIIPushTab(viewer);
1639: (*pc->ops->view)(pc,viewer);
1640: PetscViewerASCIIPopTab(viewer);
1641: }
1642: if (pc->mat) {
1643: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1644: if (pc->pmat == pc->mat) {
1645: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1646: PetscViewerASCIIPushTab(viewer);
1647: MatView(pc->mat,viewer);
1648: PetscViewerASCIIPopTab(viewer);
1649: } else {
1650: if (pc->pmat) {
1651: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1652: } else {
1653: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1654: }
1655: PetscViewerASCIIPushTab(viewer);
1656: MatView(pc->mat,viewer);
1657: if (pc->pmat) {MatView(pc->pmat,viewer);}
1658: PetscViewerASCIIPopTab(viewer);
1659: }
1660: PetscViewerPopFormat(viewer);
1661: }
1662: } else if (isstring) {
1663: PCGetType(pc,&cstr);
1664: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1665: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1666: } else if (isbinary) {
1667: PetscInt classid = PC_FILE_CLASSID;
1668: MPI_Comm comm;
1669: PetscMPIInt rank;
1670: char type[256];
1672: PetscObjectGetComm((PetscObject)pc,&comm);
1673: MPI_Comm_rank(comm,&rank);
1674: if (!rank) {
1675: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1676: PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1677: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1678: }
1679: if (pc->ops->view) {
1680: (*pc->ops->view)(pc,viewer);
1681: }
1682: } else if (isdraw) {
1683: PetscDraw draw;
1684: char str[25];
1685: PetscReal x,y,bottom,h;
1686: PetscInt n;
1688: PetscViewerDrawGetDraw(viewer,0,&draw);
1689: PetscDrawGetCurrentPoint(draw,&x,&y);
1690: if (pc->mat) {
1691: MatGetSize(pc->mat,&n,NULL);
1692: PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1693: } else {
1694: PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1695: }
1696: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1697: bottom = y - h;
1698: PetscDrawPushCurrentPoint(draw,x,bottom);
1699: if (pc->ops->view) {
1700: (*pc->ops->view)(pc,viewer);
1701: }
1702: PetscDrawPopCurrentPoint(draw);
1703: #if defined(PETSC_HAVE_SAWS)
1704: } else if (issaws) {
1705: PetscMPIInt rank;
1707: PetscObjectName((PetscObject)pc);
1708: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1709: if (!((PetscObject)pc)->amsmem && !rank) {
1710: PetscObjectViewSAWs((PetscObject)pc,viewer);
1711: }
1712: if (pc->mat) {MatView(pc->mat,viewer);}
1713: if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1714: #endif
1715: }
1716: return(0);
1717: }
1719: /*@C
1720: PCRegister - Adds a method to the preconditioner package.
1722: Not collective
1724: Input Parameters:
1725: + name_solver - name of a new user-defined solver
1726: - routine_create - routine to create method context
1728: Notes:
1729: PCRegister() may be called multiple times to add several user-defined preconditioners.
1731: Sample usage:
1732: .vb
1733: PCRegister("my_solver", MySolverCreate);
1734: .ve
1736: Then, your solver can be chosen with the procedural interface via
1737: $ PCSetType(pc,"my_solver")
1738: or at runtime via the option
1739: $ -pc_type my_solver
1741: Level: advanced
1743: .keywords: PC, register
1745: .seealso: PCRegisterAll(), PCRegisterDestroy()
1746: @*/
1747: PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1748: {
1752: PetscFunctionListAdd(&PCList,sname,function);
1753: return(0);
1754: }
1756: /*@
1757: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1759: Collective on PC
1761: Input Parameter:
1762: . pc - the preconditioner object
1764: Output Parameter:
1765: . mat - the explict preconditioned operator
1767: Notes:
1768: This computation is done by applying the operators to columns of the
1769: identity matrix.
1771: Currently, this routine uses a dense matrix format when 1 processor
1772: is used and a sparse format otherwise. This routine is costly in general,
1773: and is recommended for use only with relatively small systems.
1775: Level: advanced
1777: .keywords: PC, compute, explicit, operator
1779: .seealso: KSPComputeExplicitOperator()
1781: @*/
1782: PetscErrorCode PCComputeExplicitOperator(PC pc,Mat *mat)
1783: {
1784: Vec in,out;
1786: PetscInt i,M,m,*rows,start,end;
1787: PetscMPIInt size;
1788: MPI_Comm comm;
1789: PetscScalar *array,one = 1.0;
1795: PetscObjectGetComm((PetscObject)pc,&comm);
1796: MPI_Comm_size(comm,&size);
1798: if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1799: MatCreateVecs(pc->pmat,&in,0);
1800: VecDuplicate(in,&out);
1801: VecGetOwnershipRange(in,&start,&end);
1802: VecGetSize(in,&M);
1803: VecGetLocalSize(in,&m);
1804: PetscMalloc1(m+1,&rows);
1805: for (i=0; i<m; i++) rows[i] = start + i;
1807: MatCreate(comm,mat);
1808: MatSetSizes(*mat,m,m,M,M);
1809: if (size == 1) {
1810: MatSetType(*mat,MATSEQDENSE);
1811: MatSeqDenseSetPreallocation(*mat,NULL);
1812: } else {
1813: MatSetType(*mat,MATMPIAIJ);
1814: MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1815: }
1816: MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
1818: for (i=0; i<M; i++) {
1820: VecSet(in,0.0);
1821: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1822: VecAssemblyBegin(in);
1823: VecAssemblyEnd(in);
1825: /* should fix, allowing user to choose side */
1826: PCApply(pc,in,out);
1828: VecGetArray(out,&array);
1829: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1830: VecRestoreArray(out,&array);
1832: }
1833: PetscFree(rows);
1834: VecDestroy(&out);
1835: VecDestroy(&in);
1836: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1837: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1838: return(0);
1839: }
1841: /*@
1842: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1844: Collective on PC
1846: Input Parameters:
1847: + pc - the solver context
1848: . dim - the dimension of the coordinates 1, 2, or 3
1849: - coords - the coordinates
1851: Level: intermediate
1853: Notes: coords is an array of the 3D coordinates for the nodes on
1854: the local processor. So if there are 108 equation on a processor
1855: for a displacement finite element discretization of elasticity (so
1856: that there are 36 = 108/3 nodes) then the array must have 108
1857: double precision values (ie, 3 * 36). These x y z coordinates
1858: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1859: ... , N-1.z ].
1861: .seealso: MatSetNearNullSpace()
1862: @*/
1863: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1864: {
1870: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1871: return(0);
1872: }