Extra memory usage that cannot be released from GC when reuse_symbolic=true

[ytdHuang]

Hi,

I am now designing a package that needs to solve multiple times of A \ b, but with different As. This package really benefits a lot, thanks for the hard work of the develop team.

However, I faced some memory issue when using the cache interface. Let me simplify my problem and give an example below.

I have a function called DOS which needs a sparse matrix $M$, two vectors ( $b_+$ and $b_-$ ), and a list of $\omega$. The linear solve for this problem is to solve $x$ which satisfies $$\left(M \pm i \omega I \right)x = b_\pm$$ for different $\omega$. Here, $I$ is an identity matrix.

The definition of DOS is given below:

using  LinearSolve, SparseArrays, LinearAlgebra
import LinearSolve: set_A

function DOS(
        M::SparseMatrixCSC{ComplexF64, Int64}, 
        bp::Vector{ComplexF64},
        bm::Vector{ComplexF64},
        ωlist;
        solver,
        SOLVEROptions...
    )
    
    dim, = size(M)
    II = sparse(I, dim, dim)  # identity matrix
    result = Vector{Float64}(undef, length(ωlist)) # a vector which is used to store the results
    
    # solve for the first ω
    j = 1
    iωI = 1im * ωlist[j] * II
    
    ## (M - iωI) \ bm
    prob_m = init(LinearProblem(M - iωI, bm), solver, SOLVEROptions...)
    sol_m = solve(prob_m)
    
    ## (M + iωI) \ bp
    prob_p = init(LinearProblem(M + iωI, bp), solver, SOLVEROptions...)
    sol_p = solve(prob_p)
    
    result[j] = real(sum(sol_p.u) - sum(sol_m.u))
    
    # solve for the rest ω
    @inbounds for ω in ωlist[2:end]
        j += 1
        iωI = 1im * ω * II

        # solve (M - iωI) \ bm
        # if the sparsity doesn't change, use the cache from previous solution
        # else, initialize a new linear problem
        try 
            sol_m = solve(set_A(sol_m.cache, M - iωI))
        catch e
            if isa(e, ArgumentError)
                prob_m = init(LinearProblem(M - iωI, bm), solver, SOLVEROptions...)
                sol_m  = solve(prob_m)
            else
                throw(e)
            end
        end

        # solve (M + iωI) \ bp
        # if the sparsity doesn't change, use the cache from previous solution
        # else, initialize a new linear problem
        try
            sol_p = solve(set_A(sol_p.cache, M + iωI))
        catch e
            if isa(e, ArgumentError)
                prob_p = init(LinearProblem(M + iωI, b_p), solver, SOLVEROptions...)
                sol_p = solve(prob_p)
            else
                throw(e)
            end
        end
        result[j] = real(sum(sol_p.u) - sum(sol_m.u))
    end
    
    return result
end

The reason why I'm using the exception handling is because the sparsity pattern might be different when the $\omega$ changes (especially for $\omega = 0$).

Now, the memory issue pops out in the following testing:

julia> dim = 10000;
julia> density = 0.005;
julia> M  = sprand(ComplexF64, dim, dim, density);
julia> bp = rand(ComplexF64, dim);
julia> bm = rand(ComplexF64, dim);
julia> ωlist = -2:1:2; # [-2, -1, 0, 1, 2];

julia> # memory usage: 631 MB
julia> @time result = DOS(M, bp, bm, ωlist; solver=UMFPACKFactorization(;reuse_symbolic=false));
211.988220 seconds (5.71 M allocations: 62.507 GiB, 0.70% gc time, 2.70% compilation time)
julia> # memory usage: 13.02 GB
julia> GC.gc()
julia> # memory usage: 862 MB
julia> @time result = DOS(M, bp, bm, ωlist; solver=UMFPACKFactorization(;reuse_symbolic=true));
188.753707 seconds (876 allocations: 61.881 GiB, 0.11% gc time)
julia> # memory usage: 7.97 GB
julia> GC.gc()
julia> # memory usage: 3.29 GB

The first time I execute DOS, I set reuse_symbolic=false. For the second time, I set reuse_symbolic=true, it indeed speed up a little bit, and the number of allocations decreases a lot. However, when I clean the garbage collector after the second execution, there are some extra garbage (approximately 2.4 GB) that can not be cleaned and this value (the extra memory usage) dramatically increases when I increase the dimension of the matrix.

Did I missed anything when using the cache interface, or is there a better method to clean the garbage collection in the function so that it wouldn't produce so many un-releasable memory ?