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: {
 16:   PetscFree(pc->data);
 17:   return 0;
 18: }

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

 25:   MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,jac->fshift,jac->its,jac->lits,y);
 26:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 27:   return 0;
 28: }

 30: static PetscErrorCode PCApplyTranspose_SOR(PC pc,Vec x,Vec y)
 31: {
 32:   PC_SOR         *jac = (PC_SOR*)pc->data;
 33:   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
 34:   PetscBool      set,sym;

 36:   MatIsSymmetricKnown(pc->pmat,&set,&sym);
 38:   MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,jac->fshift,jac->its,jac->lits,y);
 39:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 40:   return 0;
 41: }

 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;
 46:   MatSORType     stype = jac->sym;

 48:   PetscInfo(pc,"Warning, convergence critera ignored, using %D iterations\n",its);
 49:   if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
 50:   MatSOR(pc->pmat,b,jac->omega,stype,jac->fshift,its*jac->its,jac->lits,y);
 51:   MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
 52:   *outits = its;
 53:   *reason = PCRICHARDSON_CONVERGED_ITS;
 54:   return 0;
 55: }

 57: PetscErrorCode PCSetFromOptions_SOR(PetscOptionItems *PetscOptionsObject,PC pc)
 58: {
 59:   PC_SOR         *jac = (PC_SOR*)pc->data;
 60:   PetscBool      flg;

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

 83: PetscErrorCode PCView_SOR(PC pc,PetscViewer viewer)
 84: {
 85:   PC_SOR         *jac = (PC_SOR*)pc->data;
 86:   MatSORType     sym  = jac->sym;
 87:   const char     *sortype;
 88:   PetscBool      iascii;

 90:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 91:   if (iascii) {
 92:     if (sym & SOR_ZERO_INITIAL_GUESS) PetscViewerASCIIPrintf(viewer,"  zero initial guess\n");
 93:     if (sym == SOR_APPLY_UPPER)                                              sortype = "apply_upper";
 94:     else if (sym == SOR_APPLY_LOWER)                                         sortype = "apply_lower";
 95:     else if (sym & SOR_EISENSTAT)                                            sortype = "Eisenstat";
 96:     else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP)             sortype = "symmetric";
 97:     else if (sym & SOR_BACKWARD_SWEEP)                                       sortype = "backward";
 98:     else if (sym & SOR_FORWARD_SWEEP)                                        sortype = "forward";
 99:     else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
100:     else if (sym & SOR_LOCAL_FORWARD_SWEEP)                                  sortype = "local_forward";
101:     else if (sym & SOR_LOCAL_BACKWARD_SWEEP)                                 sortype = "local_backward";
102:     else                                                                     sortype = "unknown";
103:     PetscViewerASCIIPrintf(viewer,"  type = %s, iterations = %D, local iterations = %D, omega = %g\n",sortype,jac->its,jac->lits,(double)jac->omega);
104:   }
105:   return 0;
106: }

108: /* ------------------------------------------------------------------------------*/
109: static PetscErrorCode  PCSORSetSymmetric_SOR(PC pc,MatSORType flag)
110: {
111:   PC_SOR *jac = (PC_SOR*)pc->data;

113:   jac->sym = flag;
114:   return 0;
115: }

117: static PetscErrorCode  PCSORSetOmega_SOR(PC pc,PetscReal omega)
118: {
119:   PC_SOR *jac = (PC_SOR*)pc->data;

122:   jac->omega = omega;
123:   return 0;
124: }

126: static PetscErrorCode  PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
127: {
128:   PC_SOR *jac = (PC_SOR*)pc->data;

130:   jac->its  = its;
131:   jac->lits = lits;
132:   return 0;
133: }

135: static PetscErrorCode  PCSORGetSymmetric_SOR(PC pc,MatSORType *flag)
136: {
137:   PC_SOR *jac = (PC_SOR*)pc->data;

139:   *flag = jac->sym;
140:   return 0;
141: }

143: static PetscErrorCode  PCSORGetOmega_SOR(PC pc,PetscReal *omega)
144: {
145:   PC_SOR *jac = (PC_SOR*)pc->data;

147:   *omega = jac->omega;
148:   return 0;
149: }

151: static PetscErrorCode  PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
152: {
153:   PC_SOR *jac = (PC_SOR*)pc->data;

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

160: /* ------------------------------------------------------------------------------*/
161: /*@
162:    PCSORGetSymmetric - Gets the form the SOR preconditioner is using;   backward, or forward relaxation.  The local variants perform SOR on
163:    each processor.  By default forward relaxation is used.

165:    Logically Collective on PC

167:    Input Parameter:
168: .  pc - the preconditioner context

170:    Output Parameter:
171: .  flag - one of the following
172: .vb
173:     SOR_FORWARD_SWEEP
174:     SOR_BACKWARD_SWEEP
175:     SOR_SYMMETRIC_SWEEP
176:     SOR_LOCAL_FORWARD_SWEEP
177:     SOR_LOCAL_BACKWARD_SWEEP
178:     SOR_LOCAL_SYMMETRIC_SWEEP
179: .ve

181:    Options Database Keys:
182: +  -pc_sor_symmetric - Activates symmetric version
183: .  -pc_sor_backward - Activates backward version
184: .  -pc_sor_local_forward - Activates local forward version
185: .  -pc_sor_local_symmetric - Activates local symmetric version
186: -  -pc_sor_local_backward - Activates local backward version

188:    Notes:
189:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
190:    which can be chosen with the option
191: .  -pc_type eisenstat - Activates Eisenstat trick

193:    Level: intermediate

195: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
196: @*/
197: PetscErrorCode  PCSORGetSymmetric(PC pc,MatSORType *flag)
198: {
200:   PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
201:   return 0;
202: }

204: /*@
205:    PCSORGetOmega - Gets the SOR relaxation coefficient, omega
206:    (where omega = 1.0 by default).

208:    Logically Collective on PC

210:    Input Parameter:
211: .  pc - the preconditioner context

213:    Output Parameter:
214: .  omega - relaxation coefficient (0 < omega < 2).

216:    Options Database Key:
217: .  -pc_sor_omega <omega> - Sets omega

219:    Level: intermediate

221: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
222: @*/
223: PetscErrorCode  PCSORGetOmega(PC pc,PetscReal *omega)
224: {
226:   PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
227:   return 0;
228: }

230: /*@
231:    PCSORGetIterations - Gets the number of inner iterations to
232:    be used by the SOR preconditioner. The default is 1.

234:    Logically Collective on PC

236:    Input Parameter:
237: .  pc - the preconditioner context

239:    Output Parameters:
240: +  lits - number of local iterations, smoothings over just variables on processor
241: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

243:    Options Database Key:
244: +  -pc_sor_its <its> - Sets number of iterations
245: -  -pc_sor_lits <lits> - Sets number of local iterations

247:    Level: intermediate

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

252: .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
253: @*/
254: PetscErrorCode  PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
255: {
257:   PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
258:   return 0;
259: }

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

266:    Logically Collective on PC

268:    Input Parameters:
269: +  pc - the preconditioner context
270: -  flag - one of the following
271: .vb
272:     SOR_FORWARD_SWEEP
273:     SOR_BACKWARD_SWEEP
274:     SOR_SYMMETRIC_SWEEP
275:     SOR_LOCAL_FORWARD_SWEEP
276:     SOR_LOCAL_BACKWARD_SWEEP
277:     SOR_LOCAL_SYMMETRIC_SWEEP
278: .ve

280:    Options Database Keys:
281: +  -pc_sor_symmetric - Activates symmetric version
282: .  -pc_sor_backward - Activates backward version
283: .  -pc_sor_local_forward - Activates local forward version
284: .  -pc_sor_local_symmetric - Activates local symmetric version
285: -  -pc_sor_local_backward - Activates local backward version

287:    Notes:
288:    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
289:    which can be chosen with the option
290: .  -pc_type eisenstat - Activates Eisenstat trick

292:    Level: intermediate

294: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
295: @*/
296: PetscErrorCode  PCSORSetSymmetric(PC pc,MatSORType flag)
297: {
300:   PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
301:   return 0;
302: }

304: /*@
305:    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
306:    (where omega = 1.0 by default).

308:    Logically Collective on PC

310:    Input Parameters:
311: +  pc - the preconditioner context
312: -  omega - relaxation coefficient (0 < omega < 2).

314:    Options Database Key:
315: .  -pc_sor_omega <omega> - Sets omega

317:    Level: intermediate

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

322: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), MatSetOption()
323: @*/
324: PetscErrorCode  PCSORSetOmega(PC pc,PetscReal omega)
325: {
328:   PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
329:   return 0;
330: }

332: /*@
333:    PCSORSetIterations - Sets the number of inner iterations to
334:    be used by the SOR preconditioner. The default is 1.

336:    Logically Collective on PC

338:    Input Parameters:
339: +  pc - the preconditioner context
340: .  lits - number of local iterations, smoothings over just variables on processor
341: -  its - number of parallel iterations to use; each parallel iteration has lits local iterations

343:    Options Database Key:
344: +  -pc_sor_its <its> - Sets number of iterations
345: -  -pc_sor_lits <lits> - Sets number of local iterations

347:    Level: intermediate

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

352: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
353: @*/
354: PetscErrorCode  PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
355: {
358:   PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
359:   return 0;
360: }

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

365:    Options Database Keys:
366: +  -pc_sor_symmetric - Activates symmetric version
367: .  -pc_sor_backward - Activates backward version
368: .  -pc_sor_forward - Activates forward version
369: .  -pc_sor_local_forward - Activates local forward version
370: .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
371: .  -pc_sor_local_backward - Activates local backward version
372: .  -pc_sor_omega <omega> - Sets omega
373: .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
374: .  -pc_sor_its <its> - Sets number of iterations   (default 1)
375: -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)

377:    Level: beginner

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

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

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

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

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

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

399: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
400:            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT, MatSetOption()
401: M*/

403: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
404: {
405:   PC_SOR         *jac;

407:   PetscNewLog(pc,&jac);

409:   pc->ops->apply           = PCApply_SOR;
410:   pc->ops->applytranspose  = PCApplyTranspose_SOR;
411:   pc->ops->applyrichardson = PCApplyRichardson_SOR;
412:   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
413:   pc->ops->setup           = NULL;
414:   pc->ops->view            = PCView_SOR;
415:   pc->ops->destroy         = PCDestroy_SOR;
416:   pc->data                 = (void*)jac;
417:   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
418:   jac->omega               = 1.0;
419:   jac->fshift              = 0.0;
420:   jac->its                 = 1;
421:   jac->lits                = 1;

423:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
424:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
425:   PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
426:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
427:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
428:   PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
429:   return 0;
430: }