1 /*
2 This file implements a Jacobi preconditioner in PETSc as part of PC.
3 You can use this as a starting point for implementing your own
4 preconditioner that is not provided with PETSc. (You might also consider
5 just using PCSHELL)
6
7 The following basic routines are required for each preconditioner.
8 PCCreate_XXX() - Creates a preconditioner context
9 PCSetFromOptions_XXX() - Sets runtime options
10 PCApply_XXX() - Applies the preconditioner
11 PCDestroy_XXX() - Destroys the preconditioner context
12 where the suffix "_XXX" denotes a particular implementation, in
13 this case we use _Jacobi (e.g., PCCreate_Jacobi, PCApply_Jacobi).
14 These routines are actually called via the common user interface
15 routines PCCreate(), PCSetFromOptions(), PCApply(), and PCDestroy(),
16 so the application code interface remains identical for all
17 preconditioners.
18
19 Another key routine is:
20 PCSetUp_XXX() - Prepares for the use of a preconditioner
21 by setting data structures and options. The interface routine PCSetUp()
22 is not usually called directly by the user, but instead is called by
23 PCApply() if necessary.
24
25 Additional basic routines are:
26 PCView_XXX() - Prints details of runtime options that
27 have actually been used.
28 These are called by application codes via the interface routines
29 PCView().
30
31 The various types of solvers (preconditioners, Krylov subspace methods,
32 nonlinear solvers, timesteppers) are all organized similarly, so the
33 above description applies to these categories also. One exception is
34 that the analogues of PCApply() for these components are KSPSolve(),
35 SNESSolve(), and TSSolve().
36
37 Additional optional functionality unique to preconditioners is left and
38 right symmetric preconditioner application via PCApplySymmetricLeft()
39 and PCApplySymmetricRight(). The Jacobi implementation is
40 PCApplySymmetricLeftOrRight_Jacobi().
41 */
42
43 /*
44 Include files needed for the Jacobi preconditioner:
45 pcimpl.h - private include file intended for use by all preconditioners
46 */
47
48 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
49
50 const char *const PCJacobiTypes[] = {"DIAGONAL", "ROWL1", "ROWMAX", "ROWSUM", "PCJacobiType", "PC_JACOBI_", NULL};
51
52 /*
53 Private context (data structure) for the Jacobi preconditioner.
54 */
55 typedef struct {
56 Vec diag; /* vector containing the reciprocals of the diagonal elements of the matrix used to construct the preconditioner */
57 Vec diagsqrt; /* vector containing the reciprocals of the square roots of
58 the diagonal elements of the matrix used to compute the preconditioner (used
59 only for symmetric preconditioner application) */
60 PCJacobiType type;
61 PetscBool useabs; /* use the absolute values of the diagonal entries */
62 PetscBool fixdiag; /* fix zero diagonal terms */
63 PetscReal scale; /* for scaling rowl1 off-diagonals */
64 } PC_Jacobi;
65
66 static PetscErrorCode PCReset_Jacobi(PC);
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 PCJacobiType old_type;
72
73 PetscFunctionBegin;
74 PetscCall(PCJacobiGetType(pc, &old_type));
75 if (old_type == type) PetscFunctionReturn(PETSC_SUCCESS);
76 PetscCall(PCReset_Jacobi(pc));
77 j->type = type;
78 PetscFunctionReturn(PETSC_SUCCESS);
79 }
80
PCJacobiGetUseAbs_Jacobi(PC pc,PetscBool * flg)81 static PetscErrorCode PCJacobiGetUseAbs_Jacobi(PC pc, PetscBool *flg)
82 {
83 PC_Jacobi *j = (PC_Jacobi *)pc->data;
84
85 PetscFunctionBegin;
86 *flg = j->useabs;
87 PetscFunctionReturn(PETSC_SUCCESS);
88 }
89
PCJacobiSetUseAbs_Jacobi(PC pc,PetscBool flg)90 static PetscErrorCode PCJacobiSetUseAbs_Jacobi(PC pc, PetscBool flg)
91 {
92 PC_Jacobi *j = (PC_Jacobi *)pc->data;
93
94 PetscFunctionBegin;
95 j->useabs = flg;
96 PetscFunctionReturn(PETSC_SUCCESS);
97 }
98
PCJacobiGetType_Jacobi(PC pc,PCJacobiType * type)99 static PetscErrorCode PCJacobiGetType_Jacobi(PC pc, PCJacobiType *type)
100 {
101 PC_Jacobi *j = (PC_Jacobi *)pc->data;
102
103 PetscFunctionBegin;
104 *type = j->type;
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
PCJacobiSetRowl1Scale_Jacobi(PC pc,PetscReal flg)108 static PetscErrorCode PCJacobiSetRowl1Scale_Jacobi(PC pc, PetscReal flg)
109 {
110 PC_Jacobi *j = (PC_Jacobi *)pc->data;
111
112 PetscFunctionBegin;
113 j->scale = flg;
114 PetscFunctionReturn(PETSC_SUCCESS);
115 }
116
PCJacobiGetRowl1Scale_Jacobi(PC pc,PetscReal * flg)117 static PetscErrorCode PCJacobiGetRowl1Scale_Jacobi(PC pc, PetscReal *flg)
118 {
119 PC_Jacobi *j = (PC_Jacobi *)pc->data;
120
121 PetscFunctionBegin;
122 *flg = j->scale;
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125
PCJacobiSetFixDiagonal_Jacobi(PC pc,PetscBool flg)126 static PetscErrorCode PCJacobiSetFixDiagonal_Jacobi(PC pc, PetscBool flg)
127 {
128 PC_Jacobi *j = (PC_Jacobi *)pc->data;
129
130 PetscFunctionBegin;
131 j->fixdiag = flg;
132 PetscFunctionReturn(PETSC_SUCCESS);
133 }
134
PCJacobiGetFixDiagonal_Jacobi(PC pc,PetscBool * flg)135 static PetscErrorCode PCJacobiGetFixDiagonal_Jacobi(PC pc, PetscBool *flg)
136 {
137 PC_Jacobi *j = (PC_Jacobi *)pc->data;
138
139 PetscFunctionBegin;
140 *flg = j->fixdiag;
141 PetscFunctionReturn(PETSC_SUCCESS);
142 }
143
PCJacobiGetDiagonal_Jacobi(PC pc,Vec diag,Vec diagsqrt)144 static PetscErrorCode PCJacobiGetDiagonal_Jacobi(PC pc, Vec diag, Vec diagsqrt)
145 {
146 PC_Jacobi *j = (PC_Jacobi *)pc->data;
147 MPI_Comm comm = PetscObjectComm((PetscObject)pc);
148
149 PetscFunctionBegin;
150 PetscCheck(j->diag || j->diagsqrt, comm, PETSC_ERR_ARG_WRONGSTATE, "Jacobi diagonal has not been created yet. Use PCApply to force creation");
151 PetscCheck(!diag || (diag && j->diag), comm, PETSC_ERR_ARG_WRONGSTATE, "Jacobi diagonal not available. Check if PC is non-symmetric");
152 PetscCheck(!diagsqrt || (diagsqrt && j->diagsqrt), comm, PETSC_ERR_ARG_WRONGSTATE, "Jacobi diagonal squareroot not available. Check if PC is symmetric");
153
154 if (diag) PetscCall(VecCopy(j->diag, diag));
155 if (diagsqrt) PetscCall(VecCopy(j->diagsqrt, diagsqrt));
156 PetscFunctionReturn(PETSC_SUCCESS);
157 }
158
159 /*
160 PCSetUp_Jacobi - Prepares for the use of the Jacobi preconditioner
161 by setting data structures and options.
162
163 Input Parameter:
164 . pc - the preconditioner context
165
166 Application Interface Routine: PCSetUp()
167
168 Note:
169 The interface routine PCSetUp() is not usually called directly by
170 the user, but instead is called by PCApply() if necessary.
171 */
PCSetUp_Jacobi(PC pc)172 static PetscErrorCode PCSetUp_Jacobi(PC pc)
173 {
174 PC_Jacobi *jac = (PC_Jacobi *)pc->data;
175 Vec diag, diagsqrt;
176 PetscInt n, i;
177 PetscBool zeroflag = PETSC_FALSE, negflag = PETSC_FALSE;
178
179 PetscFunctionBegin;
180 /*
181 For most preconditioners the code would begin here something like
182
183 if (!pc->setupcalled) { allocate space the first time this is ever called
184 PetscCall(MatCreateVecs(pc->mat,&jac->diag));
185 }
186
187 But for this preconditioner we want to support use of both the matrix' diagonal
188 elements (for left or right preconditioning) and square root of diagonal elements
189 (for symmetric preconditioning). Hence we do not allocate space here, since we
190 don't know at this point which will be needed (diag and/or diagsqrt) until the user
191 applies the preconditioner, and we don't want to allocate BOTH unless we need
192 them both. Thus, the diag and diagsqrt are allocated in PCSetUp_Jacobi_NonSymmetric()
193 and PCSetUp_Jacobi_Symmetric(), respectively.
194 */
195
196 /*
197 Here we set up the preconditioner; that is, we copy the diagonal values from
198 the matrix and put them into a format to make them quick to apply as a preconditioner.
199 */
200 diag = jac->diag;
201 diagsqrt = jac->diagsqrt;
202
203 if (diag) {
204 PetscBool isset, isspd;
205
206 PetscCall(VecLockReadPop(diag));
207 switch (jac->type) {
208 case PC_JACOBI_DIAGONAL:
209 PetscCall(MatGetDiagonal(pc->pmat, diag));
210 break;
211 case PC_JACOBI_ROWMAX:
212 PetscCall(MatGetRowMaxAbs(pc->pmat, diag, NULL));
213 break;
214 case PC_JACOBI_ROWL1:
215 PetscCall(MatGetRowSumAbs(pc->pmat, diag));
216 // fix negative rows (eg, negative definite) -- this could be done for all, not needed for userowmax
217 PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
218 if (jac->fixdiag && (!isset || !isspd)) {
219 PetscScalar *x2;
220 const PetscScalar *x;
221 Vec true_diag;
222 PetscCall(VecDuplicate(diag, &true_diag));
223 PetscCall(MatGetDiagonal(pc->pmat, true_diag));
224 PetscCall(VecGetLocalSize(diag, &n));
225 PetscCall(VecGetArrayWrite(diag, &x2));
226 PetscCall(VecGetArrayRead(true_diag, &x)); // to make more general -todo
227 for (i = 0; i < n; i++) {
228 if (PetscRealPart(x[i]) < 0.0) {
229 x2[i] = -x2[i]; // flip sign to keep DA > 0
230 negflag = PETSC_TRUE;
231 }
232 }
233 PetscCall(VecRestoreArrayRead(true_diag, &x));
234 PetscCall(VecRestoreArrayWrite(diag, &x2));
235 PetscCheck(!jac->useabs || !negflag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Jacobi use_abs and l1 not compatible with negative diagonal");
236 PetscCall(VecDestroy(&true_diag));
237 }
238 if (jac->scale != 1.0) {
239 Vec true_diag;
240 PetscCall(VecDuplicate(diag, &true_diag));
241 PetscCall(MatGetDiagonal(pc->pmat, true_diag));
242 PetscCall(VecAXPY(diag, -1, true_diag)); // subtract off diag
243 PetscCall(VecScale(diag, jac->scale)); // scale off-diag
244 PetscCall(VecAXPY(diag, 1, true_diag)); // add diag back in
245 PetscCall(VecDestroy(&true_diag));
246 }
247 break;
248 case PC_JACOBI_ROWSUM:
249 PetscCall(MatGetRowSum(pc->pmat, diag));
250 break;
251 }
252 PetscCall(VecReciprocal(diag));
253 if (jac->useabs) PetscCall(VecAbs(diag));
254 PetscCall(MatIsSPDKnown(pc->pmat, &isset, &isspd));
255 if (jac->fixdiag && (!isset || !isspd)) {
256 PetscScalar *x;
257 PetscCall(VecGetLocalSize(diag, &n));
258 PetscCall(VecGetArray(diag, &x));
259 for (i = 0; i < n; i++) {
260 if (x[i] == 0.0) {
261 x[i] = 1.0;
262 zeroflag = PETSC_TRUE;
263 }
264 }
265 PetscCall(VecRestoreArray(diag, &x));
266 }
267 PetscCall(VecLockReadPush(diag));
268 }
269 if (diagsqrt) {
270 PetscScalar *x;
271
272 PetscCall(VecLockReadPop(diagsqrt));
273 switch (jac->type) {
274 case PC_JACOBI_DIAGONAL:
275 PetscCall(MatGetDiagonal(pc->pmat, diagsqrt));
276 break;
277 case PC_JACOBI_ROWMAX:
278 PetscCall(MatGetRowMaxAbs(pc->pmat, diagsqrt, NULL));
279 break;
280 case PC_JACOBI_ROWL1:
281 PetscCall(MatGetRowSumAbs(pc->pmat, diagsqrt));
282 break;
283 case PC_JACOBI_ROWSUM:
284 PetscCall(MatGetRowSum(pc->pmat, diagsqrt));
285 break;
286 }
287 PetscCall(VecGetLocalSize(diagsqrt, &n));
288 PetscCall(VecGetArray(diagsqrt, &x));
289 for (i = 0; i < n; i++) {
290 if (PetscRealPart(x[i]) < 0.0) x[i] = 1.0 / PetscSqrtReal(PetscAbsScalar(-x[i]));
291 else if (PetscRealPart(x[i]) > 0.0) x[i] = 1.0 / PetscSqrtReal(PetscAbsScalar(x[i]));
292 else {
293 x[i] = 1.0;
294 zeroflag = PETSC_TRUE;
295 }
296 }
297 PetscCall(VecRestoreArray(diagsqrt, &x));
298 PetscCall(VecLockReadPush(diagsqrt));
299 }
300 if (zeroflag) PetscCall(PetscInfo(pc, "Zero detected in diagonal of matrix, using 1 at those locations\n"));
301 PetscFunctionReturn(PETSC_SUCCESS);
302 }
303
304 /*
305 PCSetUp_Jacobi_Symmetric - Allocates the vector needed to store the
306 inverse of the square root of the diagonal entries of the matrix. This
307 is used for symmetric application of the Jacobi preconditioner.
308
309 Input Parameter:
310 . pc - the preconditioner context
311 */
PCSetUp_Jacobi_Symmetric(PC pc)312 static PetscErrorCode PCSetUp_Jacobi_Symmetric(PC pc)
313 {
314 PC_Jacobi *jac = (PC_Jacobi *)pc->data;
315
316 PetscFunctionBegin;
317 PetscCall(MatCreateVecs(pc->pmat, &jac->diagsqrt, NULL));
318 PetscCall(VecLockReadPush(jac->diagsqrt));
319 PetscCall(PCSetUp_Jacobi(pc));
320 PetscFunctionReturn(PETSC_SUCCESS);
321 }
322
323 /*
324 PCSetUp_Jacobi_NonSymmetric - Allocates the vector needed to store the
325 inverse of the diagonal entries of the matrix. This is used for left of
326 right application of the Jacobi preconditioner.
327
328 Input Parameter:
329 . pc - the preconditioner context
330 */
PCSetUp_Jacobi_NonSymmetric(PC pc)331 static PetscErrorCode PCSetUp_Jacobi_NonSymmetric(PC pc)
332 {
333 PC_Jacobi *jac = (PC_Jacobi *)pc->data;
334
335 PetscFunctionBegin;
336 PetscCall(MatCreateVecs(pc->pmat, &jac->diag, NULL));
337 PetscCall(VecLockReadPush(jac->diag));
338 PetscCall(PCSetUp_Jacobi(pc));
339 PetscFunctionReturn(PETSC_SUCCESS);
340 }
341
342 /*
343 PCApply_Jacobi - Applies the Jacobi preconditioner to a vector.
344
345 Input Parameters:
346 . pc - the preconditioner context
347 . x - input vector
348
349 Output Parameter:
350 . y - output vector
351
352 Application Interface Routine: PCApply()
353 */
PCApply_Jacobi(PC pc,Vec x,Vec y)354 static PetscErrorCode PCApply_Jacobi(PC pc, Vec x, Vec y)
355 {
356 PC_Jacobi *jac = (PC_Jacobi *)pc->data;
357
358 PetscFunctionBegin;
359 if (!jac->diag) PetscCall(PCSetUp_Jacobi_NonSymmetric(pc));
360 PetscCall(VecPointwiseMult(y, x, jac->diag));
361 PetscFunctionReturn(PETSC_SUCCESS);
362 }
363
364 /*
365 PCApplySymmetricLeftOrRight_Jacobi - Applies the left or right part of a
366 symmetric preconditioner to a vector.
367
368 Input Parameters:
369 . pc - the preconditioner context
370 . x - input vector
371
372 Output Parameter:
373 . y - output vector
374
375 Application Interface Routines: PCApplySymmetricLeft(), PCApplySymmetricRight()
376 */
PCApplySymmetricLeftOrRight_Jacobi(PC pc,Vec x,Vec y)377 static PetscErrorCode PCApplySymmetricLeftOrRight_Jacobi(PC pc, Vec x, Vec y)
378 {
379 PC_Jacobi *jac = (PC_Jacobi *)pc->data;
380
381 PetscFunctionBegin;
382 if (!jac->diagsqrt) PetscCall(PCSetUp_Jacobi_Symmetric(pc));
383 PetscCall(VecPointwiseMult(y, x, jac->diagsqrt));
384 PetscFunctionReturn(PETSC_SUCCESS);
385 }
386
PCReset_Jacobi(PC pc)387 static PetscErrorCode PCReset_Jacobi(PC pc)
388 {
389 PC_Jacobi *jac = (PC_Jacobi *)pc->data;
390
391 PetscFunctionBegin;
392 if (jac->diag) PetscCall(VecLockReadPop(jac->diag));
393 if (jac->diagsqrt) PetscCall(VecLockReadPop(jac->diagsqrt));
394 PetscCall(VecDestroy(&jac->diag));
395 PetscCall(VecDestroy(&jac->diagsqrt));
396 PetscFunctionReturn(PETSC_SUCCESS);
397 }
398
399 /*
400 PCDestroy_Jacobi - Destroys the private context for the Jacobi preconditioner
401 that was created with PCCreate_Jacobi().
402
403 Input Parameter:
404 . pc - the preconditioner context
405
406 Application Interface Routine: PCDestroy()
407 */
PCDestroy_Jacobi(PC pc)408 static PetscErrorCode PCDestroy_Jacobi(PC pc)
409 {
410 PetscFunctionBegin;
411 PetscCall(PCReset_Jacobi(pc));
412 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetType_C", NULL));
413 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetType_C", NULL));
414 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetUseAbs_C", NULL));
415 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetUseAbs_C", NULL));
416 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetRowl1Scale_C", NULL));
417 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetRowl1Scale_C", NULL));
418 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetFixDiagonal_C", NULL));
419 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetFixDiagonal_C", NULL));
420 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetDiagonal_C", NULL));
421
422 /*
423 Free the private data structure that was hanging off the PC
424 */
425 PetscCall(PetscFree(pc->data));
426 PetscFunctionReturn(PETSC_SUCCESS);
427 }
428
PCSetFromOptions_Jacobi(PC pc,PetscOptionItems PetscOptionsObject)429 static PetscErrorCode PCSetFromOptions_Jacobi(PC pc, PetscOptionItems PetscOptionsObject)
430 {
431 PC_Jacobi *jac = (PC_Jacobi *)pc->data;
432 PetscBool flg;
433 PCJacobiType deflt, type;
434
435 PetscFunctionBegin;
436 PetscCall(PCJacobiGetType(pc, &deflt));
437 PetscOptionsHeadBegin(PetscOptionsObject, "Jacobi options");
438 PetscCall(PetscOptionsEnum("-pc_jacobi_type", "How to construct diagonal matrix", "PCJacobiSetType", PCJacobiTypes, (PetscEnum)deflt, (PetscEnum *)&type, &flg));
439 if (flg) PetscCall(PCJacobiSetType(pc, type));
440 PetscCall(PetscOptionsBool("-pc_jacobi_abs", "Use absolute values of diagonal entries", "PCJacobiSetUseAbs", jac->useabs, &jac->useabs, NULL));
441 PetscCall(PetscOptionsBool("-pc_jacobi_fixdiagonal", "Fix null terms on diagonal", "PCJacobiSetFixDiagonal", jac->fixdiag, &jac->fixdiag, NULL));
442 PetscCall(PetscOptionsRangeReal("-pc_jacobi_rowl1_scale", "scaling of off-diagonal elements for rowl1", "PCJacobiSetRowl1Scale", jac->scale, &jac->scale, NULL, 0.0, 1.0));
443 PetscOptionsHeadEnd();
444 PetscFunctionReturn(PETSC_SUCCESS);
445 }
446
PCView_Jacobi(PC pc,PetscViewer viewer)447 static PetscErrorCode PCView_Jacobi(PC pc, PetscViewer viewer)
448 {
449 PC_Jacobi *jac = (PC_Jacobi *)pc->data;
450 PetscBool isascii;
451
452 PetscFunctionBegin;
453 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
454 if (isascii) {
455 PCJacobiType type;
456 PetscBool useAbs, fixdiag;
457 PetscViewerFormat format;
458 PetscReal scale;
459
460 PetscCall(PCJacobiGetType(pc, &type));
461 PetscCall(PCJacobiGetUseAbs(pc, &useAbs));
462 PetscCall(PCJacobiGetFixDiagonal(pc, &fixdiag));
463 PetscCall(PCJacobiGetRowl1Scale(pc, &scale));
464 if (type == PC_JACOBI_ROWL1)
465 PetscCall(PetscViewerASCIIPrintf(viewer, " type %s%s%s (l1-norm off-diagonal scaling %e)\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "", !fixdiag ? ", not checking null diagonal entries" : "", (double)scale));
466 else PetscCall(PetscViewerASCIIPrintf(viewer, " type %s%s%s\n", PCJacobiTypes[type], useAbs ? ", using absolute value of entries" : "", !fixdiag ? ", not checking null diagonal entries" : ""));
467 PetscCall(PetscViewerGetFormat(viewer, &format));
468 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL && jac->diag) PetscCall(VecView(jac->diag, viewer));
469 }
470 PetscFunctionReturn(PETSC_SUCCESS);
471 }
472
473 /*
474 PCCreate_Jacobi - Creates a Jacobi preconditioner context, PC_Jacobi,
475 and sets this as the private data within the generic preconditioning
476 context, PC, that was created within PCCreate().
477
478 Input Parameter:
479 . pc - the preconditioner context
480
481 Application Interface Routine: PCCreate()
482 */
483
484 /*MC
485 PCJACOBI - Jacobi (i.e. diagonal scaling preconditioning)
486
487 Options Database Keys:
488 + -pc_jacobi_type <diagonal,rowl1,rowmax,rowsum> - approach for forming the preconditioner
489 . -pc_jacobi_abs - use the absolute value of the diagonal entry
490 . -pc_jacobi_rowl1_scale - scaling of off-diagonal terms
491 - -pc_jacobi_fixdiag - fix for zero diagonal terms by placing 1.0 in those locations
492
493 Level: beginner
494
495 Notes:
496 By using `KSPSetPCSide`(ksp,`PC_SYMMETRIC`) or -ksp_pc_side symmetric
497 can scale each side of the matrix by the square root of the diagonal entries.
498
499 Zero entries along the diagonal are replaced with the value 1.0
500
501 See `PCPBJACOBI` for fixed-size point block, `PCVPBJACOBI` for variable-sized point block, and `PCBJACOBI` for large size blocks
502
503 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
504 `PCJacobiSetType()`, `PCJacobiSetUseAbs()`, `PCJacobiGetUseAbs()`, `PCASM`,
505 `PCJacobiSetFixDiagonal()`, `PCJacobiGetFixDiagonal()`
506 `PCJacobiSetType()`, `PCJacobiSetUseAbs()`, `PCJacobiGetUseAbs()`, `PCPBJACOBI`, `PCBJACOBI`, `PCVPBJACOBI`
507 M*/
508
PCCreate_Jacobi(PC pc)509 PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC pc)
510 {
511 PC_Jacobi *jac;
512
513 PetscFunctionBegin;
514 /*
515 Creates the private data structure for this preconditioner and
516 attach it to the PC object.
517 */
518 PetscCall(PetscNew(&jac));
519 pc->data = (void *)jac;
520
521 /*
522 Initialize the pointers to vectors to ZERO; these will be used to store
523 diagonal entries of the matrix for fast preconditioner application.
524 */
525 jac->diag = NULL;
526 jac->diagsqrt = NULL;
527 jac->type = PC_JACOBI_DIAGONAL;
528 jac->useabs = PETSC_FALSE;
529 jac->fixdiag = PETSC_TRUE;
530 jac->scale = 1.0;
531
532 /*
533 Set the pointers for the functions that are provided above.
534 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
535 are called, they will automatically call these functions. Note we
536 choose not to provide a couple of these functions since they are
537 not needed.
538 */
539 pc->ops->apply = PCApply_Jacobi;
540 pc->ops->applytranspose = PCApply_Jacobi;
541 pc->ops->setup = PCSetUp_Jacobi;
542 pc->ops->reset = PCReset_Jacobi;
543 pc->ops->destroy = PCDestroy_Jacobi;
544 pc->ops->setfromoptions = PCSetFromOptions_Jacobi;
545 pc->ops->view = PCView_Jacobi;
546 pc->ops->applyrichardson = NULL;
547 pc->ops->applysymmetricleft = PCApplySymmetricLeftOrRight_Jacobi;
548 pc->ops->applysymmetricright = PCApplySymmetricLeftOrRight_Jacobi;
549
550 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetType_C", PCJacobiSetType_Jacobi));
551 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetType_C", PCJacobiGetType_Jacobi));
552 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetRowl1Scale_C", PCJacobiSetRowl1Scale_Jacobi));
553 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetRowl1Scale_C", PCJacobiGetRowl1Scale_Jacobi));
554 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetUseAbs_C", PCJacobiSetUseAbs_Jacobi));
555 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetUseAbs_C", PCJacobiGetUseAbs_Jacobi));
556 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiSetFixDiagonal_C", PCJacobiSetFixDiagonal_Jacobi));
557 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetFixDiagonal_C", PCJacobiGetFixDiagonal_Jacobi));
558 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCJacobiGetDiagonal_C", PCJacobiGetDiagonal_Jacobi));
559 PetscFunctionReturn(PETSC_SUCCESS);
560 }
561
562 /*@
563 PCJacobiSetUseAbs - Causes the Jacobi preconditioner `PCJACOBI` to use the
564 absolute values of the diagonal divisors in the preconditioner
565
566 Logically Collective
567
568 Input Parameters:
569 + pc - the preconditioner context
570 - flg - whether to use absolute values or not
571
572 Options Database Key:
573 . -pc_jacobi_abs <bool> - use absolute values
574
575 Note:
576 This takes affect at the next construction of the preconditioner
577
578 Level: intermediate
579
580 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiGetUseAbs()`
581 @*/
PCJacobiSetUseAbs(PC pc,PetscBool flg)582 PetscErrorCode PCJacobiSetUseAbs(PC pc, PetscBool flg)
583 {
584 PetscFunctionBegin;
585 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
586 PetscTryMethod(pc, "PCJacobiSetUseAbs_C", (PC, PetscBool), (pc, flg));
587 PetscFunctionReturn(PETSC_SUCCESS);
588 }
589
590 /*@
591 PCJacobiGetUseAbs - Determines if the Jacobi preconditioner `PCJACOBI` uses the
592 absolute values of the diagonal divisors in the preconditioner
593
594 Logically Collective
595
596 Input Parameter:
597 . pc - the preconditioner context
598
599 Output Parameter:
600 . flg - whether to use absolute values or not
601
602 Level: intermediate
603
604 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiSetUseAbs()`, `PCJacobiGetType()`
605 @*/
PCJacobiGetUseAbs(PC pc,PetscBool * flg)606 PetscErrorCode PCJacobiGetUseAbs(PC pc, PetscBool *flg)
607 {
608 PetscFunctionBegin;
609 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
610 PetscUseMethod(pc, "PCJacobiGetUseAbs_C", (PC, PetscBool *), (pc, flg));
611 PetscFunctionReturn(PETSC_SUCCESS);
612 }
613
614 /*@
615 PCJacobiSetRowl1Scale - Set scaling of off-diagonal of operator when computing l1 row norms, eg,
616 Remark 6.1 in "Multigrid Smoothers for Ultraparallel Computing", Baker et al, with 0.5 scaling
617
618 Logically Collective
619
620 Input Parameters:
621 + pc - the preconditioner context
622 - scale - scaling
623
624 Options Database Key:
625 . -pc_jacobi_rowl1_scale <real> - use absolute values
626
627 Level: intermediate
628
629 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiGetRowl1Scale()`
630 @*/
PCJacobiSetRowl1Scale(PC pc,PetscReal scale)631 PetscErrorCode PCJacobiSetRowl1Scale(PC pc, PetscReal scale)
632 {
633 PetscFunctionBegin;
634 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
635 PetscTryMethod(pc, "PCJacobiSetRowl1Scale_C", (PC, PetscReal), (pc, scale));
636 PetscFunctionReturn(PETSC_SUCCESS);
637 }
638
639 /*@
640 PCJacobiGetRowl1Scale - Get scaling of off-diagonal elements summed into l1-norm diagonal
641
642 Logically Collective
643
644 Input Parameter:
645 . pc - the preconditioner context
646
647 Output Parameter:
648 . scale - scaling
649
650 Level: intermediate
651
652 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiSetRowl1Scale()`, `PCJacobiGetType()`
653 @*/
PCJacobiGetRowl1Scale(PC pc,PetscReal * scale)654 PetscErrorCode PCJacobiGetRowl1Scale(PC pc, PetscReal *scale)
655 {
656 PetscFunctionBegin;
657 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
658 PetscUseMethod(pc, "PCJacobiGetRowl1Scale_C", (PC, PetscReal *), (pc, scale));
659 PetscFunctionReturn(PETSC_SUCCESS);
660 }
661
662 /*@
663 PCJacobiSetFixDiagonal - Check for zero values on the diagonal and replace them with 1.0
664
665 Logically Collective
666
667 Input Parameters:
668 + pc - the preconditioner context
669 - flg - the boolean flag
670
671 Options Database Key:
672 . -pc_jacobi_fixdiagonal <bool> - check for zero values on the diagonal
673
674 Note:
675 This takes affect at the next construction of the preconditioner
676
677 Level: intermediate
678
679 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiGetFixDiagonal()`, `PCJacobiSetUseAbs()`
680 @*/
PCJacobiSetFixDiagonal(PC pc,PetscBool flg)681 PetscErrorCode PCJacobiSetFixDiagonal(PC pc, PetscBool flg)
682 {
683 PetscFunctionBegin;
684 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
685 PetscTryMethod(pc, "PCJacobiSetFixDiagonal_C", (PC, PetscBool), (pc, flg));
686 PetscFunctionReturn(PETSC_SUCCESS);
687 }
688
689 /*@
690 PCJacobiGetFixDiagonal - Determines if the Jacobi preconditioner `PCJACOBI` checks for zero diagonal terms
691
692 Logically Collective
693
694 Input Parameter:
695 . pc - the preconditioner context
696
697 Output Parameter:
698 . flg - the boolean flag
699
700 Options Database Key:
701 . -pc_jacobi_fixdiagonal <bool> - Fix 0 terms on diagonal by using 1
702
703 Level: intermediate
704
705 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`, `PCJacobiSetFixDiagonal()`
706 @*/
PCJacobiGetFixDiagonal(PC pc,PetscBool * flg)707 PetscErrorCode PCJacobiGetFixDiagonal(PC pc, PetscBool *flg)
708 {
709 PetscFunctionBegin;
710 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
711 PetscUseMethod(pc, "PCJacobiGetFixDiagonal_C", (PC, PetscBool *), (pc, flg));
712 PetscFunctionReturn(PETSC_SUCCESS);
713 }
714
715 /*@
716 PCJacobiGetDiagonal - Returns copy of the diagonal and/or diagonal squareroot `Vec`
717
718 Logically Collective
719
720 Input Parameter:
721 . pc - the preconditioner context
722
723 Output Parameters:
724 + diagonal - Copy of `Vec` of the inverted diagonal
725 - diagonal_sqrt - Copy of `Vec` of the inverted square root diagonal
726
727 Level: developer
728
729 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetType()`
730 @*/
PCJacobiGetDiagonal(PC pc,Vec diagonal,Vec diagonal_sqrt)731 PetscErrorCode PCJacobiGetDiagonal(PC pc, Vec diagonal, Vec diagonal_sqrt)
732 {
733 PetscFunctionBegin;
734 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
735 PetscUseMethod(pc, "PCJacobiGetDiagonal_C", (PC, Vec, Vec), (pc, diagonal, diagonal_sqrt));
736 PetscFunctionReturn(PETSC_SUCCESS);
737 }
738
739 /*@
740 PCJacobiSetType - Causes the Jacobi preconditioner to use either the diagonal, the maximum entry in each row,
741 of the sum of rows entries for the diagonal preconditioner
742
743 Logically Collective
744
745 Input Parameters:
746 + pc - the preconditioner context
747 - type - `PC_JACOBI_DIAGONAL`, `PC_JACOBI_ROWL1`, `PC_JACOBI_ROWMAX`, `PC_JACOBI_ROWSUM`
748
749 Options Database Key:
750 . -pc_jacobi_type <diagonal,rowl1,rowmax,rowsum> - the type of diagonal matrix to use for Jacobi
751
752 Level: intermediate
753
754 Developer Notes:
755 Why is there a separate function for using the absolute value?
756
757 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetUseAbs()`, `PCJacobiGetType()`
758 @*/
PCJacobiSetType(PC pc,PCJacobiType type)759 PetscErrorCode PCJacobiSetType(PC pc, PCJacobiType type)
760 {
761 PetscFunctionBegin;
762 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
763 PetscTryMethod(pc, "PCJacobiSetType_C", (PC, PCJacobiType), (pc, type));
764 PetscFunctionReturn(PETSC_SUCCESS);
765 }
766
767 /*@
768 PCJacobiGetType - Gets how the diagonal matrix is produced for the preconditioner
769
770 Not Collective
771
772 Input Parameter:
773 . pc - the preconditioner context
774
775 Output Parameter:
776 . type - `PC_JACOBI_DIAGONAL`, `PC_JACOBI_ROWL1`, `PC_JACOBI_ROWMAX`, `PC_JACOBI_ROWSUM`
777
778 Level: intermediate
779
780 .seealso: [](ch_ksp), `PCJACOBI`, `PCJacobiSetUseAbs()`, `PCJacobiSetType()`
781 @*/
PCJacobiGetType(PC pc,PCJacobiType * type)782 PetscErrorCode PCJacobiGetType(PC pc, PCJacobiType *type)
783 {
784 PetscFunctionBegin;
785 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
786 PetscUseMethod(pc, "PCJacobiGetType_C", (PC, PCJacobiType *), (pc, type));
787 PetscFunctionReturn(PETSC_SUCCESS);
788 }
789