Hi all,

I am using Julia as a framework for development and testing of discretization methods in computational fluid dynamics and as a part of it, have to solve linear system of equations for velocities and pressure pretty often. Due to the fact I am eventually aiming at dense grids giving large linear systems, iterative solvers from the package `IterativeSolvers.jl`

seemed like the optimum choice for me.

However, at one point I noticed that the program I developed gives different results every time. No matter if I run it from terminal or from REPL. The differences are small, but still big enough to put that splinter in your mind that you might have failed to initialized all the variables in your code. After a thorough analysis, I realized that the source of randomness in my code is the `bicgstabl`

solver from `IterativeSolvers.jl`

Please have a look at an MWE:

```
using IterativeSolvers, SparseArrays, LinearAlgebra
n = 32
A = spdiagm(-1 => fill(-1.0, n - 1),
0 => fill( 3.0, n),
1 => fill(-2.0, n - 1))
Id = sparse(1.0I, n, n)
A = kron(A, Id) + kron(Id, A)
A = kron(A, Id) + kron(Id, A)
b = A * ones(n^3)
x = bicgstabl(A, b, 2, reltol=1.0e-5)
println("relative error = ", norm(b - A * x) / norm(b))
```

The output I get reads:

```
relative error = 9.393585200552226e-6
relative error = 2.364031042840345e-6
relative error = 9.663172610869453e-6
relative error = 3.716206357164069e-6
```

The relative error I get is invariantly smaller than the tolerance I set - which is good, but it is also different on every run - which is quite bad because it casts a shadow of a doubt on the code you developed around the solver.

A way forward in my case would be to change `bicgstabl`

from `IterativeSolvers.jl`

with the default `\`

operator from `SparseArrays.jl`

, which redirects to `umfpack`

solver and gives the same residual on every run. Another path would be to squeeze the `bicgstabl`

tolerance down until is below the number of significant digits in residuals I print, and hope the errors will not accumulate more than that during unsteady simulations of coupled non-linear equations I am solving.

But, really, is this expected behavior of an iterative solver? Living with an approximation - yes, but living with a random outcome on each run - no! It is certainly not something I have seen ever before in a linear solver and I would assume it is an indication that some arrays are not initialized properly in the package.

Any thoughts on that?

*Cheers*