JuliaLinearAlgebra / IterativeSolvers.jl

Iterative algorithms for solving linear systems, eigensystems, and singular value problems
MIT License
399 stars 105 forks source link

[cg] Default maxiter off-by one? #257

Open fredrikekre opened 4 years ago

fredrikekre commented 4 years ago

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.

KristofferC commented 4 years ago

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.