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