Actual source code: precon.c
petsc-3.3-p7 2013-05-11
2: /*
3: The PC (preconditioner) interface routines, callable by users.
4: */
5: #include <petsc-private/pcimpl.h> /*I "petscksp.h" I*/
7: /* Logging support */
8: PetscClassId PC_CLASSID;
9: PetscLogEvent PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
10: PetscLogEvent PC_ApplySymmetricRight, PC_ModifySubMatrices, PC_ApplyOnBlocks, PC_ApplyTransposeOnBlocks, PC_ApplyOnMproc;
14: PetscErrorCode PCGetDefaultType_Private(PC pc,const char* type[])
15: {
17: PetscMPIInt size;
18: PetscBool flg1,flg2,set,flg3;
21: MPI_Comm_size(((PetscObject)pc)->comm,&size);
22: if (pc->pmat) {
23: PetscErrorCode (*f)(Mat,MatReuse,Mat*);
24: PetscObjectQueryFunction((PetscObject)pc->pmat,"MatGetDiagonalBlock_C",(void (**)(void))&f);
25: if (size == 1) {
26: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ICC,&flg1);
27: MatGetFactorAvailable(pc->pmat,"petsc",MAT_FACTOR_ILU,&flg2);
28: MatIsSymmetricKnown(pc->pmat,&set,&flg3);
29: if (flg1 && (!flg2 || (set && flg3))) {
30: *type = PCICC;
31: } else if (flg2) {
32: *type = PCILU;
33: } else if (f) { /* likely is a parallel matrix run on one processor */
34: *type = PCBJACOBI;
35: } else {
36: *type = PCNONE;
37: }
38: } else {
39: if (f) {
40: *type = PCBJACOBI;
41: } else {
42: *type = PCNONE;
43: }
44: }
45: } else {
46: if (size == 1) {
47: *type = PCILU;
48: } else {
49: *type = PCBJACOBI;
50: }
51: }
52: return(0);
53: }
57: /*@
58: PCReset - Resets a PC context to the pcsetupcalled = 0 state and removes any allocated Vecs and Mats
60: Collective on PC
62: Input Parameter:
63: . pc - the preconditioner context
65: Level: developer
67: 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
69: .keywords: PC, destroy
71: .seealso: PCCreate(), PCSetUp()
72: @*/
73: PetscErrorCode PCReset(PC pc)
74: {
79: if (pc->ops->reset) {
80: (*pc->ops->reset)(pc);
81: }
82: VecDestroy(&pc->diagonalscaleright);
83: VecDestroy(&pc->diagonalscaleleft);
84: MatDestroy(&pc->pmat);
85: MatDestroy(&pc->mat);
86: pc->setupcalled = 0;
87: return(0);
88: }
92: /*@
93: PCDestroy - Destroys PC context that was created with PCCreate().
95: Collective on PC
97: Input Parameter:
98: . pc - the preconditioner context
100: Level: developer
102: .keywords: PC, destroy
104: .seealso: PCCreate(), PCSetUp()
105: @*/
106: PetscErrorCode PCDestroy(PC *pc)
107: {
111: if (!*pc) return(0);
113: if (--((PetscObject)(*pc))->refct > 0) {*pc = 0; return(0);}
115: PCReset(*pc);
117: /* if memory was published with AMS then destroy it */
118: PetscObjectDepublish((*pc));
119: if ((*pc)->ops->destroy) {(*(*pc)->ops->destroy)((*pc));}
120: DMDestroy(&(*pc)->dm);
121: PetscHeaderDestroy(pc);
122: return(0);
123: }
127: /*@C
128: PCGetDiagonalScale - Indicates if the preconditioner applies an additional left and right
129: scaling as needed by certain time-stepping codes.
131: Logically Collective on PC
133: Input Parameter:
134: . pc - the preconditioner context
136: Output Parameter:
137: . flag - PETSC_TRUE if it applies the scaling
139: Level: developer
141: Notes: If this returns PETSC_TRUE then the system solved via the Krylov method is
142: $ D M A D^{-1} y = D M b for left preconditioning or
143: $ D A M D^{-1} z = D b for right preconditioning
145: .keywords: PC
147: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCSetDiagonalScale()
148: @*/
149: PetscErrorCode PCGetDiagonalScale(PC pc,PetscBool *flag)
150: {
154: *flag = pc->diagonalscale;
155: return(0);
156: }
160: /*@
161: PCSetDiagonalScale - Indicates the left scaling to use to apply an additional left and right
162: scaling as needed by certain time-stepping codes.
164: Logically Collective on PC
166: Input Parameters:
167: + pc - the preconditioner context
168: - s - scaling vector
170: Level: intermediate
172: Notes: The system solved via the Krylov method is
173: $ D M A D^{-1} y = D M b for left preconditioning or
174: $ D A M D^{-1} z = D b for right preconditioning
176: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
178: .keywords: PC
180: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCGetDiagonalScale()
181: @*/
182: PetscErrorCode PCSetDiagonalScale(PC pc,Vec s)
183: {
189: pc->diagonalscale = PETSC_TRUE;
190: PetscObjectReference((PetscObject)s);
191: VecDestroy(&pc->diagonalscaleleft);
192: pc->diagonalscaleleft = s;
193: VecDuplicate(s,&pc->diagonalscaleright);
194: VecCopy(s,pc->diagonalscaleright);
195: VecReciprocal(pc->diagonalscaleright);
196: return(0);
197: }
201: /*@
202: PCDiagonalScaleLeft - Scales a vector by the left scaling as needed by certain time-stepping codes.
204: Logically Collective on PC
206: Input Parameters:
207: + pc - the preconditioner context
208: . in - input vector
209: + out - scaled vector (maybe the same as in)
211: Level: intermediate
213: Notes: The system solved via the Krylov method is
214: $ D M A D^{-1} y = D M b for left preconditioning or
215: $ D A M D^{-1} z = D b for right preconditioning
217: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
219: If diagonal scaling is turned off and in is not out then in is copied to out
221: .keywords: PC
223: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
224: @*/
225: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
226: {
233: if (pc->diagonalscale) {
234: VecPointwiseMult(out,pc->diagonalscaleleft,in);
235: } else if (in != out) {
236: VecCopy(in,out);
237: }
238: return(0);
239: }
243: /*@
244: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
246: Logically Collective on PC
248: Input Parameters:
249: + pc - the preconditioner context
250: . in - input vector
251: + out - scaled vector (maybe the same as in)
253: Level: intermediate
255: Notes: The system solved via the Krylov method is
256: $ D M A D^{-1} y = D M b for left preconditioning or
257: $ D A M D^{-1} z = D b for right preconditioning
259: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
261: If diagonal scaling is turned off and in is not out then in is copied to out
263: .keywords: PC
265: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
266: @*/
267: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
268: {
275: if (pc->diagonalscale) {
276: VecPointwiseMult(out,pc->diagonalscaleright,in);
277: } else if (in != out) {
278: VecCopy(in,out);
279: }
280: return(0);
281: }
283: #if 0
286: static PetscErrorCode PCPublish_Petsc(PetscObject obj)
287: {
289: return(0);
290: }
291: #endif
295: /*@
296: PCCreate - Creates a preconditioner context.
298: Collective on MPI_Comm
300: Input Parameter:
301: . comm - MPI communicator
303: Output Parameter:
304: . pc - location to put the preconditioner context
306: Notes:
307: The default preconditioner for sparse matrices is PCILU or PCICC with 0 fill on one process and block Jacobi with PCILU or ICC
308: in parallel. For dense matrices it is always PCNONE.
310: Level: developer
312: .keywords: PC, create, context
314: .seealso: PCSetUp(), PCApply(), PCDestroy()
315: @*/
316: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
317: {
318: PC pc;
323: *newpc = 0;
324: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
325: PCInitializePackage(PETSC_NULL);
326: #endif
328: PetscHeaderCreate(pc,_p_PC,struct _PCOps,PC_CLASSID,-1,"PC","Preconditioner","PC",comm,PCDestroy,PCView);
330: pc->mat = 0;
331: pc->pmat = 0;
332: pc->setupcalled = 0;
333: pc->setfromoptionscalled = 0;
334: pc->data = 0;
335: pc->diagonalscale = PETSC_FALSE;
336: pc->diagonalscaleleft = 0;
337: pc->diagonalscaleright = 0;
338: pc->reuse = 0;
340: pc->modifysubmatrices = 0;
341: pc->modifysubmatricesP = 0;
342: *newpc = pc;
343: return(0);
345: }
347: /* -------------------------------------------------------------------------------*/
351: /*@
352: PCApply - Applies the preconditioner to a vector.
354: Collective on PC and Vec
356: Input Parameters:
357: + pc - the preconditioner context
358: - x - input vector
360: Output Parameter:
361: . y - output vector
363: Level: developer
365: .keywords: PC, apply
367: .seealso: PCApplyTranspose(), PCApplyBAorAB()
368: @*/
369: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
370: {
377: if (x == y) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
378: VecValidValues(x,2,PETSC_TRUE);
379: if (pc->setupcalled < 2) {
380: PCSetUp(pc);
381: }
382: if (!pc->ops->apply) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"PC does not have apply");
383: PetscLogEventBegin(PC_Apply,pc,x,y,0);
384: (*pc->ops->apply)(pc,x,y);
385: PetscLogEventEnd(PC_Apply,pc,x,y,0);
386: VecValidValues(y,3,PETSC_FALSE);
387: return(0);
388: }
392: /*@
393: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
395: Collective on PC and Vec
397: Input Parameters:
398: + pc - the preconditioner context
399: - x - input vector
401: Output Parameter:
402: . y - output vector
404: Notes:
405: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
407: Level: developer
409: .keywords: PC, apply, symmetric, left
411: .seealso: PCApply(), PCApplySymmetricRight()
412: @*/
413: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
414: {
421: if (x == y) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
422: VecValidValues(x,2,PETSC_TRUE);
423: if (pc->setupcalled < 2) {
424: PCSetUp(pc);
425: }
426: if (!pc->ops->applysymmetricleft) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"PC does not have left symmetric apply");
427: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
428: (*pc->ops->applysymmetricleft)(pc,x,y);
429: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
430: VecValidValues(y,3,PETSC_FALSE);
431: return(0);
432: }
436: /*@
437: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
439: Collective on PC and Vec
441: Input Parameters:
442: + pc - the preconditioner context
443: - x - input vector
445: Output Parameter:
446: . y - output vector
448: Level: developer
450: Notes:
451: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
453: .keywords: PC, apply, symmetric, right
455: .seealso: PCApply(), PCApplySymmetricLeft()
456: @*/
457: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
458: {
465: if (x == y) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
466: VecValidValues(x,2,PETSC_TRUE);
467: if (pc->setupcalled < 2) {
468: PCSetUp(pc);
469: }
470: if (!pc->ops->applysymmetricright) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"PC does not have left symmetric apply");
471: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
472: (*pc->ops->applysymmetricright)(pc,x,y);
473: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
474: VecValidValues(y,3,PETSC_FALSE);
475: return(0);
476: }
480: /*@
481: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
483: Collective on PC and Vec
485: Input Parameters:
486: + pc - the preconditioner context
487: - x - input vector
489: Output Parameter:
490: . y - output vector
492: Notes: For complex numbers this applies the non-Hermitian transpose.
494: Developer Notes: We need to implement a PCApplyHermitianTranspose()
496: Level: developer
498: .keywords: PC, apply, transpose
500: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCApplyTransposeExists()
501: @*/
502: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
503: {
510: if (x == y) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
511: VecValidValues(x,2,PETSC_TRUE);
512: if (pc->setupcalled < 2) {
513: PCSetUp(pc);
514: }
515: if (!pc->ops->applytranspose) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"PC does not have apply transpose");
516: PetscLogEventBegin(PC_Apply,pc,x,y,0);
517: (*pc->ops->applytranspose)(pc,x,y);
518: PetscLogEventEnd(PC_Apply,pc,x,y,0);
519: VecValidValues(y,3,PETSC_FALSE);
520: return(0);
521: }
525: /*@
526: PCApplyTransposeExists - Test whether the preconditioner has a transpose apply operation
528: Collective on PC and Vec
530: Input Parameters:
531: . pc - the preconditioner context
533: Output Parameter:
534: . flg - PETSC_TRUE if a transpose operation is defined
536: Level: developer
538: .keywords: PC, apply, transpose
540: .seealso: PCApplyTranspose()
541: @*/
542: PetscErrorCode PCApplyTransposeExists(PC pc,PetscBool *flg)
543: {
547: if (pc->ops->applytranspose) *flg = PETSC_TRUE;
548: else *flg = PETSC_FALSE;
549: return(0);
550: }
554: /*@
555: PCApplyBAorAB - Applies the preconditioner and operator to a vector. y = B*A*x or y = A*B*x.
557: Collective on PC and Vec
559: Input Parameters:
560: + pc - the preconditioner context
561: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
562: . x - input vector
563: - work - work vector
565: Output Parameter:
566: . y - output vector
568: Level: developer
570: 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
571: specific KSPSolve() method must also be written to handle the post-solve "correction" for the diagonal scaling.
573: .keywords: PC, apply, operator
575: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
576: @*/
577: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
578: {
586: if (x == y) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
587: VecValidValues(x,3,PETSC_TRUE);
588: if (side != PC_LEFT && side != PC_SYMMETRIC && side != PC_RIGHT) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Side must be right, left, or symmetric");
589: if (pc->diagonalscale && side == PC_SYMMETRIC) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
591: if (pc->setupcalled < 2) {
592: PCSetUp(pc);
593: }
595: if (pc->diagonalscale) {
596: if (pc->ops->applyBA) {
597: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
598: VecDuplicate(x,&work2);
599: PCDiagonalScaleRight(pc,x,work2);
600: (*pc->ops->applyBA)(pc,side,work2,y,work);
601: PCDiagonalScaleLeft(pc,y,y);
602: VecDestroy(&work2);
603: } else if (side == PC_RIGHT) {
604: PCDiagonalScaleRight(pc,x,y);
605: PCApply(pc,y,work);
606: MatMult(pc->mat,work,y);
607: PCDiagonalScaleLeft(pc,y,y);
608: } else if (side == PC_LEFT) {
609: PCDiagonalScaleRight(pc,x,y);
610: MatMult(pc->mat,y,work);
611: PCApply(pc,work,y);
612: PCDiagonalScaleLeft(pc,y,y);
613: } else if (side == PC_SYMMETRIC) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
614: } else {
615: if (pc->ops->applyBA) {
616: (*pc->ops->applyBA)(pc,side,x,y,work);
617: } else if (side == PC_RIGHT) {
618: PCApply(pc,x,work);
619: MatMult(pc->mat,work,y);
620: } else if (side == PC_LEFT) {
621: MatMult(pc->mat,x,work);
622: PCApply(pc,work,y);
623: } else if (side == PC_SYMMETRIC) {
624: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
625: PCApplySymmetricRight(pc,x,work);
626: MatMult(pc->mat,work,y);
627: VecCopy(y,work);
628: PCApplySymmetricLeft(pc,work,y);
629: }
630: }
631: VecValidValues(y,4,PETSC_FALSE);
632: return(0);
633: }
637: /*@
638: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
639: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
640: NOT tr(B*A) = tr(A)*tr(B).
642: Collective on PC and Vec
644: Input Parameters:
645: + pc - the preconditioner context
646: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
647: . x - input vector
648: - work - work vector
650: Output Parameter:
651: . y - output vector
654: 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
655: defined by B'. This is why this has the funny form that it computes tr(B) * tr(A)
656:
657: Level: developer
659: .keywords: PC, apply, operator, transpose
661: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
662: @*/
663: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
664: {
672: if (x == y) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_IDN,"x and y must be different vectors");
673: VecValidValues(x,3,PETSC_TRUE);
674: if (pc->ops->applyBAtranspose) {
675: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
676: VecValidValues(y,4,PETSC_FALSE);
677: return(0);
678: }
679: if (side != PC_LEFT && side != PC_RIGHT) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
681: if (pc->setupcalled < 2) {
682: PCSetUp(pc);
683: }
685: if (side == PC_RIGHT) {
686: PCApplyTranspose(pc,x,work);
687: MatMultTranspose(pc->mat,work,y);
688: } else if (side == PC_LEFT) {
689: MatMultTranspose(pc->mat,x,work);
690: PCApplyTranspose(pc,work,y);
691: }
692: /* add support for PC_SYMMETRIC */
693: VecValidValues(y,4,PETSC_FALSE);
694: return(0);
695: }
697: /* -------------------------------------------------------------------------------*/
701: /*@
702: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
703: built-in fast application of Richardson's method.
705: Not Collective
707: Input Parameter:
708: . pc - the preconditioner
710: Output Parameter:
711: . exists - PETSC_TRUE or PETSC_FALSE
713: Level: developer
715: .keywords: PC, apply, Richardson, exists
717: .seealso: PCApplyRichardson()
718: @*/
719: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscBool *exists)
720: {
724: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
725: else *exists = PETSC_FALSE;
726: return(0);
727: }
731: /*@
732: PCApplyRichardson - Applies several steps of Richardson iteration with
733: the particular preconditioner. This routine is usually used by the
734: Krylov solvers and not the application code directly.
736: Collective on PC
738: Input Parameters:
739: + pc - the preconditioner context
740: . b - the right hand side
741: . w - one work vector
742: . rtol - relative decrease in residual norm convergence criteria
743: . abstol - absolute residual norm convergence criteria
744: . dtol - divergence residual norm increase criteria
745: . its - the number of iterations to apply.
746: - guesszero - if the input x contains nonzero initial guess
748: Output Parameter:
749: + outits - number of iterations actually used (for SOR this always equals its)
750: . reason - the reason the apply terminated
751: - y - the solution (also contains initial guess if guesszero is PETSC_FALSE
753: Notes:
754: Most preconditioners do not support this function. Use the command
755: PCApplyRichardsonExists() to determine if one does.
757: Except for the multigrid PC this routine ignores the convergence tolerances
758: and always runs for the number of iterations
759:
760: Level: developer
762: .keywords: PC, apply, Richardson
764: .seealso: PCApplyRichardsonExists()
765: @*/
766: PetscErrorCode PCApplyRichardson(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
767: {
775: if (b == y) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_IDN,"b and y must be different vectors");
776: if (pc->setupcalled < 2) {
777: PCSetUp(pc);
778: }
779: if (!pc->ops->applyrichardson) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"PC does not have apply richardson");
780: (*pc->ops->applyrichardson)(pc,b,y,w,rtol,abstol,dtol,its,guesszero,outits,reason);
781: return(0);
782: }
784: /*
785: a setupcall of 0 indicates never setup,
786: 1 needs to be resetup,
787: 2 does not need any changes.
788: */
791: /*@
792: PCSetUp - Prepares for the use of a preconditioner.
794: Collective on PC
796: Input Parameter:
797: . pc - the preconditioner context
799: Level: developer
801: .keywords: PC, setup
803: .seealso: PCCreate(), PCApply(), PCDestroy()
804: @*/
805: PetscErrorCode PCSetUp(PC pc)
806: {
808: const char *def;
812: if (!pc->mat) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");
814: if (pc->setupcalled > 1) {
815: PetscInfo(pc,"Setting PC with identical preconditioner\n");
816: return(0);
817: } else if (!pc->setupcalled) {
818: PetscInfo(pc,"Setting up new PC\n");
819: } else if (pc->flag == SAME_NONZERO_PATTERN) {
820: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
821: } else {
822: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
823: }
825: if (!((PetscObject)pc)->type_name) {
826: PCGetDefaultType_Private(pc,&def);
827: PCSetType(pc,def);
828: }
830: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
831: if (pc->ops->setup) {
832: (*pc->ops->setup)(pc);
833: }
834: pc->setupcalled = 2;
835: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
836: return(0);
837: }
841: /*@
842: PCSetUpOnBlocks - Sets up the preconditioner for each block in
843: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
844: methods.
846: Collective on PC
848: Input Parameters:
849: . pc - the preconditioner context
851: Level: developer
853: .keywords: PC, setup, blocks
855: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
856: @*/
857: PetscErrorCode PCSetUpOnBlocks(PC pc)
858: {
863: if (!pc->ops->setuponblocks) return(0);
864: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
865: (*pc->ops->setuponblocks)(pc);
866: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
867: return(0);
868: }
872: /*@C
873: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
874: submatrices that arise within certain subdomain-based preconditioners.
875: The basic submatrices are extracted from the preconditioner matrix as
876: usual; the user can then alter these (for example, to set different boundary
877: conditions for each submatrix) before they are used for the local solves.
879: Logically Collective on PC
881: Input Parameters:
882: + pc - the preconditioner context
883: . func - routine for modifying the submatrices
884: - ctx - optional user-defined context (may be null)
886: Calling sequence of func:
887: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
889: . row - an array of index sets that contain the global row numbers
890: that comprise each local submatrix
891: . col - an array of index sets that contain the global column numbers
892: that comprise each local submatrix
893: . submat - array of local submatrices
894: - ctx - optional user-defined context for private data for the
895: user-defined func routine (may be null)
897: Notes:
898: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
899: KSPSolve().
901: A routine set by PCSetModifySubMatrices() is currently called within
902: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
903: preconditioners. All other preconditioners ignore this routine.
905: Level: advanced
907: .keywords: PC, set, modify, submatrices
909: .seealso: PCModifySubMatrices(), PCASMGetSubMatrices()
910: @*/
911: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
912: {
915: pc->modifysubmatrices = func;
916: pc->modifysubmatricesP = ctx;
917: return(0);
918: }
922: /*@C
923: PCModifySubMatrices - Calls an optional user-defined routine within
924: certain preconditioners if one has been set with PCSetModifySubMarices().
926: Collective on PC
928: Input Parameters:
929: + pc - the preconditioner context
930: . nsub - the number of local submatrices
931: . row - an array of index sets that contain the global row numbers
932: that comprise each local submatrix
933: . col - an array of index sets that contain the global column numbers
934: that comprise each local submatrix
935: . submat - array of local submatrices
936: - ctx - optional user-defined context for private data for the
937: user-defined routine (may be null)
939: Output Parameter:
940: . submat - array of local submatrices (the entries of which may
941: have been modified)
943: Notes:
944: The user should NOT generally call this routine, as it will
945: automatically be called within certain preconditioners (currently
946: block Jacobi, additive Schwarz) if set.
948: The basic submatrices are extracted from the preconditioner matrix
949: as usual; the user can then alter these (for example, to set different
950: boundary conditions for each submatrix) before they are used for the
951: local solves.
953: Level: developer
955: .keywords: PC, modify, submatrices
957: .seealso: PCSetModifySubMatrices()
958: @*/
959: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
960: {
965: if (!pc->modifysubmatrices) return(0);
966: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
967: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
968: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
969: return(0);
970: }
974: /*@
975: PCSetOperators - Sets the matrix associated with the linear system and
976: a (possibly) different one associated with the preconditioner.
978: Logically Collective on PC and Mat
980: Input Parameters:
981: + pc - the preconditioner context
982: . Amat - the matrix associated with the linear system
983: . Pmat - the matrix to be used in constructing the preconditioner, usually
984: the same as Amat.
985: - flag - flag indicating information about the preconditioner matrix structure
986: during successive linear solves. This flag is ignored the first time a
987: linear system is solved, and thus is irrelevant when solving just one linear
988: system.
990: Notes:
991: The flag can be used to eliminate unnecessary work in the preconditioner
992: during the repeated solution of linear systems of the same size. The
993: available options are
994: + SAME_PRECONDITIONER -
995: Pmat is identical during successive linear solves.
996: This option is intended for folks who are using
997: different Amat and Pmat matrices and wish to reuse the
998: same preconditioner matrix. For example, this option
999: saves work by not recomputing incomplete factorization
1000: for ILU/ICC preconditioners.
1001: . SAME_NONZERO_PATTERN -
1002: Pmat has the same nonzero structure during
1003: successive linear solves.
1004: - DIFFERENT_NONZERO_PATTERN -
1005: Pmat does not have the same nonzero structure.
1007: Passing a PETSC_NULL for Amat or Pmat removes the matrix that is currently used.
1009: If you wish to replace either Amat or Pmat but leave the other one untouched then
1010: first call KSPGetOperators() to get the one you wish to keep, call PetscObjectReference()
1011: on it and then pass it back in in your call to KSPSetOperators().
1013: Caution:
1014: If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion
1015: and does not check the structure of the matrix. If you erroneously
1016: claim that the structure is the same when it actually is not, the new
1017: preconditioner will not function correctly. Thus, use this optimization
1018: feature carefully!
1020: If in doubt about whether your preconditioner matrix has changed
1021: structure or not, use the flag DIFFERENT_NONZERO_PATTERN.
1023: More Notes about Repeated Solution of Linear Systems:
1024: PETSc does NOT reset the matrix entries of either Amat or Pmat
1025: to zero after a linear solve; the user is completely responsible for
1026: matrix assembly. See the routine MatZeroEntries() if desiring to
1027: zero all elements of a matrix.
1029: Level: intermediate
1031: .keywords: PC, set, operators, matrix, linear system
1033: .seealso: PCGetOperators(), MatZeroEntries()
1034: @*/
1035: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat,MatStructure flag)
1036: {
1038: PetscInt m1,n1,m2,n2;
1046: if (pc->setupcalled && Amat && Pmat) {
1047: MatGetLocalSize(Amat,&m1,&n1);
1048: MatGetLocalSize(pc->mat,&m2,&n2);
1049: 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);
1050: MatGetLocalSize(Pmat,&m1,&n1);
1051: MatGetLocalSize(pc->pmat,&m2,&n2);
1052: 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);
1053: }
1055: /* reference first in case the matrices are the same */
1056: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1057: MatDestroy(&pc->mat);
1058: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1059: MatDestroy(&pc->pmat);
1060: pc->mat = Amat;
1061: pc->pmat = Pmat;
1063: if (pc->setupcalled == 2 && flag != SAME_PRECONDITIONER) {
1064: pc->setupcalled = 1;
1065: }
1066: pc->flag = flag;
1067: return(0);
1068: }
1072: /*@C
1073: PCGetOperators - Gets the matrix associated with the linear system and
1074: possibly a different one associated with the preconditioner.
1076: Not collective, though parallel Mats are returned if the PC is parallel
1078: Input Parameter:
1079: . pc - the preconditioner context
1081: Output Parameters:
1082: + mat - the matrix associated with the linear system
1083: . pmat - matrix associated with the preconditioner, usually the same
1084: as mat.
1085: - flag - flag indicating information about the preconditioner
1086: matrix structure. See PCSetOperators() for details.
1088: Level: intermediate
1090: Notes: Does not increase the reference count of the matrices, so you should not destroy them
1092: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1093: are created in PC and returned to the user. In this case, if both operators
1094: mat and pmat are requested, two DIFFERENT operators will be returned. If
1095: only one is requested both operators in the PC will be the same (i.e. as
1096: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1097: The user must set the sizes of the returned matrices and their type etc just
1098: as if the user created them with MatCreate(). For example,
1100: $ KSP/PCGetOperators(ksp/pc,&mat,PETSC_NULL,PETSC_NULL); is equivalent to
1101: $ set size, type, etc of mat
1103: $ MatCreate(comm,&mat);
1104: $ KSP/PCSetOperators(ksp/pc,mat,mat,SAME_NONZERO_PATTERN);
1105: $ PetscObjectDereference((PetscObject)mat);
1106: $ set size, type, etc of mat
1108: and
1110: $ KSP/PCGetOperators(ksp/pc,&mat,&pmat,PETSC_NULL); is equivalent to
1111: $ set size, type, etc of mat and pmat
1113: $ MatCreate(comm,&mat);
1114: $ MatCreate(comm,&pmat);
1115: $ KSP/PCSetOperators(ksp/pc,mat,pmat,SAME_NONZERO_PATTERN);
1116: $ PetscObjectDereference((PetscObject)mat);
1117: $ PetscObjectDereference((PetscObject)pmat);
1118: $ set size, type, etc of mat and pmat
1120: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1121: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1122: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1123: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1124: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1125: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1126: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1127: it can be created for you?
1128:
1130: .keywords: PC, get, operators, matrix, linear system
1132: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1133: @*/
1134: PetscErrorCode PCGetOperators(PC pc,Mat *mat,Mat *pmat,MatStructure *flag)
1135: {
1140: if (mat) {
1141: if (!pc->mat) {
1142: if (pc->pmat && !pmat) { /* pmat has been set, but user did not request it, so use for mat */
1143: pc->mat = pc->pmat;
1144: PetscObjectReference((PetscObject)pc->mat);
1145: } else { /* both mat and pmat are empty */
1146: MatCreate(((PetscObject)pc)->comm,&pc->mat);
1147: if (!pmat) { /* user did NOT request pmat, so make same as mat */
1148: pc->pmat = pc->mat;
1149: PetscObjectReference((PetscObject)pc->pmat);
1150: }
1151: }
1152: }
1153: *mat = pc->mat;
1154: }
1155: if (pmat) {
1156: if (!pc->pmat) {
1157: if (pc->mat && !mat) { /* mat has been set but was not requested, so use for pmat */
1158: pc->pmat = pc->mat;
1159: PetscObjectReference((PetscObject)pc->pmat);
1160: } else {
1161: MatCreate(((PetscObject)pc)->comm,&pc->pmat);
1162: if (!mat) { /* user did NOT request mat, so make same as pmat */
1163: pc->mat = pc->pmat;
1164: PetscObjectReference((PetscObject)pc->mat);
1165: }
1166: }
1167: }
1168: *pmat = pc->pmat;
1169: }
1170: if (flag) *flag = pc->flag;
1171: return(0);
1172: }
1176: /*@C
1177: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1178: possibly a different one associated with the preconditioner have been set in the PC.
1180: Not collective, though the results on all processes should be the same
1182: Input Parameter:
1183: . pc - the preconditioner context
1185: Output Parameters:
1186: + mat - the matrix associated with the linear system was set
1187: - pmat - matrix associated with the preconditioner was set, usually the same
1189: Level: intermediate
1191: .keywords: PC, get, operators, matrix, linear system
1193: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1194: @*/
1195: PetscErrorCode PCGetOperatorsSet(PC pc,PetscBool *mat,PetscBool *pmat)
1196: {
1199: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1200: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1201: return(0);
1202: }
1206: /*@
1207: PCFactorGetMatrix - Gets the factored matrix from the
1208: preconditioner context. This routine is valid only for the LU,
1209: incomplete LU, Cholesky, and incomplete Cholesky methods.
1211: Not Collective on PC though Mat is parallel if PC is parallel
1213: Input Parameters:
1214: . pc - the preconditioner context
1216: Output parameters:
1217: . mat - the factored matrix
1219: Level: advanced
1221: Notes: Does not increase the reference count for the matrix so DO NOT destroy it
1223: .keywords: PC, get, factored, matrix
1224: @*/
1225: PetscErrorCode PCFactorGetMatrix(PC pc,Mat *mat)
1226: {
1232: if (pc->ops->getfactoredmatrix) {
1233: (*pc->ops->getfactoredmatrix)(pc,mat);
1234: } else {
1235: SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"PC type does not support getting factor matrix");
1236: }
1237: return(0);
1238: }
1242: /*@C
1243: PCSetOptionsPrefix - Sets the prefix used for searching for all
1244: PC options in the database.
1246: Logically Collective on PC
1248: Input Parameters:
1249: + pc - the preconditioner context
1250: - prefix - the prefix string to prepend to all PC option requests
1252: Notes:
1253: A hyphen (-) must NOT be given at the beginning of the prefix name.
1254: The first character of all runtime options is AUTOMATICALLY the
1255: hyphen.
1257: Level: advanced
1259: .keywords: PC, set, options, prefix, database
1261: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1262: @*/
1263: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1264: {
1269: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1270: return(0);
1271: }
1275: /*@C
1276: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1277: PC options in the database.
1279: Logically Collective on PC
1281: Input Parameters:
1282: + pc - the preconditioner context
1283: - prefix - the prefix string to prepend to all PC option requests
1285: Notes:
1286: A hyphen (-) must NOT be given at the beginning of the prefix name.
1287: The first character of all runtime options is AUTOMATICALLY the
1288: hyphen.
1290: Level: advanced
1292: .keywords: PC, append, options, prefix, database
1294: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1295: @*/
1296: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1297: {
1302: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1303: return(0);
1304: }
1308: /*@C
1309: PCGetOptionsPrefix - Gets the prefix used for searching for all
1310: PC options in the database.
1312: Not Collective
1314: Input Parameters:
1315: . pc - the preconditioner context
1317: Output Parameters:
1318: . prefix - pointer to the prefix string used, is returned
1320: Notes: On the fortran side, the user should pass in a string 'prifix' of
1321: sufficient length to hold the prefix.
1323: Level: advanced
1325: .keywords: PC, get, options, prefix, database
1327: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1328: @*/
1329: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1330: {
1336: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1337: return(0);
1338: }
1342: /*@
1343: PCPreSolve - Optional pre-solve phase, intended for any
1344: preconditioner-specific actions that must be performed before
1345: the iterative solve itself.
1347: Collective on PC
1349: Input Parameters:
1350: + pc - the preconditioner context
1351: - ksp - the Krylov subspace context
1353: Level: developer
1355: Sample of Usage:
1356: .vb
1357: PCPreSolve(pc,ksp);
1358: KSPSolve(ksp,b,x);
1359: PCPostSolve(pc,ksp);
1360: .ve
1362: Notes:
1363: The pre-solve phase is distinct from the PCSetUp() phase.
1365: KSPSolve() calls this directly, so is rarely called by the user.
1367: .keywords: PC, pre-solve
1369: .seealso: PCPostSolve()
1370: @*/
1371: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1372: {
1374: Vec x,rhs;
1379: pc->presolvedone++;
1380: if (pc->presolvedone > 2) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Cannot embed PCPreSolve() more than twice");
1381: KSPGetSolution(ksp,&x);
1382: KSPGetRhs(ksp,&rhs);
1384: if (pc->ops->presolve) {
1385: (*pc->ops->presolve)(pc,ksp,rhs,x);
1386: }
1387: return(0);
1388: }
1392: /*@
1393: PCPostSolve - Optional post-solve phase, intended for any
1394: preconditioner-specific actions that must be performed after
1395: the iterative solve itself.
1397: Collective on PC
1399: Input Parameters:
1400: + pc - the preconditioner context
1401: - ksp - the Krylov subspace context
1403: Sample of Usage:
1404: .vb
1405: PCPreSolve(pc,ksp);
1406: KSPSolve(ksp,b,x);
1407: PCPostSolve(pc,ksp);
1408: .ve
1410: Note:
1411: KSPSolve() calls this routine directly, so it is rarely called by the user.
1413: Level: developer
1415: .keywords: PC, post-solve
1417: .seealso: PCPreSolve(), KSPSolve()
1418: @*/
1419: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1420: {
1422: Vec x,rhs;
1427: pc->presolvedone--;
1428: KSPGetSolution(ksp,&x);
1429: KSPGetRhs(ksp,&rhs);
1430: if (pc->ops->postsolve) {
1431: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1432: }
1433: return(0);
1434: }
1438: /*@C
1439: PCView - Prints the PC data structure.
1441: Collective on PC
1443: Input Parameters:
1444: + PC - the PC context
1445: - viewer - optional visualization context
1447: Note:
1448: The available visualization contexts include
1449: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1450: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1451: output where only the first processor opens
1452: the file. All other processors send their
1453: data to the first processor to print.
1455: The user can open an alternative visualization contexts with
1456: PetscViewerASCIIOpen() (output to a specified file).
1458: Level: developer
1460: .keywords: PC, view
1462: .seealso: KSPView(), PetscViewerASCIIOpen()
1463: @*/
1464: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1465: {
1466: const PCType cstr;
1467: PetscErrorCode ierr;
1468: PetscBool iascii,isstring;
1469: PetscViewerFormat format;
1473: if (!viewer) {
1474: PetscViewerASCIIGetStdout(((PetscObject)pc)->comm,&viewer);
1475: }
1479: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1480: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
1481: if (iascii) {
1482: PetscViewerGetFormat(viewer,&format);
1483: PetscObjectPrintClassNamePrefixType((PetscObject)pc,viewer,"PC Object");
1484: if (pc->ops->view) {
1485: PetscViewerASCIIPushTab(viewer);
1486: (*pc->ops->view)(pc,viewer);
1487: PetscViewerASCIIPopTab(viewer);
1488: }
1489: if (pc->mat) {
1490: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1491: if (pc->pmat == pc->mat) {
1492: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1493: PetscViewerASCIIPushTab(viewer);
1494: MatView(pc->mat,viewer);
1495: PetscViewerASCIIPopTab(viewer);
1496: } else {
1497: if (pc->pmat) {
1498: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1499: } else {
1500: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1501: }
1502: PetscViewerASCIIPushTab(viewer);
1503: MatView(pc->mat,viewer);
1504: if (pc->pmat) {MatView(pc->pmat,viewer);}
1505: PetscViewerASCIIPopTab(viewer);
1506: }
1507: PetscViewerPopFormat(viewer);
1508: }
1509: } else if (isstring) {
1510: PCGetType(pc,&cstr);
1511: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1512: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1513: } else {
1514: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported by PC",((PetscObject)viewer)->type_name);
1515: }
1516: return(0);
1517: }
1522: /*@
1523: PCSetInitialGuessNonzero - Tells the iterative solver that the
1524: initial guess is nonzero; otherwise PC assumes the initial guess
1525: is to be zero (and thus zeros it out before solving).
1527: Logically Collective on PC
1529: Input Parameters:
1530: + pc - iterative context obtained from PCCreate()
1531: - flg - PETSC_TRUE indicates the guess is non-zero, PETSC_FALSE indicates the guess is zero
1533: Level: Developer
1535: Notes:
1536: This is a weird function. Since PC's are linear operators on the right hand side they
1537: CANNOT use an initial guess. This function is for the "pass-through" preconditioners
1538: PCKSP, PCREDUNDANT and PCHMPI and causes the inner KSP object to use the nonzero
1539: initial guess. Not currently working for PCREDUNDANT, that has to be rewritten to use KSP.
1542: .keywords: PC, set, initial guess, nonzero
1544: .seealso: PCGetInitialGuessNonzero(), PCSetInitialGuessKnoll(), PCGetInitialGuessKnoll()
1545: @*/
1546: PetscErrorCode PCSetInitialGuessNonzero(PC pc,PetscBool flg)
1547: {
1550: pc->nonzero_guess = flg;
1551: return(0);
1552: }
1556: /*@C
1557: PCRegister - See PCRegisterDynamic()
1559: Level: advanced
1560: @*/
1561: PetscErrorCode PCRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(PC))
1562: {
1564: char fullname[PETSC_MAX_PATH_LEN];
1568: PetscFListConcat(path,name,fullname);
1569: PetscFListAdd(&PCList,sname,fullname,(void (*)(void))function);
1570: return(0);
1571: }
1575: /*@
1576: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1578: Collective on PC
1580: Input Parameter:
1581: . pc - the preconditioner object
1583: Output Parameter:
1584: . mat - the explict preconditioned operator
1586: Notes:
1587: This computation is done by applying the operators to columns of the
1588: identity matrix.
1590: Currently, this routine uses a dense matrix format when 1 processor
1591: is used and a sparse format otherwise. This routine is costly in general,
1592: and is recommended for use only with relatively small systems.
1594: Level: advanced
1595:
1596: .keywords: PC, compute, explicit, operator
1598: .seealso: KSPComputeExplicitOperator()
1600: @*/
1601: PetscErrorCode PCComputeExplicitOperator(PC pc,Mat *mat)
1602: {
1603: Vec in,out;
1605: PetscInt i,M,m,*rows,start,end;
1606: PetscMPIInt size;
1607: MPI_Comm comm;
1608: PetscScalar *array,one = 1.0;
1609:
1614: comm = ((PetscObject)pc)->comm;
1615: MPI_Comm_size(comm,&size);
1617: if (!pc->pmat) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1618: MatGetVecs(pc->pmat,&in,0);
1619: VecDuplicate(in,&out);
1620: VecGetOwnershipRange(in,&start,&end);
1621: VecGetSize(in,&M);
1622: VecGetLocalSize(in,&m);
1623: PetscMalloc((m+1)*sizeof(PetscInt),&rows);
1624: for (i=0; i<m; i++) {rows[i] = start + i;}
1626: MatCreate(comm,mat);
1627: MatSetSizes(*mat,m,m,M,M);
1628: if (size == 1) {
1629: MatSetType(*mat,MATSEQDENSE);
1630: MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
1631: } else {
1632: MatSetType(*mat,MATMPIAIJ);
1633: MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);
1634: }
1636: for (i=0; i<M; i++) {
1638: VecSet(in,0.0);
1639: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1640: VecAssemblyBegin(in);
1641: VecAssemblyEnd(in);
1643: /* should fix, allowing user to choose side */
1644: PCApply(pc,in,out);
1645:
1646: VecGetArray(out,&array);
1647: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1648: VecRestoreArray(out,&array);
1650: }
1651: PetscFree(rows);
1652: VecDestroy(&out);
1653: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1654: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1655: return(0);
1656: }
1660: /*@
1661: PCSetCoordinates - sets the coordinates of all the nodes on the local process
1663: Collective on PC
1665: Input Parameters:
1666: + pc - the solver context
1667: . dim - the dimension of the coordinates 1, 2, or 3
1668: - coords - the coordinates
1670: Level: intermediate
1672: Notes: coords is an array of the 3D coordinates for the nodes on
1673: the local processor. So if there are 108 equation on a processor
1674: for a displacement finite element discretization of elasticity (so
1675: that there are 36 = 108/3 nodes) then the array must have 108
1676: double precision values (ie, 3 * 36). These x y z coordinates
1677: should be ordered for nodes 0 to N-1 like so: [ 0.x, 0.y, 0.z, 1.x,
1678: ... , N-1.z ].
1680: .seealso: MatSetNearNullSpace
1681: @*/
1682: PetscErrorCode PCSetCoordinates(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords )
1683: {
1687: PetscTryMethod(pc,"PCSetCoordinates_C",(PC,PetscInt,PetscInt,PetscReal*),(pc,dim,nloc,coords));
1688: return(0);
1689: }