Actual source code: sor.c
petsc-3.8.4 2018-03-24
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," 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," 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: If used with KSPRICHARDSON and no monitors the convergence test is skipped to improve speed, thus it always iterates
448: the maximum number of iterations you've selected for KSP. It is usually used in this mode as a smoother for multigrid.
450: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
451: PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
452: M*/
454: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
455: {
457: PC_SOR *jac;
460: PetscNewLog(pc,&jac);
462: pc->ops->apply = PCApply_SOR;
463: pc->ops->applytranspose = PCApplyTranspose_SOR;
464: pc->ops->applyrichardson = PCApplyRichardson_SOR;
465: pc->ops->setfromoptions = PCSetFromOptions_SOR;
466: pc->ops->setup = 0;
467: pc->ops->view = PCView_SOR;
468: pc->ops->destroy = PCDestroy_SOR;
469: pc->data = (void*)jac;
470: jac->sym = SOR_LOCAL_SYMMETRIC_SWEEP;
471: jac->omega = 1.0;
472: jac->fshift = 0.0;
473: jac->its = 1;
474: jac->lits = 1;
476: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
477: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
478: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
479: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
480: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
481: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
482: return(0);
483: }