Avoiding memory allocation for solves with multiple right-hand sides

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!

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...

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> 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

1 Like

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.

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.

4 Likes

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

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

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:

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:

I wonder whether the original SuiteSparse CS functions (like cs_ltsolve) are available from within Julia’s SuiteSparse implementation…

I did that for my library, see here

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)

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?