Actual source code: precon.c
petsc-3.5.4 2015-05-23
2: /*
3: The PC (preconditioner) interface routines, callable by users.
4: */
5: #include <petsc-private/pcimpl.h> /*I "petscksp.h" I*/
6: #include <petscdm.h>
8: /* Logging support */
9: PetscClassId PC_CLASSID;
10: PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
11: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks, PC_ApplyOnMproc;
15: PetscErrorCode PCGetDefaultType_Private(PC pc,const char *type[])
16: {
18: PetscMPIInt size;
19: PetscBool flg1,flg2,set,flg3;
22: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
23: if (pc->pmat) {
24: PetscErrorCode (*f)(Mat,MatReuse,Mat*);
25: PetscObjectQueryFunction((PetscObject)pc->pmat,"MatGetDiagonalBlock_C",&f);
26: if (size == 1) {
27: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
28: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
29: MatIsSymmetricKnown(pc->pmat,&set,&flg3);
30: if (flg1 && (!flg2 || (set && flg3))) {
31: *type = PCICC;
32: } else if (flg2) {
33: *type = PCILU;
34: } else if (f) { /* likely is a parallel matrix run on one processor */
35: *type = PCBJACOBI;
36: } else {
37: *type = PCNONE;
38: }
39: } else {
40: if (f) {
41: *type = PCBJACOBI;
42: } else {
43: *type = PCNONE;
44: }
45: }
46: } else {
47: if (size == 1) {
48: *type = PCILU;
49: } else {
50: *type = PCBJACOBI;
51: }
52: }
53: return(0);
54: }
58: /*@
59: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
61: Collective on PC
63: Input Parameter:
64: . pc - the preconditioner context
66: Level: developer
68: Notes: This allows a PC to be reused for a different sized linear system but using the same options that have been previously set in the PC
70: .keywords: PC, destroy
72: .seealso: PCCreate(), PCSetUp()
73: @*/
74: PetscErrorCode PCReset(PC pc)
75: {
80: if (pc->ops->reset) {
81: (*pc->ops->reset)(pc);
82: }
83: VecDestroy(&pc->diagonalscaleright);
84: VecDestroy(&pc->diagonalscaleleft);
85: MatDestroy(&pc->pmat);
86: MatDestroy(&pc->mat);
88: pc->setupcalled = 0;
89: return(0);
90: }
94: /*@
95: PCDestroy - Destroys PC context that was created with PCCreate().
97: Collective on PC
99: Input Parameter:
100: . pc - the preconditioner context
102: Level: developer
104: .keywords: PC, destroy
106: .seealso: PCCreate(), PCSetUp()
107: @*/
108: PetscErrorCode PCDestroy(PC *pc)
109: {
113: if (!*pc) return(0);
115: if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; return(0);}
117: PCReset(*pc);
119: /* if memory was published with SAWs then destroy it */
120: PetscObjectSAWsViewOff((PetscObject)*pc);
121: if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
122: DMDestroy(&(*pc)->dm);
123: PetscHeaderDestroy(pc);
124: return(0);
125: }
129: /*@C
130: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
131: scaling as needed by certain time-stepping codes.
133: Logically Collective on PC
135: Input Parameter:
136: . pc - the preconditioner context
138: Output Parameter:
139: . flag - PETSC_TRUE if it applies the scaling
141: Level: developer
143: Notes: If this returns PETSC_TRUE then the system solved via the Krylov method is
144: $ D M A D^{-1} y = D M b for left preconditioning or
145: $ D A M D^{-1} z = D b for right preconditioning
147: .keywords: PC
149: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
150: @*/
151: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
152: {
156: *flag = pc->diagonalscale;
157: return(0);
158: }
162: /*@
163: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
164: scaling as needed by certain time-stepping codes.
166: Logically Collective on PC
168: Input Parameters:
169: + pc - the preconditioner context
170: - s - scaling vector
172: Level: intermediate
174: Notes: The system solved via the Krylov method is
175: $ D M A D^{-1} y = D M b for left preconditioning or
176: $ D A M D^{-1} z = D b for right preconditioning
178: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
180: .keywords: PC
182: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
183: @*/
184: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
185: {
191: pc->diagonalscale = PETSC_TRUE;
193: PetscObjectReference((PetscObject)s);
194: VecDestroy(&pc->diagonalscaleleft);
196: pc->diagonalscaleleft = s;
198: VecDuplicate(s,&pc->diagonalscaleright);
199: VecCopy(s,pc->diagonalscaleright);
200: VecReciprocal(pc->diagonalscaleright);
201: return(0);
202: }
206: /*@
207: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
209: Logically Collective on PC
211: Input Parameters:
212: + pc - the preconditioner context
213: . in - input vector
214: + out - scaled vector (maybe the same as in)
216: Level: intermediate
218: Notes: The system solved via the Krylov method is
219: $ D M A D^{-1} y = D M b for left preconditioning or
220: $ D A M D^{-1} z = D b for right preconditioning
222: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
224: If diagonal scaling is turned off and in is not out then in is copied to out
226: .keywords: PC
228: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
229: @*/
230: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
231: {
238: if (pc->diagonalscale) {
239: VecPointwiseMult(out,pc->diagonalscaleleft,in);
240: } else if (in != out) {
241: VecCopy(in,out);
242: }
243: return(0);
244: }
248: /*@
249: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
251: Logically Collective on PC
253: Input Parameters:
254: + pc - the preconditioner context
255: . in - input vector
256: + out - scaled vector (maybe the same as in)
258: Level: intermediate
260: Notes: The system solved via the Krylov method is
261: $ D M A D^{-1} y = D M b for left preconditioning or
262: $ D A M D^{-1} z = D b for right preconditioning
264: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
266: If diagonal scaling is turned off and in is not out then in is copied to out
268: .keywords: PC
270: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
271: @*/
272: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
273: {
280: if (pc->diagonalscale) {
281: VecPointwiseMult(out,pc->diagonalscaleright,in);
282: } else if (in != out) {
283: VecCopy(in,out);
284: }
285: return(0);
286: }
290: /*@
291: PCSetUseAmat - Sets a flag to indicate that when the preconditioner needs to apply (part of) the
292: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
293: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
295: Logically Collective on PC
297: Input Parameters:
298: + pc - the preconditioner context
299: - flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
301: Options Database Key:
302: . -pc_use_amat <true,false>
304: Notes:
305: For the common case in which the linear system matrix and the matrix used to construct the
306: preconditioner are identical, this routine is does nothing.
308: Level: intermediate
310: .seealso: PCGetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
311: @*/
312: PetscErrorCode PCSetUseAmat(PC pc,PetscBool flg)
313: {
316: pc->useAmat = flg;
317: return(0);
318: }
322: /*@
323: PCGetUseAmat - Gets a flag to indicate that when the preconditioner needs to apply (part of) the
324: operator during the preconditioning process it applies the Amat provided to TSSetRHSJacobian(),
325: TSSetIJacobian(), SNESSetJacobian(), KSPSetOperator() or PCSetOperator() not the Pmat.
327: Logically Collective on PC
329: Input Parameter:
330: . pc - the preconditioner context
332: Output Parameter:
333: . flg - PETSC_TRUE to use the Amat, PETSC_FALSE to use the Pmat (default is false)
335: Notes:
336: For the common case in which the linear system matrix and the matrix used to construct the
337: preconditioner are identical, this routine is does nothing.
339: Level: intermediate
341: .seealso: PCSetUseAmat(), PCBJACOBI, PGMG, PCFIELDSPLIT, PCCOMPOSITE
342: @*/
343: PetscErrorCode PCGetUseAmat(PC pc,PetscBool *flg)
344: {
347: *flg = pc->useAmat;
348: return(0);
349: }
353: /*@
354: PCCreate - Creates a preconditioner context.
356: Collective on MPI_Comm
358: Input Parameter:
359: . comm - MPI communicator
361: Output Parameter:
362: . pc - location to put the preconditioner context
364: Notes:
365: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or ICC
366: in parallel. For dense matrices it is always PCNONE.
368: Level: developer
370: .keywords: PC, create, context
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,_p_PC,struct _PCOps,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: /* -------------------------------------------------------------------------------*/
407: /*@
408: PCApply - Applies the preconditioner to a vector.
410: Collective on PC and Vec
412: Input Parameters:
413: + pc - the preconditioner context
414: - x - input vector
416: Output Parameter:
417: . y - output vector
419: Level: developer
421: .keywords: PC, apply
423: .seealso: PCApplyTranspose(), PCApplyBAorAB()
424: @*/
425: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
426: {
433: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
434: VecValidValues(x,2,PETSC_TRUE);
435: if (pc->setupcalled < 2) {
436: PCSetUp(pc);
437: }
438: if (!pc->ops->apply) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply");
439: PetscLogEventBegin(PC_Apply,pc,x,y,0);
440: (*pc->ops->apply)(pc,x,y);
441: PetscLogEventEnd(PC_Apply,pc,x,y,0);
442: VecValidValues(y,3,PETSC_FALSE);
443: return(0);
444: }
448: /*@
449: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
451: Collective on PC and Vec
453: Input Parameters:
454: + pc - the preconditioner context
455: - x - input vector
457: Output Parameter:
458: . y - output vector
460: Notes:
461: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
463: Level: developer
465: .keywords: PC, apply, symmetric, left
467: .seealso: PCApply(), PCApplySymmetricRight()
468: @*/
469: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
470: {
477: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
478: VecValidValues(x,2,PETSC_TRUE);
479: if (pc->setupcalled < 2) {
480: PCSetUp(pc);
481: }
482: if (!pc->ops->applysymmetricleft) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
483: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
484: (*pc->ops->applysymmetricleft)(pc,x,y);
485: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
486: VecValidValues(y,3,PETSC_FALSE);
487: return(0);
488: }
492: /*@
493: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
495: Collective on PC and Vec
497: Input Parameters:
498: + pc - the preconditioner context
499: - x - input vector
501: Output Parameter:
502: . y - output vector
504: Level: developer
506: Notes:
507: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
509: .keywords: PC, apply, symmetric, right
511: .seealso: PCApply(), PCApplySymmetricLeft()
512: @*/
513: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
514: {
521: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
522: VecValidValues(x,2,PETSC_TRUE);
523: if (pc->setupcalled < 2) {
524: PCSetUp(pc);
525: }
526: if (!pc->ops->applysymmetricright) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have left symmetric apply");
527: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
528: (*pc->ops->applysymmetricright)(pc,x,y);
529: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
530: VecValidValues(y,3,PETSC_FALSE);
531: return(0);
532: }
536: /*@
537: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
539: Collective on PC and Vec
541: Input Parameters:
542: + pc - the preconditioner context
543: - x - input vector
545: Output Parameter:
546: . y - output vector
548: Notes: For complex numbers this applies the non-Hermitian transpose.
550: Developer Notes: We need to implement a PCApplyHermitianTranspose()
552: Level: developer
554: .keywords: PC, apply, transpose
556: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
557: @*/
558: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
559: {
566: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
567: VecValidValues(x,2,PETSC_TRUE);
568: if (pc->setupcalled < 2) {
569: PCSetUp(pc);
570: }
571: if (!pc->ops->applytranspose) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply transpose");
572: PetscLogEventBegin(PC_Apply,pc,x,y,0);
573: (*pc->ops->applytranspose)(pc,x,y);
574: PetscLogEventEnd(PC_Apply,pc,x,y,0);
575: VecValidValues(y,3,PETSC_FALSE);
576: return(0);
577: }
581: /*@
582: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
584: Collective on PC and Vec
586: Input Parameters:
587: . pc - the preconditioner context
589: Output Parameter:
590: . flg - PETSC_TRUE if a transpose operation is defined
592: Level: developer
594: .keywords: PC, apply, transpose
596: .seealso: PCApplyTranspose()
597: @*/
598: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
599: {
603: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
604: else *flg = PETSC_FALSE;
605: return(0);
606: }
610: /*@
611: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
613: Collective on PC and Vec
615: Input Parameters:
616: + pc - the preconditioner context
617: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
618: . x - input vector
619: - work - work vector
621: Output Parameter:
622: . y - output vector
624: Level: developer
626: 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
627: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
629: .keywords: PC, apply, operator
631: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
632: @*/
633: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
634: {
642: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
643: VecValidValues(x,3,PETSC_TRUE);
644: 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");
645: if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
647: if (pc->setupcalled < 2) {
648: PCSetUp(pc);
649: }
651: if (pc->diagonalscale) {
652: if (pc->ops->applyBA) {
653: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
654: VecDuplicate(x,&work2);
655: PCDiagonalScaleRight(pc,x,work2);
656: (*pc->ops->applyBA)(pc,side,work2,y,work);
657: PCDiagonalScaleLeft(pc,y,y);
658: VecDestroy(&work2);
659: } else if (side == PC_RIGHT) {
660: PCDiagonalScaleRight(pc,x,y);
661: PCApply(pc,y,work);
662: MatMult(pc->mat,work,y);
663: PCDiagonalScaleLeft(pc,y,y);
664: } else if (side == PC_LEFT) {
665: PCDiagonalScaleRight(pc,x,y);
666: MatMult(pc->mat,y,work);
667: PCApply(pc,work,y);
668: PCDiagonalScaleLeft(pc,y,y);
669: } else if (side == PC_SYMMETRIC) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
670: } else {
671: if (pc->ops->applyBA) {
672: (*pc->ops->applyBA)(pc,side,x,y,work);
673: } else if (side == PC_RIGHT) {
674: PCApply(pc,x,work);
675: MatMult(pc->mat,work,y);
676: } else if (side == PC_LEFT) {
677: MatMult(pc->mat,x,work);
678: PCApply(pc,work,y);
679: } else if (side == PC_SYMMETRIC) {
680: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
681: PCApplySymmetricRight(pc,x,work);
682: MatMult(pc->mat,work,y);
683: VecCopy(y,work);
684: PCApplySymmetricLeft(pc,work,y);
685: }
686: }
687: VecValidValues(y,4,PETSC_FALSE);
688: return(0);
689: }
693: /*@
694: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
695: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
696: NOT tr(B*A) = tr(A)*tr(B).
698: Collective on PC and Vec
700: Input Parameters:
701: + pc - the preconditioner context
702: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
703: . x - input vector
704: - work - work vector
706: Output Parameter:
707: . y - output vector
710: 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
711: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
713: Level: developer
715: .keywords: PC, apply, operator, transpose
717: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
718: @*/
719: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
720: {
728: if (x == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"x and y must be different vectors");
729: VecValidValues(x,3,PETSC_TRUE);
730: if (pc->ops->applyBAtranspose) {
731: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
732: VecValidValues(y,4,PETSC_FALSE);
733: return(0);
734: }
735: if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
737: if (pc->setupcalled < 2) {
738: PCSetUp(pc);
739: }
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: VecValidValues(y,4,PETSC_FALSE);
750: return(0);
751: }
753: /* -------------------------------------------------------------------------------*/
757: /*@
758: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
759: built-in fast application of Richardson's method.
761: Not Collective
763: Input Parameter:
764: . pc - the preconditioner
766: Output Parameter:
767: . exists - PETSC_TRUE or PETSC_FALSE
769: Level: developer
771: .keywords: PC, apply, Richardson, exists
773: .seealso: PCApplyRichardson()
774: @*/
775: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
776: {
780: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
781: else *exists = PETSC_FALSE;
782: return(0);
783: }
787: /*@
788: PCApplyRichardson - Applies several steps of Richardson iteration with
789: the particular preconditioner. This routine is usually used by the
790: Krylov solvers and not the application code directly.
792: Collective on PC
794: Input Parameters:
795: + pc - the preconditioner context
796: . b - the right hand side
797: . w - one work vector
798: . rtol - relative decrease in residual norm convergence criteria
799: . abstol - absolute residual norm convergence criteria
800: . dtol - divergence residual norm increase criteria
801: . its - the number of iterations to apply.
802: - guesszero - if the input x contains nonzero initial guess
804: Output Parameter:
805: + outits - number of iterations actually used (for SOR this always equals its)
806: . reason - the reason the apply terminated
807: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
809: Notes:
810: Most preconditioners do not support this function. Use the command
811: PCApplyRichardsonExists() to determine if one does.
813: Except for the multigrid PC this routine ignores the convergence tolerances
814: and always runs for the number of iterations
816: Level: developer
818: .keywords: PC, apply, Richardson
820: .seealso: PCApplyRichardsonExists()
821: @*/
822: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
823: {
831: if (b == y) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_IDN,"b and y must be different vectors");
832: if (pc->setupcalled < 2) {
833: PCSetUp(pc);
834: }
835: if (!pc->ops->applyrichardson) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC does not have apply richardson");
836: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
837: return(0);
838: }
840: /*
841: a setupcall of 0 indicates never setup,
842: 1 indicates has been previously setup
843: */
846: /*@
847: PCSetUp - Prepares for the use of a preconditioner.
849: Collective on PC
851: Input Parameter:
852: . pc - the preconditioner context
854: Level: developer
856: .keywords: PC, setup
858: .seealso: PCCreate(), PCApply(), PCDestroy()
859: @*/
860: PetscErrorCode PCSetUp(PC pc)
861: {
862: PetscErrorCode ierr;
863: const char *def;
864: PetscObjectState matstate, matnonzerostate;
868: if (!pc->mat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
870: if (pc->setupcalled && pc->reusepreconditioner) {
871: PetscInfo(pc,"Leaving PC with identical preconditioner since reuse preconditioner is set\n");
872: return(0);
873: }
875: PetscObjectStateGet((PetscObject)pc->pmat,&matstate);
876: MatGetNonzeroState(pc->pmat,&matnonzerostate);
877: if (!pc->setupcalled) {
878: PetscInfo(pc,"Setting up PC for first time\n");
879: pc->flag = DIFFERENT_NONZERO_PATTERN;
880: } else if (matstate == pc->matstate) {
881: PetscInfo(pc,"Leaving PC with identical preconditioner since operator is unchanged\n");
882: return(0);
883: } else {
884: if (matnonzerostate > pc->matnonzerostate) {
885: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
886: pc->flag = DIFFERENT_NONZERO_PATTERN;
887: } else {
888: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
889: pc->flag = SAME_NONZERO_PATTERN;
890: }
891: }
892: pc->matstate = matstate;
893: pc->matnonzerostate = matnonzerostate;
895: if (!((PetscObject)pc)->type_name) {
896: PCGetDefaultType_Private(pc,&def);
897: PCSetType(pc,def);
898: }
900: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
901: if (pc->ops->setup) {
902: (*pc->ops->setup)(pc);
903: }
904: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
905: pc->setupcalled = 1;
906: return(0);
907: }
911: /*@
912: PCSetUpOnBlocks - Sets up the preconditioner for each block in
913: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
914: methods.
916: Collective on PC
918: Input Parameters:
919: . pc - the preconditioner context
921: Level: developer
923: .keywords: PC, setup, blocks
925: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
926: @*/
927: PetscErrorCode PCSetUpOnBlocks(PC pc)
928: {
933: if (!pc->ops->setuponblocks) return(0);
934: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
935: (*pc->ops->setuponblocks)(pc);
936: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
937: return(0);
938: }
942: /*@C
943: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
944: submatrices that arise within certain subdomain-based preconditioners.
945: The basic submatrices are extracted from the preconditioner matrix as
946: usual; the user can then alter these (for example, to set different boundary
947: conditions for each submatrix) before they are used for the local solves.
949: Logically Collective on PC
951: Input Parameters:
952: + pc - the preconditioner context
953: . func - routine for modifying the submatrices
954: - ctx - optional user-defined context (may be null)
956: Calling sequence of func:
957: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
959: . row - an array of index sets that contain the global row numbers
960: that comprise each local submatrix
961: . col - an array of index sets that contain the global column numbers
962: that comprise each local submatrix
963: . submat - array of local submatrices
964: - ctx - optional user-defined context for private data for the
965: user-defined func routine (may be null)
967: Notes:
968: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
969: KSPSolve().
971: A routine set by PCSetModifySubMatrices() is currently called within
972: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
973: preconditioners. All other preconditioners ignore this routine.
975: Level: advanced
977: .keywords: PC, set, modify, submatrices
979: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
980: @*/
981: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
982: {
985: pc->modifysubmatrices = func;
986: pc->modifysubmatricesP = ctx;
987: return(0);
988: }
992: /*@C
993: PCModifySubMatrices - Calls an optional user-defined routine within
994: certain preconditioners if one has been set with PCSetModifySubMarices().
996: Collective on PC
998: Input Parameters:
999: + pc - the preconditioner context
1000: . nsub - the number of local submatrices
1001: . row - an array of index sets that contain the global row numbers
1002: that comprise each local submatrix
1003: . col - an array of index sets that contain the global column numbers
1004: that comprise each local submatrix
1005: . submat - array of local submatrices
1006: - ctx - optional user-defined context for private data for the
1007: user-defined routine (may be null)
1009: Output Parameter:
1010: . submat - array of local submatrices (the entries of which may
1011: have been modified)
1013: Notes:
1014: The user should NOT generally call this routine, as it will
1015: automatically be called within certain preconditioners (currently
1016: block Jacobi, additive Schwarz) if set.
1018: The basic submatrices are extracted from the preconditioner matrix
1019: as usual; the user can then alter these (for example, to set different
1020: boundary conditions for each submatrix) before they are used for the
1021: local solves.
1023: Level: developer
1025: .keywords: PC, modify, submatrices
1027: .seealso: PCSetModifySubMatrices()
1028: @*/
1029: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
1030: {
1035: if (!pc->modifysubmatrices) return(0);
1036: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
1037: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
1038: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
1039: return(0);
1040: }
1044: /*@
1045: PCSetOperators - Sets the matrix associated with the linear system and
1046: a (possibly) different one associated with the preconditioner.
1048: Logically Collective on PC and Mat
1050: Input Parameters:
1051: + pc - the preconditioner context
1052: . Amat - the matrix that defines the linear system
1053: - Pmat - the matrix to be used in constructing the preconditioner, usually the same as Amat.
1055: Notes:
1056: Passing a NULL for Amat or Pmat removes the matrix that is currently used.
1058: If you wish to replace either Amat or Pmat but leave the other one untouched then
1059: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1060: on it and then pass it back in in your call to KSPSetOperators().
1062: More Notes about Repeated Solution of Linear Systems:
1063: PETSc does NOT reset the matrix entries of either Amat or Pmat
1064: to zero after a linear solve; the user is completely responsible for
1065: matrix assembly. See the routine MatZeroEntries() if desiring to
1066: zero all elements of a matrix.
1068: Level: intermediate
1070: .keywords: PC, set, operators, matrix, linear system
1072: .seealso: PCGetOperators(), MatZeroEntries()
1073: @*/
1074: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat)
1075: {
1076: PetscErrorCode ierr;
1077: PetscInt m1,n1,m2,n2;
1085: if (pc->setupcalled && pc->mat && pc->pmat && Amat && Pmat) {
1086: MatGetLocalSize(Amat,&m1,&n1);
1087: MatGetLocalSize(pc->mat,&m2,&n2);
1088: 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);
1089: MatGetLocalSize(Pmat,&m1,&n1);
1090: MatGetLocalSize(pc->pmat,&m2,&n2);
1091: 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);
1092: }
1094: if (Pmat != pc->pmat) {
1095: /* changing the operator that defines the preconditioner thus reneed to clear current states so new preconditioner is built */
1096: pc->matnonzerostate = -1;
1097: pc->matstate = -1;
1098: }
1100: /* reference first in case the matrices are the same */
1101: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1102: MatDestroy(&pc->mat);
1103: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1104: MatDestroy(&pc->pmat);
1105: pc->mat = Amat;
1106: pc->pmat = Pmat;
1107: return(0);
1108: }
1112: /*@
1113: PCSetReusePreconditioner - reuse the current preconditioner even if the operator in the preconditioner has changed.
1115: Logically Collective on PC
1117: Input Parameters:
1118: + pc - the preconditioner context
1119: - flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1121: .seealso: PCGetOperators(), MatZeroEntries(), PCGetReusePreconditioner(), KSPSetReusePreconditioner()
1122: @*/
1123: PetscErrorCode PCSetReusePreconditioner(PC pc,PetscBool flag)
1124: {
1127: pc->reusepreconditioner = flag;
1128: return(0);
1129: }
1133: /*@
1134: PCGetReusePreconditioner - Determines if the PC reuses the current preconditioner even if the operator in the preconditioner has changed.
1136: Not Collective
1138: Input Parameter:
1139: . pc - the preconditioner context
1141: Output Parameter:
1142: . flag - PETSC_TRUE do not compute a new preconditioner, PETSC_FALSE do compute a new preconditioner
1144: .seealso: PCGetOperators(), MatZeroEntries(), PCSetReusePreconditioner()
1145: @*/
1146: PetscErrorCode PCGetReusePreconditioner(PC pc,PetscBool *flag)
1147: {
1150: *flag = pc->reusepreconditioner;
1151: return(0);
1152: }
1156: /*@C
1157: PCGetOperators - Gets the matrix associated with the linear system and
1158: possibly a different one associated with the preconditioner.
1160: Not collective, though parallel Mats are returned if the PC is parallel
1162: Input Parameter:
1163: . pc - the preconditioner context
1165: Output Parameters:
1166: + Amat - the matrix defining the linear system
1167: - Pmat - the matrix from which the preconditioner is constructed, usually the same as Amat.
1169: Level: intermediate
1171: Notes: Does not increase the reference count of the matrices, so you should not destroy them
1173: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1174: are created in PC and returned to the user. In this case, if both operators
1175: mat and pmat are requested, two DIFFERENT operators will be returned. If
1176: only one is requested both operators in the PC will be the same (i.e. as
1177: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1178: The user must set the sizes of the returned matrices and their type etc just
1179: as if the user created them with MatCreate(). For example,
1181: $ KSP/PCGetOperators(ksp/pc,&Amat,NULL); is equivalent to
1182: $ set size, type, etc of Amat
1184: $ MatCreate(comm,&mat);
1185: $ KSP/PCSetOperators(ksp/pc,Amat,Amat);
1186: $ PetscObjectDereference((PetscObject)mat);
1187: $ set size, type, etc of Amat
1189: and
1191: $ KSP/PCGetOperators(ksp/pc,&Amat,&Pmat); is equivalent to
1192: $ set size, type, etc of Amat and Pmat
1194: $ MatCreate(comm,&Amat);
1195: $ MatCreate(comm,&Pmat);
1196: $ KSP/PCSetOperators(ksp/pc,Amat,Pmat);
1197: $ PetscObjectDereference((PetscObject)Amat);
1198: $ PetscObjectDereference((PetscObject)Pmat);
1199: $ set size, type, etc of Amat and Pmat
1201: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1202: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1203: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1204: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1205: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1206: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1207: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1208: it can be created for you?
1211: .keywords: PC, get, operators, matrix, linear system
1213: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1214: @*/
1215: PetscErrorCode PCGetOperators(PC pc,Mat *Amat,Mat *Pmat)
1216: {
1221: if (Amat) {
1222: if (!pc->mat) {
1223: if (pc->pmat && !Pmat) { /* Apmat has been set, but user did not request it, so use for Amat */
1224: pc->mat = pc->pmat;
1225: PetscObjectReference((PetscObject)pc->mat);
1226: } else { /* both Amat and Pmat are empty */
1227: MatCreate(PetscObjectComm((PetscObject)pc),&pc->mat);
1228: if (!Pmat) { /* user did NOT request Pmat, so make same as Amat */
1229: pc->pmat = pc->mat;
1230: PetscObjectReference((PetscObject)pc->pmat);
1231: }
1232: }
1233: }
1234: *Amat = pc->mat;
1235: }
1236: if (Pmat) {
1237: if (!pc->pmat) {
1238: if (pc->mat && !Amat) { /* Amat has been set but was not requested, so use for pmat */
1239: pc->pmat = pc->mat;
1240: PetscObjectReference((PetscObject)pc->pmat);
1241: } else {
1242: MatCreate(PetscObjectComm((PetscObject)pc),&pc->pmat);
1243: if (!Amat) { /* user did NOT request Amat, so make same as Pmat */
1244: pc->mat = pc->pmat;
1245: PetscObjectReference((PetscObject)pc->mat);
1246: }
1247: }
1248: }
1249: *Pmat = pc->pmat;
1250: }
1251: return(0);
1252: }
1256: /*@C
1257: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1258: possibly a different one associated with the preconditioner have been set in the PC.
1260: Not collective, though the results on all processes should be the same
1262: Input Parameter:
1263: . pc - the preconditioner context
1265: Output Parameters:
1266: + mat - the matrix associated with the linear system was set
1267: - pmat - matrix associated with the preconditioner was set, usually the same
1269: Level: intermediate
1271: .keywords: PC, get, operators, matrix, linear system
1273: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1274: @*/
1275: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1276: {
1279: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1280: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1281: return(0);
1282: }
1286: /*@
1287: PCFactorGetMatrix - Gets the factored matrix from the
1288: preconditioner context. This routine is valid only for the LU,
1289: incomplete LU, Cholesky, and incomplete Cholesky methods.
1291: Not Collective on PC though Mat is parallel if PC is parallel
1293: Input Parameters:
1294: . pc - the preconditioner context
1296: Output parameters:
1297: . mat - the factored matrix
1299: Level: advanced
1301: Notes: Does not increase the reference count for the matrix so DO NOT destroy it
1303: .keywords: PC, get, factored, matrix
1304: @*/
1305: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1306: {
1312: if (pc->ops->getfactoredmatrix) {
1313: (*pc->ops->getfactoredmatrix)(pc,mat);
1314: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1315: return(0);
1316: }
1320: /*@C
1321: PCSetOptionsPrefix - Sets the prefix used for searching for all
1322: PC options in the database.
1324: Logically Collective on PC
1326: Input Parameters:
1327: + pc - the preconditioner context
1328: - prefix - the prefix string to prepend to all PC option requests
1330: Notes:
1331: A hyphen (-) must NOT be given at the beginning of the prefix name.
1332: The first character of all runtime options is AUTOMATICALLY the
1333: hyphen.
1335: Level: advanced
1337: .keywords: PC, set, options, prefix, database
1339: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1340: @*/
1341: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1342: {
1347: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1348: return(0);
1349: }
1353: /*@C
1354: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1355: PC options in the database.
1357: Logically Collective on PC
1359: Input Parameters:
1360: + pc - the preconditioner context
1361: - prefix - the prefix string to prepend to all PC option requests
1363: Notes:
1364: A hyphen (-) must NOT be given at the beginning of the prefix name.
1365: The first character of all runtime options is AUTOMATICALLY the
1366: hyphen.
1368: Level: advanced
1370: .keywords: PC, append, options, prefix, database
1372: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1373: @*/
1374: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1375: {
1380: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1381: return(0);
1382: }
1386: /*@C
1387: PCGetOptionsPrefix - Gets the prefix used for searching for all
1388: PC options in the database.
1390: Not Collective
1392: Input Parameters:
1393: . pc - the preconditioner context
1395: Output Parameters:
1396: . prefix - pointer to the prefix string used, is returned
1398: Notes: On the fortran side, the user should pass in a string 'prifix' of
1399: sufficient length to hold the prefix.
1401: Level: advanced
1403: .keywords: PC, get, options, prefix, database
1405: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1406: @*/
1407: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1408: {
1414: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1415: return(0);
1416: }
1420: /*@
1421: PCPreSolve - Optional pre-solve phase, intended for any
1422: preconditioner-specific actions that must be performed before
1423: the iterative solve itself.
1425: Collective on PC
1427: Input Parameters:
1428: + pc - the preconditioner context
1429: - ksp - the Krylov subspace context
1431: Level: developer
1433: Sample of Usage:
1434: .vb
1435: PCPreSolve(pc,ksp);
1436: KSPSolve(ksp,b,x);
1437: PCPostSolve(pc,ksp);
1438: .ve
1440: Notes:
1441: The pre-solve phase is distinct from the PCSetUp() phase.
1443: KSPSolve() calls this directly, so is rarely called by the user.
1445: .keywords: PC, pre-solve
1447: .seealso: PCPostSolve()
1448: @*/
1449: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1450: {
1452: Vec x,rhs;
1457: pc->presolvedone++;
1458: if (pc->presolvedone > 2) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1459: KSPGetSolution(ksp,&x);
1460: KSPGetRhs(ksp,&rhs);
1462: if (pc->ops->presolve) {
1463: (*pc->ops->presolve)(pc,ksp,rhs,x);
1464: }
1465: return(0);
1466: }
1470: /*@
1471: PCPostSolve - Optional post-solve phase, intended for any
1472: preconditioner-specific actions that must be performed after
1473: the iterative solve itself.
1475: Collective on PC
1477: Input Parameters:
1478: + pc - the preconditioner context
1479: - ksp - the Krylov subspace context
1481: Sample of Usage:
1482: .vb
1483: PCPreSolve(pc,ksp);
1484: KSPSolve(ksp,b,x);
1485: PCPostSolve(pc,ksp);
1486: .ve
1488: Note:
1489: KSPSolve() calls this routine directly, so it is rarely called by the user.
1491: Level: developer
1493: .keywords: PC, post-solve
1495: .seealso: PCPreSolve(), KSPSolve()
1496: @*/
1497: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1498: {
1500: Vec x,rhs;
1505: pc->presolvedone--;
1506: KSPGetSolution(ksp,&x);
1507: KSPGetRhs(ksp,&rhs);
1508: if (pc->ops->postsolve) {
1509: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1510: }
1511: return(0);
1512: }
1516: /*@C
1517: PCLoad - Loads a PC that has been stored in binary with PCView().
1519: Collective on PetscViewer
1521: Input Parameters:
1522: + newdm - the newly loaded PC, this needs to have been created with PCCreate() or
1523: some related function before a call to PCLoad().
1524: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
1526: Level: intermediate
1528: Notes:
1529: The type is determined by the data in the file, any type set into the PC before this call is ignored.
1531: Notes for advanced users:
1532: Most users should not need to know the details of the binary storage
1533: format, since PCLoad() and PCView() completely hide these details.
1534: But for anyone who's interested, the standard binary matrix storage
1535: format is
1536: .vb
1537: has not yet been determined
1538: .ve
1540: .seealso: PetscViewerBinaryOpen(), PCView(), MatLoad(), VecLoad()
1541: @*/
1542: PetscErrorCode PCLoad(PC newdm, PetscViewer viewer)
1543: {
1545: PetscBool isbinary;
1546: PetscInt classid;
1547: char type[256];
1552: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1553: if (!isbinary) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
1555: PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
1556: if (classid != PC_FILE_CLASSID) SETERRQ(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not PC next in file");
1557: PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
1558: PCSetType(newdm, type);
1559: if (newdm->ops->load) {
1560: (*newdm->ops->load)(newdm,viewer);
1561: }
1562: return(0);
1563: }
1565: #include <petscdraw.h>
1566: #if defined(PETSC_HAVE_SAWS)
1567: #include <petscviewersaws.h>
1568: #endif
1571: /*@C
1572: PCView - Prints the PC data structure.
1574: Collective on PC
1576: Input Parameters:
1577: + PC - the PC context
1578: - viewer - optional visualization context
1580: Note:
1581: The available visualization contexts include
1582: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1583: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1584: output where only the first processor opens
1585: the file. All other processors send their
1586: data to the first processor to print.
1588: The user can open an alternative visualization contexts with
1589: PetscViewerASCIIOpen() (output to a specified file).
1591: Level: developer
1593: .keywords: PC, view
1595: .seealso: KSPView(), PetscViewerASCIIOpen()
1596: @*/
1597: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1598: {
1599: PCType cstr;
1600: PetscErrorCode ierr;
1601: PetscBool iascii,isstring,isbinary,isdraw;
1602: PetscViewerFormat format;
1603: #if defined(PETSC_HAVE_SAWS)
1604: PetscBool isams;
1605: #endif
1609: if (!viewer) {
1610: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)pc),&viewer);
1611: }
1615: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1616: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1617: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
1618: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
1619: #if defined(PETSC_HAVE_SAWS)
1620: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSAWS,&isams);
1621: #endif
1623: if (iascii) {
1624: PetscViewerGetFormat(viewer,&format);
1625: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer);
1626: if (!pc->setupcalled) {
1627: PetscViewerASCIIPrintf(viewer," PC has not been set up so information may be incomplete\n");
1628: }
1629: if (pc->ops->view) {
1630: PetscViewerASCIIPushTab(viewer);
1631: (*pc->ops->view)(pc,viewer);
1632: PetscViewerASCIIPopTab(viewer);
1633: }
1634: if (pc->mat) {
1635: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1636: if (pc->pmat == pc->mat) {
1637: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1638: PetscViewerASCIIPushTab(viewer);
1639: MatView(pc->mat,viewer);
1640: PetscViewerASCIIPopTab(viewer);
1641: } else {
1642: if (pc->pmat) {
1643: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1644: } else {
1645: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1646: }
1647: PetscViewerASCIIPushTab(viewer);
1648: MatView(pc->mat,viewer);
1649: if (pc->pmat) {MatView(pc->pmat,viewer);}
1650: PetscViewerASCIIPopTab(viewer);
1651: }
1652: PetscViewerPopFormat(viewer);
1653: }
1654: } else if (isstring) {
1655: PCGetType(pc,&cstr);
1656: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1657: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1658: } else if (isbinary) {
1659: PetscInt classid = PC_FILE_CLASSID;
1660: MPI_Comm comm;
1661: PetscMPIInt rank;
1662: char type[256];
1664: PetscObjectGetComm((PetscObject)pc,&comm);
1665: MPI_Comm_rank(comm,&rank);
1666: if (!rank) {
1667: PetscViewerBinaryWrite(viewer,&classid,1,PETSC_INT,PETSC_FALSE);
1668: PetscStrncpy(type,((PetscObject)pc)->type_name,256);
1669: PetscViewerBinaryWrite(viewer,type,256,PETSC_CHAR,PETSC_FALSE);
1670: }
1671: if (pc->ops->view) {
1672: (*pc->ops->view)(pc,viewer);
1673: }
1674: } else if (isdraw) {
1675: PetscDraw draw;
1676: char str[25];
1677: PetscReal x,y,bottom,h;
1678: PetscInt n;
1680: PetscViewerDrawGetDraw(viewer,0,&draw);
1681: PetscDrawGetCurrentPoint(draw,&x,&y);
1682: if (pc->mat) {
1683: MatGetSize(pc->mat,&n,NULL);
1684: PetscSNPrintf(str,25,"PC: %s (%D)",((PetscObject)pc)->type_name,n);
1685: } else {
1686: PetscSNPrintf(str,25,"PC: %s",((PetscObject)pc)->type_name);
1687: }
1688: PetscDrawBoxedString(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
1689: bottom = y - h;
1690: PetscDrawPushCurrentPoint(draw,x,bottom);
1691: if (pc->ops->view) {
1692: (*pc->ops->view)(pc,viewer);
1693: }
1694: PetscDrawPopCurrentPoint(draw);
1695: #if defined(PETSC_HAVE_SAWS)
1696: } else if (isams) {
1697: PetscMPIInt rank;
1699: PetscObjectName((PetscObject)pc);
1700: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
1701: if (!((PetscObject)pc)->amsmem && !rank) {
1702: PetscObjectViewSAWs((PetscObject)pc,viewer);
1703: }
1704: if (pc->mat) {MatView(pc->mat,viewer);}
1705: if (pc->pmat && pc->pmat != pc->mat) {MatView(pc->pmat,viewer);}
1706: #endif
1707: }
1708: return(0);
1709: }
1714: /*@
1715: PCSetInitialGuessNonzero - Tells the iterative solver that the
1716: initial guess is nonzero; otherwise PC assumes the initial guess
1717: is to be zero (and thus zeros it out before solving).
1719: Logically Collective on PC
1721: Input Parameters:
1722: + pc - iterative context obtained from PCCreate()
1723: - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
1725: Level: Developer
1727: Notes:
1728: This is a weird function. Since PC's are linear operators on the right hand side they
1729: CANNOT use an initial guess. This function is for the "pass-through" preconditioners
1730: PCKSP and PCREDUNDANT and causes the inner KSP object to use the nonzero
1731: initial guess. Not currently working for PCREDUNDANT, that has to be rewritten to use KSP.
1734: .keywords: PC, set, initial guess, nonzero
1736: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
1737: @*/
1738: PetscErrorCode PCSetInitialGuessNonzero(PC pc,PetscBool flg)
1739: {
1742: pc->nonzero_guess = flg;
1743: return(0);
1744: }
1748: /*@C
1749: PCRegister - Adds a method to the preconditioner package.
1751: Not collective
1753: Input Parameters:
1754: + name_solver - name of a new user-defined solver
1755: - routine_create - routine to create method context
1757: Notes:
1758: PCRegister() may be called multiple times to add several user-defined preconditioners.
1760: Sample usage:
1761: .vb
1762: PCRegister("my_solver", MySolverCreate);
1763: .ve
1765: Then, your solver can be chosen with the procedural interface via
1766: $ PCSetType(pc,"my_solver")
1767: or at runtime via the option
1768: $ -pc_type my_solver
1770: Level: advanced
1772: .keywords: PC, register
1774: .seealso: PCRegisterAll(), PCRegisterDestroy()
1775: @*/
1776: PetscErrorCode PCRegister(const char sname[],PetscErrorCode (*function)(PC))
1777: {
1781: PetscFunctionListAdd(&PCList,sname,function);
1782: return(0);
1783: }
1787: /*@
1788: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1790: Collective on PC
1792: Input Parameter:
1793: . pc - the preconditioner object
1795: Output Parameter:
1796: . mat - the explict preconditioned operator
1798: Notes:
1799: This computation is done by applying the operators to columns of the
1800: identity matrix.
1802: Currently, this routine uses a dense matrix format when 1 processor
1803: is used and a sparse format otherwise. This routine is costly in general,
1804: and is recommended for use only with relatively small systems.
1806: Level: advanced
1808: .keywords: PC, compute, explicit, operator
1810: .seealso: KSPComputeExplicitOperator()
1812: @*/
1813: PetscErrorCode PCComputeExplicitOperator(PC pc,Mat *mat)
1814: {
1815: Vec in,out;
1817: PetscInt i,M,m,*rows,start,end;
1818: PetscMPIInt size;
1819: MPI_Comm comm;
1820: PetscScalar *array,one = 1.0;
1826: PetscObjectGetComm((PetscObject)pc,&comm);
1827: MPI_Comm_size(comm,&size);
1829: if (!pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1830: MatGetVecs(pc->pmat,&in,0);
1831: VecDuplicate(in,&out);
1832: VecGetOwnershipRange(in,&start,&end);
1833: VecGetSize(in,&M);
1834: VecGetLocalSize(in,&m);
1835: PetscMalloc1((m+1),&rows);
1836: for (i=0; i<m; i++) rows[i] = start + i;
1838: MatCreate(comm,mat);
1839: MatSetSizes(*mat,m,m,M,M);
1840: if (size == 1) {
1841: MatSetType(*mat,MATSEQDENSE);
1842: MatSeqDenseSetPreallocation(*mat,NULL);
1843: } else {
1844: MatSetType(*mat,MATMPIAIJ);
1845: MatMPIAIJSetPreallocation(*mat,0,NULL,0,NULL);
1846: }
1847: MatSetOption(*mat,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);
1849: for (i=0; i<M; i++) {
1851: VecSet(in,0.0);
1852: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1853: VecAssemblyBegin(in);
1854: VecAssemblyEnd(in);
1856: /* should fix, allowing user to choose side */
1857: PCApply(pc,in,out);
1859: VecGetArray(out,&array);
1860: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1861: VecRestoreArray(out,&array);
1863: }
1864: PetscFree(rows);
1865: VecDestroy(&out);
1866: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1867: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1868: return(0);
1869: }
1873: /*@
1874: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1876: Collective on PC
1878: Input Parameters:
1879: + pc - the solver context
1880: . dim - the dimension of the coordinates 1, 2, or 3
1881: - coords - the coordinates
1883: Level: intermediate
1885: Notes: coords is an array of the 3D coordinates for the nodes on
1886: the local processor. So if there are 108 equation on a processor
1887: for a displacement finite element discretization of elasticity (so
1888: that there are 36 = 108/3 nodes) then the array must have 108
1889: double precision values (ie, 3 * 36). These x y z coordinates
1890: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1891: ... , N-1.z ].
1893: .seealso: MatSetNearNullSpace
1894: @*/
1895: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1896: {
1900: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1901: return(0);
1902: }