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: {
19: PetscFree(pc->data);
20: return(0);
21: }
23: static PetscErrorCode PCApply_SOR(PC pc,Vec x,Vec y)
24: {
25: PC_SOR *jac = (PC_SOR*)pc->data;
27: PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
30: MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,jac->fshift,jac->its,jac->lits,y);
31: MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
32: return(0);
33: }
35: static PetscErrorCode PCApplyTranspose_SOR(PC pc,Vec x,Vec y)
36: {
37: PC_SOR *jac = (PC_SOR*)pc->data;
39: PetscInt flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
40: PetscBool set,sym;
43: MatIsSymmetricKnown(pc->pmat,&set,&sym);
44: 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");
45: MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,jac->fshift,jac->its,jac->lits,y);
46: MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
47: return(0);
48: }
50: 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)
51: {
52: PC_SOR *jac = (PC_SOR*)pc->data;
54: MatSORType stype = jac->sym;
57: PetscInfo1(pc,"Warning, convergence critera ignored, using %D iterations\n",its);
58: if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
59: MatSOR(pc->pmat,b,jac->omega,stype,jac->fshift,its*jac->its,jac->lits,y);
60: MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason);
61: *outits = its;
62: *reason = PCRICHARDSON_CONVERGED_ITS;
63: return(0);
64: }
66: PetscErrorCode PCSetFromOptions_SOR(PetscOptionItems *PetscOptionsObject,PC pc)
67: {
68: PC_SOR *jac = (PC_SOR*)pc->data;
70: PetscBool flg;
73: PetscOptionsHead(PetscOptionsObject,"(S)SOR options");
74: PetscOptionsReal("-pc_sor_omega","relaxation factor (0 < omega < 2)","PCSORSetOmega",jac->omega,&jac->omega,NULL);
75: PetscOptionsReal("-pc_sor_diagonal_shift","Add to the diagonal entries","",jac->fshift,&jac->fshift,NULL);
76: PetscOptionsInt("-pc_sor_its","number of inner SOR iterations","PCSORSetIterations",jac->its,&jac->its,NULL);
77: PetscOptionsInt("-pc_sor_lits","number of local inner SOR iterations","PCSORSetIterations",jac->lits,&jac->lits,NULL);
78: PetscOptionsBoolGroupBegin("-pc_sor_symmetric","SSOR, not SOR","PCSORSetSymmetric",&flg);
79: if (flg) {PCSORSetSymmetric(pc,SOR_SYMMETRIC_SWEEP);}
80: PetscOptionsBoolGroup("-pc_sor_backward","use backward sweep instead of forward","PCSORSetSymmetric",&flg);
81: if (flg) {PCSORSetSymmetric(pc,SOR_BACKWARD_SWEEP);}
82: PetscOptionsBoolGroup("-pc_sor_forward","use forward sweep","PCSORSetSymmetric",&flg);
83: if (flg) {PCSORSetSymmetric(pc,SOR_FORWARD_SWEEP);}
84: PetscOptionsBoolGroup("-pc_sor_local_symmetric","use SSOR separately on each processor","PCSORSetSymmetric",&flg);
85: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP);}
86: PetscOptionsBoolGroup("-pc_sor_local_backward","use backward sweep locally","PCSORSetSymmetric",&flg);
87: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_BACKWARD_SWEEP);}
88: PetscOptionsBoolGroupEnd("-pc_sor_local_forward","use forward sweep locally","PCSORSetSymmetric",&flg);
89: if (flg) {PCSORSetSymmetric(pc,SOR_LOCAL_FORWARD_SWEEP);}
90: PetscOptionsTail();
91: return(0);
92: }
94: PetscErrorCode PCView_SOR(PC pc,PetscViewer viewer)
95: {
96: PC_SOR *jac = (PC_SOR*)pc->data;
97: MatSORType sym = jac->sym;
98: const char *sortype;
100: PetscBool iascii;
103: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
104: if (iascii) {
105: if (sym & SOR_ZERO_INITIAL_GUESS) {PetscViewerASCIIPrintf(viewer," zero initial guess\n");}
106: if (sym == SOR_APPLY_UPPER) sortype = "apply_upper";
107: else if (sym == SOR_APPLY_LOWER) sortype = "apply_lower";
108: else if (sym & SOR_EISENSTAT) sortype = "Eisenstat";
109: else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP) sortype = "symmetric";
110: else if (sym & SOR_BACKWARD_SWEEP) sortype = "backward";
111: else if (sym & SOR_FORWARD_SWEEP) sortype = "forward";
112: else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
113: else if (sym & SOR_LOCAL_FORWARD_SWEEP) sortype = "local_forward";
114: else if (sym & SOR_LOCAL_BACKWARD_SWEEP) sortype = "local_backward";
115: else sortype = "unknown";
116: PetscViewerASCIIPrintf(viewer," type = %s, iterations = %D, local iterations = %D, omega = %g\n",sortype,jac->its,jac->lits,(double)jac->omega);
117: }
118: return(0);
119: }
122: /* ------------------------------------------------------------------------------*/
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: }
132: static PetscErrorCode PCSORSetOmega_SOR(PC pc,PetscReal omega)
133: {
134: PC_SOR *jac = (PC_SOR*)pc->data;
137: if (omega >= 2.0 || omega <= 0.0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
138: jac->omega = omega;
139: return(0);
140: }
142: static PetscErrorCode PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
143: {
144: PC_SOR *jac = (PC_SOR*)pc->data;
147: jac->its = its;
148: jac->lits = lits;
149: return(0);
150: }
152: static PetscErrorCode PCSORGetSymmetric_SOR(PC pc,MatSORType *flag)
153: {
154: PC_SOR *jac = (PC_SOR*)pc->data;
157: *flag = jac->sym;
158: return(0);
159: }
161: static PetscErrorCode PCSORGetOmega_SOR(PC pc,PetscReal *omega)
162: {
163: PC_SOR *jac = (PC_SOR*)pc->data;
166: *omega = jac->omega;
167: return(0);
168: }
170: static PetscErrorCode PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
171: {
172: PC_SOR *jac = (PC_SOR*)pc->data;
175: if (its) *its = jac->its;
176: if (lits) *lits = jac->lits;
177: return(0);
178: }
180: /* ------------------------------------------------------------------------------*/
181: /*@
182: PCSORGetSymmetric - Gets the form the SOR preconditioner is using; backward, or forward relaxation. The local variants perform SOR on
183: each processor. By default forward relaxation is used.
185: Logically Collective on PC
187: Input Parameter:
188: . pc - the preconditioner context
190: Output Parameter:
191: . flag - one of the following
192: .vb
193: SOR_FORWARD_SWEEP
194: SOR_BACKWARD_SWEEP
195: SOR_SYMMETRIC_SWEEP
196: SOR_LOCAL_FORWARD_SWEEP
197: SOR_LOCAL_BACKWARD_SWEEP
198: SOR_LOCAL_SYMMETRIC_SWEEP
199: .ve
201: Options Database Keys:
202: + -pc_sor_symmetric - Activates symmetric version
203: . -pc_sor_backward - Activates backward version
204: . -pc_sor_local_forward - Activates local forward version
205: . -pc_sor_local_symmetric - Activates local symmetric version
206: - -pc_sor_local_backward - Activates local backward version
208: Notes:
209: To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
210: which can be chosen with the option
211: . -pc_type eisenstat - Activates Eisenstat trick
213: Level: intermediate
215: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
216: @*/
217: PetscErrorCode PCSORGetSymmetric(PC pc,MatSORType *flag)
218: {
223: PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
224: return(0);
225: }
227: /*@
228: PCSORGetOmega - Gets the SOR relaxation coefficient, omega
229: (where omega = 1.0 by default).
231: Logically Collective on PC
233: Input Parameter:
234: . pc - the preconditioner context
236: Output Parameter:
237: . omega - relaxation coefficient (0 < omega < 2).
239: Options Database Key:
240: . -pc_sor_omega <omega> - Sets omega
242: Level: intermediate
244: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
245: @*/
246: PetscErrorCode PCSORGetOmega(PC pc,PetscReal *omega)
247: {
252: PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
253: return(0);
254: }
256: /*@
257: PCSORGetIterations - Gets the number of inner iterations to
258: be used by the SOR preconditioner. The default is 1.
260: Logically Collective on PC
262: Input Parameter:
263: . pc - the preconditioner context
265: Output Parameter:
266: + lits - number of local iterations, smoothings over just variables on processor
267: - its - number of parallel iterations to use; each parallel iteration has lits local iterations
269: Options Database Key:
270: + -pc_sor_its <its> - Sets number of iterations
271: - -pc_sor_lits <lits> - Sets number of local iterations
273: Level: intermediate
275: Notes:
276: When run on one processor the number of smoothings is lits*its
278: .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
279: @*/
280: PetscErrorCode PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
281: {
286: PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
287: return(0);
288: }
290: /*@
291: PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
292: backward, or forward relaxation. The local variants perform SOR on
293: each processor. By default forward relaxation is used.
295: Logically Collective on PC
297: Input Parameters:
298: + pc - the preconditioner context
299: - flag - one of the following
300: .vb
301: SOR_FORWARD_SWEEP
302: SOR_BACKWARD_SWEEP
303: SOR_SYMMETRIC_SWEEP
304: SOR_LOCAL_FORWARD_SWEEP
305: SOR_LOCAL_BACKWARD_SWEEP
306: SOR_LOCAL_SYMMETRIC_SWEEP
307: .ve
309: Options Database Keys:
310: + -pc_sor_symmetric - Activates symmetric version
311: . -pc_sor_backward - Activates backward version
312: . -pc_sor_local_forward - Activates local forward version
313: . -pc_sor_local_symmetric - Activates local symmetric version
314: - -pc_sor_local_backward - Activates local backward version
316: Notes:
317: To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
318: which can be chosen with the option
319: . -pc_type eisenstat - Activates Eisenstat trick
321: Level: intermediate
323: .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
324: @*/
325: PetscErrorCode PCSORSetSymmetric(PC pc,MatSORType flag)
326: {
332: PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
333: return(0);
334: }
336: /*@
337: PCSORSetOmega - Sets the SOR relaxation coefficient, omega
338: (where omega = 1.0 by default).
340: Logically Collective on PC
342: Input Parameters:
343: + pc - the preconditioner context
344: - omega - relaxation coefficient (0 < omega < 2).
346: Options Database Key:
347: . -pc_sor_omega <omega> - Sets omega
349: Level: intermediate
351: Note:
352: If omega != 1, you will need to set the MAT_USE_INODES option to PETSC_FALSE on the matrix.
355: .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), MatSetOption()
356: @*/
357: PetscErrorCode PCSORSetOmega(PC pc,PetscReal omega)
358: {
364: PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
365: return(0);
366: }
368: /*@
369: PCSORSetIterations - Sets the number of inner iterations to
370: be used by the SOR preconditioner. The default is 1.
372: Logically Collective on PC
374: Input Parameters:
375: + pc - the preconditioner context
376: . lits - number of local iterations, smoothings over just variables on processor
377: - its - number of parallel iterations to use; each parallel iteration has lits local iterations
379: Options Database Key:
380: + -pc_sor_its <its> - Sets number of iterations
381: - -pc_sor_lits <lits> - Sets number of local iterations
383: Level: intermediate
385: Notes:
386: When run on one processor the number of smoothings is lits*its
388: .seealso: PCSORSetOmega(), PCSORSetSymmetric()
389: @*/
390: PetscErrorCode PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
391: {
397: PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
398: return(0);
399: }
401: /*MC
402: PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning
404: Options Database Keys:
405: + -pc_sor_symmetric - Activates symmetric version
406: . -pc_sor_backward - Activates backward version
407: . -pc_sor_forward - Activates forward version
408: . -pc_sor_local_forward - Activates local forward version
409: . -pc_sor_local_symmetric - Activates local symmetric version (default version)
410: . -pc_sor_local_backward - Activates local backward version
411: . -pc_sor_omega <omega> - Sets omega
412: . -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
413: . -pc_sor_its <its> - Sets number of iterations (default 1)
414: - -pc_sor_lits <lits> - Sets number of local iterations (default 1)
416: Level: beginner
418: Notes:
419: Only implemented for the AIJ and SeqBAIJ matrix formats.
420: Not a true parallel SOR, in parallel this implementation corresponds to block
421: Jacobi with SOR on each block.
423: For AIJ matrix if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
424: zero will be used and hence the KSPSolve() will terminate with KSP_DIVERGED_NANORIF. If the option
425: KSPSetErrorIfNotConverged() or -ksp_error_if_not_converged the code will terminate as soon as it detects the
426: zero pivot.
428: For SeqBAIJ matrices this implements point-block SOR, but the omega, its, lits options are not supported.
430: For SeqBAIJ the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected
431: the computation is stopped with an error
433: If used with KSPRICHARDSON and no monitors the convergence test is skipped to improve speed, thus it always iterates
434: the maximum number of iterations you've selected for KSP. It is usually used in this mode as a smoother for multigrid.
436: If omega != 1, you will need to set the MAT_USE_INODES option to PETSC_FALSE on the matrix.
438: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
439: PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT, MatSetOption()
440: M*/
442: PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
443: {
445: PC_SOR *jac;
448: PetscNewLog(pc,&jac);
450: pc->ops->apply = PCApply_SOR;
451: pc->ops->applytranspose = PCApplyTranspose_SOR;
452: pc->ops->applyrichardson = PCApplyRichardson_SOR;
453: pc->ops->setfromoptions = PCSetFromOptions_SOR;
454: pc->ops->setup = NULL;
455: pc->ops->view = PCView_SOR;
456: pc->ops->destroy = PCDestroy_SOR;
457: pc->data = (void*)jac;
458: jac->sym = SOR_LOCAL_SYMMETRIC_SWEEP;
459: jac->omega = 1.0;
460: jac->fshift = 0.0;
461: jac->its = 1;
462: jac->lits = 1;
464: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);
465: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);
466: PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);
467: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);
468: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);
469: PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);
470: return(0);
471: }