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:   return 0;
 26: }

 28: static PetscErrorCode PCApply_Eisenstat(PC pc,Vec x,Vec y)
 29: {
 30:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
 31:   PetscBool      hasop;

 33:   if (eis->usediag) {
 34:     MatHasOperation(pc->pmat,MATOP_MULT_DIAGONAL_BLOCK,&hasop);
 35:     if (hasop) {
 36:       MatMultDiagonalBlock(pc->pmat,x,y);
 37:     } else {
 38:       VecPointwiseMult(y,x,eis->diag);
 39:     }
 40:   } else VecCopy(x,y);
 41:   return 0;
 42: }

 44: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
 45: {
 46:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
 47:   PetscBool      nonzero;

 49:   if (pc->presolvedone < 2) {
 51:     /* swap shell matrix and true matrix */
 52:     eis->A  = pc->mat;
 53:     pc->mat = eis->shell;
 54:   }

 56:   if (!eis->b[pc->presolvedone-1]) {
 57:     VecDuplicate(b,&eis->b[pc->presolvedone-1]);
 58:     PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->b[pc->presolvedone-1]);
 59:   }

 61:   /* if nonzero initial guess, modify x */
 62:   KSPGetInitialGuessNonzero(ksp,&nonzero);
 63:   if (nonzero) {
 64:     VecCopy(x,eis->b[pc->presolvedone-1]);
 65:     MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
 66:   }

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

 71:   /* modify b by (L + D/omega)^{-1} */
 72:   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:   return 0;
 74: }

 76: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
 77: {
 78:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

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

 83:   /* modify x by (U + D/omega)^{-1} */
 84:   VecCopy(x,eis->b[pc->presolvedone]);
 85:   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:   return 0;
 88: }

 90: static PetscErrorCode PCReset_Eisenstat(PC pc)
 91: {
 92:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

 94:   VecDestroy(&eis->b[0]);
 95:   VecDestroy(&eis->b[1]);
 96:   MatDestroy(&eis->shell);
 97:   VecDestroy(&eis->diag);
 98:   return 0;
 99: }

101: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
102: {
103:   PCReset_Eisenstat(pc);
104:   PetscFree(pc->data);
105:   return 0;
106: }

108: static PetscErrorCode PCSetFromOptions_Eisenstat(PetscOptionItems *PetscOptionsObject,PC pc)
109: {
110:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
111:   PetscBool      set,flg;

113:   PetscOptionsHead(PetscOptionsObject,"Eisenstat SSOR options");
114:   PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,NULL);
115:   PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatSetNoDiagonalScaling",eis->usediag ? PETSC_FALSE : PETSC_TRUE,&flg,&set);
116:   if (set) {
117:     PCEisenstatSetNoDiagonalScaling(pc,flg);
118:   }
119:   PetscOptionsTail();
120:   return 0;
121: }

123: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
124: {
125:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
126:   PetscBool      iascii;

128:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
129:   if (iascii) {
130:     PetscViewerASCIIPrintf(viewer,"  omega = %g\n",(double)eis->omega);
131:     if (eis->usediag) {
132:       PetscViewerASCIIPrintf(viewer,"  Using diagonal scaling (default)\n");
133:     } else {
134:       PetscViewerASCIIPrintf(viewer,"  Not using diagonal scaling\n");
135:     }
136:   }
137:   return 0;
138: }

140: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
141: {
142:   PetscInt       M,N,m,n;
143:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

145:   if (!pc->setupcalled) {
146:     MatGetSize(pc->mat,&M,&N);
147:     MatGetLocalSize(pc->mat,&m,&n);
148:     MatCreate(PetscObjectComm((PetscObject)pc),&eis->shell);
149:     MatSetSizes(eis->shell,m,n,M,N);
150:     MatSetType(eis->shell,MATSHELL);
151:     MatSetUp(eis->shell);
152:     MatShellSetContext(eis->shell,pc);
153:     PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->shell);
154:     MatShellSetOperation(eis->shell,MATOP_MULT,(void (*)(void))PCMult_Eisenstat);
155:   }
156:   if (!eis->usediag) return 0;
157:   if (!pc->setupcalled) {
158:     MatCreateVecs(pc->pmat,&eis->diag,NULL);
159:     PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->diag);
160:   }
161:   MatGetDiagonal(pc->pmat,eis->diag);
162:   return 0;
163: }

165: /* --------------------------------------------------------------------*/

167: static PetscErrorCode  PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
168: {
169:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

172:   eis->omega = omega;
173:   return 0;
174: }

176: static PetscErrorCode  PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc,PetscBool flg)
177: {
178:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

180:   eis->usediag = flg;
181:   return 0;
182: }

184: static PetscErrorCode  PCEisenstatGetOmega_Eisenstat(PC pc,PetscReal *omega)
185: {
186:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

188:   *omega = eis->omega;
189:   return 0;
190: }

192: static PetscErrorCode  PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc,PetscBool *flg)
193: {
194:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

196:   *flg = eis->usediag;
197:   return 0;
198: }

200: /*@
201:    PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
202:    to use with Eisenstat's trick (where omega = 1.0 by default).

204:    Logically Collective on PC

206:    Input Parameters:
207: +  pc - the preconditioner context
208: -  omega - relaxation coefficient (0 < omega < 2)

210:    Options Database Key:
211: .  -pc_eisenstat_omega <omega> - Sets omega

213:    Notes:
214:    The Eisenstat trick implementation of SSOR requires about 50% of the
215:    usual amount of floating point operations used for SSOR + Krylov method;
216:    however, the preconditioned problem must be solved with both left
217:    and right preconditioning.

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

223:    Level: intermediate

225: .seealso: PCSORSetOmega()
226: @*/
227: PetscErrorCode  PCEisenstatSetOmega(PC pc,PetscReal omega)
228: {
231:   PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
232:   return 0;
233: }

235: /*@
236:    PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner
237:    not to do additional diagonal preconditioning. For matrices with a constant
238:    along the diagonal, this may save a small amount of work.

240:    Logically Collective on PC

242:    Input Parameters:
243: +  pc - the preconditioner context
244: -  flg - PETSC_TRUE turns off diagonal scaling inside the algorithm

246:    Options Database Key:
247: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

249:    Level: intermediate

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

255: .seealso: PCEisenstatSetOmega()
256: @*/
257: PetscErrorCode  PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)
258: {
260:   PetscTryMethod(pc,"PCEisenstatSetNoDiagonalScaling_C",(PC,PetscBool),(pc,flg));
261:   return 0;
262: }

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

268:    Logically Collective on PC

270:    Input Parameter:
271: .  pc - the preconditioner context

273:    Output Parameter:
274: .  omega - relaxation coefficient (0 < omega < 2)

276:    Options Database Key:
277: .  -pc_eisenstat_omega <omega> - Sets omega

279:    Notes:
280:    The Eisenstat trick implementation of SSOR requires about 50% of the
281:    usual amount of floating point operations used for SSOR + Krylov method;
282:    however, the preconditioned problem must be solved with both left
283:    and right preconditioning.

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

289:    Level: intermediate

291: .seealso: PCSORGetOmega(), PCEisenstatSetOmega()
292: @*/
293: PetscErrorCode  PCEisenstatGetOmega(PC pc,PetscReal *omega)
294: {
296:   PetscUseMethod(pc,"PCEisenstatGetOmega_C",(PC,PetscReal*),(pc,omega));
297:   return 0;
298: }

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

305:    Logically Collective on PC

307:    Input Parameter:
308: .  pc - the preconditioner context

310:    Output Parameter:
311: .  flg - PETSC_TRUE means there is no diagonal scaling applied

313:    Options Database Key:
314: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

316:    Level: intermediate

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

322: .seealso: PCEisenstatGetOmega()
323: @*/
324: PetscErrorCode  PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool *flg)
325: {
327:   PetscUseMethod(pc,"PCEisenstatGetNoDiagonalScaling_C",(PC,PetscBool*),(pc,flg));
328:   return 0;
329: }

331: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool* change)
332: {
333:   *change = PETSC_TRUE;
334:   return 0;
335: }

337: /* --------------------------------------------------------------------*/

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

343:    Options Database Keys:
344: +  -pc_eisenstat_omega <omega> - Sets omega
345: -  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

347:    Level: beginner

349:    Notes:
350:     Only implemented for the SeqAIJ matrix format.
351:           Not a true parallel SOR, in parallel this implementation corresponds to block
352:           Jacobi with SOR on each block.

354: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
355:            PCEisenstatSetNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
356: M*/

358: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
359: {
360:   PC_Eisenstat   *eis;

362:   PetscNewLog(pc,&eis);

364:   pc->ops->apply           = PCApply_Eisenstat;
365:   pc->ops->presolve        = PCPreSolve_Eisenstat;
366:   pc->ops->postsolve       = PCPostSolve_Eisenstat;
367:   pc->ops->applyrichardson = NULL;
368:   pc->ops->setfromoptions  = PCSetFromOptions_Eisenstat;
369:   pc->ops->destroy         = PCDestroy_Eisenstat;
370:   pc->ops->reset           = PCReset_Eisenstat;
371:   pc->ops->view            = PCView_Eisenstat;
372:   pc->ops->setup           = PCSetUp_Eisenstat;

374:   pc->data     = eis;
375:   eis->omega   = 1.0;
376:   eis->b[0]    = NULL;
377:   eis->b[1]    = NULL;
378:   eis->diag    = NULL;
379:   eis->usediag = PETSC_TRUE;

381:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetOmega_C",PCEisenstatSetOmega_Eisenstat);
382:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetNoDiagonalScaling_C",PCEisenstatSetNoDiagonalScaling_Eisenstat);
383:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetOmega_C",PCEisenstatGetOmega_Eisenstat);
384:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetNoDiagonalScaling_C",PCEisenstatGetNoDiagonalScaling_Eisenstat);
385:   PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_Eisenstat);
386:   return 0;
387: }