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 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 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 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 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 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 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 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 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 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 */ 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 */ 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 */ 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 */ 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 */ 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 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 */ 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 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 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 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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 @*/ 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