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