Actual source code: sor.c

petsc-master 2015-05-31
Report Typos and Errors
  2: /*
  3:    Defines a  (S)SOR  preconditioner for any Mat implementation
  4: */
  5: #include <petsc/private/pcimpl.h>               /*I "petscpc.h" I*/

  7: typedef struct {
  8:   PetscInt   its;         /* inner iterations, number of sweeps */
  9:   PetscInt   lits;        /* local inner iterations, number of sweeps applied by the local matrix mat->A */
 10:   MatSORType sym;         /* forward, reverse, symmetric etc. */
 11:   PetscReal  omega;
 12:   PetscReal  fshift;
 13: } PC_SOR;

 17: static PetscErrorCode PCDestroy_SOR(PC pc)
 18: {

 22:   PetscFree(pc->data);
 23:   return(0);
 24: }

 28: static PetscErrorCode PCApply_SOR(PC pc,Vec x,Vec y)
 29: {
 30:   PC_SOR         *jac = (PC_SOR*)pc->data;
 32:   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
 33:   PetscReal      fshift;
 34: 
 36:   fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
 37:   MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,fshift,jac->its,jac->lits,y);
 38:   return(0);
 39: }

 43: static PetscErrorCode PCApplyRichardson_SOR(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
 44: {
 45:   PC_SOR         *jac = (PC_SOR*)pc->data;
 47:   MatSORType     stype = jac->sym;
 48:   PetscReal      fshift;

 51:   PetscInfo1(pc,"Warning, convergence critera ignored, using %D iterations\n",its);
 52:   if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
 53:   fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
 54:   MatSOR(pc->pmat,b,jac->omega,stype,fshift,its*jac->its,jac->lits,y);
 55:   *outits = its;
 56:   *reason = PCRICHARDSON_CONVERGED_ITS;
 57:   return(0);
 58: }

 62: PetscErrorCode PCSetFromOptions_SOR(PetscOptions *PetscOptionsObject,PC pc)
 63: {
 64:   PC_SOR         *jac = (PC_SOR*)pc->data;
 66:   PetscBool      flg;

 69:   PetscOptionsHead(PetscOptionsObject,"(S)SOR options");
 70:   PetscOptionsReal("-pc_sor_omega","relaxation factor (0 < omega < 2)","PCSORSetOmega",jac->omega,&jac->omega,NULL);
 71:   PetscOptionsReal("-pc_sor_diagonal_shift","Add to the diagonal entries","",jac->fshift,&jac->fshift,NULL);
 72:   PetscOptionsInt("-pc_sor_its","number of inner SOR iterations","PCSORSetIterations",jac->its,&jac->its,NULL);
 73:   PetscOptionsInt("-pc_sor_lits","number of local inner SOR iterations","PCSORSetIterations",jac->lits,&jac->lits,NULL);
 74:   PetscOptionsBoolGroupBegin("-pc_sor_symmetric","SSOR, not SOR","PCSORSetSymmetric",&flg);
 75:   if (flg) {PCSORSetSymmetric(pc,SOR_SYMMETRIC_SWEEP);}
 76:   PetscOptionsBoolGroup("-pc_sor_backward","use backward sweep instead of forward","PCSORSetSymmetric",&flg);
 77:   if (flg) {PCSORSetSymmetric(pc,SOR_BACKWARD_SWEEP);}
 78:   PetscOptionsBoolGroup("-pc_sor_forward","use forward sweep","PCSORSetSymmetric",&flg);
 79:   if (flg) {PCSORSetSymmetric(pc,SOR_FORWARD_SWEEP);}
 80:   PetscOptionsBoolGroup("-pc_sor_local_symmetric","use SSOR separately on each processor","PCSORSetSymmetric",&flg);
 81:   if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP);}
 82:   PetscOptionsBoolGroup("-pc_sor_local_backward","use backward sweep locally","PCSORSetSymmetric",&flg);
 83:   if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_BACKWARD_SWEEP);}
 84:   PetscOptionsBoolGroupEnd("-pc_sor_local_forward","use forward sweep locally","PCSORSetSymmetric",&flg);
 85:   if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_FORWARD_SWEEP);}
 86:   PetscOptionsTail();
 87:   return(0);
 88: }

 92: PetscErrorCode PCView_SOR(PC pc,PetscViewer viewer)
 93: {
 94:   PC_SOR         *jac = (PC_SOR*)pc->data;
 95:   MatSORType     sym  = jac->sym;
 96:   const char     *sortype;
 98:   PetscBool      iascii;

101:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
102:   if (iascii) {
103:     if (sym & SOR_ZERO_INITIAL_GUESS) {PetscViewerASCIIPrintf(viewer,"  SOR:  zero initial guess\n");}
104:     if (sym == SOR_APPLY_UPPER)                                              sortype = "apply_upper";
105:     else if (sym == SOR_APPLY_LOWER)                                         sortype = "apply_lower";
106:     else if (sym & SOR_EISENSTAT)                                            sortype = "Eisenstat";
107:     else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP)             sortype = "symmetric";
108:     else if (sym & SOR_BACKWARD_SWEEP)                                       sortype = "backward";
109:     else if (sym & SOR_FORWARD_SWEEP)                                        sortype = "forward";
110:     else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
111:     else if (sym & SOR_LOCAL_FORWARD_SWEEP)                                  sortype = "local_forward";
112:     else if (sym & SOR_LOCAL_BACKWARD_SWEEP)                                 sortype = "local_backward";
113:     else                                                                     sortype = "unknown";
114:     PetscViewerASCIIPrintf(viewer,"  SOR: type = %s, iterations = %D, local iterations = %D, omega = %g\n",sortype,jac->its,jac->lits,(double)jac->omega);
115:   }
116:   return(0);
117: }


120: /* ------------------------------------------------------------------------------*/
123: static PetscErrorCode  PCSORSetSymmetric_SOR(PC pc,MatSORType flag)
124: {
125:   PC_SOR *jac = (PC_SOR*)pc->data;

128:   jac->sym = flag;
129:   return(0);
130: }

134: static PetscErrorCode  PCSORSetOmega_SOR(PC pc,PetscReal omega)
135: {
136:   PC_SOR *jac = (PC_SOR*)pc->data;

139:   if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
140:   jac->omega = omega;
141:   return(0);
142: }

146: static PetscErrorCode  PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
147: {
148:   PC_SOR *jac = (PC_SOR*)pc->data;

151:   jac->its  = its;
152:   jac->lits = lits;
153:   return(0);
154: }

158: static PetscErrorCode  PCSORGetSymmetric_SOR(PC pc,MatSORType *flag)
159: {
160:   PC_SOR *jac = (PC_SOR*)pc->data;

163:   *flag = jac->sym;
164:   return(0);
165: }

169: static PetscErrorCode  PCSORGetOmega_SOR(PC pc,PetscReal *omega)
170: {
171:   PC_SOR *jac = (PC_SOR*)pc->data;

174:   *omega = jac->omega;
175:   return(0);
176: }

180: static PetscErrorCode  PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
181: {
182:   PC_SOR *jac = (PC_SOR*)pc->data;

185:   if (its)  *its = jac->its;
186:   if (lits) *lits = jac->lits;
187:   return(0);
188: }

190: /* ------------------------------------------------------------------------------*/
193: /*@
194:    PCSORGetSymmetric - Gets the form the SOR preconditioner is using;   backward, or forward relaxation.  The local variants perform SOR on
195:    each processor.  By default forward relaxation is used.

197:    Logically Collective on PC

199:    Input Parameter:
200: .  pc - the preconditioner context

202:    Output Parameter:
203: .  flag - one of the following
204: .vb
205:     SOR_FORWARD_SWEEP
206:     SOR_BACKWARD_SWEEP
207:     SOR_SYMMETRIC_SWEEP
208:     SOR_LOCAL_FORWARD_SWEEP
209:     SOR_LOCAL_BACKWARD_SWEEP
210:     SOR_LOCAL_SYMMETRIC_SWEEP
211: .ve

213:    Options Database Keys:
214: +  -pc_sor_symmetric - Activates symmetric version
215: .  -pc_sor_backward - Activates backward version
216: .  -pc_sor_local_forward - Activates local forward version
217: .  -pc_sor_local_symmetric - Activates local symmetric version
218: -  -pc_sor_local_backward - Activates local backward version

220:    Notes:
221:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
222:    which can be chosen with the option
223: .  -pc_type eisenstat - Activates Eisenstat trick

225:    Level: intermediate

227: .keywords: PC, SOR, SSOR, set, relaxation, sweep, forward, backward, symmetric

229: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
230: @*/
231: PetscErrorCode  PCSORGetSymmetric(PC pc,MatSORType *flag)
232: {

237:   PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
238:   return(0);
239: }

243: /*@
244:    PCSORGetOmega - Gets the SOR relaxation coefficient, omega
245:    (where omega = 1.0 by default).

247:    Logically Collective on PC

249:    Input Parameter:
250: .  pc - the preconditioner context

252:    Output Parameter:
253: .  omega - relaxation coefficient (0 < omega < 2).

255:    Options Database Key:
256: .  -pc_sor_omega <omega> - Sets omega

258:    Level: intermediate

260: .keywords: PC, SOR, SSOR, set, relaxation, omega

262: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
263: @*/
264: PetscErrorCode  PCSORGetOmega(PC pc,PetscReal *omega)
265: {

270:   PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
271:   return(0);
272: }

276: /*@
277:    PCSORGetIterations - Gets the number of inner iterations to
278:    be used by the SOR preconditioner. The default is 1.

280:    Logically Collective on PC

282:    Input Parameter:
283: .  pc - the preconditioner context

285:    Output Parameter:
286: +  lits - number of local iterations, smoothings over just variables on processor
287: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

289:    Options Database Key:
290: +  -pc_sor_its <its> - Sets number of iterations
291: -  -pc_sor_lits <lits> - Sets number of local iterations

293:    Level: intermediate

295:    Notes: When run on one processor the number of smoothings is lits*its

297: .keywords: PC, SOR, SSOR, set, iterations

299: .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
300: @*/
301: PetscErrorCode  PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
302: {

307:   PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
308:   return(0);
309: }

313: /*@
314:    PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
315:    backward, or forward relaxation.  The local variants perform SOR on
316:    each processor.  By default forward relaxation is used.

318:    Logically Collective on PC

320:    Input Parameters:
321: +  pc - the preconditioner context
322: -  flag - one of the following
323: .vb
324:     SOR_FORWARD_SWEEP
325:     SOR_BACKWARD_SWEEP
326:     SOR_SYMMETRIC_SWEEP
327:     SOR_LOCAL_FORWARD_SWEEP
328:     SOR_LOCAL_BACKWARD_SWEEP
329:     SOR_LOCAL_SYMMETRIC_SWEEP
330: .ve

332:    Options Database Keys:
333: +  -pc_sor_symmetric - Activates symmetric version
334: .  -pc_sor_backward - Activates backward version
335: .  -pc_sor_local_forward - Activates local forward version
336: .  -pc_sor_local_symmetric - Activates local symmetric version
337: -  -pc_sor_local_backward - Activates local backward version

339:    Notes:
340:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
341:    which can be chosen with the option
342: .  -pc_type eisenstat - Activates Eisenstat trick

344:    Level: intermediate

346: .keywords: PC, SOR, SSOR, set, relaxation, sweep, forward, backward, symmetric

348: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
349: @*/
350: PetscErrorCode  PCSORSetSymmetric(PC pc,MatSORType flag)
351: {

357:   PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
358:   return(0);
359: }

363: /*@
364:    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
365:    (where omega = 1.0 by default).

367:    Logically Collective on PC

369:    Input Parameters:
370: +  pc - the preconditioner context
371: -  omega - relaxation coefficient (0 < omega < 2).

373:    Options Database Key:
374: .  -pc_sor_omega <omega> - Sets omega

376:    Level: intermediate

378: .keywords: PC, SOR, SSOR, set, relaxation, omega

380: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega()
381: @*/
382: PetscErrorCode  PCSORSetOmega(PC pc,PetscReal omega)
383: {

389:   PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
390:   return(0);
391: }

395: /*@
396:    PCSORSetIterations - Sets the number of inner iterations to
397:    be used by the SOR preconditioner. The default is 1.

399:    Logically Collective on PC

401:    Input Parameters:
402: +  pc - the preconditioner context
403: .  lits - number of local iterations, smoothings over just variables on processor
404: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

406:    Options Database Key:
407: +  -pc_sor_its <its> - Sets number of iterations
408: -  -pc_sor_lits <lits> - Sets number of local iterations

410:    Level: intermediate

412:    Notes: When run on one processor the number of smoothings is lits*its

414: .keywords: PC, SOR, SSOR, set, iterations

416: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
417: @*/
418: PetscErrorCode  PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
419: {

425:   PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
426:   return(0);
427: }

429: /*MC
430:      PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning

432:    Options Database Keys:
433: +  -pc_sor_symmetric - Activates symmetric version
434: .  -pc_sor_backward - Activates backward version
435: .  -pc_sor_forward - Activates forward version
436: .  -pc_sor_local_forward - Activates local forward version
437: .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
438: .  -pc_sor_local_backward - Activates local backward version
439: .  -pc_sor_omega <omega> - Sets omega
440: .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
441: .  -pc_sor_its <its> - Sets number of iterations   (default 1)
442: -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)

444:    Level: beginner

446:   Concepts: SOR, preconditioners, Gauss-Seidel

448:    Notes: Only implemented for the AIJ  and SeqBAIJ matrix formats.
449:           Not a true parallel SOR, in parallel this implementation corresponds to block
450:           Jacobi with SOR on each block.

452:           For AIJ matrix if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
453:           zero will be used and hence the KSPSolve() will terminate with KSP_DIVERGED_NANORIF. If the option
454:           KSPSetErrorIfNotConverged() or -ksp_error_if_not_converged the code will terminate as soon as it detects the 
455:           zero pivot.

457:           For SeqBAIJ matrices this implements point-block SOR, but the omega, its, lits options are not supported.

459:           For SeqBAIJ the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected 
460:           the computation is stopped with an error

462:   Developer Notes: We should add support for diagonal blocks that are singular to generate a Inf and thus cause KSPSolve()
463:            to terminate with KSP_DIVERGED_NANORIF instead of stopping the program allowing
464:            a nonlinear solver/ODE integrator to recover without stopping the program as currently happens.

466: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
467:            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
468: M*/

472: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
473: {
475:   PC_SOR         *jac;

478:   PetscNewLog(pc,&jac);

480:   pc->ops->apply           = PCApply_SOR;
481:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
482:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
483:   pc->ops->setup           = 0;
484:   pc->ops->view            = PCView_SOR;
485:   pc->ops->destroy         = PCDestroy_SOR;
486:   pc->data                 = (void*)jac;
487:   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
488:   jac->omega               = 1.0;
489:   jac->fshift              = 0.0;
490:   jac->its                 = 1;
491:   jac->lits                = 1;

493:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
494:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
495:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
496:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
497:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
498:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
499:   return(0);
500: }