Actual source code: eisen.c


  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>

 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;

 17: static PetscErrorCode PCMult_Eisenstat(Mat mat, Vec b, Vec x)
 18: {
 19:   PC            pc;
 20:   PC_Eisenstat *eis;

 22:   MatShellGetContext(mat, &pc);
 23:   eis = (PC_Eisenstat *)pc->data;
 24:   MatSOR(eis->A, b, eis->omega, SOR_EISENSTAT, 0.0, 1, 1, x);
 25:   MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason);
 26:   return 0;
 27: }

 29: static PetscErrorCode PCNorm_Eisenstat(Mat mat, NormType type, PetscReal *nrm)
 30: {
 31:   PC            pc;
 32:   PC_Eisenstat *eis;

 34:   MatShellGetContext(mat, &pc);
 35:   eis = (PC_Eisenstat *)pc->data;
 36:   MatNorm(eis->A, type, nrm);
 37:   return 0;
 38: }

 40: static PetscErrorCode PCApply_Eisenstat(PC pc, Vec x, Vec y)
 41: {
 42:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
 43:   PetscBool     hasop;

 45:   if (eis->usediag) {
 46:     MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop);
 47:     if (hasop) {
 48:       MatMultDiagonalBlock(pc->pmat, x, y);
 49:     } else {
 50:       VecPointwiseMult(y, x, eis->diag);
 51:     }
 52:   } else VecCopy(x, y);
 53:   return 0;
 54: }

 56: static PetscErrorCode PCApplyTranspose_Eisenstat(PC pc, Vec x, Vec y)
 57: {
 58:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
 59:   PetscBool     hasop, set, sym;

 61:   MatIsSymmetricKnown(eis->A, &set, &sym);
 63:   if (eis->usediag) {
 64:     MatHasOperation(pc->pmat, MATOP_MULT_DIAGONAL_BLOCK, &hasop);
 65:     if (hasop) {
 66:       MatMultDiagonalBlock(pc->pmat, x, y);
 67:     } else {
 68:       VecPointwiseMult(y, x, eis->diag);
 69:     }
 70:   } else VecCopy(x, y);
 71:   return 0;
 72: }

 74: static PetscErrorCode PCPreSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
 75: {
 76:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
 77:   PetscBool     nonzero;

 79:   if (pc->presolvedone < 2) {
 81:     /* swap shell matrix and true matrix */
 82:     eis->A  = pc->mat;
 83:     pc->mat = eis->shell;
 84:   }

 86:   if (!eis->b[pc->presolvedone - 1]) { VecDuplicate(b, &eis->b[pc->presolvedone - 1]); }

 88:   /* if nonzero initial guess, modify x */
 89:   KSPGetInitialGuessNonzero(ksp, &nonzero);
 90:   if (nonzero) {
 91:     VecCopy(x, eis->b[pc->presolvedone - 1]);
 92:     MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, SOR_APPLY_UPPER, 0.0, 1, 1, x);
 93:     MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason);
 94:   }

 96:   /* save true b, other option is to swap pointers */
 97:   VecCopy(b, eis->b[pc->presolvedone - 1]);

 99:   /* modify b by (L + D/omega)^{-1} */
100:   MatSOR(eis->A, eis->b[pc->presolvedone - 1], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP), 0.0, 1, 1, b);
101:   MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason);
102:   return 0;
103: }

105: static PetscErrorCode PCPostSolve_Eisenstat(PC pc, KSP ksp, Vec b, Vec x)
106: {
107:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;

109:   /* get back true b */
110:   VecCopy(eis->b[pc->presolvedone], b);

112:   /* modify x by (U + D/omega)^{-1} */
113:   VecCopy(x, eis->b[pc->presolvedone]);
114:   MatSOR(eis->A, eis->b[pc->presolvedone], eis->omega, (MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP), 0.0, 1, 1, x);
115:   MatFactorGetError(eis->A, (MatFactorError *)&pc->failedreason);
116:   if (!pc->presolvedone) pc->mat = eis->A;
117:   return 0;
118: }

120: static PetscErrorCode PCReset_Eisenstat(PC pc)
121: {
122:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;

124:   VecDestroy(&eis->b[0]);
125:   VecDestroy(&eis->b[1]);
126:   MatDestroy(&eis->shell);
127:   VecDestroy(&eis->diag);
128:   return 0;
129: }

131: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
132: {
133:   PCReset_Eisenstat(pc);
134:   PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", NULL);
135:   PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", NULL);
136:   PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", NULL);
137:   PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", NULL);
138:   PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", NULL);
139:   PetscFree(pc->data);
140:   return 0;
141: }

143: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc, PetscOptionItems *PetscOptionsObject)
144: {
145:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
146:   PetscBool     set, flg;
147:   PetscReal     omega;

149:   PetscOptionsHeadBegin(PetscOptionsObject, "Eisenstat SSOR options");
150:   PetscOptionsReal("-pc_eisenstat_omega", "Relaxation factor 0 < omega < 2", "PCEisenstatSetOmega", eis->omega, &omega, &flg);
151:   if (flg) PCEisenstatSetOmega(pc, omega);
152:   PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling", "Do not use standard diagonal scaling", "PCEisenstatSetNoDiagonalScaling", eis->usediag ? PETSC_FALSE : PETSC_TRUE, &flg, &set);
153:   if (set) PCEisenstatSetNoDiagonalScaling(pc, flg);
154:   PetscOptionsHeadEnd();
155:   return 0;
156: }

158: static PetscErrorCode PCView_Eisenstat(PC pc, PetscViewer viewer)
159: {
160:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
161:   PetscBool     iascii;

163:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
164:   if (iascii) {
165:     PetscViewerASCIIPrintf(viewer, "  omega = %g\n", (double)eis->omega);
166:     if (eis->usediag) {
167:       PetscViewerASCIIPrintf(viewer, "  Using diagonal scaling (default)\n");
168:     } else {
169:       PetscViewerASCIIPrintf(viewer, "  Not using diagonal scaling\n");
170:     }
171:   }
172:   return 0;
173: }

175: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
176: {
177:   PetscInt      M, N, m, n;
178:   PetscBool     set, sym;
179:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;

181:   if (!pc->setupcalled) {
182:     MatGetSize(pc->mat, &M, &N);
183:     MatGetLocalSize(pc->mat, &m, &n);
184:     MatIsSymmetricKnown(pc->mat, &set, &sym);
185:     MatCreate(PetscObjectComm((PetscObject)pc), &eis->shell);
186:     MatSetSizes(eis->shell, m, n, M, N);
187:     MatSetType(eis->shell, MATSHELL);
188:     MatSetUp(eis->shell);
189:     MatShellSetContext(eis->shell, pc);
190:     MatShellSetOperation(eis->shell, MATOP_MULT, (void (*)(void))PCMult_Eisenstat);
191:     if (set && sym) MatShellSetOperation(eis->shell, MATOP_MULT_TRANSPOSE, (void (*)(void))PCMult_Eisenstat);
192:     MatShellSetOperation(eis->shell, MATOP_NORM, (void (*)(void))PCNorm_Eisenstat);
193:   }
194:   if (!eis->usediag) return 0;
195:   if (!pc->setupcalled) { MatCreateVecs(pc->pmat, &eis->diag, NULL); }
196:   MatGetDiagonal(pc->pmat, eis->diag);
197:   return 0;
198: }

200: /* --------------------------------------------------------------------*/

202: static PetscErrorCode PCEisenstatSetOmega_Eisenstat(PC pc, PetscReal omega)
203: {
204:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;

207:   eis->omega = omega;
208:   return 0;
209: }

211: static PetscErrorCode PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc, PetscBool flg)
212: {
213:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;

215:   eis->usediag = flg;
216:   return 0;
217: }

219: static PetscErrorCode PCEisenstatGetOmega_Eisenstat(PC pc, PetscReal *omega)
220: {
221:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;

223:   *omega = eis->omega;
224:   return 0;
225: }

227: static PetscErrorCode PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc, PetscBool *flg)
228: {
229:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;

231:   *flg = eis->usediag;
232:   return 0;
233: }

235: /*@
236:    PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
237:    to use with Eisenstat's trick (where omega = 1.0 by default)

239:    Logically Collective

241:    Input Parameters:
242: +  pc - the preconditioner context
243: -  omega - relaxation coefficient (0 < omega < 2)

245:    Options Database Key:
246: .  -pc_eisenstat_omega <omega> - Sets omega

248:    Notes:
249:    The Eisenstat trick implementation of SSOR requires about 50% of the
250:    usual amount of floating point operations used for SSOR + Krylov method;
251:    however, the preconditioned problem must be solved with both left
252:    and right preconditioning.

254:    To use SSOR without the Eisenstat trick, employ the `PCSOR` preconditioner,
255:    which can be chosen with the database options
256: $    -pc_type  sor  -pc_sor_symmetric

258:    Level: intermediate

260: .seealso: `PCSORSetOmega()`, `PCEISENSTAT`
261: @*/
262: PetscErrorCode PCEisenstatSetOmega(PC pc, PetscReal omega)
263: {
266:   PetscTryMethod(pc, "PCEisenstatSetOmega_C", (PC, PetscReal), (pc, omega));
267:   return 0;
268: }

270: /*@
271:    PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner, `PCEISENSTAT`
272:    not to do additional diagonal preconditioning. For matrices with a constant
273:    along the diagonal, this may save a small amount of work.

275:    Logically Collective

277:    Input Parameters:
278: +  pc - the preconditioner context
279: -  flg - `PETSC_TRUE` turns off diagonal scaling inside the algorithm

281:    Options Database Key:
282: .  -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`

284:    Level: intermediate

286:    Note:
287:      If you use the `KSPSetDiagonalScaling()` or -ksp_diagonal_scale option then you will
288:    likely want to use this routine since it will save you some unneeded flops.

290: .seealso: `PCEisenstatSetOmega()`, `PCEISENSTAT`
291: @*/
292: PetscErrorCode PCEisenstatSetNoDiagonalScaling(PC pc, PetscBool flg)
293: {
295:   PetscTryMethod(pc, "PCEisenstatSetNoDiagonalScaling_C", (PC, PetscBool), (pc, flg));
296:   return 0;
297: }

299: /*@
300:    PCEisenstatGetOmega - Gets the SSOR relaxation coefficient, omega,
301:    to use with Eisenstat's trick (where omega = 1.0 by default).

303:    Logically Collective

305:    Input Parameter:
306: .  pc - the preconditioner context

308:    Output Parameter:
309: .  omega - relaxation coefficient (0 < omega < 2)

311:    Options Database Key:
312: .  -pc_eisenstat_omega <omega> - Sets omega

314:    Notes:
315:    The Eisenstat trick implementation of SSOR requires about 50% of the
316:    usual amount of floating point operations used for SSOR + Krylov method;
317:    however, the preconditioned problem must be solved with both left
318:    and right preconditioning.

320:    To use SSOR without the Eisenstat trick, employ the PCSOR preconditioner,
321:    which can be chosen with the database options
322: $    -pc_type  sor  -pc_sor_symmetric

324:    Level: intermediate

326: .seealso: `PCEISENSTAT`, `PCSORGetOmega()`, `PCEisenstatSetOmega()`
327: @*/
328: PetscErrorCode PCEisenstatGetOmega(PC pc, PetscReal *omega)
329: {
331:   PetscUseMethod(pc, "PCEisenstatGetOmega_C", (PC, PetscReal *), (pc, omega));
332:   return 0;
333: }

335: /*@
336:    PCEisenstatGetNoDiagonalScaling - Tells if the Eisenstat preconditioner
337:    not to do additional diagonal preconditioning. For matrices with a constant
338:    along the diagonal, this may save a small amount of work.

340:    Logically Collective

342:    Input Parameter:
343: .  pc - the preconditioner context

345:    Output Parameter:
346: .  flg - `PETSC_TRUE` means there is no diagonal scaling applied

348:    Options Database Key:
349: .  -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`

351:    Level: intermediate

353:    Note:
354:      If you use the KSPSetDiagonalScaling() or -ksp_diagonal_scale option then you will
355:    likely want to use this routine since it will save you some unneeded flops.

357: .seealso: , `PCEISENSTAT`, `PCEisenstatGetOmega()`
358: @*/
359: PetscErrorCode PCEisenstatGetNoDiagonalScaling(PC pc, PetscBool *flg)
360: {
362:   PetscUseMethod(pc, "PCEisenstatGetNoDiagonalScaling_C", (PC, PetscBool *), (pc, flg));
363:   return 0;
364: }

366: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool *change)
367: {
368:   *change = PETSC_TRUE;
369:   return 0;
370: }

372: /*MC
373:      PCEISENSTAT - An implementation of SSOR (symmetric successive over relaxation, symmetric Gauss-Seidel)
374:            preconditioning that incorporates Eisenstat's trick to reduce the amount of computation needed.

376:    Options Database Keys:
377: +  -pc_eisenstat_omega <omega> - Sets omega
378: -  -pc_eisenstat_no_diagonal_scaling - Activates `PCEisenstatSetNoDiagonalScaling()`

380:    Level: beginner

382:    Notes:
383:    Only implemented for the `MATAIJ` matrix format.

385:    Not a true parallel SOR, in parallel this implementation corresponds to block Jacobi with SOR on each block.

387:    Developer Note:
388:    Since this algorithm runs the Krylov method on a transformed linear system the implementation provides `PCPreSolve()` and `PCPostSolve()`
389:    routines that `KSP` uses to set up the transformed linear system.

391: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCEisenstatGetOmega()`,
392:           `PCEisenstatSetNoDiagonalScaling()`, `PCEisenstatSetOmega()`, `PCSOR`
393: M*/

395: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
396: {
397:   PC_Eisenstat *eis;

399:   PetscNew(&eis);

401:   pc->ops->apply           = PCApply_Eisenstat;
402:   pc->ops->applytranspose  = PCApplyTranspose_Eisenstat;
403:   pc->ops->presolve        = PCPreSolve_Eisenstat;
404:   pc->ops->postsolve       = PCPostSolve_Eisenstat;
405:   pc->ops->applyrichardson = NULL;
406:   pc->ops->setfromoptions  = PCSetFromOptions_Eisenstat;
407:   pc->ops->destroy         = PCDestroy_Eisenstat;
408:   pc->ops->reset           = PCReset_Eisenstat;
409:   pc->ops->view            = PCView_Eisenstat;
410:   pc->ops->setup           = PCSetUp_Eisenstat;

412:   pc->data     = eis;
413:   eis->omega   = 1.0;
414:   eis->b[0]    = NULL;
415:   eis->b[1]    = NULL;
416:   eis->diag    = NULL;
417:   eis->usediag = PETSC_TRUE;

419:   PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetOmega_C", PCEisenstatSetOmega_Eisenstat);
420:   PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatSetNoDiagonalScaling_C", PCEisenstatSetNoDiagonalScaling_Eisenstat);
421:   PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetOmega_C", PCEisenstatGetOmega_Eisenstat);
422:   PetscObjectComposeFunction((PetscObject)pc, "PCEisenstatGetNoDiagonalScaling_C", PCEisenstatGetNoDiagonalScaling_Eisenstat);
423:   PetscObjectComposeFunction((PetscObject)pc, "PCPreSolveChangeRHS_C", PCPreSolveChangeRHS_Eisenstat);
424:   return 0;
425: }