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