NaN in bicgstabl when using preconditioners

[davidschlegel]

In some cases, using a preconditioner results in NaNs in the bicgstabl function. Tested with ilu preconditioner and DiagonalPreconditioner. Here is an example:

using SparseArrays
using IncompleteLU
using IterativeSolvers

N = Int(1e5) #
d = 1e-4  # sparseness density
drop_tol = 0.1 # absolute drop tolerance for ilu
mat = sprand(N, N, d)
x = ones(N)  # solution vector
b = mat * x

prec = ilu(mat, τ=drop_tol)
@show nnz(prec)/nnz(mat)
sol, hist = bicgstabl(mat, b, Pl = prec, log = true, verbose = true)

[deltaeecs]

So do mine!

[davidschlegel]

I now think, it's actually due to the problem and not because of the algorithm itself. There are many reasons why bicgstabl might fail. One of them is the structure of the problem, i.e. the properties of the matrix. In some cases I find that making the system diagonally dominant, for instance by doing a reverse-Cuthill-McKee reordering of indices, the solver converges when the initial problem does not converge. It could be interesting though to have such a feature of index reordering as a parameter to the solver which is known in the numerics community to work well on a large class of problems.