Parallel sparse solve and inverse

Specific Domains Numerics
1

So i need to solve a really big (and sparse) system of equation and to my knowledge there are no parallel sparse solvers for julia. I tried the pardiso for MKL but that wasn’t parallel and had performance similar to the one used by julia. What i do right now is to do an lu factorization of my lhs (lu is the factorization chosen by factorize) and then use something like

"solve Ax = b"
function par_solv!(A, b, x=similar(b))
    blocks = Threads.nthreads()
    block_size = ceil(Int, size(x, 2) / blocks)
    @threads for i in 1:blocks
            s = (i-1)*block_size+1:min(i * block_size, size(x, 2))
            ldiv!(@view(x[:, s]), A, @view(b[:, s]))
    end
    return x
end

This works pretty well, but can i do better? Sometimes i need the inverse of lhs and i give a dense identity matrix to this function. Can i do better?

I only with work with SparseMatrixCSC and standard arrays. i don’t care, nor plan to support any other format.

1 Like