Actual source code: precon.c
1: /*$Id: precon.c,v 1.216 2001/08/21 21:03:13 bsmith Exp $*/
2: /*
3: The PC (preconditioner) interface routines, callable by users.
4: */
5: #include src/sles/pc/pcimpl.h
7: /* Logging support */
8: int PC_COOKIE;
9: int PC_SetUp, PC_SetUpOnBlocks, PC_Apply, PC_ApplyCoarse, PC_ApplyMultiple, PC_ApplySymmetricLeft;
10: int PC_ApplySymmetricRight, PC_ModifySubMatrices;
12: EXTERN int SLESInitializePackage(char *);
14: #undef __FUNCT__
16: /*@C
17: PCNullSpaceAttach - attaches a null space to a preconditioner object.
18: This null space will be removed from the resulting vector whenever
19: the preconditioner is applied.
21: Collective on PC
23: Input Parameters:
24: + pc - the preconditioner context
25: - nullsp - the null space object
27: Level: developer
29: Notes:
30: Overwrites any previous null space that may have been attached
32: Users manual sections:
33: . Section 4.15 Solving Singular Systems
35: .keywords: PC, destroy, null space
37: .seealso: PCCreate(), PCSetUp(), MatNullSpaceCreate(), MatNullSpace
39: @*/
40: int PCNullSpaceAttach(PC pc,MatNullSpace nullsp)
41: {
48: if (pc->nullsp) {
49: MatNullSpaceDestroy(pc->nullsp);
50: }
51: pc->nullsp = nullsp;
52: PetscObjectReference((PetscObject)nullsp);
53: return(0);
54: }
56: #undef __FUNCT__
58: /*@C
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: int 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->nullsp) {MatNullSpaceDestroy(pc->nullsp);}
85: if (pc->diagonalscaleright) {VecDestroy(pc->diagonalscaleright);}
86: if (pc->diagonalscaleleft) {VecDestroy(pc->diagonalscaleleft);}
88: PetscLogObjectDestroy(pc);
89: PetscHeaderDestroy(pc);
90: return(0);
91: }
93: #undef __FUNCT__
95: /*@C
96: PCDiagonalScale - Indicates if the preconditioner applies an additional left and right
97: scaling as needed by certain time-stepping codes.
99: Collective on PC
101: Input Parameter:
102: . pc - the preconditioner context
104: Output Parameter:
105: . flag - PETSC_TRUE if it applies the scaling
107: Level: developer
109: Notes: If this returns PETSC_TRUE then the system solved via the Krylov method is
110: $ D M A D^{-1} y = D M b for left preconditioning or
111: $ D A M D^{-1} z = D b for right preconditioning
113: .keywords: PC
115: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCDiagonalScaleSet()
116: @*/
117: int PCDiagonalScale(PC pc,PetscTruth *flag)
118: {
121: *flag = pc->diagonalscale;
122: return(0);
123: }
125: #undef __FUNCT__
127: /*@
128: PCDiagonalScaleSet - Indicates the left scaling to use to apply an additional left and right
129: scaling as needed by certain time-stepping codes.
131: Collective on PC
133: Input Parameters:
134: + pc - the preconditioner context
135: - s - scaling vector
137: Level: intermediate
139: Notes: The system solved via the Krylov method is
140: $ D M A D^{-1} y = D M b for left preconditioning or
141: $ D A M D^{-1} z = D b for right preconditioning
143: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
145: .keywords: PC
147: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleRight(), PCDiagonalScale()
148: @*/
149: int PCDiagonalScaleSet(PC pc,Vec s)
150: {
156: pc->diagonalscale = PETSC_TRUE;
157: if (pc->diagonalscaleleft) {
158: VecDestroy(pc->diagonalscaleleft);
159: }
160: pc->diagonalscaleleft = s;
161: ierr = PetscObjectReference((PetscObject)s);
162: if (!pc->diagonalscaleright) {
163: VecDuplicate(s,&pc->diagonalscaleright);
164: }
165: VecCopy(s,pc->diagonalscaleright);
166: VecReciprocal(pc->diagonalscaleright);
167: return(0);
168: }
170: #undef __FUNCT__
172: /*@C
173: PCDiagonalScaleLeft - Indicates the left scaling to use to apply an additional left and right
174: scaling as needed by certain time-stepping codes.
176: Collective on PC
178: Input Parameters:
179: + pc - the preconditioner context
180: . in - input vector
181: + out - scaled vector (maybe the same as in)
183: Level: intermediate
185: Notes: The system solved via the Krylov method is
186: $ D M A D^{-1} y = D M b for left preconditioning or
187: $ D A M D^{-1} z = D b for right preconditioning
189: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
191: If diagonal scaling is turned off and in is not out then in is copied to out
193: .keywords: PC
195: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleSet(), PCDiagonalScaleRight(), PCDiagonalScale()
196: @*/
197: int PCDiagonalScaleLeft(PC pc,Vec in,Vec out)
198: {
205: if (pc->diagonalscale) {
206: VecPointwiseMult(pc->diagonalscaleleft,in,out);
207: } else if (in != out) {
208: VecCopy(in,out);
209: }
210: return(0);
211: }
213: #undef __FUNCT__
215: /*@C
216: PCDiagonalScaleRight - Scales a vector by the right scaling as needed by certain time-stepping codes.
218: Collective on PC
220: Input Parameters:
221: + pc - the preconditioner context
222: . in - input vector
223: + out - scaled vector (maybe the same as in)
225: Level: intermediate
227: Notes: The system solved via the Krylov method is
228: $ D M A D^{-1} y = D M b for left preconditioning or
229: $ D A M D^{-1} z = D b for right preconditioning
231: PCDiagonalScaleLeft() scales a vector by D. PCDiagonalScaleRight() scales a vector by D^{-1}.
233: If diagonal scaling is turned off and in is not out then in is copied to out
235: .keywords: PC
237: .seealso: PCCreate(), PCSetUp(), PCDiagonalScaleLeft(), PCDiagonalScaleSet(), PCDiagonalScale()
238: @*/
239: int PCDiagonalScaleRight(PC pc,Vec in,Vec out)
240: {
247: if (pc->diagonalscale) {
248: VecPointwiseMult(pc->diagonalscaleright,in,out);
249: } else if (in != out) {
250: VecCopy(in,out);
251: }
252: return(0);
253: }
255: #undef __FUNCT__
257: static int PCPublish_Petsc(PetscObject obj)
258: {
259: #if defined(PETSC_HAVE_AMS)
260: PC v = (PC) obj;
261: int ierr;
262: #endif
266: #if defined(PETSC_HAVE_AMS)
267: /* if it is already published then return */
268: if (v->amem >=0) return(0);
270: PetscObjectPublishBaseBegin(obj);
271: PetscObjectPublishBaseEnd(obj);
272: #endif
274: return(0);
275: }
277: #undef __FUNCT__
279: /*@C
280: PCCreate - Creates a preconditioner context.
282: Collective on MPI_Comm
284: Input Parameter:
285: . comm - MPI communicator
287: Output Parameter:
288: . pc - location to put the preconditioner context
290: Notes:
291: The default preconditioner on one processor is PCILU with 0 fill on more
292: then one it is PCBJACOBI with ILU() on each processor.
294: Level: developer
296: .keywords: PC, create, context
298: .seealso: PCSetUp(), PCApply(), PCDestroy()
299: @*/
300: int PCCreate(MPI_Comm comm,PC *newpc)
301: {
302: PC pc;
307: *newpc = 0;
308: #ifndef PETSC_USE_DYNAMIC_LIBRARIES
309: SLESInitializePackage(PETSC_NULL);
310: #endif
312: PetscHeaderCreate(pc,_p_PC,struct _PCOps,PC_COOKIE,-1,"PC",comm,PCDestroy,PCView);
313: PetscLogObjectCreate(pc);
314: pc->bops->publish = PCPublish_Petsc;
315: pc->vec = 0;
316: pc->mat = 0;
317: pc->setupcalled = 0;
318: pc->nullsp = 0;
319: pc->data = 0;
320: pc->diagonalscale = PETSC_FALSE;
321: pc->diagonalscaleleft = 0;
322: pc->diagonalscaleright = 0;
324: pc->ops->destroy = 0;
325: pc->ops->apply = 0;
326: pc->ops->applytranspose = 0;
327: pc->ops->applyBA = 0;
328: pc->ops->applyBAtranspose = 0;
329: pc->ops->applyrichardson = 0;
330: pc->ops->view = 0;
331: pc->ops->getfactoredmatrix = 0;
332: pc->ops->applysymmetricright = 0;
333: pc->ops->applysymmetricleft = 0;
334: pc->ops->setuponblocks = 0;
336: pc->modifysubmatrices = 0;
337: pc->modifysubmatricesP = 0;
338: *newpc = pc;
339: PetscPublishAll(pc);
340: return(0);
342: }
344: /* -------------------------------------------------------------------------------*/
346: #undef __FUNCT__
348: /*@
349: PCApply - Applies the preconditioner to a vector.
351: Collective on PC and Vec
353: Input Parameters:
354: + pc - the preconditioner context
355: - x - input vector
357: Output Parameter:
358: . y - output vector
360: Level: developer
362: .keywords: PC, apply
364: .seealso: PCApplyTranspose(), PCApplyBAorAB()
365: @*/
366: int PCApply(PC pc,Vec x,Vec y)
367: {
368: int ierr;
374: if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
376: if (pc->setupcalled < 2) {
377: PCSetUp(pc);
378: }
380: PetscLogEventBegin(PC_Apply,pc,x,y,0);
381: (*pc->ops->apply)(pc,x,y);
382: PetscLogEventEnd(PC_Apply,pc,x,y,0);
384: /* Remove null space from preconditioned vector y */
385: if (pc->nullsp) {
386: MatNullSpaceRemove(pc->nullsp,y,PETSC_NULL);
387: }
388: return(0);
389: }
391: #undef __FUNCT__
393: /*@
394: PCApplySymmetricLeft - Applies the left part of a symmetric preconditioner to a vector.
396: Collective on PC and Vec
398: Input Parameters:
399: + pc - the preconditioner context
400: - x - input vector
402: Output Parameter:
403: . y - output vector
405: Notes:
406: Currently, this routine is implemented only for PCICC and PCJACOBI preconditioners.
408: Level: developer
410: .keywords: PC, apply, symmetric, left
412: .seealso: PCApply(), PCApplySymmetricRight()
413: @*/
414: int PCApplySymmetricLeft(PC pc,Vec x,Vec y)
415: {
422: if (!pc->ops->applysymmetricleft) SETERRQ(1,"PC does not have left symmetric apply");
424: if (pc->setupcalled < 2) {
425: PCSetUp(pc);
426: }
428: PetscLogEventBegin(PC_ApplySymmetricLeft,pc,x,y,0);
429: (*pc->ops->applysymmetricleft)(pc,x,y);
430: PetscLogEventEnd(PC_ApplySymmetricLeft,pc,x,y,0);
431: return(0);
432: }
434: #undef __FUNCT__
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: int PCApplySymmetricRight(PC pc,Vec x,Vec y)
458: {
465: if (!pc->ops->applysymmetricright) SETERRQ(1,"PC does not have left symmetric apply");
467: if (pc->setupcalled < 2) {
468: PCSetUp(pc);
469: }
471: PetscLogEventBegin(PC_ApplySymmetricRight,pc,x,y,0);
472: (*pc->ops->applysymmetricright)(pc,x,y);
473: PetscLogEventEnd(PC_ApplySymmetricRight,pc,x,y,0);
474: return(0);
475: }
477: #undef __FUNCT__
479: /*@
480: PCApplyTranspose - Applies the transpose of preconditioner to a vector.
482: Collective on PC and Vec
484: Input Parameters:
485: + pc - the preconditioner context
486: - x - input vector
488: Output Parameter:
489: . y - output vector
491: Level: developer
493: .keywords: PC, apply, transpose
495: .seealso: PCApplyTranspose(), PCApplyBAorAB(), PCApplyBAorABTranspose()
496: @*/
497: int PCApplyTranspose(PC pc,Vec x,Vec y)
498: {
505: if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
506: if (!pc->ops->applytranspose) SETERRQ(PETSC_ERR_SUP," ");
508: if (pc->setupcalled < 2) {
509: PCSetUp(pc);
510: }
512: PetscLogEventBegin(PC_Apply,pc,x,y,0);
513: (*pc->ops->applytranspose)(pc,x,y);
514: PetscLogEventEnd(PC_Apply,pc,x,y,0);
515: return(0);
516: }
518: #undef __FUNCT__
520: /*@
521: PCApplyBAorAB - Applies the preconditioner and operator to a vector.
523: Collective on PC and Vec
525: Input Parameters:
526: + pc - the preconditioner context
527: . side - indicates the preconditioner side, one of PC_LEFT, PC_RIGHT, or PC_SYMMETRIC
528: . x - input vector
529: - work - work vector
531: Output Parameter:
532: . y - output vector
534: Level: developer
536: .keywords: PC, apply, operator
538: .seealso: PCApply(), PCApplyTranspose(), PCApplyBAorABTranspose()
539: @*/
540: int PCApplyBAorAB(PC pc,PCSide side,Vec x,Vec y,Vec work)
541: {
542: int ierr;
549: if (x == y) SETERRQ(PETSC_ERR_ARG_IDN,"x and y must be different vectors");
550: if (side != PC_LEFT && side != PC_SYMMETRIC && side != PC_RIGHT) {
551: SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Side must be right, left, or symmetric");
552: }
553: if (pc->diagonalscale && side == PC_SYMMETRIC) {
554: SETERRQ(1,"Cannot include diagonal scaling with symmetric preconditioner application");
555: }
557: if (pc->setupcalled < 2) {
558: PCSetUp(pc);
559: }
561: if (pc->diagonalscale) {
562: if (pc->ops->applyBA) {
563: Vec work2; /* this is expensive, but to fix requires a second work vector argument to PCApplyBAorAB() */
564: VecDuplicate(x,&work2);
565: PCDiagonalScaleRight(pc,x,work2);
566: (*pc->ops->applyBA)(pc,side,work2,y,work);
567: /* Remove null space from preconditioned vector y */
568: if (pc->nullsp) {
569: MatNullSpaceRemove(pc->nullsp,y,PETSC_NULL);
570: }
571: PCDiagonalScaleLeft(pc,y,y);
572: VecDestroy(work2);
573: } else if (side == PC_RIGHT) {
574: PCDiagonalScaleRight(pc,x,y);
575: PCApply(pc,y,work);
576: MatMult(pc->mat,work,y);
577: PCDiagonalScaleLeft(pc,y,y);
578: } else if (side == PC_LEFT) {
579: PCDiagonalScaleRight(pc,x,y);
580: MatMult(pc->mat,y,work);
581: PCApply(pc,work,y);
582: PCDiagonalScaleLeft(pc,y,y);
583: } else if (side == PC_SYMMETRIC) {
584: SETERRQ(1,"Cannot provide diagonal scaling with symmetric application of preconditioner");
585: }
586: } else {
587: if (pc->ops->applyBA) {
588: (*pc->ops->applyBA)(pc,side,x,y,work);
589: /* Remove null space from preconditioned vector y */
590: if (pc->nullsp) {
591: MatNullSpaceRemove(pc->nullsp,y,PETSC_NULL);
592: }
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: }
610: #undef __FUNCT__
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: int 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: /* -------------------------------------------------------------------------------*/
669: #undef __FUNCT__
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: int PCApplyRichardsonExists(PC pc,PetscTruth *exists)
690: {
694: if (pc->ops->applyrichardson) *exists = PETSC_TRUE;
695: else *exists = PETSC_FALSE;
696: return(0);
697: }
699: #undef __FUNCT__
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: . atol - 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: int PCApplyRichardson(PC pc,Vec x,Vec y,Vec w,PetscReal rtol,PetscReal atol, PetscReal dtol,int 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,atol,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: */
757: #undef __FUNCT__
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: int PCSetUp(PC pc)
774: {
776: PetscTruth flg;
781: if (pc->setupcalled > 1) {
782: PetscLogInfo(pc,"PCSetUp:Setting PC with identical preconditionern");
783: return(0);
784: } else if (pc->setupcalled == 0) {
785: PetscLogInfo(pc,"PCSetUp:Setting up new PCn");
786: } else if (pc->flag == SAME_NONZERO_PATTERN) {
787: PetscLogInfo(pc,"PCSetUp:Setting up PC with same nonzero patternn");
788: } else {
789: PetscLogInfo(pc,"PCSetUp:Setting up PC with different nonzero patternn");
790: }
792: PetscLogEventBegin(PC_SetUp,pc,0,0,0);
793: if (!pc->vec) {SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Vector must be set first");}
794: if (!pc->mat) {SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be set first");}
796: if (!pc->type_name) {
797: int size;
799: MPI_Comm_size(pc->comm,&size);
800: if (size == 1) {
801: PetscTypeCompare((PetscObject)pc->pmat,MATSEQSBAIJ,&flg);
802: if (flg) { /* for sbaij mat */
803: PCSetType(pc,PCICC);
804: } else {
805: PCSetType(pc,PCILU);
806: }
807: } else {
808: PCSetType(pc,PCBJACOBI);
809: }
810: }
811: if (pc->ops->setup) {
812: (*pc->ops->setup)(pc);
813: }
814: pc->setupcalled = 2;
815: if (pc->nullsp) {
816: PetscTruth test;
817: PetscOptionsHasName(pc->prefix,"-pc_test_null_space",&test);
818: if (test) {
819: MatNullSpaceTest(pc->nullsp,pc->mat);
820: }
821: }
822: PetscLogEventEnd(PC_SetUp,pc,0,0,0);
823: return(0);
824: }
826: #undef __FUNCT__
828: /*@
829: PCSetUpOnBlocks - Sets up the preconditioner for each block in
830: the block Jacobi, block Gauss-Seidel, and overlapping Schwarz
831: methods.
833: Collective on PC
835: Input Parameters:
836: . pc - the preconditioner context
838: Level: developer
840: .keywords: PC, setup, blocks
842: .seealso: PCCreate(), PCApply(), PCDestroy(), PCSetUp()
843: @*/
844: int PCSetUpOnBlocks(PC pc)
845: {
850: if (!pc->ops->setuponblocks) return(0);
851: PetscLogEventBegin(PC_SetUpOnBlocks,pc,0,0,0);
852: (*pc->ops->setuponblocks)(pc);
853: PetscLogEventEnd(PC_SetUpOnBlocks,pc,0,0,0);
854: return(0);
855: }
857: #undef __FUNCT__
859: /*@
860: PCSetModifySubMatrices - Sets a user-defined routine for modifying the
861: submatrices that arise within certain subdomain-based preconditioners.
862: The basic submatrices are extracted from the preconditioner matrix as
863: usual; the user can then alter these (for example, to set different boundary
864: conditions for each submatrix) before they are used for the local solves.
866: Collective on PC
868: Input Parameters:
869: + pc - the preconditioner context
870: . func - routine for modifying the submatrices
871: - ctx - optional user-defined context (may be null)
873: Calling sequence of func:
874: $ func (PC pc,int nsub,IS *row,IS *col,Mat *submat,void *ctx);
876: . row - an array of index sets that contain the global row numbers
877: that comprise each local submatrix
878: . col - an array of index sets that contain the global column numbers
879: that comprise each local submatrix
880: . submat - array of local submatrices
881: - ctx - optional user-defined context for private data for the
882: user-defined func routine (may be null)
884: Notes:
885: PCSetModifySubMatrices() MUST be called before SLESSetUp() and
886: SLESSolve().
888: A routine set by PCSetModifySubMatrices() is currently called within
889: the block Jacobi (PCBJACOBI) and additive Schwarz (PCASM)
890: preconditioners. All other preconditioners ignore this routine.
892: Level: advanced
894: .keywords: PC, set, modify, submatrices
896: .seealso: PCModifySubMatrices()
897: @*/
898: int PCSetModifySubMatrices(PC pc,int(*func)(PC,int,IS*,IS*,Mat*,void*),void *ctx)
899: {
902: pc->modifysubmatrices = func;
903: pc->modifysubmatricesP = ctx;
904: return(0);
905: }
907: #undef __FUNCT__
909: /*@
910: PCModifySubMatrices - Calls an optional user-defined routine within
911: certain preconditioners if one has been set with PCSetModifySubMarices().
913: Collective on PC
915: Input Parameters:
916: + pc - the preconditioner context
917: . nsub - the number of local submatrices
918: . row - an array of index sets that contain the global row numbers
919: that comprise each local submatrix
920: . col - an array of index sets that contain the global column numbers
921: that comprise each local submatrix
922: . submat - array of local submatrices
923: - ctx - optional user-defined context for private data for the
924: user-defined routine (may be null)
926: Output Parameter:
927: . submat - array of local submatrices (the entries of which may
928: have been modified)
930: Notes:
931: The user should NOT generally call this routine, as it will
932: automatically be called within certain preconditioners (currently
933: block Jacobi, additive Schwarz) if set.
935: The basic submatrices are extracted from the preconditioner matrix
936: as usual; the user can then alter these (for example, to set different
937: boundary conditions for each submatrix) before they are used for the
938: local solves.
940: Level: developer
942: .keywords: PC, modify, submatrices
944: .seealso: PCSetModifySubMatrices()
945: @*/
946: int PCModifySubMatrices(PC pc,int nsub,IS *row,IS *col,Mat *submat,void *ctx)
947: {
951: if (!pc->modifysubmatrices) return(0);
952: PetscLogEventBegin(PC_ModifySubMatrices,pc,0,0,0);
953: (*pc->modifysubmatrices)(pc,nsub,row,col,submat,ctx);
954: PetscLogEventEnd(PC_ModifySubMatrices,pc,0,0,0);
955: return(0);
956: }
958: #undef __FUNCT__
960: /*@
961: PCSetOperators - Sets the matrix associated with the linear system and
962: a (possibly) different one associated with the preconditioner.
964: Collective on PC and Mat
966: Input Parameters:
967: + pc - the preconditioner context
968: . Amat - the matrix associated with the linear system
969: . Pmat - the matrix to be used in constructing the preconditioner, usually
970: the same as Amat.
971: - flag - flag indicating information about the preconditioner matrix structure
972: during successive linear solves. This flag is ignored the first time a
973: linear system is solved, and thus is irrelevant when solving just one linear
974: system.
976: Notes:
977: The flag can be used to eliminate unnecessary work in the preconditioner
978: during the repeated solution of linear systems of the same size. The
979: available options are
980: + SAME_PRECONDITIONER -
981: Pmat is identical during successive linear solves.
982: This option is intended for folks who are using
983: different Amat and Pmat matrices and wish to reuse the
984: same preconditioner matrix. For example, this option
985: saves work by not recomputing incomplete factorization
986: for ILU/ICC preconditioners.
987: . SAME_NONZERO_PATTERN -
988: Pmat has the same nonzero structure during
989: successive linear solves.
990: - DIFFERENT_NONZERO_PATTERN -
991: Pmat does not have the same nonzero structure.
993: Caution:
994: If you specify SAME_NONZERO_PATTERN, PETSc believes your assertion
995: and does not check the structure of the matrix. If you erroneously
996: claim that the structure is the same when it actually is not, the new
997: preconditioner will not function correctly. Thus, use this optimization
998: feature carefully!
1000: If in doubt about whether your preconditioner matrix has changed
1001: structure or not, use the flag DIFFERENT_NONZERO_PATTERN.
1003: More Notes about Repeated Solution of Linear Systems:
1004: PETSc does NOT reset the matrix entries of either Amat or Pmat
1005: to zero after a linear solve; the user is completely responsible for
1006: matrix assembly. See the routine MatZeroEntries() if desiring to
1007: zero all elements of a matrix.
1009: Level: intermediate
1011: .keywords: PC, set, operators, matrix, linear system
1013: .seealso: PCGetOperators(), MatZeroEntries()
1014: @*/
1015: int PCSetOperators(PC pc,Mat Amat,Mat Pmat,MatStructure flag)
1016: {
1017: int ierr;
1018: PetscTruth isbjacobi,isrowbs;
1025: /*
1026: BlockSolve95 cannot use default BJacobi preconditioning
1027: */
1028: PetscTypeCompare((PetscObject)Amat,MATMPIROWBS,&isrowbs);
1029: if (isrowbs) {
1030: PetscTypeCompare((PetscObject)pc,PCBJACOBI,&isbjacobi);
1031: if (isbjacobi) {
1032: PCSetType(pc,PCILU);
1033: PetscLogInfo(pc,"PCSetOperators:Switching default PC to PCILU since BS95 doesn't support PCBJACOBIn");
1034: }
1035: }
1037: pc->mat = Amat;
1038: pc->pmat = Pmat;
1039: if (pc->setupcalled == 2 && flag != SAME_PRECONDITIONER) {
1040: pc->setupcalled = 1;
1041: }
1042: pc->flag = flag;
1043: return(0);
1044: }
1046: #undef __FUNCT__
1048: /*@C
1049: PCGetOperators - Gets the matrix associated with the linear system and
1050: possibly a different one associated with the preconditioner.
1052: Not collective, though parallel Mats are returned if the PC is parallel
1054: Input Parameter:
1055: . pc - the preconditioner context
1057: Output Parameters:
1058: + mat - the matrix associated with the linear system
1059: . pmat - matrix associated with the preconditioner, usually the same
1060: as mat.
1061: - flag - flag indicating information about the preconditioner
1062: matrix structure. See PCSetOperators() for details.
1064: Level: intermediate
1066: .keywords: PC, get, operators, matrix, linear system
1068: .seealso: PCSetOperators()
1069: @*/
1070: int PCGetOperators(PC pc,Mat *mat,Mat *pmat,MatStructure *flag)
1071: {
1074: if (mat) *mat = pc->mat;
1075: if (pmat) *pmat = pc->pmat;
1076: if (flag) *flag = pc->flag;
1077: return(0);
1078: }
1080: #undef __FUNCT__
1082: /*@
1083: PCSetVector - Sets a vector associated with the preconditioner.
1085: Collective on PC and Vec
1087: Input Parameters:
1088: + pc - the preconditioner context
1089: - vec - the vector
1091: Notes:
1092: The vector must be set so that the preconditioner knows what type
1093: of vector to allocate if necessary.
1095: Level: intermediate
1097: .keywords: PC, set, vector
1099: .seealso: PCGetVector()
1101: @*/
1102: int PCSetVector(PC pc,Vec vec)
1103: {
1108: pc->vec = vec;
1109: return(0);
1110: }
1112: #undef __FUNCT__
1114: /*@
1115: PCGetVector - Gets a vector associated with the preconditioner; if the
1116: vector was not get set it will return a 0 pointer.
1118: Not collective, but vector is shared by all processors that share the PC
1120: Input Parameter:
1121: . pc - the preconditioner context
1123: Output Parameter:
1124: . vec - the vector
1126: Level: intermediate
1128: .keywords: PC, get, vector
1130: .seealso: PCSetVector()
1132: @*/
1133: int PCGetVector(PC pc,Vec *vec)
1134: {
1137: *vec = pc->vec;
1138: return(0);
1139: }
1141: #undef __FUNCT__
1143: /*@C
1144: PCGetFactoredMatrix - Gets the factored matrix from the
1145: preconditioner context. This routine is valid only for the LU,
1146: incomplete LU, Cholesky, and incomplete Cholesky methods.
1148: Not Collective on PC though Mat is parallel if PC is parallel
1150: Input Parameters:
1151: . pc - the preconditioner context
1153: Output parameters:
1154: . mat - the factored matrix
1156: Level: advanced
1158: .keywords: PC, get, factored, matrix
1159: @*/
1160: int PCGetFactoredMatrix(PC pc,Mat *mat)
1161: {
1166: if (pc->ops->getfactoredmatrix) {
1167: (*pc->ops->getfactoredmatrix)(pc,mat);
1168: }
1169: return(0);
1170: }
1172: #undef __FUNCT__
1174: /*@C
1175: PCSetOptionsPrefix - Sets the prefix used for searching for all
1176: PC options in the database.
1178: Collective on PC
1180: Input Parameters:
1181: + pc - the preconditioner context
1182: - prefix - the prefix string to prepend to all PC option requests
1184: Notes:
1185: A hyphen (-) must NOT be given at the beginning of the prefix name.
1186: The first character of all runtime options is AUTOMATICALLY the
1187: hyphen.
1189: Level: advanced
1191: .keywords: PC, set, options, prefix, database
1193: .seealso: PCAppendOptionsPrefix(), PCGetOptionsPrefix()
1194: @*/
1195: int PCSetOptionsPrefix(PC pc,char *prefix)
1196: {
1201: PetscObjectSetOptionsPrefix((PetscObject)pc,prefix);
1202: return(0);
1203: }
1205: #undef __FUNCT__
1207: /*@C
1208: PCAppendOptionsPrefix - Appends to the prefix used for searching for all
1209: PC options in the database.
1211: Collective on PC
1213: Input Parameters:
1214: + pc - the preconditioner context
1215: - prefix - the prefix string to prepend to all PC option requests
1217: Notes:
1218: A hyphen (-) must NOT be given at the beginning of the prefix name.
1219: The first character of all runtime options is AUTOMATICALLY the
1220: hyphen.
1222: Level: advanced
1224: .keywords: PC, append, options, prefix, database
1226: .seealso: PCSetOptionsPrefix(), PCGetOptionsPrefix()
1227: @*/
1228: int PCAppendOptionsPrefix(PC pc,char *prefix)
1229: {
1234: PetscObjectAppendOptionsPrefix((PetscObject)pc,prefix);
1235: return(0);
1236: }
1238: #undef __FUNCT__
1240: /*@C
1241: PCGetOptionsPrefix - Gets the prefix used for searching for all
1242: PC options in the database.
1244: Not Collective
1246: Input Parameters:
1247: . pc - the preconditioner context
1249: Output Parameters:
1250: . prefix - pointer to the prefix string used, is returned
1252: Notes: On the fortran side, the user should pass in a string 'prifix' of
1253: sufficient length to hold the prefix.
1255: Level: advanced
1257: .keywords: PC, get, options, prefix, database
1259: .seealso: PCSetOptionsPrefix(), PCAppendOptionsPrefix()
1260: @*/
1261: int PCGetOptionsPrefix(PC pc,char **prefix)
1262: {
1267: PetscObjectGetOptionsPrefix((PetscObject)pc,prefix);
1268: return(0);
1269: }
1271: #undef __FUNCT__
1273: /*@
1274: PCPreSolve - Optional pre-solve phase, intended for any
1275: preconditioner-specific actions that must be performed before
1276: the iterative solve itself.
1278: Collective on PC
1280: Input Parameters:
1281: + pc - the preconditioner context
1282: - ksp - the Krylov subspace context
1284: Level: developer
1286: Sample of Usage:
1287: .vb
1288: PCPreSolve(pc,ksp);
1289: KSPSolve(ksp,its);
1290: PCPostSolve(pc,ksp);
1291: .ve
1293: Notes:
1294: The pre-solve phase is distinct from the PCSetUp() phase.
1296: SLESSolve() calls this directly, so is rarely called by the user.
1298: .keywords: PC, pre-solve
1300: .seealso: PCPostSolve()
1301: @*/
1302: int PCPreSolve(PC pc,KSP ksp)
1303: {
1305: Vec x,rhs;
1306: Mat A,B;
1311: KSPGetSolution(ksp,&x);
1312: KSPGetRhs(ksp,&rhs);
1313: /*
1314: Scale the system and have the matrices use the scaled form
1315: only if the two matrices are actually the same (and hence
1316: have the same scaling
1317: */
1318: PCGetOperators(pc,&A,&B,PETSC_NULL);
1319: if (A == B) {
1320: MatScaleSystem(pc->mat,x,rhs);
1321: MatUseScaledForm(pc->mat,PETSC_TRUE);
1322: }
1324: if (pc->ops->presolve) {
1325: (*pc->ops->presolve)(pc,ksp,x,rhs);
1326: }
1327: return(0);
1328: }
1330: #undef __FUNCT__
1332: /*@
1333: PCPostSolve - Optional post-solve phase, intended for any
1334: preconditioner-specific actions that must be performed after
1335: the iterative solve itself.
1337: Collective on PC
1339: Input Parameters:
1340: + pc - the preconditioner context
1341: - ksp - the Krylov subspace context
1343: Sample of Usage:
1344: .vb
1345: PCPreSolve(pc,ksp);
1346: KSPSolve(ksp,its);
1347: PCPostSolve(pc,ksp);
1348: .ve
1350: Note:
1351: SLESSolve() calls this routine directly, so it is rarely called by the user.
1353: Level: developer
1355: .keywords: PC, post-solve
1357: .seealso: PCPreSolve(), SLESSolve()
1358: @*/
1359: int PCPostSolve(PC pc,KSP ksp)
1360: {
1362: Vec x,rhs;
1363: Mat A,B;
1367: KSPGetSolution(ksp,&x);
1368: KSPGetRhs(ksp,&rhs);
1369: if (pc->ops->postsolve) {
1370: (*pc->ops->postsolve)(pc,ksp,x,rhs);
1371: }
1373: /*
1374: Scale the system and have the matrices use the scaled form
1375: only if the two matrices are actually the same (and hence
1376: have the same scaling
1377: */
1378: PCGetOperators(pc,&A,&B,PETSC_NULL);
1379: if (A == B) {
1380: MatUnScaleSystem(pc->mat,x,rhs);
1381: MatUseScaledForm(pc->mat,PETSC_FALSE);
1382: }
1383: return(0);
1384: }
1386: #undef __FUNCT__
1388: /*@C
1389: PCView - Prints the PC data structure.
1391: Collective on PC
1393: Input Parameters:
1394: + PC - the PC context
1395: - viewer - optional visualization context
1397: Note:
1398: The available visualization contexts include
1399: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
1400: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
1401: output where only the first processor opens
1402: the file. All other processors send their
1403: data to the first processor to print.
1405: The user can open an alternative visualization contexts with
1406: PetscViewerASCIIOpen() (output to a specified file).
1408: Level: developer
1410: .keywords: PC, view
1412: .seealso: KSPView(), PetscViewerASCIIOpen()
1413: @*/
1414: int PCView(PC pc,PetscViewer viewer)
1415: {
1416: PCType cstr;
1417: int ierr;
1418: PetscTruth mat_exists,isascii,isstring;
1419: PetscViewerFormat format;
1423: if (!viewer) viewer = PETSC_VIEWER_STDOUT_(pc->comm);
1427: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
1428: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
1429: if (isascii) {
1430: PetscViewerGetFormat(viewer,&format);
1431: if (pc->prefix) {
1432: PetscViewerASCIIPrintf(viewer,"PC Object:(%s)n",pc->prefix);
1433: } else {
1434: PetscViewerASCIIPrintf(viewer,"PC Object:n");
1435: }
1436: PCGetType(pc,&cstr);
1437: if (cstr) {
1438: PetscViewerASCIIPrintf(viewer," type: %sn",cstr);
1439: } else {
1440: PetscViewerASCIIPrintf(viewer," type: not yet setn");
1441: }
1442: if (pc->ops->view) {
1443: PetscViewerASCIIPushTab(viewer);
1444: (*pc->ops->view)(pc,viewer);
1445: PetscViewerASCIIPopTab(viewer);
1446: }
1447: PetscObjectExists((PetscObject)pc->mat,&mat_exists);
1448: if (mat_exists) {
1449: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
1450: if (pc->pmat == pc->mat) {
1451: PetscViewerASCIIPrintf(viewer," linear system matrix = precond matrix:n");
1452: PetscViewerASCIIPushTab(viewer);
1453: MatView(pc->mat,viewer);
1454: PetscViewerASCIIPopTab(viewer);
1455: } else {
1456: PetscObjectExists((PetscObject)pc->pmat,&mat_exists);
1457: if (mat_exists) {
1458: PetscViewerASCIIPrintf(viewer," linear system matrix followed by preconditioner matrix:n");
1459: } else {
1460: PetscViewerASCIIPrintf(viewer," linear system matrix:n");
1461: }
1462: PetscViewerASCIIPushTab(viewer);
1463: MatView(pc->mat,viewer);
1464: if (mat_exists) {MatView(pc->pmat,viewer);}
1465: PetscViewerASCIIPopTab(viewer);
1466: }
1467: PetscViewerPopFormat(viewer);
1468: }
1469: } else if (isstring) {
1470: PCGetType(pc,&cstr);
1471: PetscViewerStringSPrintf(viewer," %-7.7s",cstr);
1472: if (pc->ops->view) {(*pc->ops->view)(pc,viewer);}
1473: } else {
1474: SETERRQ1(1,"Viewer type %s not supported by PC",((PetscObject)viewer)->type_name);
1475: }
1476: return(0);
1477: }
1479: /*MC
1480: PCRegisterDynamic - Adds a method to the preconditioner package.
1482: Synopsis:
1483: int PCRegisterDynamic(char *name_solver,char *path,char *name_create,int (*routine_create)(PC))
1485: Not collective
1487: Input Parameters:
1488: + name_solver - name of a new user-defined solver
1489: . path - path (either absolute or relative) the library containing this solver
1490: . name_create - name of routine to create method context
1491: - routine_create - routine to create method context
1493: Notes:
1494: PCRegisterDynamic() may be called multiple times to add several user-defined preconditioners.
1496: If dynamic libraries are used, then the fourth input argument (routine_create)
1497: is ignored.
1499: Sample usage:
1500: .vb
1501: PCRegisterDynamic("my_solver","/home/username/my_lib/lib/libO/solaris/mylib",
1502: "MySolverCreate",MySolverCreate);
1503: .ve
1505: Then, your solver can be chosen with the procedural interface via
1506: $ PCSetType(pc,"my_solver")
1507: or at runtime via the option
1508: $ -pc_type my_solver
1510: Level: advanced
1512: Notes: ${PETSC_ARCH}, ${PETSC_DIR}, ${PETSC_LIB_DIR}, ${BOPT}, or ${any environmental variable}
1513: occuring in pathname will be replaced with appropriate values.
1514: If your function is not being put into a shared library then use PCRegister() instead
1516: .keywords: PC, register
1518: .seealso: PCRegisterAll(), PCRegisterDestroy()
1519: M*/
1521: #undef __FUNCT__
1523: /*@C
1524: PCRegister - See PCRegisterDynamic()
1526: Level: advanced
1527: @*/
1528: int PCRegister(char *sname,char *path,char *name,int (*function)(PC))
1529: {
1530: int ierr;
1531: char fullname[256];
1535: PetscFListConcat(path,name,fullname);
1536: PetscFListAdd(&PCList,sname,fullname,(void (*)(void))function);
1537: return(0);
1538: }
1540: #undef __FUNCT__
1542: /*@
1543: PCComputeExplicitOperator - Computes the explicit preconditioned operator.
1545: Collective on PC
1547: Input Parameter:
1548: . pc - the preconditioner object
1550: Output Parameter:
1551: . mat - the explict preconditioned operator
1553: Notes:
1554: This computation is done by applying the operators to columns of the
1555: identity matrix.
1557: Currently, this routine uses a dense matrix format when 1 processor
1558: is used and a sparse format otherwise. This routine is costly in general,
1559: and is recommended for use only with relatively small systems.
1561: Level: advanced
1562:
1563: .keywords: PC, compute, explicit, operator
1565: @*/
1566: int PCComputeExplicitOperator(PC pc,Mat *mat)
1567: {
1568: Vec in,out;
1569: int ierr,i,M,m,size,*rows,start,end;
1570: MPI_Comm comm;
1571: PetscScalar *array,zero = 0.0,one = 1.0;
1577: comm = pc->comm;
1578: MPI_Comm_size(comm,&size);
1580: PCGetVector(pc,&in);
1581: VecDuplicate(in,&out);
1582: VecGetOwnershipRange(in,&start,&end);
1583: VecGetSize(in,&M);
1584: VecGetLocalSize(in,&m);
1585: PetscMalloc((m+1)*sizeof(int),&rows);
1586: for (i=0; i<m; i++) {rows[i] = start + i;}
1588: if (size == 1) {
1589: MatCreateSeqDense(comm,M,M,PETSC_NULL,mat);
1590: } else {
1591: MatCreateMPIAIJ(comm,m,m,M,M,0,0,0,0,mat);
1592: }
1594: for (i=0; i<M; i++) {
1596: VecSet(&zero,in);
1597: VecSetValues(in,1,&i,&one,INSERT_VALUES);
1598: VecAssemblyBegin(in);
1599: VecAssemblyEnd(in);
1601: PCApply(pc,in,out);
1602:
1603: VecGetArray(out,&array);
1604: MatSetValues(*mat,m,rows,1,&i,array,INSERT_VALUES);
1605: VecRestoreArray(out,&array);
1607: }
1608: PetscFree(rows);
1609: VecDestroy(out);
1610: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
1611: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
1612: return(0);
1613: }