1 /* 2 Defines a (S)SOR preconditioner for any Mat implementation 3 */ 4 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 5 6 typedef struct { 7 PetscInt its; /* inner iterations, number of sweeps */ 8 PetscInt lits; /* local inner iterations, number of sweeps applied by the local matrix mat->A */ 9 MatSORType sym; /* forward, reverse, symmetric etc. */ 10 PetscReal omega; 11 PetscReal fshift; 12 } PC_SOR; 13 14 static PetscErrorCode PCDestroy_SOR(PC pc) 15 { 16 PetscFunctionBegin; 17 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", NULL)); 18 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", NULL)); 19 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", NULL)); 20 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", NULL)); 21 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", NULL)); 22 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", NULL)); 23 PetscCall(PetscFree(pc->data)); 24 PetscFunctionReturn(PETSC_SUCCESS); 25 } 26 27 static PetscErrorCode PCApply_SOR(PC pc, Vec x, Vec y) 28 { 29 PC_SOR *jac = (PC_SOR *)pc->data; 30 PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS; 31 32 PetscFunctionBegin; 33 PetscCall(MatSOR(pc->pmat, x, jac->omega, (MatSORType)flag, jac->fshift, jac->its, jac->lits, y)); 34 PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason)); 35 PetscFunctionReturn(PETSC_SUCCESS); 36 } 37 38 static PetscErrorCode PCApplyTranspose_SOR(PC pc, Vec x, Vec y) 39 { 40 PC_SOR *jac = (PC_SOR *)pc->data; 41 PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS; 42 PetscBool set, sym; 43 44 PetscFunctionBegin; 45 PetscCall(MatIsSymmetricKnown(pc->pmat, &set, &sym)); 46 PetscCheck(set && sym && (jac->sym == SOR_SYMMETRIC_SWEEP || jac->sym == SOR_LOCAL_SYMMETRIC_SWEEP), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Can only apply transpose of SOR if matrix is symmetric and sweep is symmetric"); 47 PetscCall(MatSOR(pc->pmat, x, jac->omega, (MatSORType)flag, jac->fshift, jac->its, jac->lits, y)); 48 PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason)); 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 static PetscErrorCode PCApplyRichardson_SOR(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason) 53 { 54 PC_SOR *jac = (PC_SOR *)pc->data; 55 MatSORType stype = jac->sym; 56 57 PetscFunctionBegin; 58 if (guesszero) stype = (MatSORType)(stype | SOR_ZERO_INITIAL_GUESS); 59 PetscCall(MatSOR(pc->pmat, b, jac->omega, stype, jac->fshift, its * jac->its, jac->lits, y)); 60 PetscCall(MatFactorGetError(pc->pmat, (MatFactorError *)&pc->failedreason)); 61 *outits = its; 62 *reason = PCRICHARDSON_CONVERGED_ITS; 63 PetscFunctionReturn(PETSC_SUCCESS); 64 } 65 66 static PetscErrorCode PCSetFromOptions_SOR(PC pc, PetscOptionItems PetscOptionsObject) 67 { 68 PC_SOR *jac = (PC_SOR *)pc->data; 69 PetscBool flg; 70 PetscReal omega; 71 72 PetscFunctionBegin; 73 PetscOptionsHeadBegin(PetscOptionsObject, "(S)SOR options"); 74 PetscCall(PetscOptionsReal("-pc_sor_omega", "relaxation factor (0 < omega < 2)", "PCSORSetOmega", jac->omega, &omega, &flg)); 75 if (flg) PetscCall(PCSORSetOmega(pc, omega)); 76 PetscCall(PetscOptionsReal("-pc_sor_diagonal_shift", "Add to the diagonal entries", "", jac->fshift, &jac->fshift, NULL)); 77 PetscCall(PetscOptionsInt("-pc_sor_its", "number of inner SOR iterations", "PCSORSetIterations", jac->its, &jac->its, NULL)); 78 PetscCall(PetscOptionsInt("-pc_sor_lits", "number of local inner SOR iterations", "PCSORSetIterations", jac->lits, &jac->lits, NULL)); 79 PetscCall(PetscOptionsBoolGroupBegin("-pc_sor_symmetric", "SSOR, not SOR", "PCSORSetSymmetric", &flg)); 80 if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_SYMMETRIC_SWEEP)); 81 PetscCall(PetscOptionsBoolGroup("-pc_sor_backward", "use backward sweep instead of forward", "PCSORSetSymmetric", &flg)); 82 if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_BACKWARD_SWEEP)); 83 PetscCall(PetscOptionsBoolGroup("-pc_sor_forward", "use forward sweep", "PCSORSetSymmetric", &flg)); 84 if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_FORWARD_SWEEP)); 85 PetscCall(PetscOptionsBoolGroup("-pc_sor_local_symmetric", "use SSOR separately on each processor", "PCSORSetSymmetric", &flg)); 86 if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_SYMMETRIC_SWEEP)); 87 PetscCall(PetscOptionsBoolGroup("-pc_sor_local_backward", "use backward sweep locally", "PCSORSetSymmetric", &flg)); 88 if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_BACKWARD_SWEEP)); 89 PetscCall(PetscOptionsBoolGroupEnd("-pc_sor_local_forward", "use forward sweep locally", "PCSORSetSymmetric", &flg)); 90 if (flg) PetscCall(PCSORSetSymmetric(pc, SOR_LOCAL_FORWARD_SWEEP)); 91 PetscOptionsHeadEnd(); 92 PetscFunctionReturn(PETSC_SUCCESS); 93 } 94 95 static PetscErrorCode PCView_SOR(PC pc, PetscViewer viewer) 96 { 97 PC_SOR *jac = (PC_SOR *)pc->data; 98 MatSORType sym = jac->sym; 99 const char *sortype; 100 PetscBool isascii; 101 102 PetscFunctionBegin; 103 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 104 if (isascii) { 105 if (sym & SOR_ZERO_INITIAL_GUESS) PetscCall(PetscViewerASCIIPrintf(viewer, " zero initial guess\n")); 106 if (sym == SOR_APPLY_UPPER) sortype = "apply_upper"; 107 else if (sym == SOR_APPLY_LOWER) sortype = "apply_lower"; 108 else if (sym & SOR_EISENSTAT) sortype = "Eisenstat"; 109 else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP) sortype = "symmetric"; 110 else if (sym & SOR_BACKWARD_SWEEP) sortype = "backward"; 111 else if (sym & SOR_FORWARD_SWEEP) sortype = "forward"; 112 else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric"; 113 else if (sym & SOR_LOCAL_FORWARD_SWEEP) sortype = "local_forward"; 114 else if (sym & SOR_LOCAL_BACKWARD_SWEEP) sortype = "local_backward"; 115 else sortype = "unknown"; 116 PetscCall(PetscViewerASCIIPrintf(viewer, " type = %s, iterations = %" PetscInt_FMT ", local iterations = %" PetscInt_FMT ", omega = %g\n", sortype, jac->its, jac->lits, (double)jac->omega)); 117 } 118 PetscFunctionReturn(PETSC_SUCCESS); 119 } 120 121 static PetscErrorCode PCSORSetSymmetric_SOR(PC pc, MatSORType flag) 122 { 123 PC_SOR *jac = (PC_SOR *)pc->data; 124 125 PetscFunctionBegin; 126 jac->sym = flag; 127 PetscFunctionReturn(PETSC_SUCCESS); 128 } 129 130 static PetscErrorCode PCSORSetOmega_SOR(PC pc, PetscReal omega) 131 { 132 PC_SOR *jac = (PC_SOR *)pc->data; 133 134 PetscFunctionBegin; 135 PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range"); 136 jac->omega = omega; 137 PetscFunctionReturn(PETSC_SUCCESS); 138 } 139 140 static PetscErrorCode PCSORSetIterations_SOR(PC pc, PetscInt its, PetscInt lits) 141 { 142 PC_SOR *jac = (PC_SOR *)pc->data; 143 144 PetscFunctionBegin; 145 jac->its = its; 146 jac->lits = lits; 147 PetscFunctionReturn(PETSC_SUCCESS); 148 } 149 150 static PetscErrorCode PCSORGetSymmetric_SOR(PC pc, MatSORType *flag) 151 { 152 PC_SOR *jac = (PC_SOR *)pc->data; 153 154 PetscFunctionBegin; 155 *flag = jac->sym; 156 PetscFunctionReturn(PETSC_SUCCESS); 157 } 158 159 static PetscErrorCode PCSORGetOmega_SOR(PC pc, PetscReal *omega) 160 { 161 PC_SOR *jac = (PC_SOR *)pc->data; 162 163 PetscFunctionBegin; 164 *omega = jac->omega; 165 PetscFunctionReturn(PETSC_SUCCESS); 166 } 167 168 static PetscErrorCode PCSORGetIterations_SOR(PC pc, PetscInt *its, PetscInt *lits) 169 { 170 PC_SOR *jac = (PC_SOR *)pc->data; 171 172 PetscFunctionBegin; 173 if (its) *its = jac->its; 174 if (lits) *lits = jac->lits; 175 PetscFunctionReturn(PETSC_SUCCESS); 176 } 177 178 /*@ 179 PCSORGetSymmetric - Gets the form the SOR preconditioner is using; backward, or forward relaxation. The local variants perform SOR on 180 each processor. By default forward relaxation is used. 181 182 Logically Collective 183 184 Input Parameter: 185 . pc - the preconditioner context 186 187 Output Parameter: 188 . flag - one of the following 189 .vb 190 SOR_FORWARD_SWEEP 191 SOR_BACKWARD_SWEEP 192 SOR_SYMMETRIC_SWEEP 193 SOR_LOCAL_FORWARD_SWEEP 194 SOR_LOCAL_BACKWARD_SWEEP 195 SOR_LOCAL_SYMMETRIC_SWEEP 196 .ve 197 198 Options Database Keys: 199 + -pc_sor_symmetric - Activates symmetric version 200 . -pc_sor_backward - Activates backward version 201 . -pc_sor_local_forward - Activates local forward version 202 . -pc_sor_local_symmetric - Activates local symmetric version 203 - -pc_sor_local_backward - Activates local backward version 204 205 Note: 206 To use the Eisenstat trick with SSOR, employ the `PCEISENSTAT` preconditioner, 207 which can be chosen with the option 208 . -pc_type eisenstat - Activates Eisenstat trick 209 210 Level: intermediate 211 212 .seealso: [](ch_ksp), `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()`, `PCSORSetSymmetric()` 213 @*/ 214 PetscErrorCode PCSORGetSymmetric(PC pc, MatSORType *flag) 215 { 216 PetscFunctionBegin; 217 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 218 PetscUseMethod(pc, "PCSORGetSymmetric_C", (PC, MatSORType *), (pc, flag)); 219 PetscFunctionReturn(PETSC_SUCCESS); 220 } 221 222 /*@ 223 PCSORGetOmega - Gets the SOR relaxation coefficient, omega 224 (where omega = 1.0 by default). 225 226 Logically Collective 227 228 Input Parameter: 229 . pc - the preconditioner context 230 231 Output Parameter: 232 . omega - relaxation coefficient (0 < omega < 2). 233 234 Options Database Key: 235 . -pc_sor_omega <omega> - Sets omega 236 237 Level: intermediate 238 239 .seealso: [](ch_ksp), `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `PCSORSetOmega()` 240 @*/ 241 PetscErrorCode PCSORGetOmega(PC pc, PetscReal *omega) 242 { 243 PetscFunctionBegin; 244 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 245 PetscUseMethod(pc, "PCSORGetOmega_C", (PC, PetscReal *), (pc, omega)); 246 PetscFunctionReturn(PETSC_SUCCESS); 247 } 248 249 /*@ 250 PCSORGetIterations - Gets the number of inner iterations to 251 be used by the SOR preconditioner. The default is 1. 252 253 Logically Collective 254 255 Input Parameter: 256 . pc - the preconditioner context 257 258 Output Parameters: 259 + lits - number of local iterations, smoothings over just variables on processor 260 - its - number of parallel iterations to use; each parallel iteration has lits local iterations 261 262 Options Database Keys: 263 + -pc_sor_its <its> - Sets number of iterations 264 - -pc_sor_lits <lits> - Sets number of local iterations 265 266 Level: intermediate 267 268 Note: 269 When run on one processor the number of smoothings is lits*its 270 271 .seealso: [](ch_ksp), `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()`, `PCSORSetIterations()` 272 @*/ 273 PetscErrorCode PCSORGetIterations(PC pc, PetscInt *its, PetscInt *lits) 274 { 275 PetscFunctionBegin; 276 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 277 PetscUseMethod(pc, "PCSORGetIterations_C", (PC, PetscInt *, PetscInt *), (pc, its, lits)); 278 PetscFunctionReturn(PETSC_SUCCESS); 279 } 280 281 /*@ 282 PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR), 283 backward, or forward relaxation. The local variants perform SOR on 284 each processor. By default forward relaxation is used. 285 286 Logically Collective 287 288 Input Parameters: 289 + pc - the preconditioner context 290 - flag - one of the following 291 .vb 292 SOR_FORWARD_SWEEP 293 SOR_BACKWARD_SWEEP 294 SOR_SYMMETRIC_SWEEP 295 SOR_LOCAL_FORWARD_SWEEP 296 SOR_LOCAL_BACKWARD_SWEEP 297 SOR_LOCAL_SYMMETRIC_SWEEP 298 .ve 299 300 Options Database Keys: 301 + -pc_sor_symmetric - Activates symmetric version 302 . -pc_sor_backward - Activates backward version 303 . -pc_sor_local_forward - Activates local forward version 304 . -pc_sor_local_symmetric - Activates local symmetric version 305 - -pc_sor_local_backward - Activates local backward version 306 307 Note: 308 To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner, 309 which can be chosen with the option 310 . -pc_type eisenstat - Activates Eisenstat trick 311 312 Level: intermediate 313 314 .seealso: [](ch_ksp), `PCSOR`, `PCEisenstatSetOmega()`, `PCSORSetIterations()`, `PCSORSetOmega()` 315 @*/ 316 PetscErrorCode PCSORSetSymmetric(PC pc, MatSORType flag) 317 { 318 PetscFunctionBegin; 319 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 320 PetscValidLogicalCollectiveEnum(pc, flag, 2); 321 PetscTryMethod(pc, "PCSORSetSymmetric_C", (PC, MatSORType), (pc, flag)); 322 PetscFunctionReturn(PETSC_SUCCESS); 323 } 324 325 /*@ 326 PCSORSetOmega - Sets the SOR relaxation coefficient, omega 327 (where omega = 1.0 by default). 328 329 Logically Collective 330 331 Input Parameters: 332 + pc - the preconditioner context 333 - omega - relaxation coefficient (0 < omega < 2). 334 335 Options Database Key: 336 . -pc_sor_omega <omega> - Sets omega 337 338 Level: intermediate 339 340 Note: 341 If omega != 1, you will need to set the `MAT_USE_INODE`S option to `PETSC_FALSE` on the matrix. 342 343 .seealso: [](ch_ksp), `PCSOR`, `PCSORSetSymmetric()`, `PCSORSetIterations()`, `PCEisenstatSetOmega()`, `MatSetOption()` 344 @*/ 345 PetscErrorCode PCSORSetOmega(PC pc, PetscReal omega) 346 { 347 PetscFunctionBegin; 348 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 349 PetscValidLogicalCollectiveReal(pc, omega, 2); 350 PetscTryMethod(pc, "PCSORSetOmega_C", (PC, PetscReal), (pc, omega)); 351 PetscFunctionReturn(PETSC_SUCCESS); 352 } 353 354 /*@ 355 PCSORSetIterations - Sets the number of inner iterations to 356 be used by the SOR preconditioner. The default is 1. 357 358 Logically Collective 359 360 Input Parameters: 361 + pc - the preconditioner context 362 . lits - number of local iterations, smoothings over just variables on processor 363 - its - number of parallel iterations to use; each parallel iteration has lits local iterations 364 365 Options Database Keys: 366 + -pc_sor_its <its> - Sets number of iterations 367 - -pc_sor_lits <lits> - Sets number of local iterations 368 369 Level: intermediate 370 371 Note: 372 When run on one processor the number of smoothings is lits*its 373 374 .seealso: [](ch_ksp), `PCSOR`, `PCSORSetOmega()`, `PCSORSetSymmetric()` 375 @*/ 376 PetscErrorCode PCSORSetIterations(PC pc, PetscInt its, PetscInt lits) 377 { 378 PetscFunctionBegin; 379 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 380 PetscValidLogicalCollectiveInt(pc, its, 2); 381 PetscTryMethod(pc, "PCSORSetIterations_C", (PC, PetscInt, PetscInt), (pc, its, lits)); 382 PetscFunctionReturn(PETSC_SUCCESS); 383 } 384 385 /*MC 386 PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning 387 388 Options Database Keys: 389 + -pc_sor_symmetric - Activates symmetric version 390 . -pc_sor_backward - Activates backward version 391 . -pc_sor_forward - Activates forward version 392 . -pc_sor_local_forward - Activates local forward version 393 . -pc_sor_local_symmetric - Activates local symmetric version (default version) 394 . -pc_sor_local_backward - Activates local backward version 395 . -pc_sor_omega <omega> - Sets omega 396 . -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal 397 . -pc_sor_its <its> - Sets number of iterations (default 1) 398 - -pc_sor_lits <lits> - Sets number of local iterations (default 1) 399 400 Level: beginner 401 402 Notes: 403 Only implemented for the `MATAIJ` and `MATSEQBAIJ` matrix formats. 404 405 Not a true parallel SOR, in parallel this implementation corresponds to block 406 Jacobi with SOR on each block. 407 408 For `MATAIJ` matrices if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that 409 zero will be used and hence the `KSPSolve()` will terminate with `KSP_DIVERGED_NANORINF`. If the option 410 `KSPSetErrorIfNotConverged()` or -ksp_error_if_not_converged the code will terminate as soon as it detects the 411 zero pivot. 412 413 For `MATSEQBAIJ` matrices this implements point-block SOR, but the omega, its, lits options are not supported. 414 415 For `MATSEQBAIJ` the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected 416 the computation is stopped with an error 417 418 If used with `KSPRICHARDSON` and no monitors the convergence test is skipped to improve speed, thus it always iterates 419 the maximum number of iterations you've selected for `KSP`. It is usually used in this mode as a smoother for multigrid. 420 421 If omega != 1, you will need to set the `MAT_USE_INODES` option to `PETSC_FALSE` on the matrix. 422 423 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`, 424 `PCSORSetIterations()`, `PCSORSetSymmetric()`, `PCSORSetOmega()`, `PCEISENSTAT`, `MatSetOption()` 425 M*/ 426 427 PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc) 428 { 429 PC_SOR *jac; 430 431 PetscFunctionBegin; 432 PetscCall(PetscNew(&jac)); 433 434 pc->ops->apply = PCApply_SOR; 435 pc->ops->applytranspose = PCApplyTranspose_SOR; 436 pc->ops->applyrichardson = PCApplyRichardson_SOR; 437 pc->ops->setfromoptions = PCSetFromOptions_SOR; 438 pc->ops->setup = NULL; 439 pc->ops->view = PCView_SOR; 440 pc->ops->destroy = PCDestroy_SOR; 441 pc->data = (void *)jac; 442 jac->sym = SOR_LOCAL_SYMMETRIC_SWEEP; 443 jac->omega = 1.0; 444 jac->fshift = 0.0; 445 jac->its = 1; 446 jac->lits = 1; 447 448 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetSymmetric_C", PCSORSetSymmetric_SOR)); 449 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetOmega_C", PCSORSetOmega_SOR)); 450 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORSetIterations_C", PCSORSetIterations_SOR)); 451 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetSymmetric_C", PCSORGetSymmetric_SOR)); 452 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetOmega_C", PCSORGetOmega_SOR)); 453 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSORGetIterations_C", PCSORGetIterations_SOR)); 454 PetscFunctionReturn(PETSC_SUCCESS); 455 } 456