1
2 /* --------------------------------------------------------------------
3
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)
8
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.
20
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.
26
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().
32
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().
38
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().
43
44 -------------------------------------------------------------------- */
45
46 /*
47 Include files needed for the Jacobi preconditioner:
48 pcimpl.h - private include file intended for use by all preconditioners
49 */
50
51 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
52
53 const char *const PCJacobiTypes[] = {"DIAGONAL","ROWMAX","ROWSUM","PCJacobiType","PC_JACOBI_",NULL};
54
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;
67
PCJacobiSetType_Jacobi(PC pc,PCJacobiType type)68 static PetscErrorCode PCJacobiSetType_Jacobi(PC pc,PCJacobiType type)
69 {
70 PC_Jacobi *j = (PC_Jacobi*)pc->data;
71
72 PetscFunctionBegin;
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 PetscFunctionReturn(0);
81 }
82
PCJacobiGetType_Jacobi(PC pc,PCJacobiType * type)83 static PetscErrorCode PCJacobiGetType_Jacobi(PC pc,PCJacobiType *type)
84 {
85 PC_Jacobi *j = (PC_Jacobi*)pc->data;
86
87 PetscFunctionBegin;
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 PetscFunctionReturn(0);
96 }
97
PCJacobiSetUseAbs_Jacobi(PC pc,PetscBool flg)98 static PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc,PetscBool flg)
99 {
100 PC_Jacobi *j = (PC_Jacobi*)pc->data;
101
102 PetscFunctionBegin;
103 j->useabs = flg;
104 PetscFunctionReturn(0);
105 }
106
PCJacobiGetUseAbs_Jacobi(PC pc,PetscBool * flg)107 static PetscErrorCode PCJacobiGetUseAbs_Jacobi(PC pc,PetscBool *flg)
108 {
109 PC_Jacobi *j = (PC_Jacobi*)pc->data;
110
111 PetscFunctionBegin;
112 *flg = j->useabs;
113 PetscFunctionReturn(0);
114 }
115
116 /* -------------------------------------------------------------------------- */
117 /*
118 PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
119 by setting data structures and options.
120
121 Input Parameter:
122 . pc - the preconditioner context
123
124 Application Interface Routine: PCSetUp()
125
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 */
PCSetUp_Jacobi(PC pc)130 static PetscErrorCode PCSetUp_Jacobi(PC pc)
131 {
132 PC_Jacobi *jac = (PC_Jacobi*)pc->data;
133 Vec diag,diagsqrt;
134 PetscErrorCode ierr;
135 PetscInt n,i;
136 PetscScalar *x;
137 PetscBool zeroflag = PETSC_FALSE;
138
139 PetscFunctionBegin;
140 /*
141 For most preconditioners the code would begin here something like
142
143 if (pc->setupcalled == 0) { allocate space the first time this is ever called
144 ierr = MatCreateVecs(pc->mat,&jac->diag);CHKERRQ(ierr);
145 PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);
146 }
147
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 */
156
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;
163
164 if (diag) {
165 if (jac->userowmax) {
166 ierr = MatGetRowMaxAbs(pc->pmat,diag,NULL);CHKERRQ(ierr);
167 } else if (jac->userowsum) {
168 ierr = MatGetRowSum(pc->pmat,diag);CHKERRQ(ierr);
169 } else {
170 ierr = MatGetDiagonal(pc->pmat,diag);CHKERRQ(ierr);
171 }
172 ierr = VecReciprocal(diag);CHKERRQ(ierr);
173 ierr = VecGetLocalSize(diag,&n);CHKERRQ(ierr);
174 ierr = VecGetArray(diag,&x);CHKERRQ(ierr);
175 if (jac->useabs) {
176 for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
177 }
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 ierr = VecRestoreArray(diag,&x);CHKERRQ(ierr);
185 }
186 if (diagsqrt) {
187 if (jac->userowmax) {
188 ierr = MatGetRowMaxAbs(pc->pmat,diagsqrt,NULL);CHKERRQ(ierr);
189 } else if (jac->userowsum) {
190 ierr = MatGetRowSum(pc->pmat,diagsqrt);CHKERRQ(ierr);
191 } else {
192 ierr = MatGetDiagonal(pc->pmat,diagsqrt);CHKERRQ(ierr);
193 }
194 ierr = VecGetLocalSize(diagsqrt,&n);CHKERRQ(ierr);
195 ierr = VecGetArray(diagsqrt,&x);CHKERRQ(ierr);
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 ierr = VecRestoreArray(diagsqrt,&x);CHKERRQ(ierr);
204 }
205 if (zeroflag) {
206 ierr = PetscInfo(pc,"Zero detected in diagonal of matrix, using 1 at those locations\n");CHKERRQ(ierr);
207 }
208 PetscFunctionReturn(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.
215
216 Input Parameter:
217 . pc - the preconditioner context
218 */
PCSetUp_Jacobi_Symmetric(PC pc)219 static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
220 {
221 PetscErrorCode ierr;
222 PC_Jacobi *jac = (PC_Jacobi*)pc->data;
223
224 PetscFunctionBegin;
225 ierr = MatCreateVecs(pc->pmat,&jac->diagsqrt,NULL);CHKERRQ(ierr);
226 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diagsqrt);CHKERRQ(ierr);
227 ierr = PCSetUp_Jacobi(pc);CHKERRQ(ierr);
228 PetscFunctionReturn(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.
235
236 Input Parameter:
237 . pc - the preconditioner context
238 */
PCSetUp_Jacobi_NonSymmetric(PC pc)239 static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
240 {
241 PetscErrorCode ierr;
242 PC_Jacobi *jac = (PC_Jacobi*)pc->data;
243
244 PetscFunctionBegin;
245 ierr = MatCreateVecs(pc->pmat,&jac->diag,NULL);CHKERRQ(ierr);
246 ierr = PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->diag);CHKERRQ(ierr);
247 ierr = PCSetUp_Jacobi(pc);CHKERRQ(ierr);
248 PetscFunctionReturn(0);
249 }
250 /* -------------------------------------------------------------------------- */
251 /*
252 PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
253
254 Input Parameters:
255 . pc - the preconditioner context
256 . x - input vector
257
258 Output Parameter:
259 . y - output vector
260
261 Application Interface Routine: PCApply()
262 */
PCApply_Jacobi(PC pc,Vec x,Vec y)263 static PetscErrorCode PCApply_Jacobi(PC pc,Vec x,Vec y)
264 {
265 PC_Jacobi *jac = (PC_Jacobi*)pc->data;
266 PetscErrorCode ierr;
267
268 PetscFunctionBegin;
269 if (!jac->diag) {
270 ierr = PCSetUp_Jacobi_NonSymmetric(pc);CHKERRQ(ierr);
271 }
272 ierr = VecPointwiseMult(y,x,jac->diag);CHKERRQ(ierr);
273 PetscFunctionReturn(0);
274 }
275 /* -------------------------------------------------------------------------- */
276 /*
277 PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
278 symmetric preconditioner to a vector.
279
280 Input Parameters:
281 . pc - the preconditioner context
282 . x - input vector
283
284 Output Parameter:
285 . y - output vector
286
287 Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
288 */
PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)289 static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)
290 {
291 PetscErrorCode ierr;
292 PC_Jacobi *jac = (PC_Jacobi*)pc->data;
293
294 PetscFunctionBegin;
295 if (!jac->diagsqrt) {
296 ierr = PCSetUp_Jacobi_Symmetric(pc);CHKERRQ(ierr);
297 }
298 VecPointwiseMult(y,x,jac->diagsqrt);
299 PetscFunctionReturn(0);
300 }
301 /* -------------------------------------------------------------------------- */
PCReset_Jacobi(PC pc)302 static PetscErrorCode PCReset_Jacobi(PC pc)
303 {
304 PC_Jacobi *jac = (PC_Jacobi*)pc->data;
305 PetscErrorCode ierr;
306
307 PetscFunctionBegin;
308 ierr = VecDestroy(&jac->diag);CHKERRQ(ierr);
309 ierr = VecDestroy(&jac->diagsqrt);CHKERRQ(ierr);
310 PetscFunctionReturn(0);
311 }
312
313 /*
314 PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
315 that was created with PCCreate_Jacobi().
316
317 Input Parameter:
318 . pc - the preconditioner context
319
320 Application Interface Routine: PCDestroy()
321 */
PCDestroy_Jacobi(PC pc)322 static PetscErrorCode PCDestroy_Jacobi(PC pc)
323 {
324 PetscErrorCode ierr;
325
326 PetscFunctionBegin;
327 ierr = PCReset_Jacobi(pc);CHKERRQ(ierr);
328
329 /*
330 Free the private data structure that was hanging off the PC
331 */
332 ierr = PetscFree(pc->data);CHKERRQ(ierr);
333 PetscFunctionReturn(0);
334 }
335
PCSetFromOptions_Jacobi(PetscOptionItems * PetscOptionsObject,PC pc)336 static PetscErrorCode PCSetFromOptions_Jacobi(PetscOptionItems *PetscOptionsObject,PC pc)
337 {
338 PC_Jacobi *jac = (PC_Jacobi*)pc->data;
339 PetscErrorCode ierr;
340 PetscBool flg;
341 PCJacobiType deflt,type;
342
343 PetscFunctionBegin;
344 ierr = PCJacobiGetType(pc,&deflt);CHKERRQ(ierr);
345 ierr = PetscOptionsHead(PetscOptionsObject,"Jacobi options");CHKERRQ(ierr);
346 ierr = PetscOptionsEnum("-pc_jacobi_type","How to construct diagonal matrix","PCJacobiSetType",PCJacobiTypes,(PetscEnum)deflt,(PetscEnum*)&type,&flg);CHKERRQ(ierr);
347 if (flg) {
348 ierr = PCJacobiSetType(pc,type);CHKERRQ(ierr);
349 }
350 ierr = PetscOptionsBool("-pc_jacobi_abs","Use absolute values of diagaonal entries","PCJacobiSetUseAbs",jac->useabs,&jac->useabs,NULL);CHKERRQ(ierr);
351 ierr = PetscOptionsTail();CHKERRQ(ierr);
352 PetscFunctionReturn(0);
353 }
354
355 /* -------------------------------------------------------------------------- */
356 /*
357 PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
358 and sets this as the private data within the generic preconditioning
359 context, PC, that was created within PCCreate().
360
361 Input Parameter:
362 . pc - the preconditioner context
363
364 Application Interface Routine: PCCreate()
365 */
366
367 /*MC
368 PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)
369
370 Options Database Key:
371 + -pc_jacobi_type <diagonal,rowmax,rowsum> - approach for forming the preconditioner
372 - -pc_jacobi_abs - use the absolute value of the diagonal entry
373
374 Level: beginner
375
376 Notes:
377 By using KSPSetPCSide(ksp,PC_SYMMETRIC) or -ksp_pc_side symmetric
378 can scale each side of the matrix by the square root of the diagonal entries.
379
380 Zero entries along the diagonal are replaced with the value 1.0
381
382 See PCPBJACOBI for a point-block Jacobi preconditioner
383
384 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
385 PCJacobiSetType(), PCJacobiSetUseAbs(), PCJacobiGetUseAbs(), PCPBJACOBI
386 M*/
387
PCCreate_Jacobi(PC pc)388 PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
389 {
390 PC_Jacobi *jac;
391 PetscErrorCode ierr;
392
393 PetscFunctionBegin;
394 /*
395 Creates the private data structure for this preconditioner and
396 attach it to the PC object.
397 */
398 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr);
399 pc->data = (void*)jac;
400
401 /*
402 Initialize the pointers to vectors to ZERO; these will be used to store
403 diagonal entries of the matrix for fast preconditioner application.
404 */
405 jac->diag = NULL;
406 jac->diagsqrt = NULL;
407 jac->userowmax = PETSC_FALSE;
408 jac->userowsum = PETSC_FALSE;
409 jac->useabs = PETSC_FALSE;
410
411 /*
412 Set the pointers for the functions that are provided above.
413 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
414 are called, they will automatically call these functions. Note we
415 choose not to provide a couple of these functions since they are
416 not needed.
417 */
418 pc->ops->apply = PCApply_Jacobi;
419 pc->ops->applytranspose = PCApply_Jacobi;
420 pc->ops->setup = PCSetUp_Jacobi;
421 pc->ops->reset = PCReset_Jacobi;
422 pc->ops->destroy = PCDestroy_Jacobi;
423 pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
424 pc->ops->view = NULL;
425 pc->ops->applyrichardson = NULL;
426 pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
427 pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
428
429 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetType_C",PCJacobiSetType_Jacobi);CHKERRQ(ierr);
430 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetType_C",PCJacobiGetType_Jacobi);CHKERRQ(ierr);
431 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCJacobiSetUseAbs_C",PCJacobiSetUseAbs_Jacobi);CHKERRQ(ierr);
432 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCJacobiGetUseAbs_C",PCJacobiGetUseAbs_Jacobi);CHKERRQ(ierr);
433 PetscFunctionReturn(0);
434 }
435
436 /*@
437 PCJacobiSetUseAbs - Causes the Jacobi preconditioner to use the
438 absolute values of the digonal divisors in the preconditioner
439
440 Logically Collective on PC
441
442 Input Parameters:
443 + pc - the preconditioner context
444 - flg - whether to use absolute values or not
445
446 Options Database Key:
447 . -pc_jacobi_abs
448
449 Notes:
450 This takes affect at the next construction of the preconditioner
451
452 Level: intermediate
453
454 .seealso: PCJacobiaSetType(), PCJacobiGetUseAbs()
455
456 @*/
PCJacobiSetUseAbs(PC pc,PetscBool flg)457 PetscErrorCode PCJacobiSetUseAbs(PC pc,PetscBool flg)
458 {
459 PetscErrorCode ierr;
460
461 PetscFunctionBegin;
462 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
463 ierr = PetscTryMethod(pc,"PCJacobiSetUseAbs_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr);
464 PetscFunctionReturn(0);
465 }
466
467 /*@
468 PCJacobiGetUseAbs - Determines if the Jacobi preconditioner uses the
469 absolute values of the digonal divisors in the preconditioner
470
471 Logically Collective on PC
472
473 Input Parameter:
474 . pc - the preconditioner context
475
476 Output Parameter:
477 . flg - whether to use absolute values or not
478
479 Options Database Key:
480 . -pc_jacobi_abs
481
482 Level: intermediate
483
484 .seealso: PCJacobiaSetType(), PCJacobiSetUseAbs(), PCJacobiGetType()
485
486 @*/
PCJacobiGetUseAbs(PC pc,PetscBool * flg)487 PetscErrorCode PCJacobiGetUseAbs(PC pc,PetscBool *flg)
488 {
489 PetscErrorCode ierr;
490
491 PetscFunctionBegin;
492 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
493 ierr = PetscUseMethod(pc,"PCJacobiGetUseAbs_C",(PC,PetscBool*),(pc,flg));CHKERRQ(ierr);
494 PetscFunctionReturn(0);
495 }
496
497 /*@
498 PCJacobiSetType - Causes the Jacobi preconditioner to use either the diagonal, the maximum entry in each row,
499 of the sum of rows entries for the diagonal preconditioner
500
501 Logically Collective on PC
502
503 Input Parameters:
504 + pc - the preconditioner context
505 - type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM
506
507 Options Database Key:
508 . -pc_jacobi_type <diagonal,rowmax,rowsum>
509
510 Level: intermediate
511
512 .seealso: PCJacobiaUseAbs(), PCJacobiGetType()
513 @*/
PCJacobiSetType(PC pc,PCJacobiType type)514 PetscErrorCode PCJacobiSetType(PC pc,PCJacobiType type)
515 {
516 PetscErrorCode ierr;
517
518 PetscFunctionBegin;
519 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
520 ierr = PetscTryMethod(pc,"PCJacobiSetType_C",(PC,PCJacobiType),(pc,type));CHKERRQ(ierr);
521 PetscFunctionReturn(0);
522 }
523
524 /*@
525 PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner
526
527 Not Collective on PC
528
529 Input Parameter:
530 . pc - the preconditioner context
531
532 Output Parameter:
533 . type - PC_JACOBI_DIAGONAL, PC_JACOBI_ROWMAX, PC_JACOBI_ROWSUM
534
535 Level: intermediate
536
537 .seealso: PCJacobiaUseAbs(), PCJacobiSetType()
538 @*/
PCJacobiGetType(PC pc,PCJacobiType * type)539 PetscErrorCode PCJacobiGetType(PC pc,PCJacobiType *type)
540 {
541 PetscErrorCode ierr;
542
543 PetscFunctionBegin;
544 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
545 ierr = PetscUseMethod(pc,"PCJacobiGetType_C",(PC,PCJacobiType*),(pc,type));CHKERRQ(ierr);
546 PetscFunctionReturn(0);
547 }
548