Actual source code: jacobi.c
1: /*$Id: jacobi.c,v 1.75 2001/08/07 03:03:32 balay Exp $*/
3: /* --------------------------------------------------------------------
5: This file implements a Jacobi preconditioner for matrices that use
6: the Mat interface (various matrix formats). Actually, the only
7: matrix operation used here is MatGetDiagonal(), which extracts
8: diagonal elements of the preconditioning matrix.
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 src/sles/pc/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: } PC_Jacobi;
66: EXTERN_C_BEGIN
67: #undef __FUNCT__
69: int PCJacobiSetUseRowMax_Jacobi(PC pc)
70: {
71: PC_Jacobi *j;
74: j = (PC_Jacobi*)pc->data;
75: j->userowmax = PETSC_TRUE;
76: return(0);
77: }
78: EXTERN_C_END
80: /* -------------------------------------------------------------------------- */
81: /*
82: PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
83: by setting data structures and options.
85: Input Parameter:
86: . pc - the preconditioner context
88: Application Interface Routine: PCSetUp()
90: Notes:
91: The interface routine PCSetUp() is not usually called directly by
92: the user, but instead is called by PCApply() if necessary.
93: */
94: #undef __FUNCT__
96: static int PCSetUp_Jacobi(PC pc)
97: {
98: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
99: Vec diag,diagsqrt;
100: int ierr,n,i,zeroflag=0;
101: PetscScalar *x;
105: /*
106: For most preconditioners the code would begin here something like
108: if (pc->setupcalled == 0) { allocate space the first time this is ever called
109: VecDuplicate(pc->vec,&jac->diag);
110: PetscLogObjectParent(pc,jac->diag);
111: }
113: But for this preconditioner we want to support use of both the matrix' diagonal
114: elements (for left or right preconditioning) and square root of diagonal elements
115: (for symmetric preconditioning). Hence we do not allocate space here, since we
116: don't know at this point which will be needed (diag and/or diagsqrt) until the user
117: applies the preconditioner, and we don't want to allocate BOTH unless we need
118: them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
119: and PCSetUp_Jacobi_Symmetric(), respectively.
120: */
122: /*
123: Here we set up the preconditioner; that is, we copy the diagonal values from
124: the matrix and put them into a format to make them quick to apply as a preconditioner.
125: */
126: diag = jac->diag;
127: diagsqrt = jac->diagsqrt;
129: if (diag) {
130: if (jac->userowmax) {
131: MatGetRowMax(pc->pmat,diag);
132: } else {
133: MatGetDiagonal(pc->pmat,diag);
134: }
135: VecReciprocal(diag);
136: VecGetLocalSize(diag,&n);
137: VecGetArray(diag,&x);
138: for (i=0; i<n; i++) {
139: if (x[i] == 0.0) {
140: x[i] = 1.0;
141: zeroflag = 1;
142: }
143: }
144: VecRestoreArray(diag,&x);
145: }
146: if (diagsqrt) {
147: if (jac->userowmax) {
148: MatGetRowMax(pc->pmat,diagsqrt);
149: } else {
150: MatGetDiagonal(pc->pmat,diagsqrt);
151: }
152: VecGetLocalSize(diagsqrt,&n);
153: VecGetArray(diagsqrt,&x);
154: for (i=0; i<n; i++) {
155: if (x[i] != 0.0) x[i] = 1.0/sqrt(PetscAbsScalar(x[i]));
156: else {
157: x[i] = 1.0;
158: zeroflag = 1;
159: }
160: }
161: VecRestoreArray(diagsqrt,&x);
162: }
163: if (zeroflag) {
164: PetscLogInfo(pc,"PCSetUp_Jacobi:Zero detected in diagonal of matrix, using 1 at those locationsn");
165: }
167: return(0);
168: }
169: /* -------------------------------------------------------------------------- */
170: /*
171: PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
172: inverse of the square root of the diagonal entries of the matrix. This
173: is used for symmetric application of the Jacobi preconditioner.
175: Input Parameter:
176: . pc - the preconditioner context
177: */
178: #undef __FUNCT__
180: static int PCSetUp_Jacobi_Symmetric(PC pc)
181: {
182: int ierr;
183: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
187: VecDuplicate(pc->vec,&jac->diagsqrt);
188: PetscLogObjectParent(pc,jac->diagsqrt);
189: PCSetUp_Jacobi(pc);
190: return(0);
191: }
192: /* -------------------------------------------------------------------------- */
193: /*
194: PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
195: inverse of the diagonal entries of the matrix. This is used for left of
196: right application of the Jacobi preconditioner.
198: Input Parameter:
199: . pc - the preconditioner context
200: */
201: #undef __FUNCT__
203: static int PCSetUp_Jacobi_NonSymmetric(PC pc)
204: {
205: int ierr;
206: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
210: VecDuplicate(pc->vec,&jac->diag);
211: PetscLogObjectParent(pc,jac->diag);
212: PCSetUp_Jacobi(pc);
213: return(0);
214: }
215: /* -------------------------------------------------------------------------- */
216: /*
217: PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
219: Input Parameters:
220: . pc - the preconditioner context
221: . x - input vector
223: Output Parameter:
224: . y - output vector
226: Application Interface Routine: PCApply()
227: */
228: #undef __FUNCT__
230: static int PCApply_Jacobi(PC pc,Vec x,Vec y)
231: {
232: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
233: int ierr;
236: if (!jac->diag) {
237: PCSetUp_Jacobi_NonSymmetric(pc);
238: }
239: VecPointwiseMult(x,jac->diag,y);
240: return(0);
241: }
242: /* -------------------------------------------------------------------------- */
243: /*
244: PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
245: symmetric preconditioner to a vector.
247: Input Parameters:
248: . pc - the preconditioner context
249: . x - input vector
251: Output Parameter:
252: . y - output vector
254: Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
255: */
256: #undef __FUNCT__
258: static int PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
259: {
260: int ierr;
261: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
264: if (!jac->diagsqrt) {
265: PCSetUp_Jacobi_Symmetric(pc);
266: }
267: VecPointwiseMult(x,jac->diagsqrt,y);
268: return(0);
269: }
270: /* -------------------------------------------------------------------------- */
271: /*
272: PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
273: that was created with PCCreate_Jacobi().
275: Input Parameter:
276: . pc - the preconditioner context
278: Application Interface Routine: PCDestroy()
279: */
280: #undef __FUNCT__
282: static int PCDestroy_Jacobi(PC pc)
283: {
284: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
285: int ierr;
288: if (jac->diag) {VecDestroy(jac->diag);}
289: if (jac->diagsqrt) {VecDestroy(jac->diagsqrt);}
291: /*
292: Free the private data structure that was hanging off the PC
293: */
294: PetscFree(jac);
295: return(0);
296: }
298: #undef __FUNCT__
300: static int PCSetFromOptions_Jacobi(PC pc)
301: {
302: PC_Jacobi *jac = (PC_Jacobi*)pc->data;
303: int ierr;
306: PetscOptionsHead("Jacobi options");
307: PetscOptionsLogical("-pc_jacobi_rowmax","Use row maximums for diagonal","PCJacobiSetUseRowMax",jac->userowmax,
308: &jac->userowmax,PETSC_NULL);
309: PetscOptionsTail();
310: return(0);
311: }
313: /* -------------------------------------------------------------------------- */
314: /*
315: PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
316: and sets this as the private data within the generic preconditioning
317: context, PC, that was created within PCCreate().
319: Input Parameter:
320: . pc - the preconditioner context
322: Application Interface Routine: PCCreate()
323: */
325: /*MC
326: PCJacobi - Jacobi (i.e. diagonal scaling preconditioning)
328: Options Database Key:
329: . -pc_jacobi_rowmax - use the maximum absolute value in each row as the scaling factor,
330: rather than the diagonal
332: Level: beginner
334: Concepts: Jacobi, diagonal scaling, preconditioners
336: Notes: By using KSPSetPreconditionerSide(ksp,PC_SYMMETRIC) or -ksp_symmetric_pc you
337: can scale each side of the matrix by the squareroot of the diagonal entries.
339: Zero entries along the diagonal are replaced with the value 1.0
341: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
342: PCJacobiSetUseRowMax(),
343: M*/
345: EXTERN_C_BEGIN
346: #undef __FUNCT__
348: int PCCreate_Jacobi(PC pc)
349: {
350: PC_Jacobi *jac;
351: int ierr;
355: /*
356: Creates the private data structure for this preconditioner and
357: attach it to the PC object.
358: */
359: ierr = PetscNew(PC_Jacobi,&jac);
360: pc->data = (void*)jac;
362: /*
363: Logs the memory usage; this is not needed but allows PETSc to
364: monitor how much memory is being used for various purposes.
365: */
366: PetscLogObjectMemory(pc,sizeof(PC_Jacobi));
368: /*
369: Initialize the pointers to vectors to ZERO; these will be used to store
370: diagonal entries of the matrix for fast preconditioner application.
371: */
372: jac->diag = 0;
373: jac->diagsqrt = 0;
374: jac->userowmax = PETSC_FALSE;
376: /*
377: Set the pointers for the functions that are provided above.
378: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
379: are called, they will automatically call these functions. Note we
380: choose not to provide a couple of these functions since they are
381: not needed.
382: */
383: pc->ops->apply = PCApply_Jacobi;
384: pc->ops->applytranspose = PCApply_Jacobi;
385: pc->ops->setup = PCSetUp_Jacobi;
386: pc->ops->destroy = PCDestroy_Jacobi;
387: pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
388: pc->ops->view = 0;
389: pc->ops->applyrichardson = 0;
390: pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
391: pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
392: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCJacobiSetUseRowMax_C","PCJacobiSetUseRowMax_Jacobi",
393: PCJacobiSetUseRowMax_Jacobi);
394: return(0);
395: }
396: EXTERN_C_END
398: #undef __FUNCT__
400: /*@
401: PCJacobiSetUseRowMax - Causes the Jacobi preconditioner to use the
402: maximum entry in each row as the diagonal preconditioner, instead of
403: the diagonal entry
405: Collective on PC
407: Input Parameters:
408: . pc - the preconditioner context
411: Options Database Key:
412: . -pc_jacobi_rowmax
414: Level: intermediate
416: Concepts: Jacobi preconditioner
418: @*/
419: int PCJacobiSetUseRowMax(PC pc)
420: {
421: int ierr,(*f)(PC);
425: PetscObjectQueryFunction((PetscObject)pc,"PCJacobiSetRowMax_C",(void (**)(void))&f);
426: if (f) {
427: (*f)(pc);
428: }
429: return(0);
430: }