Actual source code: sor.c

petsc-master 2019-08-18
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;

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

186:    Logically Collective on PC

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

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

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

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

214:    Level: intermediate

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

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

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

232:    Logically Collective on PC

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

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

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

243:    Level: intermediate

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

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

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

261:    Logically Collective on PC

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

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

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

274:    Level: intermediate

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

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

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

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

296:    Logically Collective on PC

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

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

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

322:    Level: intermediate

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

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

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

341:    Logically Collective on PC

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

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

350:    Level: intermediate

352: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega()
353: @*/
354: PetscErrorCode  PCSORSetOmega(PC pc,PetscReal omega)
355: {

361:   PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
362:   return(0);
363: }

365: /*@
366:    PCSORSetIterations - Sets the number of inner iterations to
367:    be used by the SOR preconditioner. The default is 1.

369:    Logically Collective on PC

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

376:    Options Database Key:
377: +  -pc_sor_its <its> - Sets number of iterations
378: -  -pc_sor_lits <lits> - Sets number of local iterations

380:    Level: intermediate

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

385: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
386: @*/
387: PetscErrorCode  PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
388: {

394:   PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
395:   return(0);
396: }

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

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

413:    Level: beginner

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

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

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

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

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

433: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
434:            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
435: M*/

437: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
438: {
440:   PC_SOR         *jac;

443:   PetscNewLog(pc,&jac);

445:   pc->ops->apply           = PCApply_SOR;
446:   pc->ops->applytranspose  = PCApplyTranspose_SOR;
447:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
448:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
449:   pc->ops->setup           = 0;
450:   pc->ops->view            = PCView_SOR;
451:   pc->ops->destroy         = PCDestroy_SOR;
452:   pc->data                 = (void*)jac;
453:   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
454:   jac->omega               = 1.0;
455:   jac->fshift              = 0.0;
456:   jac->its                 = 1;
457:   jac->lits                = 1;

459:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
460:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
461:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
462:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
463:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
464:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
465:   return(0);
466: }