Actual source code: precon.c
1: #define PETSCKSP_DLL
3: /*
4: The PC (preconditioner) interface routines, callable by users.
5: */
6: #include private/pcimpl.h
8: /* Logging support */
9: PetscCookie PC_COOKIE = 0;
10: PetscEvent PC_SetUp = 0, PC_SetUpOnBlocks = 0, PC_Apply = 0, PC_ApplyCoarse = 0, PC_ApplyMultiple = 0, PC_ApplySymmetricLeft = 0;
11: PetscEvent PC_ApplySymmetricRight = 0, PC_ModifySubMatrices = 0;
15: PetscErrorCode PCGetDefaultType_Private(PC pc,const char* type[])
16: {
18: PetscMPIInt size;
19: PetscTruth flg1,flg2,set,flg3;
22: MPI_Comm_size(pc->comm,&size);
23: if (pc->pmat) {
24: PetscErrorCode (*f)(Mat,PetscTruth*,MatReuse,Mat*);
25: PetscObjectQueryFunction((PetscObject)pc->pmat,"MatGetDiagonalBlock_C",(void (**)(void))&f);
26: if (size == 1) {
27: MatHasOperation(pc->pmat,MATOP_ICCFACTOR_SYMBOLIC,&flg1);
28: MatHasOperation(pc->pmat,MATOP_ILUFACTOR_SYMBOLIC,&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: PCDestroy - Destroys PC context that was created with PCCreate().
61: Collective on PC
63: Input Parameter:
64: . pc - the preconditioner context
66: Level: developer
68: .keywords: PC, destroy
70: .seealso: PCCreate(), PCSetUp()
71: @*/
72: PetscErrorCode PCDestroy(PC pc)
73: {
78: if (--pc->refct > 0) return(0);
80: /* if memory was published with AMS then destroy it */
81: PetscObjectDepublish(pc);
83: if (pc->ops->destroy) { (*pc->ops->destroy)(pc);}
84: if (pc->diagonalscaleright) {VecDestroy(pc->diagonalscaleright);}
85: if (pc->diagonalscaleleft) {VecDestroy(pc->diagonalscaleleft);}
87: if (pc->pmat) {MatDestroy(pc->pmat);}
88: if (pc->mat) {MatDestroy(pc->mat);}
90: PetscHeaderDestroy(pc);
91: return(0);
92: }
96: /*@C
97: PCDiagonalScale - Indicates if the preconditioner applies an additional left and right
98: scaling as needed by certain time-stepping codes.
100: Collective on PC
102: Input Parameter:
103: . pc - the preconditioner context
105: Output Parameter:
106: . flag - PETSC_TRUE if it applies the scaling
108: Level: developer
110: Notes: If this returns PETSC_TRUE then the system solved via the Krylov method is
111: $ D M A D^{-1} y = D M b for left preconditioning or
112: $ D A M D^{-1} z = D b for right preconditioning
114: .keywords: PC
116: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCDiagonalScaleSet()
117: @*/
118: PetscErrorCode PCDiagonalScale(PC pc,PetscTruth *flag)
119: {
123: *flag = pc->diagonalscale;
124: return(0);
125: }
129: /*@
130: PCDiagonalScaleSet - Indicates the left scaling to use to apply an additional left and right
131: scaling as needed by certain time-stepping codes.
133: Collective on PC
135: Input Parameters:
136: + pc - the preconditioner context
137: - s - scaling vector
139: Level: intermediate
141: Notes: 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: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
147: .keywords: PC
149: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCDiagonalScale()
150: @*/
151: PetscErrorCode PCDiagonalScaleSet(PC pc,Vec s)
152: {
158: pc->diagonalscale = PETSC_TRUE;
159: if (pc->diagonalscaleleft) {
160: VecDestroy(pc->diagonalscaleleft);
161: }
162: pc->diagonalscaleleft = s;
163: PetscObjectReference((PetscObject)s);
164: if (!pc->diagonalscaleright) {
165: VecDuplicate(s,&pc->diagonalscaleright);
166: }
167: VecCopy(s,pc->diagonalscaleright);
168: VecReciprocal(pc->diagonalscaleright);
169: return(0);
170: }
174: /*@C
175: PCDiagonalScaleLeft - Indicates the left scaling to use to apply an additional left and right
176: scaling as needed by certain time-stepping codes.
178: Collective on PC
180: Input Parameters:
181: + pc - the preconditioner context
182: . in - input vector
183: + out - scaled vector (maybe the same as in)
185: Level: intermediate
187: Notes: The system solved via the Krylov method is
188: $ D M A D^{-1} y = D M b for left preconditioning or
189: $ D A M D^{-1} z = D b for right preconditioning
191: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
193: If diagonal scaling is turned off and in is not out then in is copied to out
195: .keywords: PC
197: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
198: @*/
199: PetscErrorCode PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
200: {
207: if (pc->diagonalscale) {
208: VecPointwiseMult(out,pc->diagonalscaleleft,in);
209: } else if (in != out) {
210: VecCopy(in,out);
211: }
212: return(0);
213: }
217: /*@C
218: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
220: Collective on PC
222: Input Parameters:
223: + pc - the preconditioner context
224: . in - input vector
225: + out - scaled vector (maybe the same as in)
227: Level: intermediate
229: Notes: The system solved via the Krylov method is
230: $ D M A D^{-1} y = D M b for left preconditioning or
231: $ D A M D^{-1} z = D b for right preconditioning
233: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
235: If diagonal scaling is turned off and in is not out then in is copied to out
237: .keywords: PC
239: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
240: @*/
241: PetscErrorCode PCDiagonalScaleRight(PC pc,Vec in,Vec out)
242: {
249: if (pc->diagonalscale) {
250: VecPointwiseMult(out,pc->diagonalscaleright,in);
251: } else if (in != out) {
252: VecCopy(in,out);
253: }
254: return(0);
255: }
259: static PetscErrorCode PCPublish_Petsc(PetscObject obj)
260: {
262: return(0);
263: }
267: /*@
268: PCCreate - Creates a preconditioner context.
270: Collective on MPI_Comm
272: Input Parameter:
273: . comm - MPI communicator
275: Output Parameter:
276: . pc - location to put the preconditioner context
278: Notes:
279: The default preconditioner on one processor is PCILU with 0 fill on more
280: then one it is PCBJACOBI with ILU() on each processor.
282: Level: developer
284: .keywords: PC, create, context
286: .seealso: PCSetUp(), PCApply(), PCDestroy()
287: @*/
288: PetscErrorCode PCCreate(MPI_Comm comm,PC *newpc)
289: {
290: PC pc;
295: *newpc = 0;
296: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
297: PCInitializePackage(PETSC_NULL);
298: #endif
300: PetscHeaderCreate(pc,_p_PC,struct _PCOps,PC_COOKIE,-1,"PC",comm,PCDestroy,PCView);
301: pc->bops->publish = PCPublish_Petsc;
302: pc->mat = 0;
303: pc->pmat = 0;
304: pc->setupcalled = 0;
305: pc->data = 0;
306: pc->diagonalscale = PETSC_FALSE;
307: pc->diagonalscaleleft = 0;
308: pc->diagonalscaleright = 0;
310: pc->ops->destroy = 0;
311: pc->ops->apply = 0;
312: pc->ops->applytranspose = 0;
313: pc->ops->applyBA = 0;
314: pc->ops->applyBAtranspose = 0;
315: pc->ops->applyrichardson = 0;
316: pc->ops->view = 0;
317: pc->ops->getfactoredmatrix = 0;
318: pc->ops->applysymmetricright = 0;
319: pc->ops->applysymmetricleft = 0;
320: pc->ops->setuponblocks = 0;
322: pc->modifysubmatrices = 0;
323: pc->modifysubmatricesP = 0;
324: *newpc = pc;
325: PetscPublishAll(pc);
326: return(0);
328: }
330: /* -------------------------------------------------------------------------------*/
334: /*@
335: PCApply - Applies the preconditioner to a vector.
337: Collective on PC and Vec
339: Input Parameters:
340: + pc - the preconditioner context
341: - x - input vector
343: Output Parameter:
344: . y - output vector
346: Level: developer
348: .keywords: PC, apply
350: .seealso: PCApplyTranspose(), PCApplyBAorAB()
351: @*/
352: PetscErrorCode PCApply(PC pc,Vec x,Vec y)
353: {
360: if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
362: if (pc->setupcalled < 2) {
363: PCSetUp(pc);
364: }
366: (*pc->ops->apply)(pc,x,y);
368: return(0);
369: }
373: /*@
374: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
376: Collective on PC and Vec
378: Input Parameters:
379: + pc - the preconditioner context
380: - x - input vector
382: Output Parameter:
383: . y - output vector
385: Notes:
386: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
388: Level: developer
390: .keywords: PC, apply, symmetric, left
392: .seealso: PCApply(), PCApplySymmetricRight()
393: @*/
394: PetscErrorCode PCApplySymmetricLeft(PC pc,Vec x,Vec y)
395: {
402: if (!pc->ops->applysymmetricleft) SETERRQ(PETSC_ERR_SUP,"PC does not have left symmetric apply");
404: if (pc->setupcalled < 2) {
405: PCSetUp(pc);
406: }
409: (*pc->ops->applysymmetricleft)(pc,x,y);
411: return(0);
412: }
416: /*@
417: PCApplySymmetricRight - Applies the right part of a symmetric preconditioner to a vector.
419: Collective on PC and Vec
421: Input Parameters:
422: + pc - the preconditioner context
423: - x - input vector
425: Output Parameter:
426: . y - output vector
428: Level: developer
430: Notes:
431: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
433: .keywords: PC, apply, symmetric, right
435: .seealso: PCApply(), PCApplySymmetricLeft()
436: @*/
437: PetscErrorCode PCApplySymmetricRight(PC pc,Vec x,Vec y)
438: {
445: if (!pc->ops->applysymmetricright) SETERRQ(PETSC_ERR_SUP,"PC does not have left symmetric apply");
447: if (pc->setupcalled < 2) {
448: PCSetUp(pc);
449: }
452: (*pc->ops->applysymmetricright)(pc,x,y);
454: return(0);
455: }
459: /*@
460: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
462: Collective on PC and Vec
464: Input Parameters:
465: + pc - the preconditioner context
466: - x - input vector
468: Output Parameter:
469: . y - output vector
471: Level: developer
473: .keywords: PC, apply, transpose
475: .seealso: PCApply(), PCApplyBAorAB(), PCApplyBAorABTranspose(), PCHasApplyTranspose()
476: @*/
477: PetscErrorCode PCApplyTranspose(PC pc,Vec x,Vec y)
478: {
485: if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
486: if (!pc->ops->applytranspose) SETERRQ(PETSC_ERR_SUP," ");
488: if (pc->setupcalled < 2) {
489: PCSetUp(pc);
490: }
493: (*pc->ops->applytranspose)(pc,x,y);
495: return(0);
496: }
500: /*@
501: PCHasApplyTranspose - Test whether the preconditioner has a transpose apply operation
503: Collective on PC and Vec
505: Input Parameters:
506: . pc - the preconditioner context
508: Output Parameter:
509: . flg - PETSC_TRUE if a transpose operation is defined
511: Level: developer
513: .keywords: PC, apply, transpose
515: .seealso: PCApplyTranspose()
516: @*/
517: PetscErrorCode PCHasApplyTranspose(PC pc,PetscTruth *flg)
518: {
522: *flg = (PetscTruth) (pc->ops->applytranspose == 0);
523: return(0);
524: }
528: /*@
529: PCApplyBAorAB - Applies the preconditioner and operator to a vector.
531: Collective on PC and Vec
533: Input Parameters:
534: + pc - the preconditioner context
535: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
536: . x - input vector
537: - work - work vector
539: Output Parameter:
540: . y - output vector
542: Level: developer
544: .keywords: PC, apply, operator
546: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
547: @*/
548: PetscErrorCode PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
549: {
557: if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
558: if (side != PC_LEFT && side != PC_SYMMETRIC && side != PC_RIGHT) {
559: SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Side must be right, left, or symmetric");
560: }
561: if (pc->diagonalscale && side == PC_SYMMETRIC) {
562: SETERRQ(PETSC_ERR_SUP,"Cannot include diagonal scaling with symmetric preconditioner application");
563: }
565: if (pc->setupcalled < 2) {
566: PCSetUp(pc);
567: }
569: if (pc->diagonalscale) {
570: if (pc->ops->applyBA) {
571: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
572: VecDuplicate(x,&work2);
573: PCDiagonalScaleRight(pc,x,work2);
574: (*pc->ops->applyBA)(pc,side,work2,y,work);
575: PCDiagonalScaleLeft(pc,y,y);
576: VecDestroy(work2);
577: } else if (side == PC_RIGHT) {
578: PCDiagonalScaleRight(pc,x,y);
579: PCApply(pc,y,work);
580: MatMult(pc->mat,work,y);
581: PCDiagonalScaleLeft(pc,y,y);
582: } else if (side == PC_LEFT) {
583: PCDiagonalScaleRight(pc,x,y);
584: MatMult(pc->mat,y,work);
585: PCApply(pc,work,y);
586: PCDiagonalScaleLeft(pc,y,y);
587: } else if (side == PC_SYMMETRIC) {
588: SETERRQ(PETSC_ERR_SUP,"Cannot provide diagonal scaling with symmetric application of preconditioner");
589: }
590: } else {
591: if (pc->ops->applyBA) {
592: (*pc->ops->applyBA)(pc,side,x,y,work);
593: } else if (side == PC_RIGHT) {
594: PCApply(pc,x,work);
595: MatMult(pc->mat,work,y);
596: } else if (side == PC_LEFT) {
597: MatMult(pc->mat,x,work);
598: PCApply(pc,work,y);
599: } else if (side == PC_SYMMETRIC) {
600: /* There's an extra copy here; maybe should provide 2 work vectors instead? */
601: PCApplySymmetricRight(pc,x,work);
602: MatMult(pc->mat,work,y);
603: VecCopy(y,work);
604: PCApplySymmetricLeft(pc,work,y);
605: }
606: }
607: return(0);
608: }
612: /*@
613: PCApplyBAorABTranspose - Applies the transpose of the preconditioner
614: and operator to a vector. That is, applies tr(B) * tr(A) with left preconditioning,
615: not tr(B*A) = tr(A)*tr(B).
617: Collective on PC and Vec
619: Input Parameters:
620: + pc - the preconditioner context
621: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
622: . x - input vector
623: - work - work vector
625: Output Parameter:
626: . y - output vector
628: Level: developer
630: .keywords: PC, apply, operator, transpose
632: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorAB()
633: @*/
634: PetscErrorCode PCApplyBAorABTranspose(PC pc,PCSide side,Vec x,Vec y,Vec work)
635: {
643: if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
644: if (pc->ops->applyBAtranspose) {
645: (*pc->ops->applyBAtranspose)(pc,side,x,y,work);
646: return(0);
647: }
648: if (side != PC_LEFT && side != PC_RIGHT) {
649: SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Side must be right or left");
650: }
652: if (pc->setupcalled < 2) {
653: PCSetUp(pc);
654: }
656: if (side == PC_RIGHT) {
657: PCApplyTranspose(pc,x,work);
658: MatMultTranspose(pc->mat,work,y);
659: } else if (side == PC_LEFT) {
660: MatMultTranspose(pc->mat,x,work);
661: PCApplyTranspose(pc,work,y);
662: }
663: /* add support for PC_SYMMETRIC */
664: return(0); /* actually will never get here */
665: }
667: /* -------------------------------------------------------------------------------*/
671: /*@
672: PCApplyRichardsonExists - Determines whether a particular preconditioner has a
673: built-in fast application of Richardson's method.
675: Not Collective
677: Input Parameter:
678: . pc - the preconditioner
680: Output Parameter:
681: . exists - PETSC_TRUE or PETSC_FALSE
683: Level: developer
685: .keywords: PC, apply, Richardson, exists
687: .seealso: PCApplyRichardson()
688: @*/
689: PetscErrorCode PCApplyRichardsonExists(PC pc,PetscTruth *exists)
690: {
694: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
695: else *exists = PETSC_FALSE;
696: return(0);
697: }
701: /*@
702: PCApplyRichardson - Applies several steps of Richardson iteration with
703: the particular preconditioner. This routine is usually used by the
704: Krylov solvers and not the application code directly.
706: Collective on PC
708: Input Parameters:
709: + pc - the preconditioner context
710: . x - the initial guess
711: . w - one work vector
712: . rtol - relative decrease in residual norm convergence criteria
713: . abstol - absolute residual norm convergence criteria
714: . dtol - divergence residual norm increase criteria
715: - its - the number of iterations to apply.
717: Output Parameter:
718: . y - the solution
720: Notes:
721: Most preconditioners do not support this function. Use the command
722: PCApplyRichardsonExists() to determine if one does.
724: Except for the multigrid PC this routine ignores the convergence tolerances
725: and always runs for the number of iterations
726:
727: Level: developer
729: .keywords: PC, apply, Richardson
731: .seealso: PCApplyRichardsonExists()
732: @*/
733: PetscErrorCode PCApplyRichardson(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its)
734: {
742: if (!pc->ops->applyrichardson) SETERRQ(PETSC_ERR_SUP," ");
744: if (pc->setupcalled < 2) {
745: PCSetUp(pc);
746: }
748: (*pc->ops->applyrichardson)(pc,x,y,w,rtol,abstol,dtol,its);
749: return(0);
750: }
752: /*
753: a setupcall of 0 indicates never setup,
754: 1 needs to be resetup,
755: 2 does not need any changes.
756: */
759: /*@
760: PCSetUp - Prepares for the use of a preconditioner.
762: Collective on PC
764: Input Parameter:
765: . pc - the preconditioner context
767: Level: developer
769: .keywords: PC, setup
771: .seealso: PCCreate(), PCApply(), PCDestroy()
772: @*/
773: PetscErrorCode PCSetUp(PC pc)
774: {
776: const char *def;
781: if (pc->setupcalled > 1) {
782: PetscInfo(pc,"Setting PC with identical preconditioner\n");
783: return(0);
784: } else if (!pc->setupcalled) {
785: PetscInfo(pc,"Setting up new PC\n");
786: } else if (pc->flag == SAME_NONZERO_PATTERN) {
787: PetscInfo(pc,"Setting up PC with same nonzero pattern\n");
788: } else {
789: PetscInfo(pc,"Setting up PC with different nonzero pattern\n");
790: }
793: if (!pc->mat) {SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");}
795: if (!pc->type_name) {
796: PCGetDefaultType_Private(pc,&def);
797: PCSetType(pc,def);
798: }
800: if (pc->ops->setup) {
801: (*pc->ops->setup)(pc);
802: }
803: pc->setupcalled = 2;
805: return(0);
806: }
810: /*@
811: PCSetUpOnBlocks - Sets up the preconditioner for each block in
812: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
813: methods.
815: Collective on PC
817: Input Parameters:
818: . pc - the preconditioner context
820: Level: developer
822: .keywords: PC, setup, blocks
824: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
825: @*/
826: PetscErrorCode PCSetUpOnBlocks(PC pc)
827: {
832: if (!pc->ops->setuponblocks) return(0);
834: (*pc->ops->setuponblocks)(pc);
836: return(0);
837: }
841: /*@C
842: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
843: submatrices that arise within certain subdomain-based preconditioners.
844: The basic submatrices are extracted from the preconditioner matrix as
845: usual; the user can then alter these (for example, to set different boundary
846: conditions for each submatrix) before they are used for the local solves.
848: Collective on PC
850: Input Parameters:
851: + pc - the preconditioner context
852: . func - routine for modifying the submatrices
853: - ctx - optional user-defined context (may be null)
855: Calling sequence of func:
856: $ func (PC pc,PetscInt nsub,IS *row,IS *col,Mat *submat,void *ctx);
858: . row - an array of index sets that contain the global row numbers
859: that comprise each local submatrix
860: . col - an array of index sets that contain the global column numbers
861: that comprise each local submatrix
862: . submat - array of local submatrices
863: - ctx - optional user-defined context for private data for the
864: user-defined func routine (may be null)
866: Notes:
867: PCSetModifySubMatrices() MUST be called before KSPSetUp() and
868: KSPSolve().
870: A routine set by PCSetModifySubMatrices() is currently called within
871: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
872: preconditioners. All other preconditioners ignore this routine.
874: Level: advanced
876: .keywords: PC, set, modify, submatrices
878: .seealso: PCModifySubMatrices()
879: @*/
880: PetscErrorCode PCSetModifySubMatrices(PC pc,PetscErrorCode (*func)(PC,PetscInt,const IS[],const IS[],Mat[],void*),void *ctx)
881: {
884: pc->modifysubmatrices = func;
885: pc->modifysubmatricesP = ctx;
886: return(0);
887: }
891: /*@C
892: PCModifySubMatrices - Calls an optional user-defined routine within
893: certain preconditioners if one has been set with PCSetModifySubMarices().
895: Collective on PC
897: Input Parameters:
898: + pc - the preconditioner context
899: . nsub - the number of local submatrices
900: . row - an array of index sets that contain the global row numbers
901: that comprise each local submatrix
902: . col - an array of index sets that contain the global column numbers
903: that comprise each local submatrix
904: . submat - array of local submatrices
905: - ctx - optional user-defined context for private data for the
906: user-defined routine (may be null)
908: Output Parameter:
909: . submat - array of local submatrices (the entries of which may
910: have been modified)
912: Notes:
913: The user should NOT generally call this routine, as it will
914: automatically be called within certain preconditioners (currently
915: block Jacobi, additive Schwarz) if set.
917: The basic submatrices are extracted from the preconditioner matrix
918: as usual; the user can then alter these (for example, to set different
919: boundary conditions for each submatrix) before they are used for the
920: local solves.
922: Level: developer
924: .keywords: PC, modify, submatrices
926: .seealso: PCSetModifySubMatrices()
927: @*/
928: PetscErrorCode PCModifySubMatrices(PC pc,PetscInt nsub,const IS row[],const IS col[],Mat submat[],void *ctx)
929: {
933: if (!pc->modifysubmatrices) return(0);
935: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
937: return(0);
938: }
942: /*@
943: PCSetOperators - Sets the matrix associated with the linear system and
944: a (possibly) different one associated with the preconditioner.
946: Collective on PC and Mat
948: Input Parameters:
949: + pc - the preconditioner context
950: . Amat - the matrix associated with the linear system
951: . Pmat - the matrix to be used in constructing the preconditioner, usually
952: the same as Amat.
953: - flag - flag indicating information about the preconditioner matrix structure
954: during successive linear solves. This flag is ignored the first time a
955: linear system is solved, and thus is irrelevant when solving just one linear
956: system.
958: Notes:
959: The flag can be used to eliminate unnecessary work in the preconditioner
960: during the repeated solution of linear systems of the same size. The
961: available options are
962: + SAME_PRECONDITIONER -
963: Pmat is identical during successive linear solves.
964: This option is intended for folks who are using
965: different Amat and Pmat matrices and wish to reuse the
966: same preconditioner matrix. For example, this option
967: saves work by not recomputing incomplete factorization
968: for ILU/ICC preconditioners.
969: . SAME_NONZERO_PATTERN -
970: Pmat has the same nonzero structure during
971: successive linear solves.
972: - DIFFERENT_NONZERO_PATTERN -
973: Pmat does not have the same nonzero structure.
975: Caution:
976: If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion
977: and does not check the structure of the matrix. If you erroneously
978: claim that the structure is the same when it actually is not, the new
979: preconditioner will not function correctly. Thus, use this optimization
980: feature carefully!
982: If in doubt about whether your preconditioner matrix has changed
983: structure or not, use the flag DIFFERENT_NONZERO_PATTERN.
985: More Notes about Repeated Solution of Linear Systems:
986: PETSc does NOT reset the matrix entries of either Amat or Pmat
987: to zero after a linear solve; the user is completely responsible for
988: matrix assembly. See the routine MatZeroEntries() if desiring to
989: zero all elements of a matrix.
991: Level: intermediate
993: .keywords: PC, set, operators, matrix, linear system
995: .seealso: PCGetOperators(), MatZeroEntries()
996: @*/
997: PetscErrorCode PCSetOperators(PC pc,Mat Amat,Mat Pmat,MatStructure flag)
998: {
1000: PetscTruth isbjacobi,isrowbs;
1009: /*
1010: BlockSolve95 cannot use default BJacobi preconditioning
1011: */
1012: if (Amat) {
1013: PetscTypeCompare((PetscObject)Amat,MATMPIROWBS,&isrowbs);
1014: if (isrowbs) {
1015: PetscTypeCompare((PetscObject)pc,PCBJACOBI,&isbjacobi);
1016: if (isbjacobi) {
1017: PCSetType(pc,PCILU);
1018: PetscInfo(pc,"Switching default PC to PCILU since BS95 doesn't support PCBJACOBI\n");
1019: }
1020: }
1021: }
1023: /* reference first in case the matrices are the same */
1024: if (Amat) {PetscObjectReference((PetscObject)Amat);}
1025: if (pc->mat) {MatDestroy(pc->mat);}
1026: if (Pmat) {PetscObjectReference((PetscObject)Pmat);}
1027: if (pc->pmat) {MatDestroy(pc->pmat);}
1028: pc->mat = Amat;
1029: pc->pmat = Pmat;
1031: if (pc->setupcalled == 2 && flag != SAME_PRECONDITIONER) {
1032: pc->setupcalled = 1;
1033: }
1034: pc->flag = flag;
1035: return(0);
1036: }
1040: /*@C
1041: PCGetOperators - Gets the matrix associated with the linear system and
1042: possibly a different one associated with the preconditioner.
1044: Not collective, though parallel Mats are returned if the PC is parallel
1046: Input Parameter:
1047: . pc - the preconditioner context
1049: Output Parameters:
1050: + mat - the matrix associated with the linear system
1051: . pmat - matrix associated with the preconditioner, usually the same
1052: as mat.
1053: - flag - flag indicating information about the preconditioner
1054: matrix structure. See PCSetOperators() for details.
1056: Level: intermediate
1058: Alternative usage: If the operators have NOT been set with KSP/PCSetOperators() then the operators
1059: are created in PC and returned to the user. In this case, if both operators
1060: mat and pmat are requested, two DIFFERENT operators will be returned. If
1061: only one is requested both operators in the PC will be the same (i.e. as
1062: if one had called KSP/PCSetOperators() with the same argument for both Mats).
1063: The user must set the sizes of the returned matrices and their type etc just
1064: as if the user created them with MatCreate(). For example,
1066: $ KSP/PCGetOperators(ksp/pc,&mat,PETSC_NULL,PETSC_NULL); is equivalent to
1067: $ set size, type, etc of mat
1069: $ MatCreate(comm,&mat);
1070: $ KSP/PCSetOperators(ksp/pc,mat,mat,SAME_NONZERO_PATTERN);
1071: $ PetscObjectDereference((PetscObject)mat);
1072: $ set size, type, etc of mat
1074: and
1076: $ KSP/PCGetOperators(ksp/pc,&mat,&pmat,PETSC_NULL); is equivalent to
1077: $ set size, type, etc of mat and pmat
1079: $ MatCreate(comm,&mat);
1080: $ MatCreate(comm,&pmat);
1081: $ KSP/PCSetOperators(ksp/pc,mat,pmat,SAME_NONZERO_PATTERN);
1082: $ PetscObjectDereference((PetscObject)mat);
1083: $ PetscObjectDereference((PetscObject)pmat);
1084: $ set size, type, etc of mat and pmat
1086: The rational for this support is so that when creating a TS, SNES, or KSP the hierarchy
1087: of underlying objects (i.e. SNES, KSP, PC, Mat) and their livespans can be completely
1088: managed by the top most level object (i.e. the TS, SNES, or KSP). Another way to look
1089: at this is when you create a SNES you do not NEED to create a KSP and attach it to
1090: the SNES object (the SNES object manages it for you). Similarly when you create a KSP
1091: you do not need to attach a PC to it (the KSP object manages the PC object for you).
1092: Thus, why should YOU have to create the Mat and attach it to the SNES/KSP/PC, when
1093: it can be created for you?
1094:
1096: .keywords: PC, get, operators, matrix, linear system
1098: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperatorsSet()
1099: @*/
1100: PetscErrorCode PCGetOperators(PC pc,Mat *mat,Mat *pmat,MatStructure *flag)
1101: {
1106: if (mat) {
1107: if (!pc->mat) {
1108: MatCreate(pc->comm,&pc->mat);
1109: if (!pc->pmat && !pmat) { /* user did NOT request pmat, so make same as mat */
1110: pc->pmat = pc->mat;
1111: PetscObjectReference((PetscObject)pc->pmat);
1112: }
1113: }
1114: *mat = pc->mat;
1115: }
1116: if (pmat) {
1117: if (!pc->pmat) {
1118: MatCreate(pc->comm,&pc->mat);
1119: if (!pc->mat && !mat) { /* user did NOT request mat, so make same as pmat */
1120: pc->mat = pc->pmat;
1121: PetscObjectReference((PetscObject)pc->mat);
1122: }
1123: }
1124: *pmat = pc->pmat;
1125: }
1126: if (flag) *flag = pc->flag;
1127: return(0);
1128: }
1132: /*@C
1133: PCGetOperatorsSet - Determines if the matrix associated with the linear system and
1134: possibly a different one associated with the preconditioner have been set in the PC.
1136: Not collective, though the results on all processes should be the same
1138: Input Parameter:
1139: . pc - the preconditioner context
1141: Output Parameters:
1142: + mat - the matrix associated with the linear system was set
1143: - pmat - matrix associated with the preconditioner was set, usually the same
1145: Level: intermediate
1147: .keywords: PC, get, operators, matrix, linear system
1149: .seealso: PCSetOperators(), KSPGetOperators(), KSPSetOperators(), PCGetOperators()
1150: @*/
1151: PetscErrorCode PCGetOperatorsSet(PC pc,PetscTruth *mat,PetscTruth *pmat)
1152: {
1155: if (mat) *mat = (pc->mat) ? PETSC_TRUE : PETSC_FALSE;
1156: if (pmat) *pmat = (pc->pmat) ? PETSC_TRUE : PETSC_FALSE;
1157: return(0);
1158: }
1162: /*@
1163: PCGetFactoredMatrix - Gets the factored matrix from the
1164: preconditioner context. This routine is valid only for the LU,
1165: incomplete LU, Cholesky, and incomplete Cholesky methods.
1167: Not Collective on PC though Mat is parallel if PC is parallel
1169: Input Parameters:
1170: . pc - the preconditioner context
1172: Output parameters:
1173: . mat - the factored matrix
1175: Level: advanced
1177: Notes: Does not increase the reference count for the matrix so DO NOT destroy it
1179: .keywords: PC, get, factored, matrix
1180: @*/
1181: PetscErrorCode PCGetFactoredMatrix(PC pc,Mat *mat)
1182: {
1188: if (pc->ops->getfactoredmatrix) {
1189: (*pc->ops->getfactoredmatrix)(pc,mat);
1190: }
1191: return(0);
1192: }
1196: /*@C
1197: PCSetOptionsPrefix - Sets the prefix used for searching for all
1198: PC options in the database.
1200: Collective on PC
1202: Input Parameters:
1203: + pc - the preconditioner context
1204: - prefix - the prefix string to prepend to all PC option requests
1206: Notes:
1207: A hyphen (-) must NOT be given at the beginning of the prefix name.
1208: The first character of all runtime options is AUTOMATICALLY the
1209: hyphen.
1211: Level: advanced
1213: .keywords: PC, set, options, prefix, database
1215: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1216: @*/
1217: PetscErrorCode PCSetOptionsPrefix(PC pc,const char prefix[])
1218: {
1223: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1224: return(0);
1225: }
1229: /*@C
1230: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1231: PC options in the database.
1233: Collective on PC
1235: Input Parameters:
1236: + pc - the preconditioner context
1237: - prefix - the prefix string to prepend to all PC option requests
1239: Notes:
1240: A hyphen (-) must NOT be given at the beginning of the prefix name.
1241: The first character of all runtime options is AUTOMATICALLY the
1242: hyphen.
1244: Level: advanced
1246: .keywords: PC, append, options, prefix, database
1248: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1249: @*/
1250: PetscErrorCode PCAppendOptionsPrefix(PC pc,const char prefix[])
1251: {
1256: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1257: return(0);
1258: }
1262: /*@C
1263: PCGetOptionsPrefix - Gets the prefix used for searching for all
1264: PC options in the database.
1266: Not Collective
1268: Input Parameters:
1269: . pc - the preconditioner context
1271: Output Parameters:
1272: . prefix - pointer to the prefix string used, is returned
1274: Notes: On the fortran side, the user should pass in a string 'prifix' of
1275: sufficient length to hold the prefix.
1277: Level: advanced
1279: .keywords: PC, get, options, prefix, database
1281: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1282: @*/
1283: PetscErrorCode PCGetOptionsPrefix(PC pc,const char *prefix[])
1284: {
1290: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1291: return(0);
1292: }
1296: /*@
1297: PCPreSolve - Optional pre-solve phase, intended for any
1298: preconditioner-specific actions that must be performed before
1299: the iterative solve itself.
1301: Collective on PC
1303: Input Parameters:
1304: + pc - the preconditioner context
1305: - ksp - the Krylov subspace context
1307: Level: developer
1309: Sample of Usage:
1310: .vb
1311: PCPreSolve(pc,ksp);
1312: KSPSolve(ksp,b,x);
1313: PCPostSolve(pc,ksp);
1314: .ve
1316: Notes:
1317: The pre-solve phase is distinct from the PCSetUp() phase.
1319: KSPSolve() calls this directly, so is rarely called by the user.
1321: .keywords: PC, pre-solve
1323: .seealso: PCPostSolve()
1324: @*/
1325: PetscErrorCode PCPreSolve(PC pc,KSP ksp)
1326: {
1328: Vec x,rhs;
1329: Mat A,B;
1334: KSPGetSolution(ksp,&x);
1335: KSPGetRhs(ksp,&rhs);
1336: /*
1337: Scale the system and have the matrices use the scaled form
1338: only if the two matrices are actually the same (and hence
1339: have the same scaling
1340: */
1341: PCGetOperators(pc,&A,&B,PETSC_NULL);
1342: if (A == B) {
1343: MatScaleSystem(pc->mat,rhs,x);
1344: MatUseScaledForm(pc->mat,PETSC_TRUE);
1345: }
1347: if (pc->ops->presolve) {
1348: (*pc->ops->presolve)(pc,ksp,rhs,x);
1349: }
1350: return(0);
1351: }
1355: /*@
1356: PCPostSolve - Optional post-solve phase, intended for any
1357: preconditioner-specific actions that must be performed after
1358: the iterative solve itself.
1360: Collective on PC
1362: Input Parameters:
1363: + pc - the preconditioner context
1364: - ksp - the Krylov subspace context
1366: Sample of Usage:
1367: .vb
1368: PCPreSolve(pc,ksp);
1369: KSPSolve(ksp,b,x);
1370: PCPostSolve(pc,ksp);
1371: .ve
1373: Note:
1374: KSPSolve() calls this routine directly, so it is rarely called by the user.
1376: Level: developer
1378: .keywords: PC, post-solve
1380: .seealso: PCPreSolve(), KSPSolve()
1381: @*/
1382: PetscErrorCode PCPostSolve(PC pc,KSP ksp)
1383: {
1385: Vec x,rhs;
1386: Mat A,B;
1391: KSPGetSolution(ksp,&x);
1392: KSPGetRhs(ksp,&rhs);
1393: if (pc->ops->postsolve) {
1394: (*pc->ops->postsolve)(pc,ksp,rhs,x);
1395: }
1397: /*
1398: Scale the system and have the matrices use the scaled form
1399: only if the two matrices are actually the same (and hence
1400: have the same scaling
1401: */
1402: PCGetOperators(pc,&A,&B,PETSC_NULL);
1403: if (A == B) {
1404: MatUnScaleSystem(pc->mat,rhs,x);
1405: MatUseScaledForm(pc->mat,PETSC_FALSE);
1406: }
1407: return(0);
1408: }
1412: /*@C
1413: PCView - Prints the PC data structure.
1415: Collective on PC
1417: Input Parameters:
1418: + PC - the PC context
1419: - viewer - optional visualization context
1421: Note:
1422: The available visualization contexts include
1423: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1424: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1425: output where only the first processor opens
1426: the file. All other processors send their
1427: data to the first processor to print.
1429: The user can open an alternative visualization contexts with
1430: PetscViewerASCIIOpen() (output to a specified file).
1432: Level: developer
1434: .keywords: PC, view
1436: .seealso: KSPView(), PetscViewerASCIIOpen()
1437: @*/
1438: PetscErrorCode PCView(PC pc,PetscViewer viewer)
1439: {
1440: PCType cstr;
1441: PetscErrorCode ierr;
1442: PetscTruth mat_exists,iascii,isstring;
1443: PetscViewerFormat format;
1447: if (!viewer) {
1448: PetscViewerASCIIGetStdout(pc->comm,&viewer);
1449: }
1453: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
1454: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
1455: if (iascii) {
1456: PetscViewerGetFormat(viewer,&format);
1457: if (pc->prefix) {
1458: PetscViewerASCIIPrintf(viewer,"PC Object:(%s)\n",pc->prefix);
1459: } else {
1460: PetscViewerASCIIPrintf(viewer,"PC Object:\n");
1461: }
1462: PCGetType(pc,&cstr);
1463: if (cstr) {
1464: PetscViewerASCIIPrintf(viewer," type: %s\n",cstr);
1465: } else {
1466: PetscViewerASCIIPrintf(viewer," type: not yet set\n");
1467: }
1468: if (pc->ops->view) {
1469: PetscViewerASCIIPushTab(viewer);
1470: (*pc->ops->view)(pc,viewer);
1471: PetscViewerASCIIPopTab(viewer);
1472: }
1473: PetscObjectExists((PetscObject)pc->mat,&mat_exists);
1474: if (mat_exists) {
1475: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1476: if (pc->pmat == pc->mat) {
1477: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:\n");
1478: PetscViewerASCIIPushTab(viewer);
1479: MatView(pc->mat,viewer);
1480: PetscViewerASCIIPopTab(viewer);
1481: } else {
1482: PetscObjectExists((PetscObject)pc->pmat,&mat_exists);
1483: if (mat_exists) {
1484: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:\n");
1485: } else {
1486: PetscViewerASCIIPrintf(viewer," linear system matrix:\n");
1487: }
1488: PetscViewerASCIIPushTab(viewer);
1489: MatView(pc->mat,viewer);
1490: if (mat_exists) {MatView(pc->pmat,viewer);}
1491: PetscViewerASCIIPopTab(viewer);
1492: }
1493: PetscViewerPopFormat(viewer);
1494: }
1495: } else if (isstring) {
1496: PCGetType(pc,&cstr);
1497: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1498: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1499: } else {
1500: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by PC",((PetscObject)viewer)->type_name);
1501: }
1502: return(0);
1503: }
1507: /*@C
1508: PCRegister - See PCRegisterDynamic()
1510: Level: advanced
1511: @*/
1512: PetscErrorCode PCRegister(const char sname[],const char path[],const char name[],PetscErrorCode (*function)(PC))
1513: {
1515: char fullname[PETSC_MAX_PATH_LEN];
1519: PetscFListConcat(path,name,fullname);
1520: PetscFListAdd(&PCList,sname,fullname,(void (*)(void))function);
1521: return(0);
1522: }
1526: /*@
1527: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1529: Collective on PC
1531: Input Parameter:
1532: . pc - the preconditioner object
1534: Output Parameter:
1535: . mat - the explict preconditioned operator
1537: Notes:
1538: This computation is done by applying the operators to columns of the
1539: identity matrix.
1541: Currently, this routine uses a dense matrix format when 1 processor
1542: is used and a sparse format otherwise. This routine is costly in general,
1543: and is recommended for use only with relatively small systems.
1545: Level: advanced
1546:
1547: .keywords: PC, compute, explicit, operator
1549: @*/
1550: PetscErrorCode PCComputeExplicitOperator(PC pc,Mat *mat)
1551: {
1552: Vec in,out;
1554: PetscInt i,M,m,*rows,start,end;
1555: PetscMPIInt size;
1556: MPI_Comm comm;
1557: PetscScalar *array,one = 1.0;
1558:
1563: comm = pc->comm;
1564: MPI_Comm_size(comm,&size);
1566: if (!pc->pmat) SETERRQ(PETSC_ERR_ORDER,"You must call KSPSetOperators() or PCSetOperators() before this call");
1567: MatGetVecs(pc->pmat,&in,0);
1568: VecDuplicate(in,&out);
1569: VecGetOwnershipRange(in,&start,&end);
1570: VecGetSize(in,&M);
1571: VecGetLocalSize(in,&m);
1572: PetscMalloc((m+1)*sizeof(PetscInt),&rows);
1573: for (i=0; i<m; i++) {rows[i] = start + i;}
1575: MatCreate(comm,mat);
1576: MatSetSizes(*mat,m,m,M,M);
1577: if (size == 1) {
1578: MatSetType(*mat,MATSEQDENSE);
1579: MatSeqDenseSetPreallocation(*mat,PETSC_NULL);
1580: } else {
1581: MatSetType(*mat,MATMPIAIJ);
1582: MatMPIAIJSetPreallocation(*mat,0,PETSC_NULL,0,PETSC_NULL);
1583: }
1585: for (i=0; i<M; i++) {
1587: VecSet(in,0.0);
1588: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1589: VecAssemblyBegin(in);
1590: VecAssemblyEnd(in);
1592: /* should fix, allowing user to choose side */
1593: PCApply(pc,in,out);
1594:
1595: VecGetArray(out,&array);
1596: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1597: VecRestoreArray(out,&array);
1599: }
1600: PetscFree(rows);
1601: VecDestroy(out);
1602: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1603: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1604: return(0);
1605: }