Blocking refers to a class of techniques that aim to break down the original problem into small chunks, so that all data needed to process each chunk fits in the cache. This way that data can be optimally reused before the next chunk is processed.

Blocking can be performed to optimize reuse of adjacent data, spatial blocking, or to optimize for reuse of the same data, temporal blocking, or both.

An example of a problem where spatial blocking can be used is creating a transposed copy of a matrix:

void transpose(double dst[SIZE][SIZE], double src[SIZE][SIZE]) {
  int i, j;
  for (i = 0; i < SIZE; i++)
    for (j = 0; j < SIZE; j++)
      dst[j][i] = src[i][j];
}
        

We see that the src matrix is read one row at a time and that the dst matrix is written one column at a time. If the matrix is large enough there won't be enough cache lines for all the rows accessed in the dst matrix in one iteration of the outer loop, and we get a very inefficient accesses pattern as described in Section 5.2.1, “Inefficient Loop Nesting”.

We could change the nesting of the loops as suggested there:

void transpose(double dst[SIZE][SIZE], double src[SIZE][SIZE]) {
  int i, j;
  for (j = 0; j < SIZE; j++)
    for (i = 0; i < SIZE; i++)
      dst[j][i] = src[i][j];
}
        

However, while this fixed the problem with the dst matrix the src is instead read one column at a time, so we have only moved the problem.

What we can do is to only copy a few columns from the src matrix at a time, so that we know there are enough cache lines to hold the corresponding rows in the dst matrix. Each such set of columns becomes a block.

void transpose(double dst[SIZE][SIZE], double src[SIZE][SIZE]) {
  int jb, i, j;
  for (jb = 0; jb < SIZE; jb += BLOCK)
    for (i = 0; i < SIZE; i++)
      for (j = jb; j < jb + BLOCK && j < SIZE; j++)
        dst[j][i] = src[i][j];
}
        

The variable jb now keeps track of which columns are being copied from the source matrix in each block. In the first iteration of the outermost loop the columns 0 to BLOCK-1 are copied, in the next iteration the columns BLOCK to 2*BLOCK-1 are copied, and so on.

By choosing a small enough value for the block size BLOCK we can make sure that the corresponding rows written to the dst matrix in each iteration of the middle loop fit in the cache.

Temporal blocking is sometimes possible for iterated algorithms. Data produced during one iteration is used in the next one. By performing several iterations for a sub-problem, the previous iteration values will remain in the cache when performing the next iteration.

The cookbook way of performing blocking is to focus on reuse of a target structure, and work outwards through loop levels until reaching a loop level that causes all elements of the target structure to be accessed for each iteration. Inner loops to this outer loop level are candidates to be broken up. Consider the following trivial example:

for (int iter = 0; iter < ITER; iter++) // "outer" loop
  for (int j = 0; j < SIZE; j++)
    func(a[j]);
        

Each iteration of the iter loop causes each element of a to be accessed. Hence we declare the iter loop to be the outer loop. Loop levels interior to this loop, namely for j in this case, are candidates to be blocked. Upon blocking, one new loop level is placed outside the outer loop and the original loop is kept in its the original location, with the restriction that it works on only a subset of the data for each iteration of the new outermost loop:

for (int jj = 0; jj < SIZE; jj += BLOCK)
  for (int iter = 0; iter < ITER; iter++) // "outer" loop
    for (int j = jj; j < jj + BLOCK; j++)
      func(a[j]);
        

Blocking will cause calculations to occur in a different order than in the original code, and sometimes this is fine. Many numerical algorithms can be reformulated to maintain accumulators to hold partial results between processing of different blocks. Some algorithms, however, will not lend themselves to blocking that easily.

One famous example is blocking a matrix multiplication. The unoptimized code looks like:

for (int i = 0; i < SIZE; i++)
  for (int j = 0; j < SIZE; j++)
    c[i][j] = 0;

for (int i = 0; i < SIZE; i++) // "outer" loop
  for (int j = 0; j < SIZE; j++)
    for (int k = 0; k < SIZE; k++)
      c[i][j] += a[i][k] * b[k][j];
        

The array b is traversed in the wrong way, and we should try to block with respect to it. Note that the for i loop level is the outer loop which causes all b elements to be revisited. for j and for k are therefore candidate loops to block. This can be done in different ways, for instance only splitting the k loop:

for (int i = 0; i < SIZE; i++)
  for (int j = 0; j < SIZE; j++)
    c[i][j] = 0;

for (int kk = 0; kk < SIZE; kk += BLOCK)
  for (int i = 0; i < SIZE; i++) // "outer" loop
    for (int j = 0; j < SIZE; j++)
      for (int k = kk; k < kk + BLOCK && k < SIZE; k++)
        c[i][j] += a[i][k] * b[k][j];
        

Or you could even split all the loops, possibly with different blocking factors:

for (int i = 0; i < SIZE; i++)
  for (int j = 0; j < SIZE; j++)
    c[i][j] = 0;

for (int ii = 0; ii < SIZE; ii += BLOCK_I)
  for (int kk = 0; kk < SIZE; kk += BLOCK_K)
    for (int jj = 0; jj < SIZE; jj += BLOCK_J)
      for (int i = ii; i < ii + BLOCK_I && i < SIZE; i++)
        for (int k = kk; k < kk + BLOCK_K && k < SIZE; k++)
          for (int j = jj; j < jj + BLOCK_J && j < SIZE; j++)
            c[i][j] += a[i][k] * b[k][j];
        

When blocking, special care needs to be taken to look for data dependencies. If one iteration depends on values produced in an earlier iteration, it follows that one must give special treatment to the boundary zones of the blocks to preserve the application semantics.

The blocking factor should be chosen to make the active problem size fit inside the target cache size. Different loop levels can be blocked to fit different cache levels.

When two loops use the same data it is beneficial to reduce the amount of data accessed between the loops. Cache lines fetched into the cache by the first loop may then still be in the cache when the second loop is run, so that it doesn't have to fetch them. The less data touched between the two loops, the better. Ultimately, if it is possible to completely merge the loop bodies, then all potential misses from the second loop will disappear.

It is not always possible to move the operations between loops. Special care must be taken not to move any data accesses past a write accesses to the same address, since that will change the meaning of the program.

Freja will identify such accesses that may preclude loop fusion. However, there are cases where it may fail, for example, when a dependence is carried in a processor register and never stored to memory, or when Freja's sparse sampling of information missed some dependency.

You therefore must always verify that the loops really are fusible, even if Freja suggests fusion.

Depending on the data dependencies at hand, there are several variations of how loop fusion can be achieved:

Note that since there are four different directions to move the instructions, the code motion barriers are evaluated independently for each potential motion direction.

Consider this example:

double vector[SIZE];
double vector2[SIZE];
int i;

for (i = 0; i < SIZE; i++)
  vector[i] = i;

f(vector2);

for (i = 0; i < SIZE; i++)
  vector[i] *= vector2[i];
        

These two loops both use all elements in vector. The first loop always misses in the cache, assuming that vector has not been accessed before.

If vector is smaller than the available cache, then it might still be in the cache when the second loop starts, but that depends on how much data is accessed in the function f. One thing to consider is whether the call to f can be moved out of the way. We can't tell from looking at this snippet whether it would be possible, since the function body is not revealed.

Assuming that f only reads the elements in vector2 it would be possible to move f before the first loop, and thus increase the likelihood that data remains in the cache between the loops.

double vector[SIZE];
double vector2[SIZE];
int i;

f(vector2);

for (i = 0; i < SIZE; i++)
  vector[i] = i;

for (i = 0; i < SIZE; i++)
  vector[i] *= vector2[i];
        

If vector is larger than the available cache the second loop will still suffer from cache misses when accessing vector, even though f has been moved away from between the loops. We should therefore try to merge the two loops.

double vector[SIZE];
double vector2[SIZE];
int i;

f(vector2);

for (i = 0; i < SIZE; i++) {
  vector[i] = i;
  vector[i] *= vector2[i];
}
        

Now Freja may indicate another fusion possibility between a potential loop inside f and the one remaining loop. In that case the programmer needs to decide whether it is worth duplicating the contents of f to merge these loops, sacrificing legibility for performance.

A variation of loop fusion occurs when there are multiple possible fusions targeting the same loop. The loop body may need to be fused multiple times in order to maintain proper program semantics. Consider this example:

for (i = 0; i < SIZE; i++) {
  vector[i] = i;
}

if (do_this) {
  for (i = 0; i < SIZE; i++) {
    vector[i]++;
  }
} else {
  for (i = 0; i < SIZE; i++) {
    vector[i]--;
  }
}
	

In this case, the do_this test is independent of the loops, and could be moved to before the whole section. But then the first loop needs to be duplicated in each of the then and else clauses. Then there are two loop fusion opportunities, which can be resolved like this:

if (do_this) {
  for (i = 0; i < SIZE; i++) {
    vector[i] = i; /* duplicated line */
    vector[i]++;
  }
} else {
  for (i = 0; i < SIZE; i++) {
    vector[i] = i; /* duplicated line */
    vector[i]--;
  }
}
	

The Expert Loop Fusion advice expresses this analysis as one step, by identifying all other candidate fusion targets. In principle, there can be multiple "first" loops and multiple "second" loops at the same time, and all of them are identified at the same time.