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