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