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

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?

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.