IterativeSolvers.jl
IterativeSolvers.jl copied to clipboard
`bicgstabl` returns NaNs for identity matrix
Julia 1.3, IterativeSolvers v0.8.1
julia> bicgstabl([1. 0; 0 1], [1.,2])
2-element Array{Float64,1}:
NaN
NaN
That's a bit of an edge case when solving a probably overdetermined system internally. If you want to cover this case, please run bicgstabl through Debugger.jl to see where the first division by zero occurs. I would not worry about it though, since it will not happen in real life.
I mean, it happened to me in real life ;)
Anyway, I don't think I'm qualified to help further with this and indeed its an edge case, but just wanted to bring this to your attention.
Ah ok, that makes it more interesting. Do you remember whether it happened when an exact solution was found?
Basically what happens in bicgstabl is that multiple steps of gmres are done at once, and a small low-dimensional system is inverted, so my guess would be that as soon as a true solution is found, this will break down.
I mean that this was the exact thing that happened to me, so all I have is the above (reduced) example. Although its obviously trivial, in my real case the matrix is the output of some other computation, which at various times does output an identity matrix (which then subsequently gives NaNs).
I just ran into this too. And noted that matlab handles this at least
>> A = eye(5,5);
>> b = [1; 2; 3; 4; 5];
>> bicgstab(A, b)
bicgstab converged at iteration 0.5 to a solution with relative residual 0.
ans =
1
2
3
4
5
>> bicgstabl(A, b)
bicgstabl converged at iteration 0.2 to a solution with relative residual 0.
ans =
1
2
3
4
5