I am performing several (thousands) matrix multiplications of an NxN sparse (~1-2%) matrix, let's call it B, with an NxM dense matrix, let's call it A (where M<N). N is large, as is M; on the order of several thousands.
Now, usually, matrix multiplications and most other matrix operations are implicitly parallelized in Matlab, i.e. they make use of multiple threads automatically. This appears NOT to be the case if either of the matrices are sparse (see e.g. this StackOverflow discussion and this largely unanswered MathWorks thread). This is a rather unhappy surprise to me. We can verify this by the following code:
clc; clear all;
N = 5000;
M = 3000;
A = randn(N,M);
B = sprand(N,N,0.015);
Bf = full(B);
m(i) = 2^(i-1);
y = B*A;
walltime(i) = toc;
speedup(i) = walltime(1)/walltime(i);
This produces the following output, which illustrates that there is no difference between using 1, 2, and 4 threads:
If, on the other hand, I replace B by its dense form, refered to as Bf above, I get significant speedup:
So, my question: is there any way at all to access a parallelized/threaded version of matrix operations for sparse matrices without converting them to dense form? Alternatively, is this something one might expect would be implemented in future versions of Matlab? I found some suggestion involving .mex files here, but it seems the links are dead and not very well documented/no feedback?
It seems to be a rather severe restriction of implicit parallelism functionality, since sparse matrices are abound in computationally heavy problems, and hyperthreaded functionality highly desirable in these cases.
I appreciate anyone's consideration of the this issue at any level - thanks a bunch in advance! -Thomas