Open ChrisRackauckas opened 4 years ago
10s vs 0.5s
I have another issue pending to be submitted to DiffEqBase
which could explain this. Somehow, I got the feeling that the tolerance tol
is not well handled in function (f::LinSolveIterativeSolvers)(x,A,b,update_matrix=false; Pl=nothing, Pr=nothing, tol=nothing, kwargs...)
. You can see that if you enable the log in IterativeSolvers, hence https://github.com/JuliaDiffEq/DiffEqBase.jl/issues/366
Yes, i printed the tolerance for a small case and it is weirdly high:
example:
"""
OrdinaryDiffEq/test/interface/linear_nonlinear_tests.jl
"""
n = 10
A = Tridiagonal(-ones(n-1),2ones(n),-ones(n-1))
rn = (du, u, p, t) -> begin
mul!(du, A, u)
end
u0 = rand(n)
prob = ODEProblem(ODEFunction(rn, jac_prototype=JacVecOperator{Float64}(rn, u0; autodiff=false)), u0, (0, 10.))
sol = solve(prob, QNDF(autodiff=false, linsolve=LinSolveGMRES()));
with DiffEqBase/linear_nonlinear.jl
modified as follows:
function (f::LinSolveIterativeSolvers)(x,A,b,update_matrix=false; Pl=nothing, Pr=nothing, reltol=eps(eltype(x)), kwargs...)
if f.iterable === nothing
Pl = ComposePreconditioner(f.Pl, Pl, true)
Pr = ComposePreconditioner(f.Pr, Pr, false)
reltol = checkreltol(reltol)
f.iterable = f.generate_iterator(x,A,b,f.args...;
initially_zero=true,restart=5,
maxiter=20,abstol=1e-12,reltol=reltol,
Pl=Pl,Pr=Pr,
kwargs...,f.kwargs...)
end
x .= false
iter = f.iterable
purge_history!(iter, x, b)
for residual in iter
println("GMRES step $(iter.k), residual: $(iter.residual.current), tol=$(iter.tol)")
end
return nothing
end
julia> include("krylov_test.jl")
GMRES step 1, residual: 0.004468827765402508, tol=5.481711392696291
GMRES step 1, residual: 0.0006976753003995138, tol=5.481711392696291
GMRES step 1, residual: 0.00010892502571461855, tol=5.481711392696291
GMRES step 1, residual: 0.0067410564826430016, tol=5.481711392696291
GMRES step 1, residual: 0.0031643581369913687, tol=5.481711392696291
GMRES step 1, residual: 0.0009081523038929753, tol=5.481711392696291
GMRES step 1, residual: 0.00013766281617396835, tol=5.481711392696291
GMRES step 1, residual: 0.0034366200233135642, tol=5.481711392696291
GMRES step 1, residual: 0.0018549989980046486, tol=5.481711392696291
GMRES step 1, residual: 0.0005839668067025519, tol=5.481711392696291
GMRES step 1, residual: 0.005929950997054515, tol=5.481711392696291
GMRES step 1, residual: 0.16739523946435644, tol=5.481711392696291
GMRES step 1, residual: 0.09269193011431896, tol=5.481711392696291
GMRES step 1, residual: 0.0733150942936091, tol=5.481711392696291
GMRES step 1, residual: 0.07245774505656238, tol=5.481711392696291
GMRES step 1, residual: 0.3740037683520435, tol=5.481711392696291
GMRES step 1, residual: 1.571981217219508, tol=5.481711392696291
Why iter.tol is 5.4817113??
Otherwise, the behavior is actually normal. The matrix is badly conditioned. See:
using BifurcationKit
n = 10
A = Tridiagonal(-ones(n-1),2ones(n),-ones(n-1)) |> Array
rhs = rand(n)
ls = GMRESIterativeSolvers(verbose = true, log = true, N = n)
sol = ls(A, rhs)
gives
=== gmres ===
rest iter resnorm
1 1 1.65e+00
1 2 1.51e+00
1 3 1.21e+00
1 4 9.02e-01
1 5 5.33e-01
1 6 3.43e-01
1 7 1.66e-01
1 8 1.11e-01
1 9 4.21e-02
1 10 5.86e-15
Better example for iterative methods:
n = 100
A = I + sprand(n,n,1/n)
rhs = rand(n)
ls = GMRESIterativeSolvers(verbose = true, log = true, N = n)
sol = ls(A, rhs)
gives
=== gmres ===
rest iter resnorm
1 1 2.34e+00
1 2 1.34e+00
1 3 5.37e-01
1 4 1.56e-01
1 5 5.98e-02
1 6 1.88e-02
1 7 8.50e-03
1 8 1.80e-03
1 9 7.18e-04
1 10 1.51e-04
1 11 4.92e-05
1 12 4.90e-06
1 13 2.41e-06
1 14 4.88e-07
1 15 4.78e-08
I apologize to use BifurcationKit for the linear solver but I am very accustomed to the interface which is unified for several packages. It should be similar to LinSolveIterativeSolvers except I do not use the iterator form of
gmres
@YingboMa