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