Actual source code: eisen.c

petsc-3.7.3 2016-08-01
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>           /*I "petscpc.h" I*/

 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;


 20: static PetscErrorCode PCMult_Eisenstat(Mat mat,Vec b,Vec x)
 21: {
 23:   PC             pc;
 24:   PC_Eisenstat   *eis;

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

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

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

 55: static PetscErrorCode PCPreSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
 56: {
 57:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
 58:   PetscBool      nonzero;

 62:   if (pc->presolvedone < 2) {
 63:     if (pc->mat != pc->pmat) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot have different mat and pmat");
 64:     /* swap shell matrix and true matrix */
 65:     eis->A  = pc->mat;
 66:     pc->mat = eis->shell;
 67:   }

 69:   if (!eis->b[pc->presolvedone-1]) {
 70:     VecDuplicate(b,&eis->b[pc->presolvedone-1]);
 71:     PetscLogObjectParent((PetscObject)pc,(PetscObject)eis->b[pc->presolvedone-1]);
 72:   }

 74:   /* if nonzero initial guess, modify x */
 75:   KSPGetInitialGuessNonzero(ksp,&nonzero);
 76:   if (nonzero) {
 77:     VecCopy(x,eis->b[pc->presolvedone-1]);
 78:     MatSOR(eis->A,eis->b[pc->presolvedone-1],eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
 79:   }

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

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

 91: static PetscErrorCode PCPostSolve_Eisenstat(PC pc,KSP ksp,Vec b,Vec x)
 92: {
 93:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

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

100:   /* modify x by (U + D/omega)^{-1} */
101:   VecCopy(x,eis->b[pc->presolvedone]);
102:   MatSOR(eis->A,eis->b[pc->presolvedone],eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_LOCAL_BACKWARD_SWEEP),0.0,1,1,x);
103:   if (!pc->presolvedone) pc->mat = eis->A;
104:   return(0);
105: }

109: static PetscErrorCode PCReset_Eisenstat(PC pc)
110: {
111:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;

115:   VecDestroy(&eis->b[0]);
116:   VecDestroy(&eis->b[1]);
117:   MatDestroy(&eis->shell);
118:   VecDestroy(&eis->diag);
119:   return(0);
120: }

124: static PetscErrorCode PCDestroy_Eisenstat(PC pc)
125: {

129:   PCReset_Eisenstat(pc);
130:   PetscFree(pc->data);
131:   return(0);
132: }

136: static PetscErrorCode PCSetFromOptions_Eisenstat(PetscOptionItems *PetscOptionsObject,PC pc)
137: {
138:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
140:   PetscBool      set,flg;

143:   PetscOptionsHead(PetscOptionsObject,"Eisenstat SSOR options");
144:   PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,NULL);
145:   PetscOptionsBool("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatSetNoDiagonalScaling",eis->usediag ? PETSC_FALSE : PETSC_TRUE,&flg,&set);
146:   if (set) {
147:     PCEisenstatSetNoDiagonalScaling(pc,flg);
148:   }
149:   PetscOptionsTail();
150:   return(0);
151: }

155: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
156: {
157:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
159:   PetscBool      iascii;

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

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

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

203: /* --------------------------------------------------------------------*/

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

212:   if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
213:   eis->omega = omega;
214:   return(0);
215: }

219: static PetscErrorCode  PCEisenstatSetNoDiagonalScaling_Eisenstat(PC pc,PetscBool flg)
220: {
221:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

224:   eis->usediag = flg;
225:   return(0);
226: }

230: static PetscErrorCode  PCEisenstatGetOmega_Eisenstat(PC pc,PetscReal *omega)
231: {
232:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

235:   *omega = eis->omega;
236:   return(0);
237: }

241: static PetscErrorCode  PCEisenstatGetNoDiagonalScaling_Eisenstat(PC pc,PetscBool *flg)
242: {
243:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

246:   *flg = eis->usediag;
247:   return(0);
248: }

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

256:    Logically Collective on PC

258:    Input Parameters:
259: +  pc - the preconditioner context
260: -  omega - relaxation coefficient (0 < omega < 2)

262:    Options Database Key:
263: .  -pc_eisenstat_omega <omega> - Sets omega

265:    Notes:
266:    The Eisenstat trick implementation of SSOR requires about 50% of the
267:    usual amount of floating point operations used for SSOR + Krylov method;
268:    however, the preconditioned problem must be solved with both left
269:    and right preconditioning.

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

275:    Level: intermediate

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

279: .seealso: PCSORSetOmega()
280: @*/
281: PetscErrorCode  PCEisenstatSetOmega(PC pc,PetscReal omega)
282: {

288:   PetscTryMethod(pc,"PCEisenstatSetOmega_C",(PC,PetscReal),(pc,omega));
289:   return(0);
290: }

294: /*@
295:    PCEisenstatSetNoDiagonalScaling - Causes the Eisenstat preconditioner
296:    not to do additional diagonal preconditioning. For matrices with a constant
297:    along the diagonal, this may save a small amount of work.

299:    Logically Collective on PC

301:    Input Parameters:
302: +  pc - the preconditioner context
303: -  flg - PETSC_TRUE turns off diagonal scaling inside the algorithm  

305:    Options Database Key:
306: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

308:    Level: intermediate

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

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

316: .seealso: PCEisenstatSetOmega()
317: @*/
318: PetscErrorCode  PCEisenstatSetNoDiagonalScaling(PC pc,PetscBool flg)
319: {

324:   PetscTryMethod(pc,"PCEisenstatSetNoDiagonalScaling_C",(PC,PetscBool),(pc,flg));
325:   return(0);
326: }

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

334:    Logically Collective on PC

336:    Input Parameter:
337: .  pc - the preconditioner context

339:    Output Parameter:
340: .  omega - relaxation coefficient (0 < omega < 2)

342:    Options Database Key:
343: .  -pc_eisenstat_omega <omega> - Sets omega

345:    Notes:
346:    The Eisenstat trick implementation of SSOR requires about 50% of the
347:    usual amount of floating point operations used for SSOR + Krylov method;
348:    however, the preconditioned problem must be solved with both left
349:    and right preconditioning.

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

355:    Level: intermediate

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

359: .seealso: PCSORGetOmega(), PCEisenstatSetOmega()
360: @*/
361: PetscErrorCode  PCEisenstatGetOmega(PC pc,PetscReal *omega)
362: {

367:   PetscUseMethod(pc,"PCEisenstatGetOmega_C",(PC,PetscReal*),(pc,omega));
368:   return(0);
369: }

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

378:    Logically Collective on PC

380:    Input Parameter:
381: .  pc - the preconditioner context

383:    Output Parameter:
384: .  flg - PETSC_TRUE means there is no diagonal scaling applied

386:    Options Database Key:
387: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

389:    Level: intermediate

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

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

397: .seealso: PCEisenstatGetOmega()
398: @*/
399: PetscErrorCode  PCEisenstatGetNoDiagonalScaling(PC pc,PetscBool *flg)
400: {

405:   PetscUseMethod(pc,"PCEisenstatGetNoDiagonalScaling_C",(PC,PetscBool*),(pc,flg));
406:   return(0);
407: }

409: /* --------------------------------------------------------------------*/

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

415:    Options Database Keys:
416: +  -pc_eisenstat_omega <omega> - Sets omega
417: -  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatSetNoDiagonalScaling()

419:    Level: beginner

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

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

427: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
428:            PCEisenstatSetNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
429: M*/

433: PETSC_EXTERN PetscErrorCode PCCreate_Eisenstat(PC pc)
434: {
436:   PC_Eisenstat   *eis;

439:   PetscNewLog(pc,&eis);

441:   pc->ops->apply           = PCApply_Eisenstat;
442:   pc->ops->presolve        = PCPreSolve_Eisenstat;
443:   pc->ops->postsolve       = PCPostSolve_Eisenstat;
444:   pc->ops->applyrichardson = 0;
445:   pc->ops->setfromoptions  = PCSetFromOptions_Eisenstat;
446:   pc->ops->destroy         = PCDestroy_Eisenstat;
447:   pc->ops->reset           = PCReset_Eisenstat;
448:   pc->ops->view            = PCView_Eisenstat;
449:   pc->ops->setup           = PCSetUp_Eisenstat;

451:   pc->data     = (void*)eis;
452:   eis->omega   = 1.0;
453:   eis->b[0]    = 0;
454:   eis->b[1]    = 0;
455:   eis->diag    = 0;
456:   eis->usediag = PETSC_TRUE;

458:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetOmega_C",PCEisenstatSetOmega_Eisenstat);
459:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatSetNoDiagonalScaling_C",PCEisenstatSetNoDiagonalScaling_Eisenstat);
460:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetOmega_C",PCEisenstatGetOmega_Eisenstat);
461:   PetscObjectComposeFunction((PetscObject)pc,"PCEisenstatGetNoDiagonalScaling_C",PCEisenstatGetNoDiagonalScaling_Eisenstat);
462:   return(0);
463: }