[petsc-users] meaning of constants in classical iterations

[Patrick Sanan]

On Fri, Jan 20, 2017 at 2:13 AM, Ed Bueler <elbueler at alaska.edu> wrote:
> Dear PETSc --
>
> In my humble opinion, the answers to the following rather basic questions
> are missing from the petsc users manual.  At least in part, I am not certain
> what the answers are.
>
>
> Question 1:   What formula does the preconditioned Richardson iteration
> (-ksp_type richardson -pc_type X) satisfy and where does the scale
> (-ksp_richardson_scale alpha) go?
>
> Answer?:  In the left-preconditioned form I would guess it is
>
> (*)         x_{k+1} = x_k + alpha M_L^{-1} (b - A x_k),
>
> where M_L is the matrix implied by -pc_type X.  If so this makes option
> combination
>
> -ksp_type richardson -pc_type jacobi [-ksp_richardson_scale 1]
>
> into the classical Jacobi iteration
>
> x_{k+1} = x_k + D^{-1} (b - A x_k) = D^{-1} (b - (L+U)x_k)
>
> where A = D+L+U and D is the diagonal of A.
>
> I am pretty sure this answer is right, but equation (*) should probably be
> somewhere in the manual!
It is at least on the KSPRICHARDSON man page (which I have added to my
list of things to eventually clean up), albeit without explicit
mention of the left preconditioning.
>
>
> Question 2.  What formula is the (non-symmetric) SOR satisfying, and where
> do the constants -pc_sor_omega and -pc_sor_diagonal_shift delta go?
>
> Answer?:  My understanding is first that the classical Jacobi, Gauss-Seidel,
> and SOR iterations are all regarded by PETSc developers as
> left-preconditioned Richardson iterations, so names "jacobi" and "sor"
> describe PC objects and these classical iterations are all KSP type
> "richardson".  That is, they are all instances of (*).
>
> (RANT:  There is no clear statement of this excellent philosophy in the
> users manual, and it would help a lot of students of PETSc!  While some
> references share it, reading standard literature, including relevant wiki
> pages, is not healthful on this topic.  Because the literature does not
> necessarily parallel PETSc thinking, it requires constant mental
> translation.  This is bad when it comes to setting parameters at runtime!)
This would be very good for one of the proposed introductory
tutorials, about KSP/PC, helping people translate algorithms into
PETSc options. I've added some notes on this to my todo list as well.
>
> Supposing the above philosophy, and that A = D+L+U with the usual meanings,
> then I guess SOR is
>
> (**)    x_{k+1} = x_k + alpha (omega^{-1} (D+delta) + L)^{-1} (b - A x_k),
>
> In case alpha=1, omega=1, delta=0 then (**) becomes Gauss-Seidel:
>
> -ksp_type richardson -pc_type sor [-ksp_richardson_scale 1 -pc_sor_omega 1
> -pc_sor_diagonal_shift 0]
>
> Especially because it involves three user-controlled constants, equation
> (**) could usefully be stated somewhere in the manual ... supposing it is
> correct.

Would you have any objection to this sort of thing being on the man
pages (and/or in tutorials), instead of in the manual? Especially with
the current state of affairs, with the  man pages google-able and the
manual not, the man pages might be a better place for these specifics
- hopefully it's clear to readers of the manual that they can click
links to be taken to man pages (if they have web access), which can
allow for finer details without interrupting the flow of the manual
too much or making it even longer.
>
>
> Question 3.  What is the probability that a PETSc developer will address the
> above two questions by telling me that classical iterations are lame?
>
> Answer:  Not the point!  These iterations are used blockwise (e.g. ASM) and
> as smoothers in many proper PETSc applications.  Clarity is helpful!
>
>
> Thanks for the great tool, as usual!
>
> Ed
>
>
> --
> Ed Bueler
> Dept of Math and Stat and Geophysical Institute
> University of Alaska Fairbanks
> Fairbanks, AK 99775-6660
> 301C Chapman and 410D Elvey
> 907 474-7693 and 907 474-7199  (fax 907 474-5394)