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