# Avoiding memory allocation for solves with multiple right-hand sides

**URL:** <https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879>\
**Category:** Performance\
**Tags:** linearalgebra, memory-allocation\
**Created:** [September 8, 2021, 3:36pm UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879 "2021-09-08T15:36:34Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [September 8, 2021, 3:36pm UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879/1 "2021-09-08T15:36:34Z")

</div>

Hi everyone,

I am new to Julia and I try to understand why a solve using the backslash operator allocates memory although the RHS vector is allocated, the solution vector is allocated and the factorisation is precomputed (here Cholesky factorisation). Is there a way to avoid such behaviour?

Below is a minimum example.

Cheers and thanks for the hints!

```julia
using LinearAlgebra, SparseArrays

function Successive_solves!( x, Kc, b, nsolves )
    for i=1:nsolves
        x .= Kc\b 
    end 
    return
end

function ThreePointStencil( N, dx )
    # Assemble matrix (not the part of interest)
    cC = 2.0/dx^2 * ones(Float64, N); cC[1] = 3.0/dx^2; cC[end] = 3.0/dx^2
    cW = -1.0/dx^2 * ones(Float64, N); cW[1] = 0.0
    cE = -1.0/dx^2 * ones(Float64, N); cE[end] = 0.0
    iC = 1:N
    iW = ones(Int64, N); iW[2:end-0] .= iC[1:end-1]
    iE = ones(Int64, N); iE[1:end-1] .= iC[2:end-0]
    I = [iC[:]; iC[:]; iC[:]]
    J = [iC[:]; iW[:]; iE[:]]
    V = [cC[:]; cW[:]; cE[:]]
    K = sparse( I, J, V )
    droptol!(K, 1e-9)
    return K
end

# Input
N = 1000000
dx = 1

# Sparse matrix
K = ThreePointStencil( N, dx )

# Allocate arrays
b = rand(N) # Right-hand-side array 
x = zeros(Float64, N) # Solution array

# Factor
Kc = cholesky(K)

# One solve
@time x .= Kc\b # why does this allocate memory although the factorisation done and the vectors (x, b) are allocated?

# Successive solves
nsolves = 10
@time Successive_solves!( x, Kc, b, nsolves ) # then, memory allocation depends on the number of solves...

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [September 8, 2021, 4:59pm UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879/2 "2021-09-08T16:59:46Z")

</div>

It allocates because the operation `Kc\b` allocates, on the right-hand side. You either have to use an in-place version, or if the matrices are small, use static arrays:

```julia
julia> M = rand(3,3); b = rand(3);

julia> @btime x = $M \ $b;
  346.107 ns (3 allocations: 384 bytes)

julia> using StaticArrays

julia> M = @SMatrix(rand(3,3)); b = @SVector(rand(3));

julia> @btime x = $M \ $b;
  5.228 ns (0 allocations: 0 bytes)

```

Probably this is what you are looking for: [Non-allocating matrix inversion - #7 by stevengj](https://discourse.julialang.org/t/non-allocating-matrix-inversion/62264/7)

---

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [September 8, 2021, 5:51pm UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879/3 "2021-09-08T17:51:01Z")

</div>

Thanks a lot for the insights!

I guess I would really need an in-place version of backslash rather than using StaticArrays since I would like to work with large systems.  
So, is there an in-place version of backslash existing or does one really need to call low-level LAPACK function like in the stevengj post?

Initially I was looking for an in-place version of cs\_ltsolve() and cs\_lsolve() from SuiteSparse but I’ve not yet found how to call them within Julia. These functions are really nice to use in combination with Cholmod.

---

<div class="post-metadata">

**Author:** ![Vasily\_Pisarev](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vasily_pisarev/32/7929_2.png) [@Vasily\_Pisarev](https://discourse.julialang.org/u/Vasily_Pisarev)\
**Post date:** [September 8, 2021, 5:57pm UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879/4 "2021-09-08T17:57:22Z")

</div>

> [@Gravlax](#):
>
> So, is there an in-place version of backslash existing or does one really need to call low-level LAPACK function like in the stevengj post?

`LinearAlgebra.ldiv!` is what you need.

> `ldiv!(Y, A, B) -> Y`  
> Compute `A \ B` in-place and store the result in `Y`, returning the result.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 8, 2021, 7:06pm UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879/5 "2021-09-08T19:06:27Z")

</div>

I would catch the LU of Kc and use it with `ldiv!`. A bit like

```julia
function Successive_solves!( x, Kc, b, nsolves )
Kclu = lu(Kc)
    for i=1:nsolves
        ldiv!(x, Kclu, b) 
    end 
    return
end

```

---

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [September 8, 2021, 7:27pm UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879/6 "2021-09-08T19:27:00Z")

</div>

Thank you all!

OK, get it ldiv!() to work with an LU factorisation, cool! However, it does not seem to work with a Cholesky factorisation as it returns:  
ERROR: LoadError: MethodError: no method matching rdiv!(::Vector{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64}, ::Vector{Float64})

Actually, if one inspires from the example given in the doc of ldiv!() and changes A to a symmetric positive-definite matrix, it seems that ldiv!() can work with Cholesky if A is dense but not when it is sparse, see below:

```julia
A = [2.0 -1.0 0.0; -1.0 2.0 -1.0; 0.0 -1.0 2.0];
Ac = cholesky(A) # works
# Ac = cholesky(sparse(A)) # does not work
X = [1; 2.5; 3];
Y = zero(X);
ldiv!(Y, Ac, X);

```

Any hints?

EDIT: Actually it looks like there is no trivial was to combine ldiv!() with a sparse Cholesky factorisation:  
[https://github.com/JuliaLang/SuiteSparse.jl/issues/19](https://github.com/JuliaLang/SuiteSparse.jl/issues/19)  
I wonder whether the original SuiteSparse CS functions (like [cs\_ltsolve](https://github.com/fieldtrip/fieldtrip/blob/master/external/dmlt/external/gpstuff/SuiteSparse/CSparse/Source/cs_ltsolve.c)) are available from within Julia’s SuiteSparse implementation…

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [September 9, 2021, 8:11am UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879/7 "2021-09-09T08:11:17Z")

</div>

I did that for my library, see [here](https://rveltz.github.io/BifurcationKit.jl/dev/tutorialsSH3d/#d-Swift-Hohenberg-equation,-Finite-differences)

```julia
Pr = cholesky(L1)
using SuiteSparse
# we need this "hack" to be able to use Pr as a preconditioner.
LinearAlgebra.ldiv!(P::SuiteSparse.CHOLMOD.Factor{Float64}, v) = -(P \ v)

# rtol must be small enough to pass the Fold points and to get precise eigenvalues
# we know that the jacobian is symmetric so we tell the solver
ls = GMRESKrylovKit(verbose = 0, rtol = 1e-9, maxiter = 150, ishermitian = true, Pl = Pr)

```

---

<div class="post-metadata">

**Author:** ![Gravlax](https://avatars.discourse-cdn.com/v4/letter/g/dbc845/32.png) [@Gravlax](https://discourse.julialang.org/u/Gravlax)\
**Post date:** [September 10, 2021, 8:57am UTC](https://discourse.julialang.org/t/avoiding-memory-allocation-for-solves-with-multiple-right-hand-sides/67879/9 "2021-09-10T08:57:39Z")

</div>

Thanks! If I get it right you actually overlad the ldiv!() function such that it can work with Cholesky factors. Is that right?

Would there also be a hack for the 3 argument version of ldiv!() such that the entries of v are not replaced by the entries of the solution vector?