Parallel Matrix Multiplication

I am looking for a short tutorial/example explaining how we do matrix multiplication in parallel. Note that, the size of matrix is currently 200 * 200 (40000 elements). Thanks.

Answer 1

[James Tursa]

The best way to do matrix multiply in MATLAB is to use the * operator, as you normally would. This will call highly optimized BLAS routines that have parallel algorithms in the background as appropriate. If you have sparse matrices, the * operator will call specialized sparse matrix multiply routines. You will be very hard pressed to do better on your own.

[Walter Roberson]

An implication to James' response is that if you have a lightly populated array, you should be sure to use sparse() so that it can do the optimizations when it multiplies.

[Mohammad]

James, Kindly attach a .m file that demonstrates this.

[James Tursa]

What is it exactly that you want demonstrated? That * is multi-threaded? How to do sparse matrix operations? Or what?

E.g., here is a matrix multiply on a quad core system:

>> a = rand(5000);
>> b = rand(5000);
>> tic;a*b;toc
Elapsed time is 2.307906 seconds.

Here is the same matrix multiply on the same quad core system when starting MATLAB with the -singleCompThread option:

>> a = rand(5000);
>> b = rand(5000);
>> tic;a*b;toc
Elapsed time is 9.796487 seconds.

So the matrix multiply operation * for full matrices is obviously multi-threaded in the background.

And here is a demonstration of taking advantage of sparse matrix multiply:

>> a(a>.01) = 0;
>> sa = sparse(a);
>> b(b>.01) = 0;
>> sb = sparse(b);
>> tic;a*b;toc
Elapsed time is 2.252949 seconds.
>> tic;sa*sb;toc
Elapsed time is 0.433186 seconds.

[Joel Lynch]

"If you have sparse matrices, the * operator will call specialized sparse matrix multiply routines. You will be very hard pressed to do better on your own."

Unfortunately, this statement is highly misleading, if not incorrect. At least as of R2023, sparse matrix multiplication on CPU's are limited to a single thread. James's tests don't show this, but it's easy enough to see:

N = 20000;
density = 0.2;
A = sprand(N,N, density);
b = rand(N, 1);
Nmax_threads = maxNumCompThreads('automatic')

Nmax_threads = 2

timeit(@() A*b)

ans = 0.0510

maxNumCompThreads(1);
timeit(@() A*b)

ans = 0.0506

It's incredibly suprising, as MKL has supported multithreaded sparse matrix math for quite a while. Until this gets updated, the best best for a large sparse matrix problem is to use gpuArray, even a commercial card can do significantly better.