Actual source code: sor.c

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

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

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

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

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

 30:   MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,jac->fshift,jac->its,jac->lits,y);
 31:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 32:   return(0);
 33: }

 35: static PetscErrorCode PCApplyTranspose_SOR(PC pc,Vec x,Vec y)
 36: {
 37:   PC_SOR         *jac = (PC_SOR*)pc->data;
 39:   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
 40:   PetscBool      set,sym;

 43:   MatIsSymmetricKnown(pc->pmat,&set,&sym);
 44:   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");
 45:   MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,jac->fshift,jac->its,jac->lits,y);
 46:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 47:   return(0);
 48: }

 50: 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)
 51: {
 52:   PC_SOR         *jac = (PC_SOR*)pc->data;
 54:   MatSORType     stype = jac->sym;

 57:   PetscInfo1(pc,"Warning, convergence critera ignored, using %D iterations\n",its);
 58:   if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
 59:   MatSOR(pc->pmat,b,jac->omega,stype,jac->fshift,its*jac->its,jac->lits,y);
 60:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 61:   *outits = its;
 62:   *reason = PCRICHARDSON_CONVERGED_ITS;
 63:   return(0);
 64: }

 66: PetscErrorCode PCSetFromOptions_SOR(PetscOptionItems *PetscOptionsObject,PC pc)
 67: {
 68:   PC_SOR         *jac = (PC_SOR*)pc->data;
 70:   PetscBool      flg;

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

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

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


122: /* ------------------------------------------------------------------------------*/
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: }

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

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

142: static PetscErrorCode  PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
143: {
144:   PC_SOR *jac = (PC_SOR*)pc->data;

147:   jac->its  = its;
148:   jac->lits = lits;
149:   return(0);
150: }

152: static PetscErrorCode  PCSORGetSymmetric_SOR(PC pc,MatSORType *flag)
153: {
154:   PC_SOR *jac = (PC_SOR*)pc->data;

157:   *flag = jac->sym;
158:   return(0);
159: }

161: static PetscErrorCode  PCSORGetOmega_SOR(PC pc,PetscReal *omega)
162: {
163:   PC_SOR *jac = (PC_SOR*)pc->data;

166:   *omega = jac->omega;
167:   return(0);
168: }

170: static PetscErrorCode  PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
171: {
172:   PC_SOR *jac = (PC_SOR*)pc->data;

175:   if (its)  *its = jac->its;
176:   if (lits) *lits = jac->lits;
177:   return(0);
178: }

180: /* ------------------------------------------------------------------------------*/
181: /*@
182:    PCSORGetSymmetric - Gets the form the SOR preconditioner is using;   backward, or forward relaxation.  The local variants perform SOR on
183:    each processor.  By default forward relaxation is used.

185:    Logically Collective on PC

187:    Input Parameter:
188: .  pc - the preconditioner context

190:    Output Parameter:
191: .  flag - one of the following
192: .vb
193:     SOR_FORWARD_SWEEP
194:     SOR_BACKWARD_SWEEP
195:     SOR_SYMMETRIC_SWEEP
196:     SOR_LOCAL_FORWARD_SWEEP
197:     SOR_LOCAL_BACKWARD_SWEEP
198:     SOR_LOCAL_SYMMETRIC_SWEEP
199: .ve

201:    Options Database Keys:
202: +  -pc_sor_symmetric - Activates symmetric version
203: .  -pc_sor_backward - Activates backward version
204: .  -pc_sor_local_forward - Activates local forward version
205: .  -pc_sor_local_symmetric - Activates local symmetric version
206: -  -pc_sor_local_backward - Activates local backward version

208:    Notes:
209:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
210:    which can be chosen with the option
211: .  -pc_type eisenstat - Activates Eisenstat trick

213:    Level: intermediate

215: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
216: @*/
217: PetscErrorCode  PCSORGetSymmetric(PC pc,MatSORType *flag)
218: {

223:   PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
224:   return(0);
225: }

227: /*@
228:    PCSORGetOmega - Gets the SOR relaxation coefficient, omega
229:    (where omega = 1.0 by default).

231:    Logically Collective on PC

233:    Input Parameter:
234: .  pc - the preconditioner context

236:    Output Parameter:
237: .  omega - relaxation coefficient (0 < omega < 2).

239:    Options Database Key:
240: .  -pc_sor_omega <omega> - Sets omega

242:    Level: intermediate

244: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
245: @*/
246: PetscErrorCode  PCSORGetOmega(PC pc,PetscReal *omega)
247: {

252:   PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
253:   return(0);
254: }

256: /*@
257:    PCSORGetIterations - Gets the number of inner iterations to
258:    be used by the SOR preconditioner. The default is 1.

260:    Logically Collective on PC

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

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

269:    Options Database Key:
270: +  -pc_sor_its <its> - Sets number of iterations
271: -  -pc_sor_lits <lits> - Sets number of local iterations

273:    Level: intermediate

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

278: .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
279: @*/
280: PetscErrorCode  PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
281: {

286:   PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
287:   return(0);
288: }

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

295:    Logically Collective on PC

297:    Input Parameters:
298: +  pc - the preconditioner context
299: -  flag - one of the following
300: .vb
301:     SOR_FORWARD_SWEEP
302:     SOR_BACKWARD_SWEEP
303:     SOR_SYMMETRIC_SWEEP
304:     SOR_LOCAL_FORWARD_SWEEP
305:     SOR_LOCAL_BACKWARD_SWEEP
306:     SOR_LOCAL_SYMMETRIC_SWEEP
307: .ve

309:    Options Database Keys:
310: +  -pc_sor_symmetric - Activates symmetric version
311: .  -pc_sor_backward - Activates backward version
312: .  -pc_sor_local_forward - Activates local forward version
313: .  -pc_sor_local_symmetric - Activates local symmetric version
314: -  -pc_sor_local_backward - Activates local backward version

316:    Notes:
317:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
318:    which can be chosen with the option
319: .  -pc_type eisenstat - Activates Eisenstat trick

321:    Level: intermediate

323: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
324: @*/
325: PetscErrorCode  PCSORSetSymmetric(PC pc,MatSORType flag)
326: {

332:   PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
333:   return(0);
334: }

336: /*@
337:    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
338:    (where omega = 1.0 by default).

340:    Logically Collective on PC

342:    Input Parameters:
343: +  pc - the preconditioner context
344: -  omega - relaxation coefficient (0 < omega < 2).

346:    Options Database Key:
347: .  -pc_sor_omega <omega> - Sets omega

349:    Level: intermediate

351:    Note:
352:    If omega != 1, you will need to set the MAT_USE_INODES option to PETSC_FALSE on the matrix.


355: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), MatSetOption()
356: @*/
357: PetscErrorCode  PCSORSetOmega(PC pc,PetscReal omega)
358: {

364:   PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
365:   return(0);
366: }

368: /*@
369:    PCSORSetIterations - Sets the number of inner iterations to
370:    be used by the SOR preconditioner. The default is 1.

372:    Logically Collective on PC

374:    Input Parameters:
375: +  pc - the preconditioner context
376: .  lits - number of local iterations, smoothings over just variables on processor
377: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

379:    Options Database Key:
380: +  -pc_sor_its <its> - Sets number of iterations
381: -  -pc_sor_lits <lits> - Sets number of local iterations

383:    Level: intermediate

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

388: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
389: @*/
390: PetscErrorCode  PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
391: {

397:   PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
398:   return(0);
399: }

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

404:    Options Database Keys:
405: +  -pc_sor_symmetric - Activates symmetric version
406: .  -pc_sor_backward - Activates backward version
407: .  -pc_sor_forward - Activates forward version
408: .  -pc_sor_local_forward - Activates local forward version
409: .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
410: .  -pc_sor_local_backward - Activates local backward version
411: .  -pc_sor_omega <omega> - Sets omega
412: .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
413: .  -pc_sor_its <its> - Sets number of iterations   (default 1)
414: -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)

416:    Level: beginner

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

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

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

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

433:           If used with KSPRICHARDSON and no monitors the convergence test is skipped to improve speed, thus it always iterates
434:           the maximum number of iterations you've selected for KSP. It is usually used in this mode as a smoother for multigrid.

436:           If omega != 1, you will need to set the MAT_USE_INODES option to PETSC_FALSE on the matrix.

438: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
439:            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT, MatSetOption()
440: M*/

442: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
443: {
445:   PC_SOR         *jac;

448:   PetscNewLog(pc,&jac);

450:   pc->ops->apply           = PCApply_SOR;
451:   pc->ops->applytranspose  = PCApplyTranspose_SOR;
452:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
453:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
454:   pc->ops->setup           = NULL;
455:   pc->ops->view            = PCView_SOR;
456:   pc->ops->destroy         = PCDestroy_SOR;
457:   pc->data                 = (void*)jac;
458:   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
459:   jac->omega               = 1.0;
460:   jac->fshift              = 0.0;
461:   jac->its                 = 1;
462:   jac->lits                = 1;

464:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
465:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
466:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
467:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
468:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
469:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
470:   return(0);
471: }