Actual source code: jacobi.c
1: #define PETSCKSP_DLL
3: /* --------------------------------------------------------------------
5: This file implements a Jacobi preconditioner in PETSc as part of PC.
6: You can use this as a starting point for implementing your own
7: preconditioner that is not provided with PETSc. (You might also consider
8: just using PCSHELL)
10: The following basic routines are required for each preconditioner.
11: PCCreate_XXX() - Creates a preconditioner context
12: PCSetFromOptions_XXX() - Sets runtime options
13: PCApply_XXX() - Applies the preconditioner
14: PCDestroy_XXX() - Destroys the preconditioner context
15: where the suffix "_XXX" denotes a particular implementation, in
16: this case we use _Jacobi (e.g., PCCreate_Jacobi, PCApply_Jacobi).
17: These routines are actually called via the common user interface
18: routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(),
19: so the application code interface remains identical for all
20: preconditioners.
22: Another key routine is:
23: PCSetUp_XXX() - Prepares for the use of a preconditioner
24: by setting data structures and options. The interface routine PCSetUp()
25: is not usually called directly by the user, but instead is called by
26: PCApply() if necessary.
28: Additional basic routines are:
29: PCView_XXX() - Prints details of runtime options that
30: have actually been used.
31: These are called by application codes via the interface routines
32: PCView().
34: The various types of solvers (preconditioners, Krylov subspace methods,
35: nonlinear solvers, timesteppers) are all organized similarly, so the
36: above description applies to these categories also. One exception is
37: that the analogues of PCApply() for these components are KSPSolve(),
38: SNESSolve(), and TSSolve().
40: Additional optional functionality unique to preconditioners is left and
41: right symmetric preconditioner application via PCApplySymmetricLeft()
42: and PCApplySymmetricRight(). The Jacobi implementation is
43: PCApplySymmetricLeftOrRight_Jacobi().
45: -------------------------------------------------------------------- */
47: /*
48: Include files needed for the Jacobi preconditioner:
49: pcimpl.h - private include file intended for use by all preconditioners
50: */
52: #include private/pcimpl.h
54: /*
55: Private context (data structure) for the Jacobi preconditioner.
56: */
57: typedef struct {
58: Vec diag; /* vector containing the reciprocals of the diagonal elements
59: 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: PetscTruth userowmax;
64: PetscTruth useabs; /* use the absolute values of the diagonal entries */
65: } PC_Jacobi;
70: PetscErrorCode PCJacobiSetUseRowMax_Jacobi(PC pc)
71: {
72: PC_Jacobi *j;
75: j = (PC_Jacobi*)pc->data;
76: j->userowmax = PETSC_TRUE;
77: return(0);
78: }
84: PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc)
85: {
86: PC_Jacobi *j;
89: j = (PC_Jacobi*)pc->data;
90: j->useabs = PETSC_TRUE;
91: return(0);
92: }
95: /* -------------------------------------------------------------------------- */
96: /*
97: PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
98: by setting data structures and options.
100: Input Parameter:
101: . pc - the preconditioner context
103: Application Interface Routine: PCSetUp()
105: Notes:
106: The interface routine PCSetUp() is not usually called directly by
107: the user, but instead is called by PCApply() if necessary.
108: */
111: static PetscErrorCode PCSetUp_Jacobi(PC pc)
112: {
113: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
114: Vec diag,diagsqrt;
116: PetscInt n,i;
117: PetscScalar *x;
118: PetscTruth zeroflag = PETSC_FALSE;
121: /*
122: For most preconditioners the code would begin here something like
124: if (pc->setupcalled == 0) { allocate space the first time this is ever called
125: MatGetVecs(pc->mat,&jac->diag);
126: PetscLogObjectParent(pc,jac->diag);
127: }
129: But for this preconditioner we want to support use of both the matrix' diagonal
130: elements (for left or right preconditioning) and square root of diagonal elements
131: (for symmetric preconditioning). Hence we do not allocate space here, since we
132: don't know at this point which will be needed (diag and/or diagsqrt) until the user
133: applies the preconditioner, and we don't want to allocate BOTH unless we need
134: them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
135: and PCSetUp_Jacobi_Symmetric(), respectively.
136: */
138: /*
139: Here we set up the preconditioner; that is, we copy the diagonal values from
140: the matrix and put them into a format to make them quick to apply as a preconditioner.
141: */
142: diag = jac->diag;
143: diagsqrt = jac->diagsqrt;
145: if (diag) {
146: if (jac->userowmax) {
147: MatGetRowMax(pc->pmat,diag);
148: } else {
149: MatGetDiagonal(pc->pmat,diag);
150: }
151: VecReciprocal(diag);
152: VecGetLocalSize(diag,&n);
153: VecGetArray(diag,&x);
154: if (jac->useabs) {
155: for (i=0; i<n; i++) {
156: x[i] = PetscAbsScalar(x[i]);
157: }
158: }
159: for (i=0; i<n; i++) {
160: if (x[i] == 0.0) {
161: x[i] = 1.0;
162: zeroflag = PETSC_TRUE;
163: }
164: }
165: VecRestoreArray(diag,&x);
166: }
167: if (diagsqrt) {
168: if (jac->userowmax) {
169: MatGetRowMax(pc->pmat,diagsqrt);
170: } else {
171: MatGetDiagonal(pc->pmat,diagsqrt);
172: }
173: VecGetLocalSize(diagsqrt,&n);
174: VecGetArray(diagsqrt,&x);
175: for (i=0; i<n; i++) {
176: if (x[i] != 0.0) x[i] = 1.0/sqrt(PetscAbsScalar(x[i]));
177: else {
178: x[i] = 1.0;
179: zeroflag = PETSC_TRUE;
180: }
181: }
182: VecRestoreArray(diagsqrt,&x);
183: }
184: if (zeroflag) {
185: PetscInfo(pc,"Zero detected in diagonal of matrix, using 1 at those locations\n");
186: }
187: return(0);
188: }
189: /* -------------------------------------------------------------------------- */
190: /*
191: PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
192: inverse of the square root of the diagonal entries of the matrix. This
193: is used for symmetric application of the Jacobi preconditioner.
195: Input Parameter:
196: . pc - the preconditioner context
197: */
200: static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
201: {
203: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
206: MatGetVecs(pc->pmat,&jac->diagsqrt,0);
207: PetscLogObjectParent(pc,jac->diagsqrt);
208: PCSetUp_Jacobi(pc);
209: return(0);
210: }
211: /* -------------------------------------------------------------------------- */
212: /*
213: PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
214: inverse of the diagonal entries of the matrix. This is used for left of
215: right application of the Jacobi preconditioner.
217: Input Parameter:
218: . pc - the preconditioner context
219: */
222: static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
223: {
225: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
228: MatGetVecs(pc->pmat,&jac->diag,0);
229: PetscLogObjectParent(pc,jac->diag);
230: PCSetUp_Jacobi(pc);
231: return(0);
232: }
233: /* -------------------------------------------------------------------------- */
234: /*
235: PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
237: Input Parameters:
238: . pc - the preconditioner context
239: . x - input vector
241: Output Parameter:
242: . y - output vector
244: Application Interface Routine: PCApply()
245: */
248: static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
249: {
250: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
254: if (!jac->diag) {
255: PCSetUp_Jacobi_NonSymmetric(pc);
256: }
257: VecPointwiseMult(y,x,jac->diag);
258: return(0);
259: }
260: /* -------------------------------------------------------------------------- */
261: /*
262: PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
263: symmetric preconditioner to a vector.
265: Input Parameters:
266: . pc - the preconditioner context
267: . x - input vector
269: Output Parameter:
270: . y - output vector
272: Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
273: */
276: static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
277: {
279: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
282: if (!jac->diagsqrt) {
283: PCSetUp_Jacobi_Symmetric(pc);
284: }
285: VecPointwiseMult(y,x,jac->diagsqrt);
286: return(0);
287: }
288: /* -------------------------------------------------------------------------- */
289: /*
290: PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
291: that was created with PCCreate_Jacobi().
293: Input Parameter:
294: . pc - the preconditioner context
296: Application Interface Routine: PCDestroy()
297: */
300: static PetscErrorCode PCDestroy_Jacobi(PC pc)
301: {
302: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
306: if (jac->diag) {VecDestroy(jac->diag);}
307: if (jac->diagsqrt) {VecDestroy(jac->diagsqrt);}
309: /*
310: Free the private data structure that was hanging off the PC
311: */
312: PetscFree(jac);
313: return(0);
314: }
318: static PetscErrorCode PCSetFromOptions_Jacobi(PC pc)
319: {
320: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
324: PetscOptionsHead("Jacobi options");
325: PetscOptionsTruth("-pc_jacobi_rowmax","Use row maximums for diagonal","PCJacobiSetUseRowMax",jac->userowmax,
326: &jac->userowmax,PETSC_NULL);
327: PetscOptionsTruth("-pc_jacobi_abs","Use absolute values of diagaonal entries","PCJacobiSetUseAbs",jac->useabs,
328: &jac->useabs,PETSC_NULL);
329: PetscOptionsTail();
330: return(0);
331: }
333: /* -------------------------------------------------------------------------- */
334: /*
335: PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
336: and sets this as the private data within the generic preconditioning
337: context, PC, that was created within PCCreate().
339: Input Parameter:
340: . pc - the preconditioner context
342: Application Interface Routine: PCCreate()
343: */
345: /*MC
346: PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)
348: Options Database Key:
349: + -pc_jacobi_rowmax - use the maximum absolute value in each row as the scaling factor,
350: rather than the diagonal
351: - -pc_jacobi_abs - use the absolute value of the diagaonl entry
353: Level: beginner
355: Concepts: Jacobi, diagonal scaling, preconditioners
357: Notes: By using KSPSetPreconditionerSide(ksp,PC_SYMMETRIC) or -ksp_symmetric_pc you
358: can scale each side of the matrix by the squareroot of the diagonal entries.
360: Zero entries along the diagonal are replaced with the value 1.0
362: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
363: PCJacobiSetUseRowMax(), PCJacobiSetUseAbs()
364: M*/
369: PetscErrorCode PCCreate_Jacobi(PC pc)
370: {
371: PC_Jacobi *jac;
375: /*
376: Creates the private data structure for this preconditioner and
377: attach it to the PC object.
378: */
379: PetscNew(PC_Jacobi,&jac);
380: pc->data = (void*)jac;
382: /*
383: Logs the memory usage; this is not needed but allows PETSc to
384: monitor how much memory is being used for various purposes.
385: */
386: PetscLogObjectMemory(pc,sizeof(PC_Jacobi));
388: /*
389: Initialize the pointers to vectors to ZERO; these will be used to store
390: diagonal entries of the matrix for fast preconditioner application.
391: */
392: jac->diag = 0;
393: jac->diagsqrt = 0;
394: jac->userowmax = PETSC_FALSE;
395: jac->useabs = PETSC_FALSE;
397: /*
398: Set the pointers for the functions that are provided above.
399: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
400: are called, they will automatically call these functions. Note we
401: choose not to provide a couple of these functions since they are
402: not needed.
403: */
404: pc->ops->apply = PCApply_Jacobi;
405: pc->ops->applytranspose = PCApply_Jacobi;
406: pc->ops->setup = PCSetUp_Jacobi;
407: pc->ops->destroy = PCDestroy_Jacobi;
408: pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
409: pc->ops->view = 0;
410: pc->ops->applyrichardson = 0;
411: pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
412: pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
413: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseRowMax_C","PCJacobiSetUseRowMax_Jacobi",PCJacobiSetUseRowMax_Jacobi);
414: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseAbs_C","PCJacobiSetUseAbs_Jacobi",PCJacobiSetUseAbs_Jacobi);
415: return(0);
416: }
422: /*@
423: PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
424: absolute value of the diagonal to for the preconditioner
426: Collective on PC
428: Input Parameters:
429: . pc - the preconditioner context
432: Options Database Key:
433: . -pc_jacobi_abs
435: Level: intermediate
437: Concepts: Jacobi preconditioner
439: .seealso: PCJacobiaUseRowMax()
441: @*/
442: PetscErrorCode PCJacobiSetUseAbs(PC pc)
443: {
444: PetscErrorCode ierr,(*f)(PC);
448: PetscObjectQueryFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",(void (**)(void))&f);
449: if (f) {
450: (*f)(pc);
451: }
452: return(0);
453: }
457: /*@
458: PCJacobiSetUseRowMax - Causes the Jacobi preconditioner to use the
459: maximum entry in each row as the diagonal preconditioner, instead of
460: the diagonal entry
462: Collective on PC
464: Input Parameters:
465: . pc - the preconditioner context
468: Options Database Key:
469: . -pc_jacobi_rowmax
471: Level: intermediate
473: Concepts: Jacobi preconditioner
475: .seealso: PCJacobiaUseAbs()
476: @*/
477: PetscErrorCode PCJacobiSetUseRowMax(PC pc)
478: {
479: PetscErrorCode ierr,(*f)(PC);
483: PetscObjectQueryFunction((PetscObject)pc,"PCJacobiSetUseRowMax_C",(void (**)(void))&f);
484: if (f) {
485: (*f)(pc);
486: }
487: return(0);
488: }