xref: /petsc/src/ksp/pc/impls/sor/sor.c (revision c3b5f7ba6bc5ce25a01a67bb37ba5d34b02bbbd7)
1 /*
2    Defines a  (S)SOR  preconditioner for any Mat implementation
3 */
4 #include <petsc/private/pcimpl.h>               /*I "petscpc.h" I*/
5 
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;
13 
14 static PetscErrorCode PCDestroy_SOR(PC pc)
15 {
16   PetscFunctionBegin;
17   PetscCall(PetscFree(pc->data));
18   PetscFunctionReturn(0);
19 }
20 
21 static PetscErrorCode PCApply_SOR(PC pc,Vec x,Vec y)
22 {
23   PC_SOR         *jac = (PC_SOR*)pc->data;
24   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
25 
26   PetscFunctionBegin;
27   PetscCall(MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,jac->fshift,jac->its,jac->lits,y));
28   PetscCall(MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason));
29   PetscFunctionReturn(0);
30 }
31 
32 static PetscErrorCode PCApplyTranspose_SOR(PC pc,Vec x,Vec y)
33 {
34   PC_SOR         *jac = (PC_SOR*)pc->data;
35   PetscInt       flag = jac->sym | SOR_ZERO_INITIAL_GUESS;
36   PetscBool      set,sym;
37 
38   PetscFunctionBegin;
39   PetscCall(MatIsSymmetricKnown(pc->pmat,&set,&sym));
40   PetscCheck(set && sym && (jac->sym == SOR_SYMMETRIC_SWEEP || jac->sym == SOR_LOCAL_SYMMETRIC_SWEEP),PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Can only apply transpose of SOR if matrix is symmetric and sweep is symmetric");
41   PetscCall(MatSOR(pc->pmat,x,jac->omega,(MatSORType)flag,jac->fshift,jac->its,jac->lits,y));
42   PetscCall(MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason));
43   PetscFunctionReturn(0);
44 }
45 
46 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)
47 {
48   PC_SOR         *jac = (PC_SOR*)pc->data;
49   MatSORType     stype = jac->sym;
50 
51   PetscFunctionBegin;
52   PetscCall(PetscInfo(pc,"Warning, convergence critera ignored, using %D iterations\n",its));
53   if (guesszero) stype = (MatSORType) (stype | SOR_ZERO_INITIAL_GUESS);
54   PetscCall(MatSOR(pc->pmat,b,jac->omega,stype,jac->fshift,its*jac->its,jac->lits,y));
55   PetscCall(MatFactorGetError(pc->pmat,(MatFactorError*)&pc->failedreason));
56   *outits = its;
57   *reason = PCRICHARDSON_CONVERGED_ITS;
58   PetscFunctionReturn(0);
59 }
60 
61 PetscErrorCode PCSetFromOptions_SOR(PetscOptionItems *PetscOptionsObject,PC pc)
62 {
63   PC_SOR         *jac = (PC_SOR*)pc->data;
64   PetscBool      flg;
65 
66   PetscFunctionBegin;
67   PetscCall(PetscOptionsHead(PetscOptionsObject,"(S)SOR options"));
68   PetscCall(PetscOptionsReal("-pc_sor_omega","relaxation factor (0 < omega < 2)","PCSORSetOmega",jac->omega,&jac->omega,NULL));
69   PetscCall(PetscOptionsReal("-pc_sor_diagonal_shift","Add to the diagonal entries","",jac->fshift,&jac->fshift,NULL));
70   PetscCall(PetscOptionsInt("-pc_sor_its","number of inner SOR iterations","PCSORSetIterations",jac->its,&jac->its,NULL));
71   PetscCall(PetscOptionsInt("-pc_sor_lits","number of local inner SOR iterations","PCSORSetIterations",jac->lits,&jac->lits,NULL));
72   PetscCall(PetscOptionsBoolGroupBegin("-pc_sor_symmetric","SSOR, not SOR","PCSORSetSymmetric",&flg));
73   if (flg) PetscCall(PCSORSetSymmetric(pc,SOR_SYMMETRIC_SWEEP));
74   PetscCall(PetscOptionsBoolGroup("-pc_sor_backward","use backward sweep instead of forward","PCSORSetSymmetric",&flg));
75   if (flg) PetscCall(PCSORSetSymmetric(pc,SOR_BACKWARD_SWEEP));
76   PetscCall(PetscOptionsBoolGroup("-pc_sor_forward","use forward sweep","PCSORSetSymmetric",&flg));
77   if (flg) PetscCall(PCSORSetSymmetric(pc,SOR_FORWARD_SWEEP));
78   PetscCall(PetscOptionsBoolGroup("-pc_sor_local_symmetric","use SSOR separately on each processor","PCSORSetSymmetric",&flg));
79   if (flg) PetscCall(PCSORSetSymmetric(pc,SOR_LOCAL_SYMMETRIC_SWEEP));
80   PetscCall(PetscOptionsBoolGroup("-pc_sor_local_backward","use backward sweep locally","PCSORSetSymmetric",&flg));
81   if (flg) PetscCall(PCSORSetSymmetric(pc,SOR_LOCAL_BACKWARD_SWEEP));
82   PetscCall(PetscOptionsBoolGroupEnd("-pc_sor_local_forward","use forward sweep locally","PCSORSetSymmetric",&flg));
83   if (flg) PetscCall(PCSORSetSymmetric(pc,SOR_LOCAL_FORWARD_SWEEP));
84   PetscCall(PetscOptionsTail());
85   PetscFunctionReturn(0);
86 }
87 
88 PetscErrorCode PCView_SOR(PC pc,PetscViewer viewer)
89 {
90   PC_SOR         *jac = (PC_SOR*)pc->data;
91   MatSORType     sym  = jac->sym;
92   const char     *sortype;
93   PetscBool      iascii;
94 
95   PetscFunctionBegin;
96   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii));
97   if (iascii) {
98     if (sym & SOR_ZERO_INITIAL_GUESS) PetscCall(PetscViewerASCIIPrintf(viewer,"  zero initial guess\n"));
99     if (sym == SOR_APPLY_UPPER)                                              sortype = "apply_upper";
100     else if (sym == SOR_APPLY_LOWER)                                         sortype = "apply_lower";
101     else if (sym & SOR_EISENSTAT)                                            sortype = "Eisenstat";
102     else if ((sym & SOR_SYMMETRIC_SWEEP) == SOR_SYMMETRIC_SWEEP)             sortype = "symmetric";
103     else if (sym & SOR_BACKWARD_SWEEP)                                       sortype = "backward";
104     else if (sym & SOR_FORWARD_SWEEP)                                        sortype = "forward";
105     else if ((sym & SOR_LOCAL_SYMMETRIC_SWEEP) == SOR_LOCAL_SYMMETRIC_SWEEP) sortype = "local_symmetric";
106     else if (sym & SOR_LOCAL_FORWARD_SWEEP)                                  sortype = "local_forward";
107     else if (sym & SOR_LOCAL_BACKWARD_SWEEP)                                 sortype = "local_backward";
108     else                                                                     sortype = "unknown";
109     PetscCall(PetscViewerASCIIPrintf(viewer,"  type = %s, iterations = %D, local iterations = %D, omega = %g\n",sortype,jac->its,jac->lits,(double)jac->omega));
110   }
111   PetscFunctionReturn(0);
112 }
113 
114 /* ------------------------------------------------------------------------------*/
115 static PetscErrorCode  PCSORSetSymmetric_SOR(PC pc,MatSORType flag)
116 {
117   PC_SOR *jac = (PC_SOR*)pc->data;
118 
119   PetscFunctionBegin;
120   jac->sym = flag;
121   PetscFunctionReturn(0);
122 }
123 
124 static PetscErrorCode  PCSORSetOmega_SOR(PC pc,PetscReal omega)
125 {
126   PC_SOR *jac = (PC_SOR*)pc->data;
127 
128   PetscFunctionBegin;
129   PetscCheckFalse(omega >= 2.0 || omega <= 0.0,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relaxation out of range");
130   jac->omega = omega;
131   PetscFunctionReturn(0);
132 }
133 
134 static PetscErrorCode  PCSORSetIterations_SOR(PC pc,PetscInt its,PetscInt lits)
135 {
136   PC_SOR *jac = (PC_SOR*)pc->data;
137 
138   PetscFunctionBegin;
139   jac->its  = its;
140   jac->lits = lits;
141   PetscFunctionReturn(0);
142 }
143 
144 static PetscErrorCode  PCSORGetSymmetric_SOR(PC pc,MatSORType *flag)
145 {
146   PC_SOR *jac = (PC_SOR*)pc->data;
147 
148   PetscFunctionBegin;
149   *flag = jac->sym;
150   PetscFunctionReturn(0);
151 }
152 
153 static PetscErrorCode  PCSORGetOmega_SOR(PC pc,PetscReal *omega)
154 {
155   PC_SOR *jac = (PC_SOR*)pc->data;
156 
157   PetscFunctionBegin;
158   *omega = jac->omega;
159   PetscFunctionReturn(0);
160 }
161 
162 static PetscErrorCode  PCSORGetIterations_SOR(PC pc,PetscInt *its,PetscInt *lits)
163 {
164   PC_SOR *jac = (PC_SOR*)pc->data;
165 
166   PetscFunctionBegin;
167   if (its)  *its = jac->its;
168   if (lits) *lits = jac->lits;
169   PetscFunctionReturn(0);
170 }
171 
172 /* ------------------------------------------------------------------------------*/
173 /*@
174    PCSORGetSymmetric - Gets the form the SOR preconditioner is using;   backward, or forward relaxation.  The local variants perform SOR on
175    each processor.  By default forward relaxation is used.
176 
177    Logically Collective on PC
178 
179    Input Parameter:
180 .  pc - the preconditioner context
181 
182    Output Parameter:
183 .  flag - one of the following
184 .vb
185     SOR_FORWARD_SWEEP
186     SOR_BACKWARD_SWEEP
187     SOR_SYMMETRIC_SWEEP
188     SOR_LOCAL_FORWARD_SWEEP
189     SOR_LOCAL_BACKWARD_SWEEP
190     SOR_LOCAL_SYMMETRIC_SWEEP
191 .ve
192 
193    Options Database Keys:
194 +  -pc_sor_symmetric - Activates symmetric version
195 .  -pc_sor_backward - Activates backward version
196 .  -pc_sor_local_forward - Activates local forward version
197 .  -pc_sor_local_symmetric - Activates local symmetric version
198 -  -pc_sor_local_backward - Activates local backward version
199 
200    Notes:
201    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
202    which can be chosen with the option
203 .  -pc_type eisenstat - Activates Eisenstat trick
204 
205    Level: intermediate
206 
207 .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega(), PCSORSetSymmetric()
208 @*/
209 PetscErrorCode  PCSORGetSymmetric(PC pc,MatSORType *flag)
210 {
211   PetscFunctionBegin;
212   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
213   PetscUseMethod(pc,"PCSORGetSymmetric_C",(PC,MatSORType*),(pc,flag));
214   PetscFunctionReturn(0);
215 }
216 
217 /*@
218    PCSORGetOmega - Gets the SOR relaxation coefficient, omega
219    (where omega = 1.0 by default).
220 
221    Logically Collective on PC
222 
223    Input Parameter:
224 .  pc - the preconditioner context
225 
226    Output Parameter:
227 .  omega - relaxation coefficient (0 < omega < 2).
228 
229    Options Database Key:
230 .  -pc_sor_omega <omega> - Sets omega
231 
232    Level: intermediate
233 
234 .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), PCSORSetOmega()
235 @*/
236 PetscErrorCode  PCSORGetOmega(PC pc,PetscReal *omega)
237 {
238   PetscFunctionBegin;
239   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
240   PetscUseMethod(pc,"PCSORGetOmega_C",(PC,PetscReal*),(pc,omega));
241   PetscFunctionReturn(0);
242 }
243 
244 /*@
245    PCSORGetIterations - Gets the number of inner iterations to
246    be used by the SOR preconditioner. The default is 1.
247 
248    Logically Collective on PC
249 
250    Input Parameter:
251 .  pc - the preconditioner context
252 
253    Output Parameters:
254 +  lits - number of local iterations, smoothings over just variables on processor
255 -  its - number of parallel iterations to use; each parallel iteration has lits local iterations
256 
257    Options Database Key:
258 +  -pc_sor_its <its> - Sets number of iterations
259 -  -pc_sor_lits <lits> - Sets number of local iterations
260 
261    Level: intermediate
262 
263    Notes:
264     When run on one processor the number of smoothings is lits*its
265 
266 .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
267 @*/
268 PetscErrorCode  PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
269 {
270   PetscFunctionBegin;
271   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
272   PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));
273   PetscFunctionReturn(0);
274 }
275 
276 /*@
277    PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
278    backward, or forward relaxation.  The local variants perform SOR on
279    each processor.  By default forward relaxation is used.
280 
281    Logically Collective on PC
282 
283    Input Parameters:
284 +  pc - the preconditioner context
285 -  flag - one of the following
286 .vb
287     SOR_FORWARD_SWEEP
288     SOR_BACKWARD_SWEEP
289     SOR_SYMMETRIC_SWEEP
290     SOR_LOCAL_FORWARD_SWEEP
291     SOR_LOCAL_BACKWARD_SWEEP
292     SOR_LOCAL_SYMMETRIC_SWEEP
293 .ve
294 
295    Options Database Keys:
296 +  -pc_sor_symmetric - Activates symmetric version
297 .  -pc_sor_backward - Activates backward version
298 .  -pc_sor_local_forward - Activates local forward version
299 .  -pc_sor_local_symmetric - Activates local symmetric version
300 -  -pc_sor_local_backward - Activates local backward version
301 
302    Notes:
303    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
304    which can be chosen with the option
305 .  -pc_type eisenstat - Activates Eisenstat trick
306 
307    Level: intermediate
308 
309 .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
310 @*/
311 PetscErrorCode  PCSORSetSymmetric(PC pc,MatSORType flag)
312 {
313   PetscFunctionBegin;
314   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
315   PetscValidLogicalCollectiveEnum(pc,flag,2);
316   PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));
317   PetscFunctionReturn(0);
318 }
319 
320 /*@
321    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
322    (where omega = 1.0 by default).
323 
324    Logically Collective on PC
325 
326    Input Parameters:
327 +  pc - the preconditioner context
328 -  omega - relaxation coefficient (0 < omega < 2).
329 
330    Options Database Key:
331 .  -pc_sor_omega <omega> - Sets omega
332 
333    Level: intermediate
334 
335    Note:
336    If omega != 1, you will need to set the MAT_USE_INODES option to PETSC_FALSE on the matrix.
337 
338 .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega(), MatSetOption()
339 @*/
340 PetscErrorCode  PCSORSetOmega(PC pc,PetscReal omega)
341 {
342   PetscFunctionBegin;
343   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
344   PetscValidLogicalCollectiveReal(pc,omega,2);
345   PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));
346   PetscFunctionReturn(0);
347 }
348 
349 /*@
350    PCSORSetIterations - Sets the number of inner iterations to
351    be used by the SOR preconditioner. The default is 1.
352 
353    Logically Collective on PC
354 
355    Input Parameters:
356 +  pc - the preconditioner context
357 .  lits - number of local iterations, smoothings over just variables on processor
358 -  its - number of parallel iterations to use; each parallel iteration has lits local iterations
359 
360    Options Database Key:
361 +  -pc_sor_its <its> - Sets number of iterations
362 -  -pc_sor_lits <lits> - Sets number of local iterations
363 
364    Level: intermediate
365 
366    Notes:
367     When run on one processor the number of smoothings is lits*its
368 
369 .seealso: PCSORSetOmega(), PCSORSetSymmetric()
370 @*/
371 PetscErrorCode  PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
372 {
373   PetscFunctionBegin;
374   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
375   PetscValidLogicalCollectiveInt(pc,its,2);
376   PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));
377   PetscFunctionReturn(0);
378 }
379 
380 /*MC
381      PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning
382 
383    Options Database Keys:
384 +  -pc_sor_symmetric - Activates symmetric version
385 .  -pc_sor_backward - Activates backward version
386 .  -pc_sor_forward - Activates forward version
387 .  -pc_sor_local_forward - Activates local forward version
388 .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
389 .  -pc_sor_local_backward - Activates local backward version
390 .  -pc_sor_omega <omega> - Sets omega
391 .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
392 .  -pc_sor_its <its> - Sets number of iterations   (default 1)
393 -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)
394 
395    Level: beginner
396 
397    Notes:
398     Only implemented for the AIJ  and SeqBAIJ matrix formats.
399           Not a true parallel SOR, in parallel this implementation corresponds to block
400           Jacobi with SOR on each block.
401 
402           For AIJ matrix if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
403           zero will be used and hence the KSPSolve() will terminate with KSP_DIVERGED_NANORIF. If the option
404           KSPSetErrorIfNotConverged() or -ksp_error_if_not_converged the code will terminate as soon as it detects the
405           zero pivot.
406 
407           For SeqBAIJ matrices this implements point-block SOR, but the omega, its, lits options are not supported.
408 
409           For SeqBAIJ the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected
410           the computation is stopped with an error
411 
412           If used with KSPRICHARDSON and no monitors the convergence test is skipped to improve speed, thus it always iterates
413           the maximum number of iterations you've selected for KSP. It is usually used in this mode as a smoother for multigrid.
414 
415           If omega != 1, you will need to set the MAT_USE_INODES option to PETSC_FALSE on the matrix.
416 
417 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
418            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT, MatSetOption()
419 M*/
420 
421 PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
422 {
423   PC_SOR         *jac;
424 
425   PetscFunctionBegin;
426   PetscCall(PetscNewLog(pc,&jac));
427 
428   pc->ops->apply           = PCApply_SOR;
429   pc->ops->applytranspose  = PCApplyTranspose_SOR;
430   pc->ops->applyrichardson = PCApplyRichardson_SOR;
431   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
432   pc->ops->setup           = NULL;
433   pc->ops->view            = PCView_SOR;
434   pc->ops->destroy         = PCDestroy_SOR;
435   pc->data                 = (void*)jac;
436   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
437   jac->omega               = 1.0;
438   jac->fshift              = 0.0;
439   jac->its                 = 1;
440   jac->lits                = 1;
441 
442   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR));
443   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR));
444   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR));
445   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR));
446   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR));
447   PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR));
448   PetscFunctionReturn(0);
449 }
450