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.