Actual source code: jacobi.c
2: /* --------------------------------------------------------------------
4: This file implements a Jacobi preconditioner in PETSc as part of PC.
5: You can use this as a starting point for implementing your own
6: preconditioner that is not provided with PETSc. (You might also consider
7: just using PCSHELL)
9: The following basic routines are required for each preconditioner.
10: PCCreate_XXX() - Creates a preconditioner context
11: PCSetFromOptions_XXX() - Sets runtime options
12: PCApply_XXX() - Applies the preconditioner
13: PCDestroy_XXX() - Destroys the preconditioner context
14: where the suffix "_XXX" denotes a particular implementation, in
15: this case we use _Jacobi (e.g., PCCreate_Jacobi, PCApply_Jacobi).
16: These routines are actually called via the common user interface
17: routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(),
18: so the application code interface remains identical for all
19: preconditioners.
21: Another key routine is:
22: PCSetUp_XXX() - Prepares for the use of a preconditioner
23: by setting data structures and options. The interface routine PCSetUp()
24: is not usually called directly by the user, but instead is called by
25: PCApply() if necessary.
27: Additional basic routines are:
28: PCView_XXX() - Prints details of runtime options that
29: have actually been used.
30: These are called by application codes via the interface routines
31: PCView().
33: The various types of solvers (preconditioners, Krylov subspace methods,
34: nonlinear solvers, timesteppers) are all organized similarly, so the
35: above description applies to these categories also. One exception is
36: that the analogues of PCApply() for these components are KSPSolve(),
37: SNESSolve(), and TSSolve().
39: Additional optional functionality unique to preconditioners is left and
40: right symmetric preconditioner application via PCApplySymmetricLeft()
41: and PCApplySymmetricRight(). The Jacobi implementation is
42: PCApplySymmetricLeftOrRight_Jacobi().
44: -------------------------------------------------------------------- */
46: /*
47: Include files needed for the Jacobi preconditioner:
48: pcimpl.h - private include file intended for use by all preconditioners
49: */
51: #include <petsc/private/pcimpl.h>
53: const char *const PCJacobiTypes[] = {"DIAGONAL","ROWMAX","ROWSUM","PCJacobiType","PC_JACOBI_",NULL};
55: /*
56: Private context (data structure) for the Jacobi preconditioner.
57: */
58: typedef struct {
59: Vec diag; /* vector containing the reciprocals of the diagonal elements of the preconditioner matrix */
60: Vec diagsqrt; /* vector containing the reciprocals of the square roots of
61: the diagonal elements of the preconditioner matrix (used
62: only for symmetric preconditioner application) */
63: PetscBool userowmax; /* set with PCJacobiSetType() */
64: PetscBool userowsum;
65: PetscBool useabs; /* use the absolute values of the diagonal entries */
66: PetscBool fixdiag; /* fix zero diagonal terms */
67: } PC_Jacobi;
69: static PetscErrorCode PCJacobiSetType_Jacobi(PC pc,PCJacobiType type)
70: {
71: PC_Jacobi *j = (PC_Jacobi*)pc->data;
74: j->userowmax = PETSC_FALSE;
75: j->userowsum = PETSC_FALSE;
76: if (type == PC_JACOBI_ROWMAX) {
77: j->userowmax = PETSC_TRUE;
78: } else if (type == PC_JACOBI_ROWSUM) {
79: j->userowsum = PETSC_TRUE;
80: }
81: return(0);
82: }
84: static PetscErrorCode PCJacobiGetType_Jacobi(PC pc,PCJacobiType *type)
85: {
86: PC_Jacobi *j = (PC_Jacobi*)pc->data;
89: if (j->userowmax) {
90: *type = PC_JACOBI_ROWMAX;
91: } else if (j->userowsum) {
92: *type = PC_JACOBI_ROWSUM;
93: } else {
94: *type = PC_JACOBI_DIAGONAL;
95: }
96: return(0);
97: }
99: static PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc,PetscBool flg)
100: {
101: PC_Jacobi *j = (PC_Jacobi*)pc->data;
104: j->useabs = flg;
105: return(0);
106: }
108: static PetscErrorCode PCJacobiGetUseAbs_Jacobi(PC pc,PetscBool *flg)
109: {
110: PC_Jacobi *j = (PC_Jacobi*)pc->data;
113: *flg = j->useabs;
114: return(0);
115: }
117: static PetscErrorCode PCJacobiSetFixDiagonal_Jacobi(PC pc,PetscBool flg)
118: {
119: PC_Jacobi *j = (PC_Jacobi*)pc->data;
122: j->fixdiag = flg;
123: return(0);
124: }
126: static PetscErrorCode PCJacobiGetFixDiagonal_Jacobi(PC pc,PetscBool *flg)
127: {
128: PC_Jacobi *j = (PC_Jacobi*)pc->data;
131: *flg = j->fixdiag;
132: return(0);
133: }
135: /* -------------------------------------------------------------------------- */
136: /*
137: PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
138: by setting data structures and options.
140: Input Parameter:
141: . pc - the preconditioner context
143: Application Interface Routine: PCSetUp()
145: Notes:
146: The interface routine PCSetUp() is not usually called directly by
147: the user, but instead is called by PCApply() if necessary.
148: */
149: static PetscErrorCode PCSetUp_Jacobi(PC pc)
150: {
151: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
152: Vec diag,diagsqrt;
154: PetscInt n,i;
155: PetscScalar *x;
156: PetscBool zeroflag = PETSC_FALSE;
159: /*
160: For most preconditioners the code would begin here something like
162: if (pc->setupcalled == 0) { allocate space the first time this is ever called
163: MatCreateVecs(pc->mat,&jac->diag);
164: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
165: }
167: But for this preconditioner we want to support use of both the matrix' diagonal
168: elements (for left or right preconditioning) and square root of diagonal elements
169: (for symmetric preconditioning). Hence we do not allocate space here, since we
170: don't know at this point which will be needed (diag and/or diagsqrt) until the user
171: applies the preconditioner, and we don't want to allocate BOTH unless we need
172: them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
173: and PCSetUp_Jacobi_Symmetric(), respectively.
174: */
176: /*
177: Here we set up the preconditioner; that is, we copy the diagonal values from
178: the matrix and put them into a format to make them quick to apply as a preconditioner.
179: */
180: diag = jac->diag;
181: diagsqrt = jac->diagsqrt;
183: if (diag) {
184: PetscBool isspd;
186: if (jac->userowmax) {
187: MatGetRowMaxAbs(pc->pmat,diag,NULL);
188: } else if (jac->userowsum) {
189: MatGetRowSum(pc->pmat,diag);
190: } else {
191: MatGetDiagonal(pc->pmat,diag);
192: }
193: VecReciprocal(diag);
194: if (jac->useabs) {
195: VecAbs(diag);
196: }
197: MatGetOption(pc->pmat,MAT_SPD,&isspd);
198: if (jac->fixdiag && !isspd) {
199: VecGetLocalSize(diag,&n);
200: VecGetArray(diag,&x);
201: for (i=0; i<n; i++) {
202: if (x[i] == 0.0) {
203: x[i] = 1.0;
204: zeroflag = PETSC_TRUE;
205: }
206: }
207: VecRestoreArray(diag,&x);
208: }
209: }
210: if (diagsqrt) {
211: if (jac->userowmax) {
212: MatGetRowMaxAbs(pc->pmat,diagsqrt,NULL);
213: } else if (jac->userowsum) {
214: MatGetRowSum(pc->pmat,diagsqrt);
215: } else {
216: MatGetDiagonal(pc->pmat,diagsqrt);
217: }
218: VecGetLocalSize(diagsqrt,&n);
219: VecGetArray(diagsqrt,&x);
220: for (i=0; i<n; i++) {
221: if (x[i] != 0.0) x[i] = 1.0/PetscSqrtReal(PetscAbsScalar(x[i]));
222: else {
223: x[i] = 1.0;
224: zeroflag = PETSC_TRUE;
225: }
226: }
227: VecRestoreArray(diagsqrt,&x);
228: }
229: if (zeroflag) {
230: PetscInfo(pc,"Zero detected in diagonal of matrix, using 1 at those locations\n");
231: }
232: return(0);
233: }
234: /* -------------------------------------------------------------------------- */
235: /*
236: PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
237: inverse of the square root of the diagonal entries of the matrix. This
238: is used for symmetric application of the Jacobi preconditioner.
240: Input Parameter:
241: . pc - the preconditioner context
242: */
243: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
244: {
246: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
249: MatCreateVecs(pc->pmat,&jac->diagsqrt,NULL);
250: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diagsqrt);
251: PCSetUp_Jacobi(pc);
252: return(0);
253: }
254: /* -------------------------------------------------------------------------- */
255: /*
256: PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
257: inverse of the diagonal entries of the matrix. This is used for left of
258: right application of the Jacobi preconditioner.
260: Input Parameter:
261: . pc - the preconditioner context
262: */
263: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
264: {
266: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
269: MatCreateVecs(pc->pmat,&jac->diag,NULL);
270: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
271: PCSetUp_Jacobi(pc);
272: return(0);
273: }
274: /* -------------------------------------------------------------------------- */
275: /*
276: PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
278: Input Parameters:
279: . pc - the preconditioner context
280: . x - input vector
282: Output Parameter:
283: . y - output vector
285: Application Interface Routine: PCApply()
286: */
287: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
288: {
289: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
293: if (!jac->diag) {
294: PCSetUp_Jacobi_NonSymmetric(pc);
295: }
296: VecPointwiseMult(y,x,jac->diag);
297: return(0);
298: }
299: /* -------------------------------------------------------------------------- */
300: /*
301: PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
302: symmetric preconditioner to a vector.
304: Input Parameters:
305: . pc - the preconditioner context
306: . x - input vector
308: Output Parameter:
309: . y - output vector
311: Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
312: */
313: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
314: {
316: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
319: if (!jac->diagsqrt) {
320: PCSetUp_Jacobi_Symmetric(pc);
321: }
322: VecPointwiseMult(y,x,jac->diagsqrt);
323: return(0);
324: }
326: /* -------------------------------------------------------------------------- */
327: static PetscErrorCode PCReset_Jacobi(PC pc)
328: {
329: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
333: VecDestroy(&jac->diag);
334: VecDestroy(&jac->diagsqrt);
335: return(0);
336: }
338: /*
339: PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
340: that was created with PCCreate_Jacobi().
342: Input Parameter:
343: . pc - the preconditioner context
345: Application Interface Routine: PCDestroy()
346: */
347: static PetscErrorCode PCDestroy_Jacobi(PC pc)
348: {
352: PCReset_Jacobi(pc);
353: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",NULL);
354: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",NULL);
355: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",NULL);
356: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",NULL);
357: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetFixDiagonal_C",NULL);
358: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetFixDiagonal_C",NULL);
360: /*
361: Free the private data structure that was hanging off the PC
362: */
363: PetscFree(pc->data);
364: return(0);
365: }
367: static PetscErrorCode PCSetFromOptions_Jacobi(PetscOptionItems *PetscOptionsObject,PC pc)
368: {
369: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
371: PetscBool flg;
372: PCJacobiType deflt,type;
375: PCJacobiGetType(pc,&deflt);
376: PetscOptionsHead(PetscOptionsObject,"Jacobi options");
377: PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCJacobiSetType",PCJacobiTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);
378: if (flg) {
379: PCJacobiSetType(pc,type);
380: }
381: PetscOptionsBool("-pc_jacobi_abs","Use absolute values of diagonal entries","PCJacobiSetUseAbs",jac->useabs,&jac->useabs,NULL);
382: PetscOptionsBool("-pc_jacobi_fixdiagonal","Fix null terms on diagonal","PCJacobiSetFixDiagonal",jac->fixdiag,&jac->fixdiag,NULL);
383: PetscOptionsTail();
384: return(0);
385: }
387: static PetscErrorCode PCView_Jacobi(PC pc, PetscViewer viewer)
388: {
389: PC_Jacobi *jac = (PC_Jacobi *) pc->data;
390: PetscBool iascii;
394: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
395: if (iascii) {
396: PCJacobiType type;
397: PetscBool useAbs,fixdiag;
398: PetscViewerFormat format;
400: PCJacobiGetType(pc, &type);
401: PCJacobiGetUseAbs(pc, &useAbs);
402: PCJacobiGetFixDiagonal(pc, &fixdiag);
403: PetscViewerASCIIPrintf(viewer, " type %s%s%s\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "", !fixdiag ? ", not checking null diagonal entries" : "");
404: PetscViewerGetFormat(viewer, &format);
405: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
406: VecView(jac->diag, viewer);
407: }
408: }
409: return(0);
410: }
412: /* -------------------------------------------------------------------------- */
413: /*
414: PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
415: and sets this as the private data within the generic preconditioning
416: context, PC, that was created within PCCreate().
418: Input Parameter:
419: . pc - the preconditioner context
421: Application Interface Routine: PCCreate()
422: */
424: /*MC
425: PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)
427: Options Database Key:
428: + -pc_jacobi_type <diagonal,rowmax,rowsum> - approach for forming the preconditioner
429: . -pc_jacobi_abs - use the absolute value of the diagonal entry
430: - -pc_jacobi_fixdiag - fix for zero diagonal terms
432: Level: beginner
434: Notes:
435: By using KSPSetPCSide(ksp,PC_SYMMETRIC) or -ksp_pc_side symmetric
436: can scale each side of the matrix by the square root of the diagonal entries.
438: Zero entries along the diagonal are replaced with the value 1.0
440: See PCPBJACOBI for a point-block Jacobi preconditioner
442: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
443: PCJacobiSetType(), PCJacobiSetUseAbs(), PCJacobiGetUseAbs(),
444: PCJacobiSetFixDiagonal(), PCJacobiGetFixDiagonal(), PCPBJACOBI
445: M*/
447: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
448: {
449: PC_Jacobi *jac;
453: /*
454: Creates the private data structure for this preconditioner and
455: attach it to the PC object.
456: */
457: PetscNewLog(pc,&jac);
458: pc->data = (void*)jac;
460: /*
461: Initialize the pointers to vectors to ZERO; these will be used to store
462: diagonal entries of the matrix for fast preconditioner application.
463: */
464: jac->diag = NULL;
465: jac->diagsqrt = NULL;
466: jac->userowmax = PETSC_FALSE;
467: jac->userowsum = PETSC_FALSE;
468: jac->useabs = PETSC_FALSE;
469: jac->fixdiag = PETSC_TRUE;
471: /*
472: Set the pointers for the functions that are provided above.
473: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
474: are called, they will automatically call these functions. Note we
475: choose not to provide a couple of these functions since they are
476: not needed.
477: */
478: pc->ops->apply = PCApply_Jacobi;
479: pc->ops->applytranspose = PCApply_Jacobi;
480: pc->ops->setup = PCSetUp_Jacobi;
481: pc->ops->reset = PCReset_Jacobi;
482: pc->ops->destroy = PCDestroy_Jacobi;
483: pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
484: pc->ops->view = PCView_Jacobi;
485: pc->ops->applyrichardson = NULL;
486: pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
487: pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
489: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",PCJacobiSetType_Jacobi);
490: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",PCJacobiGetType_Jacobi);
491: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",PCJacobiSetUseAbs_Jacobi);
492: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",PCJacobiGetUseAbs_Jacobi);
493: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetFixDiagonal_C",PCJacobiSetFixDiagonal_Jacobi);
494: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetFixDiagonal_C",PCJacobiGetFixDiagonal_Jacobi);
495: return(0);
496: }
498: /*@
499: PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
500: absolute values of the diagonal divisors in the preconditioner
502: Logically Collective on PC
504: Input Parameters:
505: + pc - the preconditioner context
506: - flg - whether to use absolute values or not
508: Options Database Key:
509: . -pc_jacobi_abs
511: Notes:
512: This takes affect at the next construction of the preconditioner
514: Level: intermediate
516: .seealso: PCJacobiaSetType(), PCJacobiGetUseAbs()
518: @*/
519: PetscErrorCode PCJacobiSetUseAbs(PC pc,PetscBool flg)
520: {
525: PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC,PetscBool),(pc,flg));
526: return(0);
527: }
529: /*@
530: PCJacobiGetUseAbs - Determines if the Jacobi preconditioner uses the
531: absolute values of the diagonal divisors in the preconditioner
533: Logically Collective on PC
535: Input Parameter:
536: . pc - the preconditioner context
538: Output Parameter:
539: . flg - whether to use absolute values or not
541: Options Database Key:
542: . -pc_jacobi_abs
544: Level: intermediate
546: .seealso: PCJacobiaSetType(), PCJacobiSetUseAbs(), PCJacobiGetType()
548: @*/
549: PetscErrorCode PCJacobiGetUseAbs(PC pc,PetscBool *flg)
550: {
555: PetscUseMethod(pc,"PCJacobiGetUseAbs_C",(PC,PetscBool*),(pc,flg));
556: return(0);
557: }
559: /*@
560: PCJacobiSetFixDiagonal - Do not check for zero values on diagonal
562: Logically Collective on PC
564: Input Parameters:
565: + pc - the preconditioner context
566: - flg - the boolean flag
568: Options Database Key:
569: . -pc_jacobi_fixdiagonal
571: Notes:
572: This takes affect at the next construction of the preconditioner
574: Level: intermediate
576: .seealso: PCJacobiSetType(), PCJacobiGetFixDiagonal()
578: @*/
579: PetscErrorCode PCJacobiSetFixDiagonal(PC pc,PetscBool flg)
580: {
585: PetscTryMethod(pc,"PCJacobiSetFixDiagonal_C",(PC,PetscBool),(pc,flg));
586: return(0);
587: }
589: /*@
590: PCJacobiGetFixDiagonal - Determines if the Jacobi preconditioner checks for zero diagonal terms
592: Logically Collective on PC
594: Input Parameter:
595: . pc - the preconditioner context
597: Output Parameter:
598: . flg - the boolean flag
600: Options Database Key:
601: . -pc_jacobi_fixdiagonal
603: Level: intermediate
605: .seealso: PCJacobiSetType(), PCJacobiSetFixDiagonal()
607: @*/
608: PetscErrorCode PCJacobiGetFixDiagonal(PC pc,PetscBool *flg)
609: {
614: PetscUseMethod(pc,"PCJacobiGetFixDiagonal_C",(PC,PetscBool*),(pc,flg));
615: return(0);
616: }
618: /*@
619: PCJacobiSetType - Causes the Jacobi preconditioner to use either the diagonal, the maximum entry in each row,
620: of the sum of rows entries for the diagonal preconditioner
622: Logically Collective on PC
624: Input Parameters:
625: + pc - the preconditioner context
626: - type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM
628: Options Database Key:
629: . -pc_jacobi_type <diagonal,rowmax,rowsum>
631: Level: intermediate
633: .seealso: PCJacobiaUseAbs(), PCJacobiGetType()
634: @*/
635: PetscErrorCode PCJacobiSetType(PC pc,PCJacobiType type)
636: {
641: PetscTryMethod(pc,"PCJacobiSetType_C",(PC,PCJacobiType),(pc,type));
642: return(0);
643: }
645: /*@
646: PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner
648: Not Collective on PC
650: Input Parameter:
651: . pc - the preconditioner context
653: Output Parameter:
654: . type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM
656: Level: intermediate
658: .seealso: PCJacobiaUseAbs(), PCJacobiSetType()
659: @*/
660: PetscErrorCode PCJacobiGetType(PC pc,PCJacobiType *type)
661: {
666: PetscUseMethod(pc,"PCJacobiGetType_C",(PC,PCJacobiType*),(pc,type));
667: return(0);
668: }