How to efficiently update a Sparse Matrix (Diagonal + Sparse) in-place inside a hot loop?

General Usage
, ,
1

Hi everyone,

I am working on an iterative algorithm where I need to construct and solve a linear system inside a “hot loop”.

Currently, my code allocates a new sparse matrix A at every iteration because of the addition between a Diagonal matrix and a SparseMatrix.

something like:

using LinearAlgebra, SparseArrays,  BenchmarkTools

function mwe(y::Vector{T},  max_iter::Int=10) where T
    m = length(y)
    
    D = spdiagm(0 => -2.0 * ones(T, m), 1 => ones(T, m-1), -1 => ones(T, m-1))
    AA = D' * D
    
    # Buffers
    w = ones(T, m)
    tmp_rhs = similar(y)
    result = similar(y)
    residuals = similar(y)

    for _ in 1:max_iter
        
        A = Diagonal(w) * AA

       
        tmp_rhs .= w .* y
        result .= A \ tmp_rhs

       # not used in mwe
        residuals .= y .- result
        w .= 1.0 ./ (1.0 .+ abs.(residuals)) 
    end
    
    return result
end

ideally, i would like to update results in-place.
How can I do that?

2

Just to clarify, inside your code I see a multiplication between diagonal and sparse, not an addition. Which is it?