IterativeSolvers.jl
IterativeSolvers.jl copied to clipboard
[cg] Default maxiter off-by one?
Example:
julia> using LinearAlgebra, IterativeSolvers, Random
julia> Random.seed!(1234);
A = rand(10,10) + I; A = A'A;
b = rand(10);
julia> x = zeros(10); cg!(x, A, b; verbose = true, maxiter=100);
1 7.77e-01
2 5.46e-01
3 3.72e-01
4 2.06e-01
5 1.58e-01
6 3.01e-01
7 1.97e-01
8 1.28e-01
9 1.38e-02
10 2.28e-01
11 3.37e-13
Looks like the 11th step is needed.
And explicitly:
julia> x = zeros(10); cg!(x, A, b); norm(A*x-b)
0.10052838755750378
julia> x = zeros(10); cg!(x, A, b; maxiter=11); norm(A*x-b)
2.972481305376441e-12
There seems to be some counting error since after 10 iterations we should have the exact solution.