I need to get a sum of approx 10000 matrices, that are constructed by multiplication with sparse matrices.
using LinearAlgebra, SparseArrays
n = 101
p = 101
m = 25
#indK = sprand(m, n^2*p^2, 0.2) EDIT on density
indK = sprand(m, n^2*p^2, 0.01)
Sm = rand(m, m)
o = rand(n*p)
pom = zeros(n*p, n*p)
pom+= indK[:, (j-1)*n*p+1:j*n*p]'*Sm*indK[:, (j-1)*n*p+1:j*n*p]*o[j]
This is extremely slow, so I tranformed the line in for loop into
mul!(pom, indK[:, (j-1)*n*p+1:j*n*p]'*Sm, indK[:, (j-1)*n*p+1:j*n*p], o[j], 1)
but this is still too slow for me. This is part of another loop, indK does not change, but Sm and o do.
Any ideas on how to improve it or is this the best I can get?
Row-wise access is bad for CSC matrices. A lot of the numbers in the columns will likely be zero.
Note also that
.2 is not sparse. Generally for sparse methods to be effective you need sparsity of
.01 or so or lower. Also, using
@views will make this much faster.
20% sparsity is probably too large to be worthwhile using sparse matrices. Does your real problem have more sparsity?
The matrix in my original problem really is sparse, but its construction is too long to put it here.
Thanks for the note on sparsity, the matrix in my original problem is sparse, I will change the density.
@views pom += o[j]*indK[:, (j-1)*n*p+1:j*n*p]'*Sm*indK[:, (j-1)*n*p+1:j*n*p]
is way slower than the version with
mul!, even when I convert the sparse matrix into dense.
Yeah, I was thinking about creating an array with matrices instead of one huge matrix, but it does not bring any improvement.
In that case, you lose the sparse data structure when you multiply by
Sm, unless you use
Sm = sparse(rand(m,m))… And
pom is also not a sparse matrix.
Omg, you are right, how could I miss that. Together with taking
indK as an array of sparse matrices, this solves my problem. Thanks a lot.