Actual source code: sor.c

petsc-master 2016-12-08
Report Typos and Errors

  2: /*
  3:    Defines a  (S)SOR  preconditioner for any Mat implementation
  4: */
  5:  #include <petsc/private/pcimpl.h>

  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;

 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:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 39:   return(0);
 40: }

 44: static PetscErrorCode PCApplyTranspose_SOR(PC pc,Vec x,Vec y)
 45: {
 46:   PC_SOR         *jac = (PC_SOR*)pc->data;
 48:   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
 49:   PetscReal      fshift;
 50:   PetscBool      set,sym;

 53:   MatIsSymmetricKnown(pc->pmat,&set,&sym);
 54:   if (!set || !sym || (jac->sym != SOR_SYMMETRIC_SWEEP && jac->sym != SOR_LOCAL_SYMMETRIC_SWEEP)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Can only apply transpose of SOR if matrix is symmetric and sweep is symmetric");
 55:   fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
 56:   MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,fshift,jac->its,jac->lits,y);
 57:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 58:   return(0);
 59: }

 63: 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)
 64: {
 65:   PC_SOR         *jac = (PC_SOR*)pc->data;
 67:   MatSORType     stype = jac->sym;
 68:   PetscReal      fshift;

 71:   PetscInfo1(pc,"Warning, convergence critera ignored, using %D iterations\n",its);
 72:   if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
 73:   fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
 74:   MatSOR(pc->pmat,b,jac->omega,stype,fshift,its*jac->its,jac->lits,y);
 75:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 76:   *outits = its;
 77:   *reason = PCRICHARDSON_CONVERGED_ITS;
 78:   return(0);
 79: }

 83: PetscErrorCode PCSetFromOptions_SOR(PetscOptionItems *PetscOptionsObject,PC pc)
 84: {
 85:   PC_SOR         *jac = (PC_SOR*)pc->data;
 87:   PetscBool      flg;

 90:   PetscOptionsHead(PetscOptionsObject,"(S)SOR options");
 91:   PetscOptionsReal("-pc_sor_omega","relaxation factor (0 < omega < 2)","PCSORSetOmega",jac->omega,&jac->omega,NULL);
 92:   PetscOptionsReal("-pc_sor_diagonal_shift","Add to the diagonal entries","",jac->fshift,&jac->fshift,NULL);
 93:   PetscOptionsInt("-pc_sor_its","number of inner SOR iterations","PCSORSetIterations",jac->its,&jac->its,NULL);
 94:   PetscOptionsInt("-pc_sor_lits","number of local inner SOR iterations","PCSORSetIterations",jac->lits,&jac->lits,NULL);
 95:   PetscOptionsBoolGroupBegin("-pc_sor_symmetric","SSOR, not SOR","PCSORSetSymmetric",&flg);
 96:   if (flg) {PCSORSetSymmetric(pc,SOR_SYMMETRIC_SWEEP);}
 97:   PetscOptionsBoolGroup("-pc_sor_backward","use backward sweep instead of forward","PCSORSetSymmetric",&flg);
 98:   if (flg) {PCSORSetSymmetric(pc,SOR_BACKWARD_SWEEP);}
 99:   PetscOptionsBoolGroup("-pc_sor_forward","use forward sweep","PCSORSetSymmetric",&flg);
100:   if (flg) {PCSORSetSymmetric(pc,SOR_FORWARD_SWEEP);}
101:   PetscOptionsBoolGroup("-pc_sor_local_symmetric","use SSOR separately on each processor","PCSORSetSymmetric",&flg);
102:   if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP);}
103:   PetscOptionsBoolGroup("-pc_sor_local_backward","use backward sweep locally","PCSORSetSymmetric",&flg);
104:   if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_BACKWARD_SWEEP);}
105:   PetscOptionsBoolGroupEnd("-pc_sor_local_forward","use forward sweep locally","PCSORSetSymmetric",&flg);
106:   if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_FORWARD_SWEEP);}
107:   PetscOptionsTail();
108:   return(0);
109: }

113: PetscErrorCode PCView_SOR(PC pc,PetscViewer viewer)
114: {
115:   PC_SOR         *jac = (PC_SOR*)pc->data;
116:   MatSORType     sym  = jac->sym;
117:   const char     *sortype;
119:   PetscBool      iascii;

122:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
123:   if (iascii) {
124:     if (sym & SOR_ZERO_INITIAL_GUESS) {PetscViewerASCIIPrintf(viewer,"  SOR:  zero initial guess\n");}
125:     if (sym == SOR_APPLY_UPPER)                                              sortype = "apply_upper";
126:     else if (sym == SOR_APPLY_LOWER)                                         sortype = "apply_lower";
127:     else if (sym & SOR_EISENSTAT)                                            sortype = "Eisenstat";
128:     else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP)             sortype = "symmetric";
129:     else if (sym & SOR_BACKWARD_SWEEP)                                       sortype = "backward";
130:     else if (sym & SOR_FORWARD_SWEEP)                                        sortype = "forward";
131:     else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
132:     else if (sym & SOR_LOCAL_FORWARD_SWEEP)                                  sortype = "local_forward";
133:     else if (sym & SOR_LOCAL_BACKWARD_SWEEP)                                 sortype = "local_backward";
134:     else                                                                     sortype = "unknown";
135:     PetscViewerASCIIPrintf(viewer,"  SOR: type = %s, iterations = %D, local iterations = %D, omega = %g\n",sortype,jac->its,jac->lits,(double)jac->omega);
136:   }
137:   return(0);
138: }


141: /* ------------------------------------------------------------------------------*/
144: static PetscErrorCode  PCSORSetSymmetric_SOR(PC pc,MatSORType flag)
145: {
146:   PC_SOR *jac = (PC_SOR*)pc->data;

149:   jac->sym = flag;
150:   return(0);
151: }

155: static PetscErrorCode  PCSORSetOmega_SOR(PC pc,PetscReal omega)
156: {
157:   PC_SOR *jac = (PC_SOR*)pc->data;

160:   if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
161:   jac->omega = omega;
162:   return(0);
163: }

167: static PetscErrorCode  PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
168: {
169:   PC_SOR *jac = (PC_SOR*)pc->data;

172:   jac->its  = its;
173:   jac->lits = lits;
174:   return(0);
175: }

179: static PetscErrorCode  PCSORGetSymmetric_SOR(PC pc,MatSORType *flag)
180: {
181:   PC_SOR *jac = (PC_SOR*)pc->data;

184:   *flag = jac->sym;
185:   return(0);
186: }

190: static PetscErrorCode  PCSORGetOmega_SOR(PC pc,PetscReal *omega)
191: {
192:   PC_SOR *jac = (PC_SOR*)pc->data;

195:   *omega = jac->omega;
196:   return(0);
197: }

201: static PetscErrorCode  PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
202: {
203:   PC_SOR *jac = (PC_SOR*)pc->data;

206:   if (its)  *its = jac->its;
207:   if (lits) *lits = jac->lits;
208:   return(0);
209: }

211: /* ------------------------------------------------------------------------------*/
214: /*@
215:    PCSORGetSymmetric - Gets the form the SOR preconditioner is using;   backward, or forward relaxation.  The local variants perform SOR on
216:    each processor.  By default forward relaxation is used.

218:    Logically Collective on PC

220:    Input Parameter:
221: .  pc - the preconditioner context

223:    Output Parameter:
224: .  flag - one of the following
225: .vb
226:     SOR_FORWARD_SWEEP
227:     SOR_BACKWARD_SWEEP
228:     SOR_SYMMETRIC_SWEEP
229:     SOR_LOCAL_FORWARD_SWEEP
230:     SOR_LOCAL_BACKWARD_SWEEP
231:     SOR_LOCAL_SYMMETRIC_SWEEP
232: .ve

234:    Options Database Keys:
235: +  -pc_sor_symmetric - Activates symmetric version
236: .  -pc_sor_backward - Activates backward version
237: .  -pc_sor_local_forward - Activates local forward version
238: .  -pc_sor_local_symmetric - Activates local symmetric version
239: -  -pc_sor_local_backward - Activates local backward version

241:    Notes:
242:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
243:    which can be chosen with the option
244: .  -pc_type eisenstat - Activates Eisenstat trick

246:    Level: intermediate

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

250: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
251: @*/
252: PetscErrorCode  PCSORGetSymmetric(PC pc,MatSORType *flag)
253: {

258:   PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
259:   return(0);
260: }

264: /*@
265:    PCSORGetOmega - Gets the SOR relaxation coefficient, omega
266:    (where omega = 1.0 by default).

268:    Logically Collective on PC

270:    Input Parameter:
271: .  pc - the preconditioner context

273:    Output Parameter:
274: .  omega - relaxation coefficient (0 < omega < 2).

276:    Options Database Key:
277: .  -pc_sor_omega <omega> - Sets omega

279:    Level: intermediate

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

283: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
284: @*/
285: PetscErrorCode  PCSORGetOmega(PC pc,PetscReal *omega)
286: {

291:   PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
292:   return(0);
293: }

297: /*@
298:    PCSORGetIterations - Gets the number of inner iterations to
299:    be used by the SOR preconditioner. The default is 1.

301:    Logically Collective on PC

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

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

310:    Options Database Key:
311: +  -pc_sor_its <its> - Sets number of iterations
312: -  -pc_sor_lits <lits> - Sets number of local iterations

314:    Level: intermediate

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

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

320: .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
321: @*/
322: PetscErrorCode  PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
323: {

328:   PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
329:   return(0);
330: }

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

339:    Logically Collective on PC

341:    Input Parameters:
342: +  pc - the preconditioner context
343: -  flag - one of the following
344: .vb
345:     SOR_FORWARD_SWEEP
346:     SOR_BACKWARD_SWEEP
347:     SOR_SYMMETRIC_SWEEP
348:     SOR_LOCAL_FORWARD_SWEEP
349:     SOR_LOCAL_BACKWARD_SWEEP
350:     SOR_LOCAL_SYMMETRIC_SWEEP
351: .ve

353:    Options Database Keys:
354: +  -pc_sor_symmetric - Activates symmetric version
355: .  -pc_sor_backward - Activates backward version
356: .  -pc_sor_local_forward - Activates local forward version
357: .  -pc_sor_local_symmetric - Activates local symmetric version
358: -  -pc_sor_local_backward - Activates local backward version

360:    Notes:
361:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
362:    which can be chosen with the option
363: .  -pc_type eisenstat - Activates Eisenstat trick

365:    Level: intermediate

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

369: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
370: @*/
371: PetscErrorCode  PCSORSetSymmetric(PC pc,MatSORType flag)
372: {

378:   PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
379:   return(0);
380: }

384: /*@
385:    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
386:    (where omega = 1.0 by default).

388:    Logically Collective on PC

390:    Input Parameters:
391: +  pc - the preconditioner context
392: -  omega - relaxation coefficient (0 < omega < 2).

394:    Options Database Key:
395: .  -pc_sor_omega <omega> - Sets omega

397:    Level: intermediate

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

401: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega()
402: @*/
403: PetscErrorCode  PCSORSetOmega(PC pc,PetscReal omega)
404: {

410:   PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
411:   return(0);
412: }

416: /*@
417:    PCSORSetIterations - Sets the number of inner iterations to
418:    be used by the SOR preconditioner. The default is 1.

420:    Logically Collective on PC

422:    Input Parameters:
423: +  pc - the preconditioner context
424: .  lits - number of local iterations, smoothings over just variables on processor
425: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

427:    Options Database Key:
428: +  -pc_sor_its <its> - Sets number of iterations
429: -  -pc_sor_lits <lits> - Sets number of local iterations

431:    Level: intermediate

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

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

437: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
438: @*/
439: PetscErrorCode  PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
440: {

446:   PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
447:   return(0);
448: }

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

453:    Options Database Keys:
454: +  -pc_sor_symmetric - Activates symmetric version
455: .  -pc_sor_backward - Activates backward version
456: .  -pc_sor_forward - Activates forward version
457: .  -pc_sor_local_forward - Activates local forward version
458: .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
459: .  -pc_sor_local_backward - Activates local backward version
460: .  -pc_sor_omega <omega> - Sets omega
461: .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
462: .  -pc_sor_its <its> - Sets number of iterations   (default 1)
463: -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)

465:    Level: beginner

467:   Concepts: SOR, preconditioners, Gauss-Seidel

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

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

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

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

483: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
484:            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
485: M*/

489: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
490: {
492:   PC_SOR         *jac;

495:   PetscNewLog(pc,&jac);

497:   pc->ops->apply           = PCApply_SOR;
498:   pc->ops->applytranspose  = PCApplyTranspose_SOR;
499:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
500:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
501:   pc->ops->setup           = 0;
502:   pc->ops->view            = PCView_SOR;
503:   pc->ops->destroy         = PCDestroy_SOR;
504:   pc->data                 = (void*)jac;
505:   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
506:   jac->omega               = 1.0;
507:   jac->fshift              = 0.0;
508:   jac->its                 = 1;
509:   jac->lits                = 1;

511:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
512:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
513:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
514:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
515:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
516:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
517:   return(0);
518: }