xref: /petsc/src/ksp/pc/impls/sor/sor.c (revision 65f8aed5f7eaa1e2ef2ddeffe666264e0669c876)
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:
281     When run on one processor the number of smoothings is lits*its
282 
283 .keywords: PC, SOR, SSOR, set, iterations
284 
285 .seealso: PCSORSetOmega(), PCSORSetSymmetric(), PCSORSetIterations()
286 @*/
287 PetscErrorCode  PCSORGetIterations(PC pc,PetscInt *its,PetscInt *lits)
288 {
289   PetscErrorCode ierr;
290 
291   PetscFunctionBegin;
292   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
293   ierr = PetscUseMethod(pc,"PCSORGetIterations_C",(PC,PetscInt*,PetscInt*),(pc,its,lits));CHKERRQ(ierr);
294   PetscFunctionReturn(0);
295 }
296 
297 /*@
298    PCSORSetSymmetric - Sets the SOR preconditioner to use symmetric (SSOR),
299    backward, or forward relaxation.  The local variants perform SOR on
300    each processor.  By default forward relaxation is used.
301 
302    Logically Collective on PC
303 
304    Input Parameters:
305 +  pc - the preconditioner context
306 -  flag - one of the following
307 .vb
308     SOR_FORWARD_SWEEP
309     SOR_BACKWARD_SWEEP
310     SOR_SYMMETRIC_SWEEP
311     SOR_LOCAL_FORWARD_SWEEP
312     SOR_LOCAL_BACKWARD_SWEEP
313     SOR_LOCAL_SYMMETRIC_SWEEP
314 .ve
315 
316    Options Database Keys:
317 +  -pc_sor_symmetric - Activates symmetric version
318 .  -pc_sor_backward - Activates backward version
319 .  -pc_sor_local_forward - Activates local forward version
320 .  -pc_sor_local_symmetric - Activates local symmetric version
321 -  -pc_sor_local_backward - Activates local backward version
322 
323    Notes:
324    To use the Eisenstat trick with SSOR, employ the PCEISENSTAT preconditioner,
325    which can be chosen with the option
326 .  -pc_type eisenstat - Activates Eisenstat trick
327 
328    Level: intermediate
329 
330 .keywords: PC, SOR, SSOR, set, relaxation, sweep, forward, backward, symmetric
331 
332 .seealso: PCEisenstatSetOmega(), PCSORSetIterations(), PCSORSetOmega()
333 @*/
334 PetscErrorCode  PCSORSetSymmetric(PC pc,MatSORType flag)
335 {
336   PetscErrorCode ierr;
337 
338   PetscFunctionBegin;
339   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
340   PetscValidLogicalCollectiveEnum(pc,flag,2);
341   ierr = PetscTryMethod(pc,"PCSORSetSymmetric_C",(PC,MatSORType),(pc,flag));CHKERRQ(ierr);
342   PetscFunctionReturn(0);
343 }
344 
345 /*@
346    PCSORSetOmega - Sets the SOR relaxation coefficient, omega
347    (where omega = 1.0 by default).
348 
349    Logically Collective on PC
350 
351    Input Parameters:
352 +  pc - the preconditioner context
353 -  omega - relaxation coefficient (0 < omega < 2).
354 
355    Options Database Key:
356 .  -pc_sor_omega <omega> - Sets omega
357 
358    Level: intermediate
359 
360 .keywords: PC, SOR, SSOR, set, relaxation, omega
361 
362 .seealso: PCSORSetSymmetric(), PCSORSetIterations(), PCEisenstatSetOmega()
363 @*/
364 PetscErrorCode  PCSORSetOmega(PC pc,PetscReal omega)
365 {
366   PetscErrorCode ierr;
367 
368   PetscFunctionBegin;
369   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
370   PetscValidLogicalCollectiveReal(pc,omega,2);
371   ierr = PetscTryMethod(pc,"PCSORSetOmega_C",(PC,PetscReal),(pc,omega));CHKERRQ(ierr);
372   PetscFunctionReturn(0);
373 }
374 
375 /*@
376    PCSORSetIterations - Sets the number of inner iterations to
377    be used by the SOR preconditioner. The default is 1.
378 
379    Logically Collective on PC
380 
381    Input Parameters:
382 +  pc - the preconditioner context
383 .  lits - number of local iterations, smoothings over just variables on processor
384 -  its - number of parallel iterations to use; each parallel iteration has lits local iterations
385 
386    Options Database Key:
387 +  -pc_sor_its <its> - Sets number of iterations
388 -  -pc_sor_lits <lits> - Sets number of local iterations
389 
390    Level: intermediate
391 
392    Notes:
393     When run on one processor the number of smoothings is lits*its
394 
395 .keywords: PC, SOR, SSOR, set, iterations
396 
397 .seealso: PCSORSetOmega(), PCSORSetSymmetric()
398 @*/
399 PetscErrorCode  PCSORSetIterations(PC pc,PetscInt its,PetscInt lits)
400 {
401   PetscErrorCode ierr;
402 
403   PetscFunctionBegin;
404   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
405   PetscValidLogicalCollectiveInt(pc,its,2);
406   ierr = PetscTryMethod(pc,"PCSORSetIterations_C",(PC,PetscInt,PetscInt),(pc,its,lits));CHKERRQ(ierr);
407   PetscFunctionReturn(0);
408 }
409 
410 /*MC
411      PCSOR - (S)SOR (successive over relaxation, Gauss-Seidel) preconditioning
412 
413    Options Database Keys:
414 +  -pc_sor_symmetric - Activates symmetric version
415 .  -pc_sor_backward - Activates backward version
416 .  -pc_sor_forward - Activates forward version
417 .  -pc_sor_local_forward - Activates local forward version
418 .  -pc_sor_local_symmetric - Activates local symmetric version  (default version)
419 .  -pc_sor_local_backward - Activates local backward version
420 .  -pc_sor_omega <omega> - Sets omega
421 .  -pc_sor_diagonal_shift <shift> - shift the diagonal entries; useful if the matrix has zeros on the diagonal
422 .  -pc_sor_its <its> - Sets number of iterations   (default 1)
423 -  -pc_sor_lits <lits> - Sets number of local iterations  (default 1)
424 
425    Level: beginner
426 
427   Concepts: SOR, preconditioners, Gauss-Seidel
428 
429    Notes:
430     Only implemented for the AIJ  and SeqBAIJ matrix formats.
431           Not a true parallel SOR, in parallel this implementation corresponds to block
432           Jacobi with SOR on each block.
433 
434           For AIJ matrix if a diagonal entry is zero (and the diagonal shift is zero) then by default the inverse of that
435           zero will be used and hence the KSPSolve() will terminate with KSP_DIVERGED_NANORIF. If the option
436           KSPSetErrorIfNotConverged() or -ksp_error_if_not_converged the code will terminate as soon as it detects the
437           zero pivot.
438 
439           For SeqBAIJ matrices this implements point-block SOR, but the omega, its, lits options are not supported.
440 
441           For SeqBAIJ the diagonal blocks are inverted using dense LU with partial pivoting. If a zero pivot is detected
442           the computation is stopped with an error
443 
444           If used with KSPRICHARDSON and no monitors the convergence test is skipped to improve speed, thus it always iterates
445           the maximum number of iterations you've selected for KSP. It is usually used in this mode as a smoother for multigrid.
446 
447 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
448            PCSORSetIterations(), PCSORSetSymmetric(), PCSORSetOmega(), PCEISENSTAT
449 M*/
450 
451 PETSC_EXTERN PetscErrorCode PCCreate_SOR(PC pc)
452 {
453   PetscErrorCode ierr;
454   PC_SOR         *jac;
455 
456   PetscFunctionBegin;
457   ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr);
458 
459   pc->ops->apply           = PCApply_SOR;
460   pc->ops->applytranspose  = PCApplyTranspose_SOR;
461   pc->ops->applyrichardson = PCApplyRichardson_SOR;
462   pc->ops->setfromoptions  = PCSetFromOptions_SOR;
463   pc->ops->setup           = 0;
464   pc->ops->view            = PCView_SOR;
465   pc->ops->destroy         = PCDestroy_SOR;
466   pc->data                 = (void*)jac;
467   jac->sym                 = SOR_LOCAL_SYMMETRIC_SWEEP;
468   jac->omega               = 1.0;
469   jac->fshift              = 0.0;
470   jac->its                 = 1;
471   jac->lits                = 1;
472 
473   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSORSetSymmetric_C",PCSORSetSymmetric_SOR);CHKERRQ(ierr);
474   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSORSetOmega_C",PCSORSetOmega_SOR);CHKERRQ(ierr);
475   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSORSetIterations_C",PCSORSetIterations_SOR);CHKERRQ(ierr);
476   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSORGetSymmetric_C",PCSORGetSymmetric_SOR);CHKERRQ(ierr);
477   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSORGetOmega_C",PCSORGetOmega_SOR);CHKERRQ(ierr);
478   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSORGetIterations_C",PCSORGetIterations_SOR);CHKERRQ(ierr);
479   PetscFunctionReturn(0);
480 }
481 
482 
483 
484 
485 
486