Actual source code: eisen.c

  1: #define PETSCKSP_DLL

  3: /*
  4:    Defines a  Eisenstat trick SSOR  preconditioner. This uses about 
  5:  %50 of the usual amount of floating point ops used for SSOR + Krylov 
  6:  method. But it requires actually solving the preconditioned problem 
  7:  with both left and right preconditioning. 
  8: */
 9:  #include private/pcimpl.h

 11: typedef struct {
 12:   Mat        shell,A;
 13:   Vec        b,diag;     /* temporary storage for true right hand side */
 14:   PetscReal  omega;
 15:   PetscTruth usediag;    /* indicates preconditioner should include diagonal scaling*/
 16: } PC_Eisenstat;


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

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

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

 42:   if (eis->usediag)  {VecPointwiseMult(y,x,eis->diag);}
 43:   else               {VecCopy(x,y);}
 44:   return(0);
 45: }

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

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

 62:   if (!eis->b) {
 63:     VecDuplicate(b,&eis->b);
 64:     PetscLogObjectParent(pc,eis->b);
 65:   }
 66: 
 67:   /* save true b, other option is to swap pointers */
 68:   VecCopy(b,eis->b);

 70:   /* if nonzero initial guess, modify x */
 71:   KSPGetInitialGuessNonzero(ksp,&nonzero);
 72:   if (nonzero) {
 73:     MatRelax(eis->A,x,eis->omega,SOR_APPLY_UPPER,0.0,1,1,x);
 74:   }

 76:   /* modify b by (L + D)^{-1} */
 77:     MatRelax(eis->A,b,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_FORWARD_SWEEP),0.0,1,1,b);
 78:   return(0);
 79: }

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

 89:     MatRelax(eis->A,x,eis->omega,(MatSORType)(SOR_ZERO_INITIAL_GUESS | SOR_BACKWARD_SWEEP),0.0,1,1,x);
 90:   pc->mat = eis->A;
 91:   /* get back true b */
 92:   VecCopy(eis->b,b);
 93:   return(0);
 94: }

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

104:   if (eis->b)     {VecDestroy(eis->b);}
105:   if (eis->shell) {MatDestroy(eis->shell);}
106:   if (eis->diag)  {VecDestroy(eis->diag);}
107:   PetscFree(eis);
108:   return(0);
109: }

113: static PetscErrorCode PCSetFromOptions_Eisenstat(PC pc)
114: {
115:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
117:   PetscTruth     flg;

120:   PetscOptionsHead("Eisenstat SSOR options");
121:     PetscOptionsReal("-pc_eisenstat_omega","Relaxation factor 0 < omega < 2","PCEisenstatSetOmega",eis->omega,&eis->omega,0);
122:     PetscOptionsName("-pc_eisenstat_no_diagonal_scaling","Do not use standard diagonal scaling","PCEisenstatNoDiagonalScaling",&flg);
123:     if (flg) {
124:       PCEisenstatNoDiagonalScaling(pc);
125:     }
126:   PetscOptionsTail();
127:   return(0);
128: }

132: static PetscErrorCode PCView_Eisenstat(PC pc,PetscViewer viewer)
133: {
134:   PC_Eisenstat   *eis = (PC_Eisenstat*)pc->data;
136:   PetscTruth     iascii;

139:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
140:   if (iascii) {
141:     PetscViewerASCIIPrintf(viewer,"Eisenstat: omega = %G\n",eis->omega);
142:     if (eis->usediag) {
143:       PetscViewerASCIIPrintf(viewer,"Eisenstat: Using diagonal scaling (default)\n");
144:     } else {
145:       PetscViewerASCIIPrintf(viewer,"Eisenstat: Not using diagonal scaling\n");
146:     }
147:   } else {
148:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for Eisenstat PC",((PetscObject)viewer)->type_name);
149:   }
150:   return(0);
151: }

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

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

181: /* --------------------------------------------------------------------*/

186: PetscErrorCode  PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
187: {
188:   PC_Eisenstat  *eis;

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

201: PetscErrorCode  PCEisenstatNoDiagonalScaling_Eisenstat(PC pc)
202: {
203:   PC_Eisenstat *eis;

206:   eis = (PC_Eisenstat*)pc->data;
207:   eis->usediag = PETSC_FALSE;
208:   return(0);
209: }

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

218:    Collective on PC

220:    Input Parameters:
221: +  pc - the preconditioner context
222: -  omega - relaxation coefficient (0 < omega < 2)

224:    Options Database Key:
225: .  -pc_eisenstat_omega <omega> - Sets omega

227:    Notes: 
228:    The Eisenstat trick implementation of SSOR requires about 50% of the
229:    usual amount of floating point operations used for SSOR + Krylov method;
230:    however, the preconditioned problem must be solved with both left 
231:    and right preconditioning.

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

237:    Level: intermediate

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

241: .seealso: PCSORSetOmega()
242: @*/
243: PetscErrorCode  PCEisenstatSetOmega(PC pc,PetscReal omega)
244: {
245:   PetscErrorCode ierr,(*f)(PC,PetscReal);

249:   PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatSetOmega_C",(void (**)(void))&f);
250:   if (f) {
251:     (*f)(pc,omega);
252:   }
253:   return(0);
254: }

258: /*@
259:    PCEisenstatNoDiagonalScaling - Causes the Eisenstat preconditioner
260:    not to do additional diagonal preconditioning. For matrices with a constant 
261:    along the diagonal, this may save a small amount of work.

263:    Collective on PC

265:    Input Parameter:
266: .  pc - the preconditioner context

268:    Options Database Key:
269: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()

271:    Level: intermediate

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

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

279: .seealso: PCEisenstatSetOmega()
280: @*/
281: PetscErrorCode  PCEisenstatNoDiagonalScaling(PC pc)
282: {
283:   PetscErrorCode ierr,(*f)(PC);

287:   PetscObjectQueryFunction((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",(void (**)(void))&f);
288:   if (f) {
289:     (*f)(pc);
290:   }
291:   return(0);
292: }

294: /* --------------------------------------------------------------------*/

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

300:    Options Database Keys:
301: +  -pc_eisenstat_omega <omega> - Sets omega
302: -  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()

304:    Level: beginner

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

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

312: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
313:            PCEisenstatNoDiagonalScaling(), PCEisenstatSetOmega(), PCSOR
314: M*/

319: PetscErrorCode  PCCreate_Eisenstat(PC pc)
320: {
322:   PC_Eisenstat   *eis;

325:   PetscNew(PC_Eisenstat,&eis);
326:   PetscLogObjectMemory(pc,sizeof(PC_Eisenstat));

328:   pc->ops->apply           = PCApply_Eisenstat;
329:   pc->ops->presolve        = PCPreSolve_Eisenstat;
330:   pc->ops->postsolve       = PCPostSolve_Eisenstat;
331:   pc->ops->applyrichardson = 0;
332:   pc->ops->setfromoptions  = PCSetFromOptions_Eisenstat;
333:   pc->ops->destroy         = PCDestroy_Eisenstat;
334:   pc->ops->view            = PCView_Eisenstat;
335:   pc->ops->setup           = PCSetUp_Eisenstat;

337:   pc->data           = (void*)eis;
338:   eis->omega         = 1.0;
339:   eis->b             = 0;
340:   eis->diag          = 0;
341:   eis->usediag       = PETSC_TRUE;

343:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatSetOmega_C","PCEisenstatSetOmega_Eisenstat",
344:                     PCEisenstatSetOmega_Eisenstat);
345:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCEisenstatNoDiagonalScaling_C",
346:                     "PCEisenstatNoDiagonalScaling_Eisenstat",
347:                     PCEisenstatNoDiagonalScaling_Eisenstat);
348:  return(0);
349: }