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 PCReset_Jacobi(PC);
71: static PetscErrorCode PCJacobiSetType_Jacobi(PC pc, PCJacobiType type)
72: {
73: PC_Jacobi *j = (PC_Jacobi *)pc->data;
74: PCJacobiType old_type;
76: PCJacobiGetType(pc, &old_type);
77: if (old_type == type) return 0;
78: PCReset_Jacobi(pc);
79: j->userowmax = PETSC_FALSE;
80: j->userowsum = PETSC_FALSE;
81: if (type == PC_JACOBI_ROWMAX) {
82: j->userowmax = PETSC_TRUE;
83: } else if (type == PC_JACOBI_ROWSUM) {
84: j->userowsum = PETSC_TRUE;
85: }
86: return 0;
87: }
89: static PetscErrorCode PCJacobiGetType_Jacobi(PC pc, PCJacobiType *type)
90: {
91: PC_Jacobi *j = (PC_Jacobi *)pc->data;
93: if (j->userowmax) {
94: *type = PC_JACOBI_ROWMAX;
95: } else if (j->userowsum) {
96: *type = PC_JACOBI_ROWSUM;
97: } else {
98: *type = PC_JACOBI_DIAGONAL;
99: }
100: return 0;
101: }
103: static PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc, PetscBool flg)
104: {
105: PC_Jacobi *j = (PC_Jacobi *)pc->data;
107: j->useabs = flg;
108: return 0;
109: }
111: static PetscErrorCode PCJacobiGetUseAbs_Jacobi(PC pc, PetscBool *flg)
112: {
113: PC_Jacobi *j = (PC_Jacobi *)pc->data;
115: *flg = j->useabs;
116: return 0;
117: }
119: static PetscErrorCode PCJacobiSetFixDiagonal_Jacobi(PC pc, PetscBool flg)
120: {
121: PC_Jacobi *j = (PC_Jacobi *)pc->data;
123: j->fixdiag = flg;
124: return 0;
125: }
127: static PetscErrorCode PCJacobiGetFixDiagonal_Jacobi(PC pc, PetscBool *flg)
128: {
129: PC_Jacobi *j = (PC_Jacobi *)pc->data;
131: *flg = j->fixdiag;
132: return 0;
133: }
135: /*
136: PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
137: by setting data structures and options.
139: Input Parameter:
140: . pc - the preconditioner context
142: Application Interface Routine: PCSetUp()
144: Note:
145: The interface routine PCSetUp() is not usually called directly by
146: the user, but instead is called by PCApply() if necessary.
147: */
148: static PetscErrorCode PCSetUp_Jacobi(PC pc)
149: {
150: PC_Jacobi *jac = (PC_Jacobi *)pc->data;
151: Vec diag, diagsqrt;
152: PetscInt n, i;
153: PetscScalar *x;
154: PetscBool zeroflag = PETSC_FALSE;
156: /*
157: For most preconditioners the code would begin here something like
159: if (pc->setupcalled == 0) { allocate space the first time this is ever called
160: MatCreateVecs(pc->mat,&jac->diag);
161: }
163: But for this preconditioner we want to support use of both the matrix' diagonal
164: elements (for left or right preconditioning) and square root of diagonal elements
165: (for symmetric preconditioning). Hence we do not allocate space here, since we
166: don't know at this point which will be needed (diag and/or diagsqrt) until the user
167: applies the preconditioner, and we don't want to allocate BOTH unless we need
168: them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
169: and PCSetUp_Jacobi_Symmetric(), respectively.
170: */
172: /*
173: Here we set up the preconditioner; that is, we copy the diagonal values from
174: the matrix and put them into a format to make them quick to apply as a preconditioner.
175: */
176: diag = jac->diag;
177: diagsqrt = jac->diagsqrt;
179: if (diag) {
180: PetscBool isset, isspd;
182: if (jac->userowmax) {
183: MatGetRowMaxAbs(pc->pmat, diag, NULL);
184: } else if (jac->userowsum) {
185: MatGetRowSum(pc->pmat, diag);
186: } else {
187: MatGetDiagonal(pc->pmat, diag);
188: }
189: VecReciprocal(diag);
190: if (jac->useabs) VecAbs(diag);
191: MatIsSPDKnown(pc->pmat, &isset, &isspd);
192: if (jac->fixdiag && (!isset || !isspd)) {
193: VecGetLocalSize(diag, &n);
194: VecGetArray(diag, &x);
195: for (i = 0; i < n; i++) {
196: if (x[i] == 0.0) {
197: x[i] = 1.0;
198: zeroflag = PETSC_TRUE;
199: }
200: }
201: VecRestoreArray(diag, &x);
202: }
203: }
204: if (diagsqrt) {
205: if (jac->userowmax) {
206: MatGetRowMaxAbs(pc->pmat, diagsqrt, NULL);
207: } else if (jac->userowsum) {
208: MatGetRowSum(pc->pmat, diagsqrt);
209: } else {
210: MatGetDiagonal(pc->pmat, diagsqrt);
211: }
212: VecGetLocalSize(diagsqrt, &n);
213: VecGetArray(diagsqrt, &x);
214: for (i = 0; i < n; i++) {
215: if (x[i] != 0.0) x[i] = 1.0 / PetscSqrtReal(PetscAbsScalar(x[i]));
216: else {
217: x[i] = 1.0;
218: zeroflag = PETSC_TRUE;
219: }
220: }
221: VecRestoreArray(diagsqrt, &x);
222: }
223: if (zeroflag) PetscInfo(pc, "Zero detected in diagonal of matrix, using 1 at those locations\n");
224: return 0;
225: }
227: /*
228: PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
229: inverse of the square root of the diagonal entries of the matrix. This
230: is used for symmetric application of the Jacobi preconditioner.
232: Input Parameter:
233: . pc - the preconditioner context
234: */
235: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
236: {
237: PC_Jacobi *jac = (PC_Jacobi *)pc->data;
239: MatCreateVecs(pc->pmat, &jac->diagsqrt, NULL);
240: PCSetUp_Jacobi(pc);
241: return 0;
242: }
244: /*
245: PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
246: inverse of the diagonal entries of the matrix. This is used for left of
247: right application of the Jacobi preconditioner.
249: Input Parameter:
250: . pc - the preconditioner context
251: */
252: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
253: {
254: PC_Jacobi *jac = (PC_Jacobi *)pc->data;
256: MatCreateVecs(pc->pmat, &jac->diag, NULL);
257: PCSetUp_Jacobi(pc);
258: return 0;
259: }
261: /*
262: PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
264: Input Parameters:
265: . pc - the preconditioner context
266: . x - input vector
268: Output Parameter:
269: . y - output vector
271: Application Interface Routine: PCApply()
272: */
273: static PetscErrorCode PCApply_Jacobi(PC pc, Vec x, Vec y)
274: {
275: PC_Jacobi *jac = (PC_Jacobi *)pc->data;
277: if (!jac->diag) PCSetUp_Jacobi_NonSymmetric(pc);
278: VecPointwiseMult(y, x, jac->diag);
279: return 0;
280: }
282: /*
283: PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
284: symmetric preconditioner to a vector.
286: Input Parameters:
287: . pc - the preconditioner context
288: . x - input vector
290: Output Parameter:
291: . y - output vector
293: Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
294: */
295: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc, Vec x, Vec y)
296: {
297: PC_Jacobi *jac = (PC_Jacobi *)pc->data;
299: if (!jac->diagsqrt) PCSetUp_Jacobi_Symmetric(pc);
300: VecPointwiseMult(y, x, jac->diagsqrt);
301: return 0;
302: }
304: static PetscErrorCode PCReset_Jacobi(PC pc)
305: {
306: PC_Jacobi *jac = (PC_Jacobi *)pc->data;
308: VecDestroy(&jac->diag);
309: VecDestroy(&jac->diagsqrt);
310: return 0;
311: }
313: /*
314: PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
315: that was created with PCCreate_Jacobi().
317: Input Parameter:
318: . pc - the preconditioner context
320: Application Interface Routine: PCDestroy()
321: */
322: static PetscErrorCode PCDestroy_Jacobi(PC pc)
323: {
324: PCReset_Jacobi(pc);
325: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetType_C", NULL);
326: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetType_C", NULL);
327: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetUseAbs_C", NULL);
328: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetUseAbs_C", NULL);
329: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetFixDiagonal_C", NULL);
330: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetFixDiagonal_C", NULL);
332: /*
333: Free the private data structure that was hanging off the PC
334: */
335: PetscFree(pc->data);
336: return 0;
337: }
339: static PetscErrorCode PCSetFromOptions_Jacobi(PC pc, PetscOptionItems *PetscOptionsObject)
340: {
341: PC_Jacobi *jac = (PC_Jacobi *)pc->data;
342: PetscBool flg;
343: PCJacobiType deflt, type;
345: PCJacobiGetType(pc, &deflt);
346: PetscOptionsHeadBegin(PetscOptionsObject, "Jacobi options");
347: PetscOptionsEnum("-pc_jacobi_type", "How to construct diagonal matrix", "PCJacobiSetType", PCJacobiTypes, (PetscEnum)deflt, (PetscEnum *)&type, &flg);
348: if (flg) PCJacobiSetType(pc, type);
349: PetscOptionsBool("-pc_jacobi_abs", "Use absolute values of diagonal entries", "PCJacobiSetUseAbs", jac->useabs, &jac->useabs, NULL);
350: PetscOptionsBool("-pc_jacobi_fixdiagonal", "Fix null terms on diagonal", "PCJacobiSetFixDiagonal", jac->fixdiag, &jac->fixdiag, NULL);
351: PetscOptionsHeadEnd();
352: return 0;
353: }
355: static PetscErrorCode PCView_Jacobi(PC pc, PetscViewer viewer)
356: {
357: PC_Jacobi *jac = (PC_Jacobi *)pc->data;
358: PetscBool iascii;
360: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
361: if (iascii) {
362: PCJacobiType type;
363: PetscBool useAbs, fixdiag;
364: PetscViewerFormat format;
366: PCJacobiGetType(pc, &type);
367: PCJacobiGetUseAbs(pc, &useAbs);
368: PCJacobiGetFixDiagonal(pc, &fixdiag);
369: PetscViewerASCIIPrintf(viewer, " type %s%s%s\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "", !fixdiag ? ", not checking null diagonal entries" : "");
370: PetscViewerGetFormat(viewer, &format);
371: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL && jac->diag) VecView(jac->diag, viewer);
372: }
373: return 0;
374: }
376: /*
377: PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
378: and sets this as the private data within the generic preconditioning
379: context, PC, that was created within PCCreate().
381: Input Parameter:
382: . pc - the preconditioner context
384: Application Interface Routine: PCCreate()
385: */
387: /*MC
388: PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)
390: Options Database Keys:
391: + -pc_jacobi_type <diagonal,rowmax,rowsum> - approach for forming the preconditioner
392: . -pc_jacobi_abs - use the absolute value of the diagonal entry
393: - -pc_jacobi_fixdiag - fix for zero diagonal terms by placing 1.0 in those locations
395: Level: beginner
397: Notes:
398: By using `KSPSetPCSide`(ksp,`PC_SYMMETRIC`) or -ksp_pc_side symmetric
399: can scale each side of the matrix by the square root of the diagonal entries.
401: Zero entries along the diagonal are replaced with the value 1.0
403: See `PCPBJACOBI` for fixed-size point block, `PCVPBJACOBI` for variable-sized point block, and `PCBJACOBI` for large size blocks
405: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
406: `PCJacobiSetType()`, `PCJacobiSetUseAbs()`, `PCJacobiGetUseAbs()`, `PCASM`,
407: `PCJacobiSetFixDiagonal()`, `PCJacobiGetFixDiagonal()`
408: `PCJacobiSetType()`, `PCJacobiSetUseAbs()`, `PCJacobiGetUseAbs()`, `PCPBJACOBI`, `PCBJACOBI`, `PCVPBJACOBI`
409: M*/
411: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
412: {
413: PC_Jacobi *jac;
415: /*
416: Creates the private data structure for this preconditioner and
417: attach it to the PC object.
418: */
419: PetscNew(&jac);
420: pc->data = (void *)jac;
422: /*
423: Initialize the pointers to vectors to ZERO; these will be used to store
424: diagonal entries of the matrix for fast preconditioner application.
425: */
426: jac->diag = NULL;
427: jac->diagsqrt = NULL;
428: jac->userowmax = PETSC_FALSE;
429: jac->userowsum = PETSC_FALSE;
430: jac->useabs = PETSC_FALSE;
431: jac->fixdiag = PETSC_TRUE;
433: /*
434: Set the pointers for the functions that are provided above.
435: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
436: are called, they will automatically call these functions. Note we
437: choose not to provide a couple of these functions since they are
438: not needed.
439: */
440: pc->ops->apply = PCApply_Jacobi;
441: pc->ops->applytranspose = PCApply_Jacobi;
442: pc->ops->setup = PCSetUp_Jacobi;
443: pc->ops->reset = PCReset_Jacobi;
444: pc->ops->destroy = PCDestroy_Jacobi;
445: pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
446: pc->ops->view = PCView_Jacobi;
447: pc->ops->applyrichardson = NULL;
448: pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
449: pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
451: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetType_C", PCJacobiSetType_Jacobi);
452: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetType_C", PCJacobiGetType_Jacobi);
453: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetUseAbs_C", PCJacobiSetUseAbs_Jacobi);
454: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetUseAbs_C", PCJacobiGetUseAbs_Jacobi);
455: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetFixDiagonal_C", PCJacobiSetFixDiagonal_Jacobi);
456: PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetFixDiagonal_C", PCJacobiGetFixDiagonal_Jacobi);
457: return 0;
458: }
460: /*@
461: PCJacobiSetUseAbs - Causes the Jacobi preconditioner `PCJACOBI` to use the
462: absolute values of the diagonal divisors in the preconditioner
464: Logically Collective
466: Input Parameters:
467: + pc - the preconditioner context
468: - flg - whether to use absolute values or not
470: Options Database Key:
471: . -pc_jacobi_abs <bool> - use absolute values
473: Note:
474: This takes affect at the next construction of the preconditioner
476: Level: intermediate
478: .seealso: `PCJACOBI`, `PCJacobiaSetType()`, `PCJacobiGetUseAbs()`
479: @*/
480: PetscErrorCode PCJacobiSetUseAbs(PC pc, PetscBool flg)
481: {
483: PetscTryMethod(pc, "PCJacobiSetUseAbs_C", (PC, PetscBool), (pc, flg));
484: return 0;
485: }
487: /*@
488: PCJacobiGetUseAbs - Determines if the Jacobi preconditioner `PCJACOBI` uses the
489: absolute values of the diagonal divisors in the preconditioner
491: Logically Collective
493: Input Parameter:
494: . pc - the preconditioner context
496: Output Parameter:
497: . flg - whether to use absolute values or not
499: Level: intermediate
501: .seealso: `PCJACOBI`, `PCJacobiaSetType()`, `PCJacobiSetUseAbs()`, `PCJacobiGetType()`
502: @*/
503: PetscErrorCode PCJacobiGetUseAbs(PC pc, PetscBool *flg)
504: {
506: PetscUseMethod(pc, "PCJacobiGetUseAbs_C", (PC, PetscBool *), (pc, flg));
507: return 0;
508: }
510: /*@
511: PCJacobiSetFixDiagonal - Check for zero values on the diagonal and replace them with 1.0
513: Logically Collective
515: Input Parameters:
516: + pc - the preconditioner context
517: - flg - the boolean flag
519: Options Database Key:
520: . -pc_jacobi_fixdiagonal <bool> - check for zero values on the diagonal
522: Note:
523: This takes affect at the next construction of the preconditioner
525: Level: intermediate
527: .seealso: `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiGetFixDiagonal()`, `PCJacobiSetUseAbs()`
528: @*/
529: PetscErrorCode PCJacobiSetFixDiagonal(PC pc, PetscBool flg)
530: {
532: PetscTryMethod(pc, "PCJacobiSetFixDiagonal_C", (PC, PetscBool), (pc, flg));
533: return 0;
534: }
536: /*@
537: PCJacobiGetFixDiagonal - Determines if the Jacobi preconditioner `PCJACOBI` checks for zero diagonal terms
539: Logically Collective
541: Input Parameter:
542: . pc - the preconditioner context
544: Output Parameter:
545: . flg - the boolean flag
547: Options Database Key:
548: . -pc_jacobi_fixdiagonal <bool> - Fix 0 terms on diagonal by using 1
550: Level: intermediate
552: .seealso: `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiSetFixDiagonal()`
553: @*/
554: PetscErrorCode PCJacobiGetFixDiagonal(PC pc, PetscBool *flg)
555: {
557: PetscUseMethod(pc, "PCJacobiGetFixDiagonal_C", (PC, PetscBool *), (pc, flg));
558: return 0;
559: }
561: /*@
562: PCJacobiSetType - Causes the Jacobi preconditioner to use either the diagonal, the maximum entry in each row,
563: of the sum of rows entries for the diagonal preconditioner
565: Logically Collective
567: Input Parameters:
568: + pc - the preconditioner context
569: - type - `PC_JACOBI_DIAGONAL`, `PC_JACOBI_ROWMAX`, `PC_JACOBI_ROWSUM`
571: Options Database Key:
572: . -pc_jacobi_type <diagonal,rowmax,rowsum> - the type of diagonal matrix to use for Jacobi
574: Level: intermediate
576: Developer Note:
577: Why is there a separate function for using the absolute value?
579: .seealso: `PCJACOBI`, `PCJacobiSetUseAbs()`, `PCJacobiGetType()`
580: @*/
581: PetscErrorCode PCJacobiSetType(PC pc, PCJacobiType type)
582: {
584: PetscTryMethod(pc, "PCJacobiSetType_C", (PC, PCJacobiType), (pc, type));
585: return 0;
586: }
588: /*@
589: PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner
591: Not Collective
593: Input Parameter:
594: . pc - the preconditioner context
596: Output Parameter:
597: . type - `PC_JACOBI_DIAGONAL`, `PC_JACOBI_ROWMAX`, `PC_JACOBI_ROWSUM`
599: Level: intermediate
601: .seealso: `PCJACOBI`, `PCJacobiaUseAbs()`, `PCJacobiSetType()`
602: @*/
603: PetscErrorCode PCJacobiGetType(PC pc, PCJacobiType *type)
604: {
606: PetscUseMethod(pc, "PCJacobiGetType_C", (PC, PCJacobiType *), (pc, type));
607: return 0;
608: }