1 2 /* 3 Defines a Eisenstat trick SSOR preconditioner. This uses about 4 %50 of the usual amount of floating point ops used for SSOR + Krylov 5 method. But it requires actually solving the preconditioned problem 6 with both left and right preconditioning. 7 */ 8 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 9 10 typedef struct { 11 Mat shell, A; 12 Vec b[2], diag; /* temporary storage for true right hand side */ 13 PetscReal omega; 14 PetscBool usediag; /* indicates preconditioner should include diagonal scaling*/ 15 } PC_Eisenstat; 16 17 static PetscErrorCode PCMult_Eisenstat(Mat mat, Vec b, Vec x) { 18 PC pc; 19 PC_Eisenstat *eis; 20 21 PetscFunctionBegin; 22 PetscCall(MatShellGetContext(mat, &pc)); 23 eis = (PC_Eisenstat *)pc->data; 24 PetscCall(MatSOR(eis->A, b, eis->omega, SOR_EISENSTAT, 0.0, 1, 1, x)); 25 PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason)); 26 PetscFunctionReturn(0); 27 } 28 29 static PetscErrorCode PCNorm_Eisenstat(Mat mat, NormType type, PetscReal *nrm) { 30 PC pc; 31 PC_Eisenstat *eis; 32 33 PetscFunctionBegin; 34 PetscCall(MatShellGetContext(mat, &pc)); 35 eis = (PC_Eisenstat *)pc->data; 36 PetscCall(MatNorm(eis->A, type, nrm)); 37 PetscFunctionReturn(0); 38 } 39 40 static PetscErrorCode PCApply_Eisenstat(PC pc, Vec x, Vec y) { 41 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 42 PetscBool hasop; 43 44 PetscFunctionBegin; 45 if (eis->usediag) { 46 PetscCall(MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop)); 47 if (hasop) { 48 PetscCall(MatMultDiagonalBlock(pc->pmat, x, y)); 49 } else { 50 PetscCall(VecPointwiseMult(y, x, eis->diag)); 51 } 52 } else PetscCall(VecCopy(x, y)); 53 PetscFunctionReturn(0); 54 } 55 56 static PetscErrorCode PCPreSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x) { 57 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 58 PetscBool nonzero; 59 60 PetscFunctionBegin; 61 if (pc->presolvedone < 2) { 62 PetscCheck(pc->mat == pc->pmat, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot have different mat and pmat"); 63 /* swap shell matrix and true matrix */ 64 eis->A = pc->mat; 65 pc->mat = eis->shell; 66 } 67 68 if (!eis->b[pc->presolvedone - 1]) { 69 PetscCall(VecDuplicate(b, &eis->b[pc->presolvedone - 1])); 70 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)eis->b[pc->presolvedone - 1])); 71 } 72 73 /* if nonzero initial guess, modify x */ 74 PetscCall(KSPGetInitialGuessNonzero(ksp, &nonzero)); 75 if (nonzero) { 76 PetscCall(VecCopy(x, eis->b[pc->presolvedone - 1])); 77 PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, SOR_APPLY_UPPER, 0.0, 1, 1, x)); 78 PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason)); 79 } 80 81 /* save true b, other option is to swap pointers */ 82 PetscCall(VecCopy(b, eis->b[pc->presolvedone - 1])); 83 84 /* modify b by (L + D/omega)^{-1} */ 85 PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP), 0.0, 1, 1, b)); 86 PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason)); 87 PetscFunctionReturn(0); 88 } 89 90 static PetscErrorCode PCPostSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x) { 91 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 92 93 PetscFunctionBegin; 94 /* get back true b */ 95 PetscCall(VecCopy(eis->b[pc->presolvedone], b)); 96 97 /* modify x by (U + D/omega)^{-1} */ 98 PetscCall(VecCopy(x, eis->b[pc->presolvedone])); 99 PetscCall(MatSOR(eis->A, eis->b[pc->presolvedone], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP), 0.0, 1, 1, x)); 100 PetscCall(MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason)); 101 if (!pc->presolvedone) pc->mat = eis->A; 102 PetscFunctionReturn(0); 103 } 104 105 static PetscErrorCode PCReset_Eisenstat(PC pc) { 106 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 107 108 PetscFunctionBegin; 109 PetscCall(VecDestroy(&eis->b[0])); 110 PetscCall(VecDestroy(&eis->b[1])); 111 PetscCall(MatDestroy(&eis->shell)); 112 PetscCall(VecDestroy(&eis->diag)); 113 PetscFunctionReturn(0); 114 } 115 116 static PetscErrorCode PCDestroy_Eisenstat(PC pc) { 117 PetscFunctionBegin; 118 PetscCall(PCReset_Eisenstat(pc)); 119 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", NULL)); 120 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", NULL)); 121 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", NULL)); 122 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", NULL)); 123 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL)); 124 PetscCall(PetscFree(pc->data)); 125 PetscFunctionReturn(0); 126 } 127 128 static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc, PetscOptionItems *PetscOptionsObject) { 129 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 130 PetscBool set, flg; 131 132 PetscFunctionBegin; 133 PetscOptionsHeadBegin(PetscOptionsObject, "Eisenstat SSOR options"); 134 PetscCall(PetscOptionsReal("-pc_eisenstat_omega", "Relaxation factor 0 < omega < 2", "PCEisenstatSetOmega", eis->omega, &eis->omega, NULL)); 135 PetscCall(PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling", "Do not use standard diagonal scaling", "PCEisenstatSetNoDiagonalScaling", eis->usediag ? PETSC_FALSE : PETSC_TRUE, &flg, &set)); 136 if (set) PetscCall(PCEisenstatSetNoDiagonalScaling(pc, flg)); 137 PetscOptionsHeadEnd(); 138 PetscFunctionReturn(0); 139 } 140 141 static PetscErrorCode PCView_Eisenstat(PC pc, PetscViewer viewer) { 142 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 143 PetscBool iascii; 144 145 PetscFunctionBegin; 146 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 147 if (iascii) { 148 PetscCall(PetscViewerASCIIPrintf(viewer, " omega = %g\n", (double)eis->omega)); 149 if (eis->usediag) { 150 PetscCall(PetscViewerASCIIPrintf(viewer, " Using diagonal scaling (default)\n")); 151 } else { 152 PetscCall(PetscViewerASCIIPrintf(viewer, " Not using diagonal scaling\n")); 153 } 154 } 155 PetscFunctionReturn(0); 156 } 157 158 static PetscErrorCode PCSetUp_Eisenstat(PC pc) { 159 PetscInt M, N, m, n; 160 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 161 162 PetscFunctionBegin; 163 if (!pc->setupcalled) { 164 PetscCall(MatGetSize(pc->mat, &M, &N)); 165 PetscCall(MatGetLocalSize(pc->mat, &m, &n)); 166 PetscCall(MatCreate(PetscObjectComm((PetscObject)pc), &eis->shell)); 167 PetscCall(MatSetSizes(eis->shell, m, n, M, N)); 168 PetscCall(MatSetType(eis->shell, MATSHELL)); 169 PetscCall(MatSetUp(eis->shell)); 170 PetscCall(MatShellSetContext(eis->shell, pc)); 171 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)eis->shell)); 172 PetscCall(MatShellSetOperation(eis->shell, MATOP_MULT, (void (*)(void))PCMult_Eisenstat)); 173 PetscCall(MatShellSetOperation(eis->shell, MATOP_NORM, (void (*)(void))PCNorm_Eisenstat)); 174 } 175 if (!eis->usediag) PetscFunctionReturn(0); 176 if (!pc->setupcalled) { 177 PetscCall(MatCreateVecs(pc->pmat, &eis->diag, NULL)); 178 PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)eis->diag)); 179 } 180 PetscCall(MatGetDiagonal(pc->pmat, eis->diag)); 181 PetscFunctionReturn(0); 182 } 183 184 /* --------------------------------------------------------------------*/ 185 186 static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc, PetscReal omega) { 187 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 188 189 PetscFunctionBegin; 190 PetscCheck(omega > 0.0 && omega < 2.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relaxation out of range"); 191 eis->omega = omega; 192 PetscFunctionReturn(0); 193 } 194 195 static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc, PetscBool flg) { 196 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 197 198 PetscFunctionBegin; 199 eis->usediag = flg; 200 PetscFunctionReturn(0); 201 } 202 203 static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc, PetscReal *omega) { 204 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 205 206 PetscFunctionBegin; 207 *omega = eis->omega; 208 PetscFunctionReturn(0); 209 } 210 211 static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc, PetscBool *flg) { 212 PC_Eisenstat *eis = (PC_Eisenstat *)pc->data; 213 214 PetscFunctionBegin; 215 *flg = eis->usediag; 216 PetscFunctionReturn(0); 217 } 218 219 /*@ 220 PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega, 221 to use with Eisenstat's trick (where omega = 1.0 by default). 222 223 Logically Collective on PC 224 225 Input Parameters: 226 + pc - the preconditioner context 227 - omega - relaxation coefficient (0 < omega < 2) 228 229 Options Database Key: 230 . -pc_eisenstat_omega <omega> - Sets omega 231 232 Notes: 233 The Eisenstat trick implementation of SSOR requires about 50% of the 234 usual amount of floating point operations used for SSOR + Krylov method; 235 however, the preconditioned problem must be solved with both left 236 and right preconditioning. 237 238 To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner, 239 which can be chosen with the database options 240 $ -pc_type sor -pc_sor_symmetric 241 242 Level: intermediate 243 244 .seealso: `PCSORSetOmega()` 245 @*/ 246 PetscErrorCode PCEisenstatSetOmega(PC pc, PetscReal omega) { 247 PetscFunctionBegin; 248 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 249 PetscValidLogicalCollectiveReal(pc, omega, 2); 250 PetscTryMethod(pc, "PCEisenstatSetOmega_C", (PC, PetscReal), (pc, omega)); 251 PetscFunctionReturn(0); 252 } 253 254 /*@ 255 PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner 256 not to do additional diagonal preconditioning. For matrices with a constant 257 along the diagonal, this may save a small amount of work. 258 259 Logically Collective on PC 260 261 Input Parameters: 262 + pc - the preconditioner context 263 - flg - PETSC_TRUE turns off diagonal scaling inside the algorithm 264 265 Options Database Key: 266 . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling() 267 268 Level: intermediate 269 270 Note: 271 If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will 272 likely want to use this routine since it will save you some unneeded flops. 273 274 .seealso: `PCEisenstatSetOmega()` 275 @*/ 276 PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg) { 277 PetscFunctionBegin; 278 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 279 PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg)); 280 PetscFunctionReturn(0); 281 } 282 283 /*@ 284 PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega, 285 to use with Eisenstat's trick (where omega = 1.0 by default). 286 287 Logically Collective on PC 288 289 Input Parameter: 290 . pc - the preconditioner context 291 292 Output Parameter: 293 . omega - relaxation coefficient (0 < omega < 2) 294 295 Options Database Key: 296 . -pc_eisenstat_omega <omega> - Sets omega 297 298 Notes: 299 The Eisenstat trick implementation of SSOR requires about 50% of the 300 usual amount of floating point operations used for SSOR + Krylov method; 301 however, the preconditioned problem must be solved with both left 302 and right preconditioning. 303 304 To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner, 305 which can be chosen with the database options 306 $ -pc_type sor -pc_sor_symmetric 307 308 Level: intermediate 309 310 .seealso: `PCSORGetOmega()`, `PCEisenstatSetOmega()` 311 @*/ 312 PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega) { 313 PetscFunctionBegin; 314 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 315 PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega)); 316 PetscFunctionReturn(0); 317 } 318 319 /*@ 320 PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner 321 not to do additional diagonal preconditioning. For matrices with a constant 322 along the diagonal, this may save a small amount of work. 323 324 Logically Collective on PC 325 326 Input Parameter: 327 . pc - the preconditioner context 328 329 Output Parameter: 330 . flg - PETSC_TRUE means there is no diagonal scaling applied 331 332 Options Database Key: 333 . -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling() 334 335 Level: intermediate 336 337 Note: 338 If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will 339 likely want to use this routine since it will save you some unneeded flops. 340 341 .seealso: `PCEisenstatGetOmega()` 342 @*/ 343 PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg) { 344 PetscFunctionBegin; 345 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 346 PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg)); 347 PetscFunctionReturn(0); 348 } 349 350 static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change) { 351 PetscFunctionBegin; 352 *change = PETSC_TRUE; 353 PetscFunctionReturn(0); 354 } 355 356 /* --------------------------------------------------------------------*/ 357 358 /*MC 359 PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel) 360 preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed. 361 362 Options Database Keys: 363 + -pc_eisenstat_omega <omega> - Sets omega 364 - -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling() 365 366 Level: beginner 367 368 Notes: 369 Only implemented for the SeqAIJ matrix format. 370 Not a true parallel SOR, in parallel this implementation corresponds to block 371 Jacobi with SOR on each block. 372 373 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 374 `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR` 375 M*/ 376 377 PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc) { 378 PC_Eisenstat *eis; 379 380 PetscFunctionBegin; 381 PetscCall(PetscNewLog(pc, &eis)); 382 383 pc->ops->apply = PCApply_Eisenstat; 384 pc->ops->presolve = PCPreSolve_Eisenstat; 385 pc->ops->postsolve = PCPostSolve_Eisenstat; 386 pc->ops->applyrichardson = NULL; 387 pc->ops->setfromoptions = PCSetFromOptions_Eisenstat; 388 pc->ops->destroy = PCDestroy_Eisenstat; 389 pc->ops->reset = PCReset_Eisenstat; 390 pc->ops->view = PCView_Eisenstat; 391 pc->ops->setup = PCSetUp_Eisenstat; 392 393 pc->data = eis; 394 eis->omega = 1.0; 395 eis->b[0] = NULL; 396 eis->b[1] = NULL; 397 eis->diag = NULL; 398 eis->usediag = PETSC_TRUE; 399 400 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", PCEisenstatSetOmega_Eisenstat)); 401 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", PCEisenstatSetNoDiagonalScaling_Eisenstat)); 402 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", PCEisenstatGetOmega_Eisenstat)); 403 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", PCEisenstatGetNoDiagonalScaling_Eisenstat)); 404 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_Eisenstat)); 405 PetscFunctionReturn(0); 406 } 407