How much can matrix multiplication algorithm be parallelized?

Question

You all may know that simple Matrix Multiplication algorithm:

for(i = 0; i < n; i++)
   for(j = 0; j < n; j++)
      for(k = 0; k < n; k++)
         C[i,j] += A[i,k] * B[k,j]
      end for
   end for
end for 

When we want to parallelize it, we can, at first, parallelize the outer loop, which is the most easiest thing to understand, since it is clear that it does not have data dependencies with other loops.

Most tutorial on the web stops here, like we can not parallelize anymore.

But it seems for me that even the second outer loop (the one in the middle):

for(j = 0; j < n; j++)

can be parallelized.

What about the inner loop?

Is thus possible to parallelize all three loops of a matrix multiplication algorithm or must we stop with parallelizing only the first two loops?

Answer 1

Ypu're starting from a completely wrong point.

The execution time of matrix multiplication does not come from the number of multiplications and additions, it's the number of uncached memory accesses that kill you. Reading a number that's not in any processor cache takes about 100 times longer than a multiplication.

So the first step is rearranging the order of operations to perform as many operations as possible using only data that is present in the fastest processor cache. That's your first step before you even think about doing things in parallel.

The next step is adding multiple sums in parallel, still in one thread. Instead of summing up C[i,j] for example you add six sums C[i, j], C[i,j+1], C[i, j+2], C[i+1, j], C[i+1,j+1], C[i+1, j+2] in parallel. This means you are limited by the throughput of operations, not the latency.

The next step is using SIMD instructions. Your processor quite likely has instructions that perform 2, 4 or 8 floating-point operations just as fast as a single floating-point operation.

Once you have done all this, especially the first part where the work is divided up into "blocks", you can easily process multiple blocks on multiple threads. What you need to be careful about is that you have enough cache memory to do so.

Answer 2

You can very well parallelize both outer loops, since no two distinct threads will write to the same memory location.

However, when parallelizing the third loop, conflicts can happen as multiple threads want to update the same location. When you want to do this, you will need atomic operations (supported on all common hardware architectures, in C++ included as std::atomic).

Note, that addition is usually as fast as one can get, so parallelizing over the inner loop is usually not because of the speed-up, but for saving memory (the threads won't need all all of the row/column of A/B to calculate their portion anymore)

Answer 3

Here is a parallel algorithm running in $O((\log n))$. As @sonneXo pointed out, each iteration of the two outer-loops is independent of all other iterations assuming a concurrent read model. So the only problem is in summing up the elements in the inner loop. Here is a divide and conquer technique that solves the problem in $O(log(n))$ using $O(n)$ threads.

So given $i, j, s, t$, we want to compute $\sum\limits_{k=s}^t a_{ik} b_{kj}$ recursively. As input we give $i, j, 1, n$ and the algorithm outputs $c_{ij}$. You can then call the algorithm in parallel on each pair $i,j$ completing the matrix.

Now we describe the recursive algorithm. Given $i, j, s, t$. if $t-s < c$ for some constant $c$, say $c = 3$, then we compute the sum by hand in constant time. Else, let $m = s + \left\lfloor \frac{t-s}{2} \right\rfloor$. We build two instances of the problem in two new threads one of them solves $i, j, s, m$ and the other solves $i, j, m+1, t$. When the recursion returns, the algorithm only sums up two numbers returned by the recursion and returns the answer.

It is easy to see, why the running time is $O(\log n)$, since in each step we divide the size of the array by two and hence after $\left \lceil \log n \right \rceil$ steps we reach a constant size.

Note that this algorithm is not very efficient in practice due to the over head of creating threads. One algorithm that is efficient in practice, is to divide the array into $\sqrt n$ many pieces, sum up each of them in parallel in $O(\sqrt n)$ time and then sum up the answers again in $O(\sqrt n)$ time resulting in a running time of $O(\sqrt n)$ assuming you can affoard $\Omega(\sqrt n)$ threads running in parallel.

Note The number of used threads can be improved into $O(n / \log n)$ by dividing the column/row into chunks of length $\Theta(\log n)$ and each thread computes the answer sequentially in each chunk in time $O(\log n)$ with at most $O(n/\log n)$ threads and then we sum up the answer in the same divide and conquer scheme as above. Since we are summing at most $O(n / \log n)$ elements we need at most so many threads.