JuliaLinearAlgebra / IterativeSolvers.jl / 34848212235

98%

push

github

Fix the Chebyshev iteration recurrence (#376)

After #367 the iteration converges, but it still does not produce the
Chebyshev residual polynomial r_k = T_k((μI - A)/ρ) r_0 / T_k(μ/ρ):

- the iteration counter starts at 0, so the `iteration == 1` branch ran
  on the second step instead of the first, resetting the direction and
  using α = 2/μ there;
- `α = 1/(μ - β)` mixed a dimensionless β with an eigenvalue; the
  Templates form is α = 1/(μ - β/α_prev);
- the β used on the second step must be half the generic value, since the
  three-term recurrence starts from T_{-1} = T_1 rather than T_0.

Rewrite `iterate` in the coupled two-term form of Sleijpen's lecture
notes (Exercise 5.11), tracking ζ_k = T_k(μ/ρ) / T_{k-1}(μ/ρ) with
ζ_{k+1} = 2μ/ρ - 1/ζ_k, β_k = 1/ζ_k² and α_k = 2/(ρ ζ_{k+1}).

On a diagonal SPD matrix with spectrum in [1, 100] the residual norm after
60 iterations is now 1.2e-4 (the exact Chebyshev value), where the
previous code gave 2.2e-1 and, before #367, diverged to 4e10. The
existing tests did not catch this because their matrix has condition
number below 3.

Add a test that checks the residual history against the exact Chebyshev
polynomial evaluated per eigenvalue, with and without a diagonal left
preconditioner, for real and complex element types. Also fix the
`maxiter` docstring, which mentioned GMRES.

Fixes #375.

6 of 6 new or added lines in 1 file covered. (100.0%)

2027 of 2069 relevant lines covered (97.97%)

7446.69 hits per line