search

SciML / DiffEqBase.jl

The lightweight Base library for shared types and functionality for defining differential equation and scientific machine learning (SciML) problems
Other
302 stars 106 forks source link

Error with preconditioning in LinSolveGMRES #365

Open rveltz opened 4 years ago

rveltz commented 4 years ago

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 !!

ChrisRackauckas commented 4 years ago

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/issues/953 might be related too. Something seems off making it slow.

rveltz commented 4 years ago

Any clue about this issue? It seems in 366, the options are not set properly either.

ChrisRackauckas commented 4 years ago

@yingboma