Actual source code: sor.c
petsc-3.6.4 2016-04-12
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: }
28: static PetscErrorCode PCApply_SOR(PC pc,Vec x,Vec y)
29: {
30: PC_SOR *jac = (PC_SOR*)pc->data;
32: PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
33: PetscReal fshift;
34:
36: fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
37: MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,fshift,jac->its,jac->lits,y);
38: return(0);
39: }
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;
47: MatSORType stype = jac->sym;
48: PetscReal fshift;
51: PetscInfo1(pc,"Warning, convergence critera ignored, using %D iterations\n",its);
52: if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
53: fshift = (jac->fshift ? jac->fshift : pc->erroriffailure ? 0.0 : -1.0);
54: MatSOR(pc->pmat,b,jac->omega,stype,fshift,its*jac->its,jac->lits,y);
55: *outits = its;
56: *reason = PCRICHARDSON_CONVERGED_ITS;
57: return(0);
58: }
62: PetscErrorCode PCSetFromOptions_SOR(PetscOptions *PetscOptionsObject,PC pc)
63: {
64: PC_SOR *jac = (PC_SOR*)pc->data;
66: PetscBool flg;
69: PetscOptionsHead(PetscOptionsObject,"(S)SOR options");
70: PetscOptionsReal("-pc_sor_omega","relaxation factor (0 < omega < 2)","PCSORSetOmega",jac->omega,&jac->omega,NULL);
71: PetscOptionsReal("-pc_sor_diagonal_shift","Add to the diagonal entries","",jac->fshift,&jac->fshift,NULL);
72: PetscOptionsInt("-pc_sor_its","number of inner SOR iterations","PCSORSetIterations",jac->its,&jac->its,NULL);
73: PetscOptionsInt("-pc_sor_lits","number of local inner SOR iterations","PCSORSetIterations",jac->lits,&jac->lits,NULL);
74: PetscOptionsBoolGroupBegin("-pc_sor_symmetric","SSOR, not SOR","PCSORSetSymmetric",&flg);
75: if (flg) {PCSORSetSymmetric(pc,SOR_SYMMETRIC_SWEEP);}
76: PetscOptionsBoolGroup("-pc_sor_backward","use backward sweep instead of forward","PCSORSetSymmetric",&flg);
77: if (flg) {PCSORSetSymmetric(pc,SOR_BACKWARD_SWEEP);}
78: PetscOptionsBoolGroup("-pc_sor_forward","use forward sweep","PCSORSetSymmetric",&flg);
79: if (flg) {PCSORSetSymmetric(pc,SOR_FORWARD_SWEEP);}
80: PetscOptionsBoolGroup("-pc_sor_local_symmetric","use SSOR separately on each processor","PCSORSetSymmetric",&flg);
81: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP);}
82: PetscOptionsBoolGroup("-pc_sor_local_backward","use backward sweep locally","PCSORSetSymmetric",&flg);
83: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_BACKWARD_SWEEP);}
84: PetscOptionsBoolGroupEnd("-pc_sor_local_forward","use forward sweep locally","PCSORSetSymmetric",&flg);
85: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_FORWARD_SWEEP);}
86: PetscOptionsTail();
87: return(0);
88: }
92: PetscErrorCode PCView_SOR(PC pc,PetscViewer viewer)
93: {
94: PC_SOR *jac = (PC_SOR*)pc->data;
95: MatSORType sym = jac->sym;
96: const char *sortype;
98: PetscBool iascii;
101: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
102: if (iascii) {
103: if (sym & SOR_ZERO_INITIAL_GUESS) {PetscViewerASCIIPrintf(viewer," SOR: zero initial guess\n");}
104: if (sym == SOR_APPLY_UPPER) sortype = "apply_upper";
105: else if (sym == SOR_APPLY_LOWER) sortype = "apply_lower";
106: else if (sym & SOR_EISENSTAT) sortype = "Eisenstat";
107: else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP) sortype = "symmetric";
108: else if (sym & SOR_BACKWARD_SWEEP) sortype = "backward";
109: else if (sym & SOR_FORWARD_SWEEP) sortype = "forward";
110: else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
111: else if (sym & SOR_LOCAL_FORWARD_SWEEP) sortype = "local_forward";
112: else if (sym & SOR_LOCAL_BACKWARD_SWEEP) sortype = "local_backward";
113: else sortype = "unknown";
114: PetscViewerASCIIPrintf(viewer," SOR: type = %s, iterations = %D, local iterations = %D, omega = %g\n",sortype,jac->its,jac->lits,(double)jac->omega);
115: }
116: return(0);
117: }
120: /* ------------------------------------------------------------------------------*/
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: }
134: static PetscErrorCode PCSORSetOmega_SOR(PC pc,PetscReal omega)
135: {
136: PC_SOR *jac = (PC_SOR*)pc->data;
139: if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
140: jac->omega = omega;
141: return(0);
142: }
146: static PetscErrorCode PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
147: {
148: PC_SOR *jac = (PC_SOR*)pc->data;
151: jac->its = its;
152: jac->lits = lits;
153: return(0);
154: }
158: static PetscErrorCode PCSORGetSymmetric_SOR(PC pc,MatSORType *flag)
159: {
160: PC_SOR *jac = (PC_SOR*)pc->data;
163: *flag = jac->sym;
164: return(0);
165: }
169: static PetscErrorCode PCSORGetOmega_SOR(PC pc,PetscReal *omega)
170: {
171: PC_SOR *jac = (PC_SOR*)pc->data;
174: *omega = jac->omega;
175: return(0);
176: }
180: static PetscErrorCode PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
181: {
182: PC_SOR *jac = (PC_SOR*)pc->data;
185: if (its) *its = jac->its;
186: if (lits) *lits = jac->lits;
187: return(0);
188: }
190: /* ------------------------------------------------------------------------------*/
193: /*@
194: PCSORGetSymmetric - Gets the form the SOR preconditioner is using; backward, or forward relaxation. The local variants perform SOR on
195: each processor. By default forward relaxation is used.
197: Logically Collective on PC
199: Input Parameter:
200: . pc - the preconditioner context
202: Output Parameter:
203: . flag - one of the following
204: .vb
205: SOR_FORWARD_SWEEP
206: SOR_BACKWARD_SWEEP
207: SOR_SYMMETRIC_SWEEP
208: SOR_LOCAL_FORWARD_SWEEP
209: SOR_LOCAL_BACKWARD_SWEEP
210: SOR_LOCAL_SYMMETRIC_SWEEP
211: .ve
213: Options Database Keys:
214: + -pc_sor_symmetric - Activates symmetric version
215: . -pc_sor_backward - Activates backward version
216: . -pc_sor_local_forward - Activates local forward version
217: . -pc_sor_local_symmetric - Activates local symmetric version
218: - -pc_sor_local_backward - Activates local backward version
220: Notes:
221: To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
222: which can be chosen with the option
223: . -pc_type eisenstat - Activates Eisenstat trick
225: Level: intermediate
227: .keywords: PC, SOR, SSOR, set, relaxation, sweep, forward, backward, symmetric
229: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
230: @*/
231: PetscErrorCode PCSORGetSymmetric(PC pc,MatSORType *flag)
232: {
237: PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
238: return(0);
239: }
243: /*@
244: PCSORGetOmega - Gets the SOR relaxation coefficient, omega
245: (where omega = 1.0 by default).
247: Logically Collective on PC
249: Input Parameter:
250: . pc - the preconditioner context
252: Output Parameter:
253: . omega - relaxation coefficient (0 < omega < 2).
255: Options Database Key:
256: . -pc_sor_omega <omega> - Sets omega
258: Level: intermediate
260: .keywords: PC, SOR, SSOR, set, relaxation, omega
262: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
263: @*/
264: PetscErrorCode PCSORGetOmega(PC pc,PetscReal *omega)
265: {
270: PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
271: return(0);
272: }
276: /*@
277: PCSORGetIterations - Gets the number of inner iterations to
278: be used by the SOR preconditioner. The default is 1.
280: Logically Collective on PC
282: Input Parameter:
283: . pc - the preconditioner context
285: Output Parameter:
286: + lits - number of local iterations, smoothings over just variables on processor
287: - its - number of parallel iterations to use; each parallel iteration has lits local iterations
289: Options Database Key:
290: + -pc_sor_its <its> - Sets number of iterations
291: - -pc_sor_lits <lits> - Sets number of local iterations
293: Level: intermediate
295: Notes: When run on one processor the number of smoothings is lits*its
297: .keywords: PC, SOR, SSOR, set, iterations
299: .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
300: @*/
301: PetscErrorCode PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
302: {
307: PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
308: return(0);
309: }
313: /*@
314: PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
315: backward, or forward relaxation. The local variants perform SOR on
316: each processor. By default forward relaxation is used.
318: Logically Collective on PC
320: Input Parameters:
321: + pc - the preconditioner context
322: - flag - one of the following
323: .vb
324: SOR_FORWARD_SWEEP
325: SOR_BACKWARD_SWEEP
326: SOR_SYMMETRIC_SWEEP
327: SOR_LOCAL_FORWARD_SWEEP
328: SOR_LOCAL_BACKWARD_SWEEP
329: SOR_LOCAL_SYMMETRIC_SWEEP
330: .ve
332: Options Database Keys:
333: + -pc_sor_symmetric - Activates symmetric version
334: . -pc_sor_backward - Activates backward version
335: . -pc_sor_local_forward - Activates local forward version
336: . -pc_sor_local_symmetric - Activates local symmetric version
337: - -pc_sor_local_backward - Activates local backward version
339: Notes:
340: To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
341: which can be chosen with the option
342: . -pc_type eisenstat - Activates Eisenstat trick
344: Level: intermediate
346: .keywords: PC, SOR, SSOR, set, relaxation, sweep, forward, backward, symmetric
348: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
349: @*/
350: PetscErrorCode PCSORSetSymmetric(PC pc,MatSORType flag)
351: {
357: PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
358: return(0);
359: }
363: /*@
364: PCSORSetOmega - Sets the SOR relaxation coefficient, omega
365: (where omega = 1.0 by default).
367: Logically Collective on PC
369: Input Parameters:
370: + pc - the preconditioner context
371: - omega - relaxation coefficient (0 < omega < 2).
373: Options Database Key:
374: . -pc_sor_omega <omega> - Sets omega
376: Level: intermediate
378: .keywords: PC, SOR, SSOR, set, relaxation, omega
380: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega()
381: @*/
382: PetscErrorCode PCSORSetOmega(PC pc,PetscReal omega)
383: {
389: PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
390: return(0);
391: }
395: /*@
396: PCSORSetIterations - Sets the number of inner iterations to
397: be used by the SOR preconditioner. The default is 1.
399: Logically Collective on PC
401: Input Parameters:
402: + pc - the preconditioner context
403: . lits - number of local iterations, smoothings over just variables on processor
404: - its - number of parallel iterations to use; each parallel iteration has lits local iterations
406: Options Database Key:
407: + -pc_sor_its <its> - Sets number of iterations
408: - -pc_sor_lits <lits> - Sets number of local iterations
410: Level: intermediate
412: Notes: When run on one processor the number of smoothings is lits*its
414: .keywords: PC, SOR, SSOR, set, iterations
416: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
417: @*/
418: PetscErrorCode PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
419: {
425: PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
426: return(0);
427: }
429: /*MC
430: PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning
432: Options Database Keys:
433: + -pc_sor_symmetric - Activates symmetric version
434: . -pc_sor_backward - Activates backward version
435: . -pc_sor_forward - Activates forward version
436: . -pc_sor_local_forward - Activates local forward version
437: . -pc_sor_local_symmetric - Activates local symmetric version (default version)
438: . -pc_sor_local_backward - Activates local backward version
439: . -pc_sor_omega <omega> - Sets omega
440: . -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
441: . -pc_sor_its <its> - Sets number of iterations (default 1)
442: - -pc_sor_lits <lits> - Sets number of local iterations (default 1)
444: Level: beginner
446: Concepts: SOR, preconditioners, Gauss-Seidel
448: Notes: Only implemented for the AIJ and SeqBAIJ matrix formats.
449: Not a true parallel SOR, in parallel this implementation corresponds to block
450: Jacobi with SOR on each block.
452: For AIJ matrix if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
453: zero will be used and hence the KSPSolve() will terminate with KSP_DIVERGED_NANORIF. If the option
454: KSPSetErrorIfNotConverged() or -ksp_error_if_not_converged the code will terminate as soon as it detects the
455: zero pivot.
457: For SeqBAIJ matrices this implements point-block SOR, but the omega, its, lits options are not supported.
459: For SeqBAIJ the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected
460: the computation is stopped with an error
462: Developer Notes: We should add support for diagonal blocks that are singular to generate a Inf and thus cause KSPSolve()
463: to terminate with KSP_DIVERGED_NANORIF instead of stopping the program allowing
464: a nonlinear solver/ODE integrator to recover without stopping the program as currently happens.
466: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
467: PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
468: M*/
472: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
473: {
475: PC_SOR *jac;
478: PetscNewLog(pc,&jac);
480: pc->ops->apply = PCApply_SOR;
481: pc->ops->applyrichardson = PCApplyRichardson_SOR;
482: pc->ops->setfromoptions = PCSetFromOptions_SOR;
483: pc->ops->setup = 0;
484: pc->ops->view = PCView_SOR;
485: pc->ops->destroy = PCDestroy_SOR;
486: pc->data = (void*)jac;
487: jac->sym = SOR_LOCAL_SYMMETRIC_SWEEP;
488: jac->omega = 1.0;
489: jac->fshift = 0.0;
490: jac->its = 1;
491: jac->lits = 1;
493: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
494: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
495: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
496: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
497: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
498: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
499: return(0);
500: }