Actual source code: sor.c
petsc-3.7.7 2017-09-25
2: /*
3: Defines a (S)SOR preconditioner for any Mat implementation
4: */
5: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
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: }
26: #include <petsc/private/matimpl.h>
29: static PetscErrorCode PCApply_SOR(PC pc,Vec x,Vec y)
30: {
31: PC_SOR *jac = (PC_SOR*)pc->data;
33: PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
34: PetscReal fshift;
35:
37: fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
38: MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,fshift,jac->its,jac->lits,y);
39: if (pc->pmat->errortype) pc->failedreason = (PCFailedReason)pc->pmat->errortype;
40: return(0);
41: }
45: 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)
46: {
47: PC_SOR *jac = (PC_SOR*)pc->data;
49: MatSORType stype = jac->sym;
50: PetscReal fshift;
53: PetscInfo1(pc,"Warning, convergence critera ignored, using %D iterations\n",its);
54: if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
55: fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
56: MatSOR(pc->pmat,b,jac->omega,stype,fshift,its*jac->its,jac->lits,y);
57: if (pc->pmat->errortype) pc->failedreason = (PCFailedReason)pc->pmat->errortype;
58: *outits = its;
59: *reason = PCRICHARDSON_CONVERGED_ITS;
60: return(0);
61: }
65: PetscErrorCode PCSetFromOptions_SOR(PetscOptionItems *PetscOptionsObject,PC pc)
66: {
67: PC_SOR *jac = (PC_SOR*)pc->data;
69: PetscBool flg;
72: PetscOptionsHead(PetscOptionsObject,"(S)SOR options");
73: PetscOptionsReal("-pc_sor_omega","relaxation factor (0 < omega < 2)","PCSORSetOmega",jac->omega,&jac->omega,NULL);
74: PetscOptionsReal("-pc_sor_diagonal_shift","Add to the diagonal entries","",jac->fshift,&jac->fshift,NULL);
75: PetscOptionsInt("-pc_sor_its","number of inner SOR iterations","PCSORSetIterations",jac->its,&jac->its,NULL);
76: PetscOptionsInt("-pc_sor_lits","number of local inner SOR iterations","PCSORSetIterations",jac->lits,&jac->lits,NULL);
77: PetscOptionsBoolGroupBegin("-pc_sor_symmetric","SSOR, not SOR","PCSORSetSymmetric",&flg);
78: if (flg) {PCSORSetSymmetric(pc,SOR_SYMMETRIC_SWEEP);}
79: PetscOptionsBoolGroup("-pc_sor_backward","use backward sweep instead of forward","PCSORSetSymmetric",&flg);
80: if (flg) {PCSORSetSymmetric(pc,SOR_BACKWARD_SWEEP);}
81: PetscOptionsBoolGroup("-pc_sor_forward","use forward sweep","PCSORSetSymmetric",&flg);
82: if (flg) {PCSORSetSymmetric(pc,SOR_FORWARD_SWEEP);}
83: PetscOptionsBoolGroup("-pc_sor_local_symmetric","use SSOR separately on each processor","PCSORSetSymmetric",&flg);
84: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP);}
85: PetscOptionsBoolGroup("-pc_sor_local_backward","use backward sweep locally","PCSORSetSymmetric",&flg);
86: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_BACKWARD_SWEEP);}
87: PetscOptionsBoolGroupEnd("-pc_sor_local_forward","use forward sweep locally","PCSORSetSymmetric",&flg);
88: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_FORWARD_SWEEP);}
89: PetscOptionsTail();
90: return(0);
91: }
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," SOR: 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," SOR: type = %s, iterations = %D, local iterations = %D, omega = %g\n",sortype,jac->its,jac->lits,(double)jac->omega);
118: }
119: return(0);
120: }
123: /* ------------------------------------------------------------------------------*/
126: static PetscErrorCode PCSORSetSymmetric_SOR(PC pc,MatSORType flag)
127: {
128: PC_SOR *jac = (PC_SOR*)pc->data;
131: jac->sym = flag;
132: return(0);
133: }
137: static PetscErrorCode PCSORSetOmega_SOR(PC pc,PetscReal omega)
138: {
139: PC_SOR *jac = (PC_SOR*)pc->data;
142: if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
143: jac->omega = omega;
144: return(0);
145: }
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: }
161: static PetscErrorCode PCSORGetSymmetric_SOR(PC pc,MatSORType *flag)
162: {
163: PC_SOR *jac = (PC_SOR*)pc->data;
166: *flag = jac->sym;
167: return(0);
168: }
172: static PetscErrorCode PCSORGetOmega_SOR(PC pc,PetscReal *omega)
173: {
174: PC_SOR *jac = (PC_SOR*)pc->data;
177: *omega = jac->omega;
178: return(0);
179: }
183: static PetscErrorCode PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
184: {
185: PC_SOR *jac = (PC_SOR*)pc->data;
188: if (its) *its = jac->its;
189: if (lits) *lits = jac->lits;
190: return(0);
191: }
193: /* ------------------------------------------------------------------------------*/
196: /*@
197: PCSORGetSymmetric - Gets the form the SOR preconditioner is using; backward, or forward relaxation. The local variants perform SOR on
198: each processor. By default forward relaxation is used.
200: Logically Collective on PC
202: Input Parameter:
203: . pc - the preconditioner context
205: Output Parameter:
206: . flag - one of the following
207: .vb
208: SOR_FORWARD_SWEEP
209: SOR_BACKWARD_SWEEP
210: SOR_SYMMETRIC_SWEEP
211: SOR_LOCAL_FORWARD_SWEEP
212: SOR_LOCAL_BACKWARD_SWEEP
213: SOR_LOCAL_SYMMETRIC_SWEEP
214: .ve
216: Options Database Keys:
217: + -pc_sor_symmetric - Activates symmetric version
218: . -pc_sor_backward - Activates backward version
219: . -pc_sor_local_forward - Activates local forward version
220: . -pc_sor_local_symmetric - Activates local symmetric version
221: - -pc_sor_local_backward - Activates local backward version
223: Notes:
224: To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
225: which can be chosen with the option
226: . -pc_type eisenstat - Activates Eisenstat trick
228: Level: intermediate
230: .keywords: PC, SOR, SSOR, set, relaxation, sweep, forward, backward, symmetric
232: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
233: @*/
234: PetscErrorCode PCSORGetSymmetric(PC pc,MatSORType *flag)
235: {
240: PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
241: return(0);
242: }
246: /*@
247: PCSORGetOmega - Gets the SOR relaxation coefficient, omega
248: (where omega = 1.0 by default).
250: Logically Collective on PC
252: Input Parameter:
253: . pc - the preconditioner context
255: Output Parameter:
256: . omega - relaxation coefficient (0 < omega < 2).
258: Options Database Key:
259: . -pc_sor_omega <omega> - Sets omega
261: Level: intermediate
263: .keywords: PC, SOR, SSOR, set, relaxation, omega
265: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
266: @*/
267: PetscErrorCode PCSORGetOmega(PC pc,PetscReal *omega)
268: {
273: PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
274: return(0);
275: }
279: /*@
280: PCSORGetIterations - Gets the number of inner iterations to
281: be used by the SOR preconditioner. The default is 1.
283: Logically Collective on PC
285: Input Parameter:
286: . pc - the preconditioner context
288: Output Parameter:
289: + lits - number of local iterations, smoothings over just variables on processor
290: - its - number of parallel iterations to use; each parallel iteration has lits local iterations
292: Options Database Key:
293: + -pc_sor_its <its> - Sets number of iterations
294: - -pc_sor_lits <lits> - Sets number of local iterations
296: Level: intermediate
298: Notes: When run on one processor the number of smoothings is lits*its
300: .keywords: PC, SOR, SSOR, set, iterations
302: .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
303: @*/
304: PetscErrorCode PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
305: {
310: PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
311: return(0);
312: }
316: /*@
317: PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
318: backward, or forward relaxation. The local variants perform SOR on
319: each processor. By default forward relaxation is used.
321: Logically Collective on PC
323: Input Parameters:
324: + pc - the preconditioner context
325: - flag - one of the following
326: .vb
327: SOR_FORWARD_SWEEP
328: SOR_BACKWARD_SWEEP
329: SOR_SYMMETRIC_SWEEP
330: SOR_LOCAL_FORWARD_SWEEP
331: SOR_LOCAL_BACKWARD_SWEEP
332: SOR_LOCAL_SYMMETRIC_SWEEP
333: .ve
335: Options Database Keys:
336: + -pc_sor_symmetric - Activates symmetric version
337: . -pc_sor_backward - Activates backward version
338: . -pc_sor_local_forward - Activates local forward version
339: . -pc_sor_local_symmetric - Activates local symmetric version
340: - -pc_sor_local_backward - Activates local backward version
342: Notes:
343: To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
344: which can be chosen with the option
345: . -pc_type eisenstat - Activates Eisenstat trick
347: Level: intermediate
349: .keywords: PC, SOR, SSOR, set, relaxation, sweep, forward, backward, symmetric
351: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
352: @*/
353: PetscErrorCode PCSORSetSymmetric(PC pc,MatSORType flag)
354: {
360: PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
361: return(0);
362: }
366: /*@
367: PCSORSetOmega - Sets the SOR relaxation coefficient, omega
368: (where omega = 1.0 by default).
370: Logically Collective on PC
372: Input Parameters:
373: + pc - the preconditioner context
374: - omega - relaxation coefficient (0 < omega < 2).
376: Options Database Key:
377: . -pc_sor_omega <omega> - Sets omega
379: Level: intermediate
381: .keywords: PC, SOR, SSOR, set, relaxation, omega
383: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega()
384: @*/
385: PetscErrorCode PCSORSetOmega(PC pc,PetscReal omega)
386: {
392: PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
393: return(0);
394: }
398: /*@
399: PCSORSetIterations - Sets the number of inner iterations to
400: be used by the SOR preconditioner. The default is 1.
402: Logically Collective on PC
404: Input Parameters:
405: + pc - the preconditioner context
406: . lits - number of local iterations, smoothings over just variables on processor
407: - its - number of parallel iterations to use; each parallel iteration has lits local iterations
409: Options Database Key:
410: + -pc_sor_its <its> - Sets number of iterations
411: - -pc_sor_lits <lits> - Sets number of local iterations
413: Level: intermediate
415: Notes: When run on one processor the number of smoothings is lits*its
417: .keywords: PC, SOR, SSOR, set, iterations
419: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
420: @*/
421: PetscErrorCode PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
422: {
428: PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
429: return(0);
430: }
432: /*MC
433: PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning
435: Options Database Keys:
436: + -pc_sor_symmetric - Activates symmetric version
437: . -pc_sor_backward - Activates backward version
438: . -pc_sor_forward - Activates forward version
439: . -pc_sor_local_forward - Activates local forward version
440: . -pc_sor_local_symmetric - Activates local symmetric version (default version)
441: . -pc_sor_local_backward - Activates local backward version
442: . -pc_sor_omega <omega> - Sets omega
443: . -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
444: . -pc_sor_its <its> - Sets number of iterations (default 1)
445: - -pc_sor_lits <lits> - Sets number of local iterations (default 1)
447: Level: beginner
449: Concepts: SOR, preconditioners, Gauss-Seidel
451: Notes: Only implemented for the AIJ and SeqBAIJ matrix formats.
452: Not a true parallel SOR, in parallel this implementation corresponds to block
453: Jacobi with SOR on each block.
455: For AIJ matrix if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
456: zero will be used and hence the KSPSolve() will terminate with KSP_DIVERGED_NANORIF. If the option
457: KSPSetErrorIfNotConverged() or -ksp_error_if_not_converged the code will terminate as soon as it detects the
458: zero pivot.
460: For SeqBAIJ matrices this implements point-block SOR, but the omega, its, lits options are not supported.
462: For SeqBAIJ the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected
463: the computation is stopped with an error
465: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
466: PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
467: M*/
471: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
472: {
474: PC_SOR *jac;
477: PetscNewLog(pc,&jac);
479: pc->ops->apply = PCApply_SOR;
480: pc->ops->applyrichardson = PCApplyRichardson_SOR;
481: pc->ops->setfromoptions = PCSetFromOptions_SOR;
482: pc->ops->setup = 0;
483: pc->ops->view = PCView_SOR;
484: pc->ops->destroy = PCDestroy_SOR;
485: pc->data = (void*)jac;
486: jac->sym = SOR_LOCAL_SYMMETRIC_SWEEP;
487: jac->omega = 1.0;
488: jac->fshift = 0.0;
489: jac->its = 1;
490: jac->lits = 1;
492: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
493: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
494: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
495: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
496: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
497: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
498: return(0);
499: }