Actual source code: eisen.c

petsc-3.8.4 2018-03-24
Report Typos and Errors

  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;


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

 25:   MatShellGetContext(mat,(void**)&pc);
 26:   eis  = (PC_Eisenstat*)pc->data;
 27:   MatSOR(eis->A,b,eis->omega,SOR_EISENSTAT,0.0,1,1,x);
 28:   return(0);
 29: }

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

 38:   if (eis->usediag) {
 39:     MatHasOperation(pc->pmat,MATOP_MULT_DIAGONAL_BLOCK,&hasop);
 40:     if (hasop) {
 41:       MatMultDiagonalBlock(pc->pmat,x,y);
 42:     } else {
 43:       VecPointwiseMult(y,x,eis->diag);
 44:     }
 45:   } else {VecCopy(x,y);}
 46:   return(0);
 47: }

 49: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
 50: {
 51:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
 52:   PetscBool      nonzero;

 56:   if (pc->presolvedone < 2) {
 57:     if (pc->mat != pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot have different mat and pmat");
 58:     /* swap shell matrix and true matrix */
 59:     eis->A  = pc->mat;
 60:     pc->mat = eis->shell;
 61:   }

 63:   if (!eis->b[pc->presolvedone-1]) {
 64:     VecDuplicate(b,&eis->b[pc->presolvedone-1]);
 65:     PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->b[pc->presolvedone-1]);
 66:   }

 68:   /* if nonzero initial guess, modify x */
 69:   KSPGetInitialGuessNonzero(ksp,&nonzero);
 70:   if (nonzero) {
 71:     VecCopy(x,eis->b[pc->presolvedone-1]);
 72:     MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
 73:   }

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

 78:   /* modify b by (L + D/omega)^{-1} */
 79:     MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_FORWARD_SWEEP),0.0,1,1,b);
 80:   return(0);
 81: }

 83: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
 84: {
 85:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

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

 92:   /* modify x by (U + D/omega)^{-1} */
 93:   VecCopy(x,eis->b[pc->presolvedone]);
 94:   MatSOR(eis->A,eis->b[pc->presolvedone],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP),0.0,1,1,x);
 95:   if (!pc->presolvedone) pc->mat = eis->A;
 96:   return(0);
 97: }

 99: static PetscErrorCode PCReset_Eisenstat(PC pc)
100: {
101:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

105:   VecDestroy(&eis->b[0]);
106:   VecDestroy(&eis->b[1]);
107:   MatDestroy(&eis->shell);
108:   VecDestroy(&eis->diag);
109:   return(0);
110: }

112: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
113: {

117:   PCReset_Eisenstat(pc);
118:   PetscFree(pc->data);
119:   return(0);
120: }

122: static PetscErrorCode PCSetFromOptions_Eisenstat(PetscOptionItems *PetscOptionsObject,PC pc)
123: {
124:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
126:   PetscBool      set,flg;

129:   PetscOptionsHead(PetscOptionsObject,"Eisenstat SSOR options");
130:   PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,NULL);
131:   PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatSetNoDiagonalScaling",eis->usediag ? PETSC_FALSE : PETSC_TRUE,&flg,&set);
132:   if (set) {
133:     PCEisenstatSetNoDiagonalScaling(pc,flg);
134:   }
135:   PetscOptionsTail();
136:   return(0);
137: }

139: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
140: {
141:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
143:   PetscBool      iascii;

146:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
147:   if (iascii) {
148:     PetscViewerASCIIPrintf(viewer,"  omega = %g\n",(double)eis->omega);
149:     if (eis->usediag) {
150:       PetscViewerASCIIPrintf(viewer,"  Using diagonal scaling (default)\n");
151:     } else {
152:       PetscViewerASCIIPrintf(viewer,"  Not using diagonal scaling\n");
153:     }
154:   }
155:   return(0);
156: }

158: static PetscErrorCode PCSetUp_Eisenstat(PC pc)
159: {
161:   PetscInt       M,N,m,n;
162:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

165:   if (!pc->setupcalled) {
166:     MatGetSize(pc->mat,&M,&N);
167:     MatGetLocalSize(pc->mat,&m,&n);
168:     MatCreate(PetscObjectComm((PetscObject)pc),&eis->shell);
169:     MatSetSizes(eis->shell,m,n,M,N);
170:     MatSetType(eis->shell,MATSHELL);
171:     MatSetUp(eis->shell);
172:     MatShellSetContext(eis->shell,(void*)pc);
173:     PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->shell);
174:     MatShellSetOperation(eis->shell,MATOP_MULT,(void (*)(void))PCMult_Eisenstat);
175:   }
176:   if (!eis->usediag) return(0);
177:   if (!pc->setupcalled) {
178:     MatCreateVecs(pc->pmat,&eis->diag,0);
179:     PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->diag);
180:   }
181:   MatGetDiagonal(pc->pmat,eis->diag);
182:   return(0);
183: }

185: /* --------------------------------------------------------------------*/

187: static PetscErrorCode  PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
188: {
189:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

192:   if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
193:   eis->omega = omega;
194:   return(0);
195: }

197: static PetscErrorCode  PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc,PetscBool flg)
198: {
199:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

202:   eis->usediag = flg;
203:   return(0);
204: }

206: static PetscErrorCode  PCEisenstatGetOmega_Eisenstat(PC pc,PetscReal *omega)
207: {
208:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

211:   *omega = eis->omega;
212:   return(0);
213: }

215: static PetscErrorCode  PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc,PetscBool *flg)
216: {
217:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

220:   *flg = eis->usediag;
221:   return(0);
222: }

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

228:    Logically Collective on PC

230:    Input Parameters:
231: +  pc - the preconditioner context
232: -  omega - relaxation coefficient (0 < omega < 2)

234:    Options Database Key:
235: .  -pc_eisenstat_omega <omega> - Sets omega

237:    Notes:
238:    The Eisenstat trick implementation of SSOR requires about 50% of the
239:    usual amount of floating point operations used for SSOR + Krylov method;
240:    however, the preconditioned problem must be solved with both left
241:    and right preconditioning.

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

247:    Level: intermediate

249: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega

251: .seealso: PCSORSetOmega()
252: @*/
253: PetscErrorCode  PCEisenstatSetOmega(PC pc,PetscReal omega)
254: {

260:   PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
261:   return(0);
262: }

264: /*@
265:    PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner
266:    not to do additional diagonal preconditioning. For matrices with a constant
267:    along the diagonal, this may save a small amount of work.

269:    Logically Collective on PC

271:    Input Parameters:
272: +  pc - the preconditioner context
273: -  flg - PETSC_TRUE turns off diagonal scaling inside the algorithm  

275:    Options Database Key:
276: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

278:    Level: intermediate

280:    Note:
281:      If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
282:    likley want to use this routine since it will save you some unneeded flops.

284: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR

286: .seealso: PCEisenstatSetOmega()
287: @*/
288: PetscErrorCode  PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)
289: {

294:   PetscTryMethod(pc,"PCEisenstatSetNoDiagonalScaling_C",(PC,PetscBool),(pc,flg));
295:   return(0);
296: }

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

302:    Logically Collective on PC

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

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

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

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

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

323:    Level: intermediate

325: .keywords: PC, Eisenstat, set, SOR, SSOR, relaxation, omega

327: .seealso: PCSORGetOmega(), PCEisenstatSetOmega()
328: @*/
329: PetscErrorCode  PCEisenstatGetOmega(PC pc,PetscReal *omega)
330: {

335:   PetscUseMethod(pc,"PCEisenstatGetOmega_C",(PC,PetscReal*),(pc,omega));
336:   return(0);
337: }

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

344:    Logically Collective on PC

346:    Input Parameter:
347: .  pc - the preconditioner context

349:    Output Parameter:
350: .  flg - PETSC_TRUE means there is no diagonal scaling applied

352:    Options Database Key:
353: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

355:    Level: intermediate

357:    Note:
358:      If you use the KPSSetDiagonalScaling() or -ksp_diagonal_scale option then you will
359:    likley want to use this routine since it will save you some unneeded flops.

361: .keywords: PC, Eisenstat, use, diagonal, scaling, SSOR

363: .seealso: PCEisenstatGetOmega()
364: @*/
365: PetscErrorCode  PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool *flg)
366: {

371:   PetscUseMethod(pc,"PCEisenstatGetNoDiagonalScaling_C",(PC,PetscBool*),(pc,flg));
372:   return(0);
373: }

375: static PetscErrorCode PCPreSolveChangeRHS_Eisenstat(PC pc, PetscBool* change)
376: {
378:   *change = PETSC_TRUE;
379:   return(0);
380: }

382: /* --------------------------------------------------------------------*/

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

388:    Options Database Keys:
389: +  -pc_eisenstat_omega <omega> - Sets omega
390: -  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

392:    Level: beginner

394:   Concepts: SOR, preconditioners, Gauss-Seidel, Eisenstat's trick

396:    Notes: Only implemented for the SeqAIJ matrix format.
397:           Not a true parallel SOR, in parallel this implementation corresponds to block
398:           Jacobi with SOR on each block.

400: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
401:            PCEisenstatSetNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
402: M*/

404: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
405: {
407:   PC_Eisenstat   *eis;

410:   PetscNewLog(pc,&eis);

412:   pc->ops->apply           = PCApply_Eisenstat;
413:   pc->ops->presolve        = PCPreSolve_Eisenstat;
414:   pc->ops->postsolve       = PCPostSolve_Eisenstat;
415:   pc->ops->applyrichardson = 0;
416:   pc->ops->setfromoptions  = PCSetFromOptions_Eisenstat;
417:   pc->ops->destroy         = PCDestroy_Eisenstat;
418:   pc->ops->reset           = PCReset_Eisenstat;
419:   pc->ops->view            = PCView_Eisenstat;
420:   pc->ops->setup           = PCSetUp_Eisenstat;

422:   pc->data     = (void*)eis;
423:   eis->omega   = 1.0;
424:   eis->b[0]    = 0;
425:   eis->b[1]    = 0;
426:   eis->diag    = 0;
427:   eis->usediag = PETSC_TRUE;

429:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetOmega_C",PCEisenstatSetOmega_Eisenstat);
430:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetNoDiagonalScaling_C",PCEisenstatSetNoDiagonalScaling_Eisenstat);
431:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetOmega_C",PCEisenstatGetOmega_Eisenstat);
432:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetNoDiagonalScaling_C",PCEisenstatGetNoDiagonalScaling_Eisenstat);
433:   PetscObjectComposeFunction((PetscObject)pc,"PCPreSolveChangeRHS_C",PCPreSolveChangeRHS_Eisenstat);
434:   return(0);
435: }