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