Actual source code: eisen.c

  1: /*$Id: eisen.c,v 1.113 2001/08/06 21:16:28 bsmith Exp $*/

  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 src/sles/pc/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;


 19: #undef __FUNCT__  
 21: static int PCMult_Eisenstat(Mat mat,Vec b,Vec x)
 22: {
 23:   int          ierr;
 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: }

 34: #undef __FUNCT__  
 36: static int PCApply_Eisenstat(PC pc,Vec x,Vec y)
 37: {
 38:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
 39:   int          ierr;

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

 47: #undef __FUNCT__  
 49: static int PCPre_Eisenstat(PC pc,KSP ksp,Vec x,Vec b)
 50: {
 51:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
 52:   PetscTruth   nonzero;
 53:   int          ierr;

 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 |
 78:                                         SOR_FORWARD_SWEEP),0.0,1,1,b);
 79:   return(0);
 80: }

 82: #undef __FUNCT__  
 84: static int PCPost_Eisenstat(PC pc,KSP ksp,Vec x,Vec b)
 85: {
 86:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
 87:   int          ierr;

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

 98: #undef __FUNCT__  
100: static int PCDestroy_Eisenstat(PC pc)
101: {
102:   PC_Eisenstat *eis = (PC_Eisenstat *)pc->data;
103:   int          ierr;

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

113: #undef __FUNCT__  
115: static int PCSetFromOptions_Eisenstat(PC pc)
116: {
117:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
118:   int        ierr;
119:   PetscTruth flg;

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

132: #undef __FUNCT__  
134: static int PCView_Eisenstat(PC pc,PetscViewer viewer)
135: {
136:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;
137:   int          ierr;
138:   PetscTruth   isascii;

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

155: #undef __FUNCT__  
157: static int PCSetUp_Eisenstat(PC pc)
158: {
159:   int          ierr,M,N,m,n;
160:   PC_Eisenstat *eis = (PC_Eisenstat*)pc->data;

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

179: /* --------------------------------------------------------------------*/

181: EXTERN_C_BEGIN
182: #undef __FUNCT__  
184: int PCEisenstatSetOmega_Eisenstat(PC pc,PetscReal omega)
185: {
186:   PC_Eisenstat  *eis;

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

196: EXTERN_C_BEGIN
197: #undef __FUNCT__  
199: int PCEisenstatNoDiagonalScaling_Eisenstat(PC pc)
200: {
201:   PC_Eisenstat *eis;

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

210: #undef __FUNCT__  
212: /*@ 
213:    PCEisenstatSetOmega - Sets the SSOR relaxation coefficient, omega,
214:    to use with Eisenstat's trick (where omega = 1.0 by default).

216:    Collective on PC

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

222:    Options Database Key:
223: .  -pc_eisenstat_omega <omega> - Sets omega

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

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

235:    Level: intermediate

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

239: .seealso: PCSORSetOmega()
240: @*/
241: int PCEisenstatSetOmega(PC pc,PetscReal omega)
242: {
243:   int ierr,(*f)(PC,PetscReal);

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

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

261:    Collective on PC

263:    Input Parameter:
264: .  pc - the preconditioner context

266:    Options Database Key:
267: .  -pc_eisenstat_no_diagonal_scaling - Activates PCEisenstatNoDiagonalScaling()

269:    Level: intermediate

271:    Note:
272:      If you use the SLESSetDiagonalScaling() or -sles_diagonal_scale option then you will
273:    likley want to use this routine since it will save you some unneeded flops.

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

277: .seealso: PCEisenstatSetOmega()
278: @*/
279: int PCEisenstatNoDiagonalScaling(PC pc)
280: {
281:   int ierr,(*f)(PC);

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

292: /* --------------------------------------------------------------------*/

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

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

302:    Level: beginner

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

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

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

314: EXTERN_C_BEGIN
315: #undef __FUNCT__  
317: int PCCreate_Eisenstat(PC pc)
318: {
319:   int          ierr;
320:   PC_Eisenstat *eis;

323:   PetscNew(PC_Eisenstat,&eis);
324:   PetscLogObjectMemory(pc,sizeof(PC_Eisenstat));

326:   pc->ops->apply           = PCApply_Eisenstat;
327:   pc->ops->presolve        = PCPre_Eisenstat;
328:   pc->ops->postsolve       = PCPost_Eisenstat;
329:   pc->ops->applyrichardson = 0;
330:   pc->ops->setfromoptions  = PCSetFromOptions_Eisenstat;
331:   pc->ops->destroy         = PCDestroy_Eisenstat;
332:   pc->ops->view            = PCView_Eisenstat;
333:   pc->ops->setup           = PCSetUp_Eisenstat;

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

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