Loop Transforms to Increase Parallelism

I have the need to parallelize as much as I can (in order to make it run fast) the following part of an algorithm:

do i=1, n
do j=1, n
A[i, j] = A[i-4, j-4] + A[i-3, j-3] + B[i, j]
end do
end do


I was thinking that the best i could do here is a loop skewing (but with what factor?) and loop interchange that would allow me to finally have a DOALL on the inner loop. Any feedack and help regarding those thoughts?

• The diagonals are indepedendent, so you can compute those in parallel. Note that this will be memory-bound, so you can only expect a speedup of you store A and B by diagonals from the start. This will give you a speedup if you stay sequential, too, since it uses the cache more effectively. – Raphael Feb 13 '18 at 12:29
• @Raphael Why not increase parallelism even more by computing (within each diagonal as you mentioned) every two sequentially and those tuples on parallel ? (something like loop unrolling) – Jason Feb 13 '18 at 14:04
• Try it out and benchmark it. My feeling is that unless it reduces the number of cache misses, it won't do much for you. Also, keep cache syncing in mind. – Raphael Feb 13 '18 at 16:07

This is an incomplete question because you have no base case; A[-3, -3] is not defined, for example. But let's assume we're looking at a window into a larger array, so those values have some consistent, distinct and independent value.

Denote those predetermined values with #, and look at the top-left corner. The values we can fill in are those three topmost rows and three leftmost columns, denoted by ·.

# # # # # # # # # # # # # #
# # # # # # # # # # # # # #
# # # # # # # # # # # # # #
# # # # # # # # # # # # # #
# # # # · · · · · · · · · ·
# # # # · · · · · · · · · ·
# # # # · · · · · · · · · ·
# # # # · · ·
# # # # · · ·
# # # # · · ·
# # # # · · ·
# # # # · · ·
# # # # · · ·
# # # # · · ·


Doing those in parallel means that you need to take n / 3 sequential steps. Since the matrix is square, exclusively going down balances the costs better, but still gives you n / 3 sequential steps with parallelism of 3n.

This is remarkably excessive for a typical CPU for significantly sized n; rather you should be streaming memory into SIMD operations as fast as possible. You can do a straightforward N-wide SIMD loop over the x-axis.

A[i:i+N, j] = A[i-4:i-4+N, j-4] + A[i-3:i-3+N, j-3] + B[i:i+N, j]


This probably won't max out your bandwidth entirely, but it will get pretty darn close and it's unlikely you can get more from a multicore algorithm, though you are free to try partitioning the work along diagonals.