# IterativeSolvers.jl not working as expected with mutating and allocation
**URL:** https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277
**Category:** General Usage
**Tags:** benchmarktools, preallocation, iterative-solvers
**Created:** [February 16, 2024, 2:40am UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277 "2024-02-16T02:40:33Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [February 16, 2024, 2:40am UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/1 "2024-02-16T02:40:33Z")
I’m trying to solve a simple linear system. However, it seems like the preallocated versions of any IterativeSolvers.jl functions are either much slower or allocating more memory.
Here’s the example:
```julia
using IterativeSolvers
using LinearMaps
function SecondOrderCentralDiffmul!(C, B)
C[1] = -2B[1] + B[2]
for i in 2:length(B)-1
C[i] = B[i-1] - 2B[i] + B[i+1]
end
C[end] = B[end-1] - 2B[end]
return C
end
nnn = 51
A = LinearMap(SecondOrderCentralDiffmul!, nnn; issymmetric=true, ismutating=true)
b = rand(nnn)
U = gmres(A, b,maxiter=10000000)
norm(A*U - b)
```
With the following tests:
`@btime gmres(A, b,maxiter=10000000);`
`697.949 μs (43 allocations: 19.88 KiB)`
And mutating:
```julia
xxx = zeros(Float64, nnn)
@btime gmres!(xxx, A, b,maxiter=10000000);
```
`13.037 s (500015 allocations: 106.82 MiB)`
However, when I run without `@btime` :
`@time gmres!(xxx, A, b,maxiter=10000000);`
`0.000785 seconds (42 allocations: 19.391 KiB)`
So first, what’s going on with `@btime` for these tests? Second, why is the allocations the same for both `gmres` and `gmres!` ?
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [February 16, 2024, 10:38am UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/2 "2024-02-16T10:38:19Z")
This is actually not related to IterativeSolvers, but to BenchmarkTools subtleties, mostly because it runs the function several times to get a better estimate of the time.
First, as a rule, you should always [interpolate global variables](https://juliaci.github.io/BenchmarkTools.jl/stable/manual/#Interpolating-values-into-benchmark-expressions) when benchmarking.
Furthermore, running the mutating version several times is biased. Indeed, for `gmres!`, the first argument is not just a scratch space: it is also the initial guess. Thus, once you have run the solver, the initial guess actually contains the optimal solution, and this seems to cause weird behavior (although I’m not sure why we see a slowdown and not a speedup).
To escape this issue, you need to
- create a [`setup` phase](https://juliaci.github.io/BenchmarkTools.jl/stable/manual/#Setup-and-teardown-phases) in your benchmark to initialize it with a new zero vector every time
- set `evals = 1` in the [benchmark parameters](https://juliaci.github.io/BenchmarkTools.jl/stable/manual/#Benchmark-Parameters) to make sure that each benchmark sample runs `gmres!` only once, thus avoiding the bias
When you do all that, you realize that `@btime` does indeed give you a more accurate (and lower) result than `@time`:
```julia
julia> @btime gmres($A, $b, maxiter=10000000);
603.156 μs (43 allocations: 19.86 KiB)
julia> @btime gmres!($xxx, $A, $b, maxiter=10000000);
12.053 s (500015 allocations: 106.82 MiB)
julia> @btime gmres!(_xxx, $A, $b, maxiter=10000000) evals=1 setup=(_xxx = zeros(nnn));
603.605 μs (43 allocations: 19.41 KiB)
julia> @time gmres(A, b, maxiter=10000000);
0.000894 seconds (43 allocations: 19.859 KiB)
```
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [February 16, 2024, 5:14pm UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/3 "2024-02-16T17:14:53Z")
@gdalle awesome!
One last thing is that I still don’t understand why the mutating `gmres!` is allocating the same amount of memory as non-mutating version?
---
### Author:  [@leespen1](https://discourse.julialang.org/u/leespen1)
#### Post date: [February 16, 2024, 6:39pm UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/4 "2024-02-16T18:39:48Z")
I would also note that mutating `gmres!` allocates memory even without using `LinearMaps`.
```julia
julia> using IterativeSolvers, BenchmarkTools
julia> function f!(x, A, b, n_runs=1)
for i in 1:n_runs
x .= 0
gmres!(x, A, b)
end
return nothing
end
f! (generic function with 2 methods)
julia> x = zeros(2)
2-element Vector{Float64}:
0.0
0.0
julia> A = [1.0 2.0; 3.0 4.0]
2×2 Matrix{Float64}:
1.0 2.0
3.0 4.0
julia> b = [5.0, 6.0]
2-element Vector{Float64}:
5.0
6.0
julia> f!(x, A, b, 1) # Run once to compile
julia> @btime f!(x, A, b, 1)
1.393 μs (15 allocations: 1.23 KiB)
julia> @btime f!(x, A, b, 10)
13.151 μs (150 allocations: 12.34 KiB)
julia> @btime f!(x, A, b, 100)
131.335 μs (1500 allocations: 123.44 KiB)
julia> @btime gmres(A, b)
1.374 μs (16 allocations: 1.31 KiB)
```
As you can see, each call to `gmres!` makes 15 allocations even for the simplest possible case of ordinary matrices and vectors. `gmres` makes 16 allocations, one more to allocate the return vector. But most of the allocations appear to be used internally by the algorithm.
It would be nice if we had an option to provide a ‘cache’ variable, so that `gmres!` didn’t have to allocate memory each time it ran.
The good news is that (at least for this small example), the time taken to allocate memory doesn’t tank performance.
```julia
julia> function allocate_memory(n_runs; alloc_size=1)
for i in 1:n_runs
a = Vector{Int64}(undef, alloc_size)
a[1] += 1 # If we don't do some operation on a, the compiler won't allocate a in the first place
end
end
allocate_memory (generic function with 2 methods)
julia> @btime allocate_memory(15, alloc_size=1)
258.078 ns (15 allocations: 960 bytes)
julia> @btime allocate_memory(15, alloc_size=1000)
1.601 μs (15 allocations: 119.06 KiB)
julia> @btime f!(x, A, b, 1) # Run once to compile
1.396 μs (15 allocations: 1.23 KiB)
julia> @btime allocate_memory(15, alloc_size=2)
263.934 ns (15 allocations: 1.17 KiB)
```
It looks like `alloc_size=2` allocates the amount of memory closest to that of `gmres!` (probably allocates 15 vectors the size of the system, or something like that). So I would estimate that in this example 0.263 microseconds of the 1.393 microsecond runtime is spent allocating memory. So `gmres!` takes maybe 25% longer than it would have if we managed to avoid allocating memory.
This could change significantly depending on the values and sizes of `A` and `b` (perhaps more allocations are performed if more iterations are required, I don’t know the internals of the function). Using `gmres!` to solve a 2x2 system `Ax=b`, where `A` and `b` are known explicitly, is not a typical use case. (backslash solves the same problem in \< 0.3 microseconds, and I think almost all of that time is spent allocating memory).
I tried a random 1000x1000 example, and gmres took ~100 ms, while I estimate the memory allocation took \< 0.1 ms.
Still, it would be nice to have a ‘cached’ option, just to be sure that memory allocation isn’t significantly affecting performance.
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [February 16, 2024, 7:13pm UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/5 "2024-02-16T19:13:43Z")
> [@leespen1](#):
>
> Still, it would be nice to have a ‘cached’ option, just to be sure that memory allocation isn’t significantly affecting performance.
I definitely agree. The fact that we can’t provide this option makes it so that for many specific cases the user will have to write their own Iterative solver in order to manage the memory and speed better.
Thanks again for the answers!
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [February 16, 2024, 9:36pm UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/6 "2024-02-16T21:36:25Z")
Actually I think there is an advanced interface where you can avoid repeated allocation. See here: [The iterator approach · IterativeSolvers.jl](https://iterativesolvers.julialinearalgebra.org/dev/iterators/#Example:-avoiding-unnecessary-initialization)
---
### Author:  [@erny123](https://discourse.julialang.org/u/erny123)
#### Post date: [February 18, 2024, 5:21pm UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/7 "2024-02-18T17:21:52Z")
@abraemer wow completely missed this section. Good catch
---
### Author:  [@leespen1](https://discourse.julialang.org/u/leespen1)
#### Post date: [March 5, 2024, 9:18pm UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/8 "2024-03-05T21:18:50Z")
You’re right, thanks!
I would note that it’s not immediately obvious how to avoid repeated allocation with gmres, since the user has to repeat some of the initialization process when setting a new value for `b`. But some of the solvers, like jacobi, are much easier (just set a new `b` and do the iteration).
Here is a function which I think is sufficient to ‘reinitialize’ the gmres iterable with a new initial guess `x` and a new right-hand-side `b`:
```julia
function update_gmres_iterable!(iterable, x, b)
iterable.b .= b
iterable.x .= x
iterable.mv_products = 0
iterable.arnoldi.H .= 0
iterable.arnoldi.V .= 0
iterable.residual.accumulator = 1
iterable.residual.current = 1
iterable.residual.nullvec .= 1
iterable.residual.β = 1
iterable.residual.current = IterativeSolvers.init!(
iterable.arnoldi, iterable.x, iterable.b, iterable.Pl, iterable.Ax,
initially_zero=false
)
iterable.residual.nullvec .= 1
IterativeSolvers.init_residual!(iterable.residual, iterable.residual.current)
iterable.β = iterable.residual.current
return nothing
end
```
And here is an example of it in action:
```julia
julia> A = [1.0 2;3 4];
julia> b = [2.0, 4];
julia> x = [-1.0, 1];
julia> gmres_iter = IterativeSolvers.gmres_iterable!(x, A, b, abstol=1e-10, reltol=1e-10);
julia> for (i, iter) in enumerate(gmres_iter)
println("iteration $i done")
end
iteration 1 done
iteration 2 done
julia> A * x
2-element Vector{Float64}:
2.0
4.0
julia> x2 = [-2.0, 2.0]; b2 = [4.0, 8.0];
julia> update_gmres_iterable!(gmres_iter, x2, b2)
julia> for (i, iter) in enumerate(gmres_iter)
println("iteration $i done")
end
iteration 1 done
iteration 2 done
julia> A*x
2-element Vector{Float64}:
4.0
8.0
```
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [March 5, 2024, 9:40pm UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/9 "2024-03-05T21:40:39Z")
That seems like useful information that could be preserved. I suggest opening a pull request to IterativeSolvers.jl to put this into the documentation as an example usage of the iterator interface 🙂
---
### Author:  [@leespen1](https://discourse.julialang.org/u/leespen1)
#### Post date: [March 5, 2024, 9:49pm UTC](https://discourse.julialang.org/t/iterativesolvers-jl-not-working-as-expected-with-mutating-and-allocation/110277/10 "2024-03-05T21:49:00Z")
Thanks, I’m working on it now.