A lot of big sparse matrices

General Usage Performance
1

Hi everyone,

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)
for j=1:n*p
     pom+= indK[:, (j-1)*n*p+1:j*n*p]'*Sm*indK[:, (j-1)*n*p+1:j*n*p]*o[j]
end

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?

2

Row-wise access is bad for CSC matrices. A lot of the numbers in the columns will likely be zero.

3

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.

2 Likes
4

20% sparsity is probably too large to be worthwhile using sparse matrices. Does your real problem have more sparsity?

1 Like
5

The matrix in my original problem really is sparse, but its construction is too long to put it here.

6

Thanks for the note on sparsity, the matrix in my original problem is sparse, I will change the density.
Using @views as

@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.

7

Yeah, I was thinking about creating an array with matrices instead of one huge matrix, but it does not bring any improvement.

8

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.

3 Likes
9

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.