DiffEqBase.jl
DiffEqBase.jl copied to clipboard
Put options to have log in Iterative solvers
When using Iterative solvers, it is very useful to have the log of convergence to attest the usefullness / efficiency of a preconditionner. The option is not present in the current file linear_nonlinear.jl.
I suggest to do the following unless you have a better way of triggering the verbose
option:
function (f::LinSolveIterativeSolvers)(x,A,b,update_matrix=false;verbose = false, Pl=nothing, Pr=nothing, tol=nothing, kwargs...)
if f.iterable === nothing
Pl = ComposePreconditioner(f.Pl, Pl, true)
Pr = ComposePreconditioner(f.Pr, Pr, false)
@warn "Romain modified here"
f.iterable = f.generate_iterator(x,A,b,f.args...;
initially_zero=true,restart=5,
maxiter=5,tol=1e-16,Pl=Pl,Pr=Pr,
f.kwargs...,kwargs...)
tol′ = get(f.kwargs, :tol, nothing)
tol′ !== nothing && (tol = tol′)
f.iterable.reltol = tol
end
x .= false
iter = f.iterable
purge_history!(iter, x, b)
verbose && @printf("=== gmres ===\n%4s\t%4s\t%7s\n","rest","iter","resnorm")
restart = get(f.kwargs, :restart, 5)
for (iteration, residual) = enumerate(iter)
verbose && @printf("%3d\t%3d\t%1.2e\n", 1 + div(iteration - 1, restart), 1 + mod(iteration - 1, restart), residual)
end
return nothing
end
and also:
function (p::DefaultLinSolve)(x,A,b,update_matrix=false;tol=nothing,verbose=false, kwargs...)
if p.iterable isa Vector && eltype(p.iterable) <: LinearAlgebra.BlasInt # `iterable` here is the pivoting vector
F = LU{eltype(A)}(A, p.iterable, zero(LinearAlgebra.BlasInt))
ldiv!(x, F, b)
return nothing
end
if update_matrix
if typeof(A) <: Matrix
blasvendor = BLAS.vendor()
if (blasvendor === :openblas || blasvendor === :openblas64) && size(A,1) <= 500 && ArrayInterface.can_setindex(x) # if the user doesn't use OpenBLAS, we assume that is a much better BLAS implementation like MKL
p.A = RecursiveFactorization.lu!(A)
else
p.A = lu!(A)
end
elseif typeof(A) <: Tridiagonal
p.A = lu!(A)
elseif typeof(A) <: Union{SymTridiagonal}
p.A = ldlt!(A)
elseif typeof(A) <: Union{Symmetric,Hermitian}
p.A = bunchkaufman!(A)
elseif typeof(A) <: SparseMatrixCSC
p.A = lu(A)
elseif ArrayInterface.isstructured(A)
p.A = factorize(A)
elseif !(typeof(A) <: AbstractDiffEqOperator)
# Most likely QR is the one that is overloaded
# Works on things like CuArrays
p.A = qr(A)
end
end
if typeof(A) <: Union{Matrix,SymTridiagonal,Tridiagonal,Symmetric,Hermitian} # No 2-arg form for SparseArrays!
x .= b
ldiv!(p.A,x)
# Missing a little bit of efficiency in a rare case
#elseif typeof(A) <: DiffEqArrayOperator
# ldiv!(x,p.A,b)
elseif ArrayInterface.isstructured(A) || A isa SparseMatrixCSC
ldiv!(x,p.A,b)
elseif typeof(A) <: AbstractDiffEqOperator
# No good starting guess, so guess zero
if p.iterable === nothing
p.iterable = IterativeSolvers.gmres_iterable!(x,A,b;initially_zero=true,restart=15,maxiter=15,tol=1e-16,kwargs...)
p.iterable.reltol = tol
end
x .= false
iter = p.iterable
purge_history!(iter, x, b)
verbose && @printf("=== gmres ===\n%4s\t%4s\t%7s\n","rest","iter","resnorm")
for residual in iter
end
else
ldiv!(x,p.A,b)
end
return nothing
end
Any chance to have this improved? The issue is that generate_iterator
does not accept the verbose
keyword. It has to be either shaved from f.kwargs
, or passed as above, or put the field verbose
in LinSolveIterativeSolvers
Also, the defaults should be handled more carefully, like (maybe)
f.iterable = f.generate_iterator(x,A,b,f.args...;
initially_zero=true,restart=get(f.kwargs, :restart, 5),
maxiter=get(f.kwargs, :maxiter, 5),tol=get(f.kwargs, :tol, 1e-16),Pl=Pl,Pr=Pr,
f.kwargs...,kwargs...)
Yeah that looks good to me.
would be great to allow for a user-defined function within the Krylov iteration.