~ $ cd research/loops-into-zones && cat README.md

Loops into zones

A graphics card is not many processors. It is two kinds of parallelism, one inside the other. Outside, blocks: many of them, loosely coupled, which should hardly talk to each other. Inside each block, threads: tightly coupled, in step, sharing a small memory that is private to the block and very fast. A program is quick on a card when its work has that same shape, and slow, whatever else is done to it, when it has not.

I had met that shape before. The NAS multi-zone benchmarks, which had scaled so well on Cyclops and over distributed memory, are built the same way: coarse zones outside that barely communicate, fine loops inside each. So the method of the second half of the thesis is one sentence: take the algorithm you have and make it multi-zone. The rest is how to do that to an algorithm nobody wrote that way, without changing what it computes.

The algorithm

Graduated assignment matches two graphs by refining a matrix of probabilities, P[a,i]: how likely it is that vertex a of one graph is vertex i of the other. Each round computes, for every pair,

Q[a,i] = sum over b, j of  P[b,j] · C[a,i,b,j]
P[a,i] = exp(β · Q[a,i])

where C says how compatible matching a with i is with matching b with j. Written the way anyone would write it first, it is two loop nests:

for a in 1 .. R:
  for i in 1 .. R:
    Q[a,i] = 0
    for b in 1 .. R:
      for j in 1 .. R:
        Q[a,i] += P[b,j] · C[a,i,b,j]

for a in 1 .. R:
  for i in 1 .. R:
    P[a,i] = exp(β · Q[a,i])

Four loops deep over the vertices, R⁴ products a round. It is correct, and it has no shape at all: every iteration may touch any part of P and of C, so on a card every thread wants both whole, and nothing fits in a block's memory. The two steps that follow are done to the short nest first, where they are easy to see, and then to the long one.

Step one: tile the loops

Replace each index by a tile and a position inside it: a = c·B + d, i = k·B + l. The loop over a becomes two, one over the tiles c and one over the positions d, and the loop over i likewise:

for a in 1 .. R:          before

for c in 0 .. R/B:        after: which tile,
  for d in 1 .. B:        and where in it
    a = c·B + d

Done to both loops of the short nest:

for c in 0 .. R/B:
  for d in 1 .. B:
    for k in 0 .. R/B:
      for l in 1 .. B:
        a = c·B + d
        i = k·B + l
        P[a,i] = exp(β · Q[a,i])

Nothing has changed: the same assignments happen, in the same order. But the matrices are now visibly cut into sub-matrices of B × B, and a sub-matrix is something that fits in a block. These are the zones.

Step two: reorder the loops

Now move the loops until the nest has the card's shape: the tile loops outside, the position loops inside.

for c in 0 .. R/B:            ← blocks
  for k in 0 .. R/B:          ← blocks
    for d in 1 .. B:          ← threads
      for l in 1 .. B:        ← threads
        P[a,i] = exp(β · Q[a,i])

The long nest gets the same two steps. Its inner indices are tiled too, b = e·B + f and j = u·B + v, and the tile loops e, u go outside the position loops f, v:

for c, k in tiles:                      ← blocks
  for d, l in positions:                ← threads
    Q[a,i] = 0
    for e, u in tiles:                  ← one sub-matrix of P at a time
      for f, v in positions:
        b = e·B + f
        j = u·B + v
        Q[a,i] += P[b,j] · C[a,i,b,j]

Now the innermost two loops walk one B × B sub-matrix of P from corner to corner before moving to the next, which is exactly what a block's small memory can hold. One thing is still in the way: C has four dimensions and no sub-matrix of it is small. So C is replaced by what it is made of — the adjacency of each graph and the compatibility of their attributes — four small factors that can each be fetched a sub-matrix at a time.

Reordering is where a programmer knows what a compiler does not: that these iterations do not depend on each other. Compilers tile and reorder loops by themselves when they can prove it is safe, and here they cannot.

Step three: say it in the margin

The loops now have the right shape; what is left is to say which is which. Two annotations, in the manner of OpenMP:

#pragma hy parallel for [into(threads)] [reduction(OP:r)]
#pragma hy parallel fetch(m : sizes : origin : indexes [: permutation])

parallel for alone spreads a loop across blocks; with into(threads), across the threads of a block. parallel fetch copies a sub-matrix of m into the block's private memory and redirects every access inside the statement to the copy, in whatever order of dimensions keeps the accesses contiguous; if the statement writes to it, it is flushed back. It comes from the peek of the stream compiler of the first half.

On a processor, into(threads) and fetch are ignored and the rest is OpenMP. The same source is the CPU version and the GPU version, which is what the whole thesis had been after: a serial program, still readable, with its parallelism written beside it.

One algorithm, two parallelisations

The method does not give one answer, and that turned out to be a result.

Large graphs

The graph does not fit in a block. Tile it, as above, and let blocks work on different sub-matrices. It scales with the size of the graph — up to 1,024 vertices, where the serial version was estimated at forty-two days and not run — and is poor on small graphs, which do not have enough tiles to fill a card.

Many small graphs

A whole graph fits in a block. So a block is a graph: the outer level is many matchings at once that never talk to each other, the inner level is one matching. That is a multi-zone program exactly, and it is the one measured against OpenMP and one thread.

The common labelling of many graphs is a different algorithm, and the method carried over: reordering to find loops that several parts of it shared and merge them into tightly coupled kernels, and fission to split what did not belong together. Of the tile sizes tried, B = 8 was the best for CUDA.

What it does not change

The result. Tiling and reordering move the same operations around; they do not approximate, drop or relax anything, and the thesis reports the parallel versions giving the same values as the serial one. That matters more than the speed: a faster algorithm that answers slightly differently is a new algorithm, and has to be validated again by the people who trusted the old one. This one does not.