Actual source code: sor.c

petsc-master 2017-04-25
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;

 15: static PetscErrorCode PCDestroy_SOR(PC pc)
 16: {

 20:   PetscFree(pc->data);
 21:   return(0);
 22: }

 24: static PetscErrorCode PCApply_SOR(PC pc,Vec x,Vec y)
 25: {
 26:   PC_SOR         *jac = (PC_SOR*)pc->data;
 28:   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
 29:   PetscReal      fshift;

 32:   fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
 33:   MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,fshift,jac->its,jac->lits,y);
 34:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 35:   return(0);
 36: }

 38: static PetscErrorCode PCApplyTranspose_SOR(PC pc,Vec x,Vec y)
 39: {
 40:   PC_SOR         *jac = (PC_SOR*)pc->data;
 42:   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
 43:   PetscReal      fshift;
 44:   PetscBool      set,sym;

 47:   MatIsSymmetricKnown(pc->pmat,&set,&sym);
 48:   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");
 49:   fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
 50:   MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,fshift,jac->its,jac->lits,y);
 51:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 52:   return(0);
 53: }

 55: 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)
 56: {
 57:   PC_SOR         *jac = (PC_SOR*)pc->data;
 59:   MatSORType     stype = jac->sym;
 60:   PetscReal      fshift;

 63:   PetscInfo1(pc,"Warning, convergence critera ignored, using %D iterations\n",its);
 64:   if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
 65:   fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
 66:   MatSOR(pc->pmat,b,jac->omega,stype,fshift,its*jac->its,jac->lits,y);
 67:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 68:   *outits = its;
 69:   *reason = PCRICHARDSON_CONVERGED_ITS;
 70:   return(0);
 71: }

 73: PetscErrorCode PCSetFromOptions_SOR(PetscOptionItems *PetscOptionsObject,PC pc)
 74: {
 75:   PC_SOR         *jac = (PC_SOR*)pc->data;
 77:   PetscBool      flg;

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

101: PetscErrorCode PCView_SOR(PC pc,PetscViewer viewer)
102: {
103:   PC_SOR         *jac = (PC_SOR*)pc->data;
104:   MatSORType     sym  = jac->sym;
105:   const char     *sortype;
107:   PetscBool      iascii;

110:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
111:   if (iascii) {
112:     if (sym & SOR_ZERO_INITIAL_GUESS) {PetscViewerASCIIPrintf(viewer,"  SOR:  zero initial guess\n");}
113:     if (sym == SOR_APPLY_UPPER)                                              sortype = "apply_upper";
114:     else if (sym == SOR_APPLY_LOWER)                                         sortype = "apply_lower";
115:     else if (sym & SOR_EISENSTAT)                                            sortype = "Eisenstat";
116:     else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP)             sortype = "symmetric";
117:     else if (sym & SOR_BACKWARD_SWEEP)                                       sortype = "backward";
118:     else if (sym & SOR_FORWARD_SWEEP)                                        sortype = "forward";
119:     else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
120:     else if (sym & SOR_LOCAL_FORWARD_SWEEP)                                  sortype = "local_forward";
121:     else if (sym & SOR_LOCAL_BACKWARD_SWEEP)                                 sortype = "local_backward";
122:     else                                                                     sortype = "unknown";
123:     PetscViewerASCIIPrintf(viewer,"  SOR: type = %s, iterations = %D, local iterations = %D, omega = %g\n",sortype,jac->its,jac->lits,(double)jac->omega);
124:   }
125:   return(0);
126: }


129: /* ------------------------------------------------------------------------------*/
130: static PetscErrorCode  PCSORSetSymmetric_SOR(PC pc,MatSORType flag)
131: {
132:   PC_SOR *jac = (PC_SOR*)pc->data;

135:   jac->sym = flag;
136:   return(0);
137: }

139: static PetscErrorCode  PCSORSetOmega_SOR(PC pc,PetscReal omega)
140: {
141:   PC_SOR *jac = (PC_SOR*)pc->data;

144:   if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
145:   jac->omega = omega;
146:   return(0);
147: }

149: static PetscErrorCode  PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
150: {
151:   PC_SOR *jac = (PC_SOR*)pc->data;

154:   jac->its  = its;
155:   jac->lits = lits;
156:   return(0);
157: }

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

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

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

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

177: static PetscErrorCode  PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
178: {
179:   PC_SOR *jac = (PC_SOR*)pc->data;

182:   if (its)  *its = jac->its;
183:   if (lits) *lits = jac->lits;
184:   return(0);
185: }

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

192:    Logically Collective on PC

194:    Input Parameter:
195: .  pc - the preconditioner context

197:    Output Parameter:
198: .  flag - one of the following
199: .vb
200:     SOR_FORWARD_SWEEP
201:     SOR_BACKWARD_SWEEP
202:     SOR_SYMMETRIC_SWEEP
203:     SOR_LOCAL_FORWARD_SWEEP
204:     SOR_LOCAL_BACKWARD_SWEEP
205:     SOR_LOCAL_SYMMETRIC_SWEEP
206: .ve

208:    Options Database Keys:
209: +  -pc_sor_symmetric - Activates symmetric version
210: .  -pc_sor_backward - Activates backward version
211: .  -pc_sor_local_forward - Activates local forward version
212: .  -pc_sor_local_symmetric - Activates local symmetric version
213: -  -pc_sor_local_backward - Activates local backward version

215:    Notes:
216:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
217:    which can be chosen with the option
218: .  -pc_type eisenstat - Activates Eisenstat trick

220:    Level: intermediate

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

224: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
225: @*/
226: PetscErrorCode  PCSORGetSymmetric(PC pc,MatSORType *flag)
227: {

232:   PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
233:   return(0);
234: }

236: /*@
237:    PCSORGetOmega - Gets the SOR relaxation coefficient, omega
238:    (where omega = 1.0 by default).

240:    Logically Collective on PC

242:    Input Parameter:
243: .  pc - the preconditioner context

245:    Output Parameter:
246: .  omega - relaxation coefficient (0 < omega < 2).

248:    Options Database Key:
249: .  -pc_sor_omega <omega> - Sets omega

251:    Level: intermediate

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

255: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
256: @*/
257: PetscErrorCode  PCSORGetOmega(PC pc,PetscReal *omega)
258: {

263:   PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
264:   return(0);
265: }

267: /*@
268:    PCSORGetIterations - Gets the number of inner iterations to
269:    be used by the SOR preconditioner. The default is 1.

271:    Logically Collective on PC

273:    Input Parameter:
274: .  pc - the preconditioner context

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

280:    Options Database Key:
281: +  -pc_sor_its <its> - Sets number of iterations
282: -  -pc_sor_lits <lits> - Sets number of local iterations

284:    Level: intermediate

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

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

290: .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
291: @*/
292: PetscErrorCode  PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
293: {

298:   PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
299:   return(0);
300: }

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

307:    Logically Collective on PC

309:    Input Parameters:
310: +  pc - the preconditioner context
311: -  flag - one of the following
312: .vb
313:     SOR_FORWARD_SWEEP
314:     SOR_BACKWARD_SWEEP
315:     SOR_SYMMETRIC_SWEEP
316:     SOR_LOCAL_FORWARD_SWEEP
317:     SOR_LOCAL_BACKWARD_SWEEP
318:     SOR_LOCAL_SYMMETRIC_SWEEP
319: .ve

321:    Options Database Keys:
322: +  -pc_sor_symmetric - Activates symmetric version
323: .  -pc_sor_backward - Activates backward version
324: .  -pc_sor_local_forward - Activates local forward version
325: .  -pc_sor_local_symmetric - Activates local symmetric version
326: -  -pc_sor_local_backward - Activates local backward version

328:    Notes:
329:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
330:    which can be chosen with the option
331: .  -pc_type eisenstat - Activates Eisenstat trick

333:    Level: intermediate

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

337: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
338: @*/
339: PetscErrorCode  PCSORSetSymmetric(PC pc,MatSORType flag)
340: {

346:   PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
347:   return(0);
348: }

350: /*@
351:    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
352:    (where omega = 1.0 by default).

354:    Logically Collective on PC

356:    Input Parameters:
357: +  pc - the preconditioner context
358: -  omega - relaxation coefficient (0 < omega < 2).

360:    Options Database Key:
361: .  -pc_sor_omega <omega> - Sets omega

363:    Level: intermediate

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

367: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega()
368: @*/
369: PetscErrorCode  PCSORSetOmega(PC pc,PetscReal omega)
370: {

376:   PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
377:   return(0);
378: }

380: /*@
381:    PCSORSetIterations - Sets the number of inner iterations to
382:    be used by the SOR preconditioner. The default is 1.

384:    Logically Collective on PC

386:    Input Parameters:
387: +  pc - the preconditioner context
388: .  lits - number of local iterations, smoothings over just variables on processor
389: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

391:    Options Database Key:
392: +  -pc_sor_its <its> - Sets number of iterations
393: -  -pc_sor_lits <lits> - Sets number of local iterations

395:    Level: intermediate

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

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

401: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
402: @*/
403: PetscErrorCode  PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
404: {

410:   PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
411:   return(0);
412: }

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

417:    Options Database Keys:
418: +  -pc_sor_symmetric - Activates symmetric version
419: .  -pc_sor_backward - Activates backward version
420: .  -pc_sor_forward - Activates forward version
421: .  -pc_sor_local_forward - Activates local forward version
422: .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
423: .  -pc_sor_local_backward - Activates local backward version
424: .  -pc_sor_omega <omega> - Sets omega
425: .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
426: .  -pc_sor_its <its> - Sets number of iterations   (default 1)
427: -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)

429:    Level: beginner

431:   Concepts: SOR, preconditioners, Gauss-Seidel

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

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

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

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

447: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
448:            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
449: M*/

451: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
452: {
454:   PC_SOR         *jac;

457:   PetscNewLog(pc,&jac);

459:   pc->ops->apply           = PCApply_SOR;
460:   pc->ops->applytranspose  = PCApplyTranspose_SOR;
461:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
462:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
463:   pc->ops->setup           = 0;
464:   pc->ops->view            = PCView_SOR;
465:   pc->ops->destroy         = PCDestroy_SOR;
466:   pc->data                 = (void*)jac;
467:   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
468:   jac->omega               = 1.0;
469:   jac->fshift              = 0.0;
470:   jac->its                 = 1;
471:   jac->lits                = 1;

473:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
474:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
475:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
476:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
477:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
478:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
479:   return(0);
480: }