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: } PC_Jacobi;
68: static PetscErrorCode PCJacobiSetType_Jacobi(PC pc,PCJacobiType type)
69: {
70: PC_Jacobi *j = (PC_Jacobi*)pc->data;
73: j->userowmax = PETSC_FALSE;
74: j->userowsum = PETSC_FALSE;
75: if (type == PC_JACOBI_ROWMAX) {
76: j->userowmax = PETSC_TRUE;
77: } else if (type == PC_JACOBI_ROWSUM) {
78: j->userowsum = PETSC_TRUE;
79: }
80: return(0);
81: }
83: static PetscErrorCode PCJacobiGetType_Jacobi(PC pc,PCJacobiType *type)
84: {
85: PC_Jacobi *j = (PC_Jacobi*)pc->data;
88: if (j->userowmax) {
89: *type = PC_JACOBI_ROWMAX;
90: } else if (j->userowsum) {
91: *type = PC_JACOBI_ROWSUM;
92: } else {
93: *type = PC_JACOBI_DIAGONAL;
94: }
95: return(0);
96: }
98: static PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc,PetscBool flg)
99: {
100: PC_Jacobi *j = (PC_Jacobi*)pc->data;
103: j->useabs = flg;
104: return(0);
105: }
107: static PetscErrorCode PCJacobiGetUseAbs_Jacobi(PC pc,PetscBool *flg)
108: {
109: PC_Jacobi *j = (PC_Jacobi*)pc->data;
112: *flg = j->useabs;
113: return(0);
114: }
116: /* -------------------------------------------------------------------------- */
117: /*
118: PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
119: by setting data structures and options.
121: Input Parameter:
122: . pc - the preconditioner context
124: Application Interface Routine: PCSetUp()
126: Notes:
127: The interface routine PCSetUp() is not usually called directly by
128: the user, but instead is called by PCApply() if necessary.
129: */
130: static PetscErrorCode PCSetUp_Jacobi(PC pc)
131: {
132: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
133: Vec diag,diagsqrt;
135: PetscInt n,i;
136: PetscScalar *x;
137: PetscBool zeroflag = PETSC_FALSE;
140: /*
141: For most preconditioners the code would begin here something like
143: if (pc->setupcalled == 0) { allocate space the first time this is ever called
144: MatCreateVecs(pc->mat,&jac->diag);
145: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
146: }
148: But for this preconditioner we want to support use of both the matrix' diagonal
149: elements (for left or right preconditioning) and square root of diagonal elements
150: (for symmetric preconditioning). Hence we do not allocate space here, since we
151: don't know at this point which will be needed (diag and/or diagsqrt) until the user
152: applies the preconditioner, and we don't want to allocate BOTH unless we need
153: them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
154: and PCSetUp_Jacobi_Symmetric(), respectively.
155: */
157: /*
158: Here we set up the preconditioner; that is, we copy the diagonal values from
159: the matrix and put them into a format to make them quick to apply as a preconditioner.
160: */
161: diag = jac->diag;
162: diagsqrt = jac->diagsqrt;
164: if (diag) {
165: if (jac->userowmax) {
166: MatGetRowMaxAbs(pc->pmat,diag,NULL);
167: } else if (jac->userowsum) {
168: MatGetRowSum(pc->pmat,diag);
169: } else {
170: MatGetDiagonal(pc->pmat,diag);
171: }
172: VecReciprocal(diag);
173: VecGetLocalSize(diag,&n);
174: if (jac->useabs) {
175: VecAbs(diag);
176: }
177: VecGetArray(diag,&x);
178: for (i=0; i<n; i++) {
179: if (x[i] == 0.0) {
180: x[i] = 1.0;
181: zeroflag = PETSC_TRUE;
182: }
183: }
184: VecRestoreArray(diag,&x);
185: }
186: if (diagsqrt) {
187: if (jac->userowmax) {
188: MatGetRowMaxAbs(pc->pmat,diagsqrt,NULL);
189: } else if (jac->userowsum) {
190: MatGetRowSum(pc->pmat,diagsqrt);
191: } else {
192: MatGetDiagonal(pc->pmat,diagsqrt);
193: }
194: VecGetLocalSize(diagsqrt,&n);
195: VecGetArray(diagsqrt,&x);
196: for (i=0; i<n; i++) {
197: if (x[i] != 0.0) x[i] = 1.0/PetscSqrtReal(PetscAbsScalar(x[i]));
198: else {
199: x[i] = 1.0;
200: zeroflag = PETSC_TRUE;
201: }
202: }
203: VecRestoreArray(diagsqrt,&x);
204: }
205: if (zeroflag) {
206: PetscInfo(pc,"Zero detected in diagonal of matrix, using 1 at those locations\n");
207: }
208: return(0);
209: }
210: /* -------------------------------------------------------------------------- */
211: /*
212: PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
213: inverse of the square root of the diagonal entries of the matrix. This
214: is used for symmetric application of the Jacobi preconditioner.
216: Input Parameter:
217: . pc - the preconditioner context
218: */
219: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
220: {
222: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
225: MatCreateVecs(pc->pmat,&jac->diagsqrt,NULL);
226: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diagsqrt);
227: PCSetUp_Jacobi(pc);
228: return(0);
229: }
230: /* -------------------------------------------------------------------------- */
231: /*
232: PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
233: inverse of the diagonal entries of the matrix. This is used for left of
234: right application of the Jacobi preconditioner.
236: Input Parameter:
237: . pc - the preconditioner context
238: */
239: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
240: {
242: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
245: MatCreateVecs(pc->pmat,&jac->diag,NULL);
246: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
247: PCSetUp_Jacobi(pc);
248: return(0);
249: }
250: /* -------------------------------------------------------------------------- */
251: /*
252: PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
254: Input Parameters:
255: . pc - the preconditioner context
256: . x - input vector
258: Output Parameter:
259: . y - output vector
261: Application Interface Routine: PCApply()
262: */
263: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
264: {
265: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
269: if (!jac->diag) {
270: PCSetUp_Jacobi_NonSymmetric(pc);
271: }
272: VecPointwiseMult(y,x,jac->diag);
273: return(0);
274: }
275: /* -------------------------------------------------------------------------- */
276: /*
277: PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
278: symmetric preconditioner to a vector.
280: Input Parameters:
281: . pc - the preconditioner context
282: . x - input vector
284: Output Parameter:
285: . y - output vector
287: Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
288: */
289: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
290: {
292: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
295: if (!jac->diagsqrt) {
296: PCSetUp_Jacobi_Symmetric(pc);
297: }
298: VecPointwiseMult(y,x,jac->diagsqrt);
299: return(0);
300: }
301: /* -------------------------------------------------------------------------- */
302: static PetscErrorCode PCReset_Jacobi(PC pc)
303: {
304: 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: {
327: PCReset_Jacobi(pc);
329: /*
330: Free the private data structure that was hanging off the PC
331: */
332: PetscFree(pc->data);
333: return(0);
334: }
336: static PetscErrorCode PCSetFromOptions_Jacobi(PetscOptionItems *PetscOptionsObject,PC pc)
337: {
338: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
340: PetscBool flg;
341: PCJacobiType deflt,type;
344: PCJacobiGetType(pc,&deflt);
345: PetscOptionsHead(PetscOptionsObject,"Jacobi options");
346: PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCJacobiSetType",PCJacobiTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);
347: if (flg) {
348: PCJacobiSetType(pc,type);
349: }
350: PetscOptionsBool("-pc_jacobi_abs","Use absolute values of diagaonal entries","PCJacobiSetUseAbs",jac->useabs,&jac->useabs,NULL);
351: PetscOptionsTail();
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;
362: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
363: if (iascii) {
364: PCJacobiType type;
365: PetscBool useAbs;
366: PetscViewerFormat format;
368: PCJacobiGetType(pc, &type);
369: PCJacobiGetUseAbs(pc, &useAbs);
370: PetscViewerASCIIPrintf(viewer, " type %s%s\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "");
371: PetscViewerGetFormat(viewer, &format);
372: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
373: VecView(jac->diag, viewer);
374: }
375: }
376: return(0);
377: }
379: /* -------------------------------------------------------------------------- */
380: /*
381: PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
382: and sets this as the private data within the generic preconditioning
383: context, PC, that was created within PCCreate().
385: Input Parameter:
386: . pc - the preconditioner context
388: Application Interface Routine: PCCreate()
389: */
391: /*MC
392: PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)
394: Options Database Key:
395: + -pc_jacobi_type <diagonal,rowmax,rowsum> - approach for forming the preconditioner
396: - -pc_jacobi_abs - use the absolute value of the diagonal entry
398: Level: beginner
400: Notes:
401: By using KSPSetPCSide(ksp,PC_SYMMETRIC) or -ksp_pc_side symmetric
402: can scale each side of the matrix by the square root of the diagonal entries.
404: Zero entries along the diagonal are replaced with the value 1.0
406: See PCPBJACOBI for a point-block Jacobi preconditioner
408: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
409: PCJacobiSetType(), PCJacobiSetUseAbs(), PCJacobiGetUseAbs(), PCPBJACOBI
410: M*/
412: PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
413: {
414: PC_Jacobi *jac;
418: /*
419: Creates the private data structure for this preconditioner and
420: attach it to the PC object.
421: */
422: PetscNewLog(pc,&jac);
423: pc->data = (void*)jac;
425: /*
426: Initialize the pointers to vectors to ZERO; these will be used to store
427: diagonal entries of the matrix for fast preconditioner application.
428: */
429: jac->diag = NULL;
430: jac->diagsqrt = NULL;
431: jac->userowmax = PETSC_FALSE;
432: jac->userowsum = PETSC_FALSE;
433: jac->useabs = PETSC_FALSE;
435: /*
436: Set the pointers for the functions that are provided above.
437: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
438: are called, they will automatically call these functions. Note we
439: choose not to provide a couple of these functions since they are
440: not needed.
441: */
442: pc->ops->apply = PCApply_Jacobi;
443: pc->ops->applytranspose = PCApply_Jacobi;
444: pc->ops->setup = PCSetUp_Jacobi;
445: pc->ops->reset = PCReset_Jacobi;
446: pc->ops->destroy = PCDestroy_Jacobi;
447: pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
448: pc->ops->view = PCView_Jacobi;
449: pc->ops->applyrichardson = NULL;
450: pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
451: pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
453: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",PCJacobiSetType_Jacobi);
454: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",PCJacobiGetType_Jacobi);
455: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",PCJacobiSetUseAbs_Jacobi);
456: PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",PCJacobiGetUseAbs_Jacobi);
457: return(0);
458: }
460: /*@
461: PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
462: absolute values of the digonal divisors in the preconditioner
464: Logically Collective on PC
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
473: Notes:
474: This takes affect at the next construction of the preconditioner
476: Level: intermediate
478: .seealso: PCJacobiaSetType(), PCJacobiGetUseAbs()
480: @*/
481: PetscErrorCode PCJacobiSetUseAbs(PC pc,PetscBool flg)
482: {
487: PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC,PetscBool),(pc,flg));
488: return(0);
489: }
491: /*@
492: PCJacobiGetUseAbs - Determines if the Jacobi preconditioner uses the
493: absolute values of the digonal divisors in the preconditioner
495: Logically Collective on PC
497: Input Parameter:
498: . pc - the preconditioner context
500: Output Parameter:
501: . flg - whether to use absolute values or not
503: Options Database Key:
504: . -pc_jacobi_abs
506: Level: intermediate
508: .seealso: PCJacobiaSetType(), PCJacobiSetUseAbs(), PCJacobiGetType()
510: @*/
511: PetscErrorCode PCJacobiGetUseAbs(PC pc,PetscBool *flg)
512: {
517: PetscUseMethod(pc,"PCJacobiGetUseAbs_C",(PC,PetscBool*),(pc,flg));
518: return(0);
519: }
521: /*@
522: PCJacobiSetType - Causes the Jacobi preconditioner to use either the diagonal, the maximum entry in each row,
523: of the sum of rows entries for the diagonal preconditioner
525: Logically Collective on PC
527: Input Parameters:
528: + pc - the preconditioner context
529: - type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM
531: Options Database Key:
532: . -pc_jacobi_type <diagonal,rowmax,rowsum>
534: Level: intermediate
536: .seealso: PCJacobiaUseAbs(), PCJacobiGetType()
537: @*/
538: PetscErrorCode PCJacobiSetType(PC pc,PCJacobiType type)
539: {
544: PetscTryMethod(pc,"PCJacobiSetType_C",(PC,PCJacobiType),(pc,type));
545: return(0);
546: }
548: /*@
549: PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner
551: Not Collective on PC
553: Input Parameter:
554: . pc - the preconditioner context
556: Output Parameter:
557: . type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM
559: Level: intermediate
561: .seealso: PCJacobiaUseAbs(), PCJacobiSetType()
562: @*/
563: PetscErrorCode PCJacobiGetType(PC pc,PCJacobiType *type)
564: {
569: PetscUseMethod(pc,"PCJacobiGetType_C",(PC,PCJacobiType*),(pc,type));
570: return(0);
571: }