DiffEqBase.jl
DiffEqBase.jl copied to clipboard
Error with preconditioning in LinSolveGMRES
I partially managed to partially make a wrap from PseudoArclLengthContinuation.jl
to DiffEq. The code is here.
There is an example where it is shown how to use the interface. It does not work very well for GMRES, see line 97 in the example ; it returns a LAPACKexception(1)
.
Possible error location
I think I tracked back the problem. To me the issue is with DiffEqBase.ComposePreconditioner{Identity,Nothing}(Identity(), nothing, false)
. Indeed, at some point LinSolveGMRES
calls expand!
in gmres.jl
from IterativeSolvers.jl
. In this function, the line ldiv!(nextV, Pr, view(arnoldi.V, :, k))
puts nextV=0
despite the fact that Pr
is the identity and view(arnoldi.V, :, k)
is not zero!!
possible fix
If I modify line 148 into:
f.iterable = f.generate_iterator(x,A,b,f.args...; initially_zero=true,restart=5, maxiter=5,tol=1e-16,Pl=Pl.P,Pr=Pr.P, f.kwargs...,kwargs...)
in order to remove the use of ComposePreconditioner
, then the error disappears. Hence, to me, there is a bug with the preconditioning !!
https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/953 might be related too. Something seems off making it slow.
Any clue about this issue? It seems in 366, the options are not set properly either.
@yingboma