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