Traversal order and loop tiling

Two nested loops that visit the same cells of a matrix can differ by a factor of seven in cache misses just because of which loop is inside. This lecture compares reading row by row and column by column, then cuts a transpose into tiles so that a block stays in the cache while it is used, and measures both on a matrix made from SPY closes.

A matrix with n rows and n columns is stored in row-major order. A loop that moves along a row, changing j, reads neighbouring values. A loop that moves down a column, changing i, jumps a whole row each step. The two loops below visit the same cells, one with j on the inside and one with i. We print the first six addresses of a 4 by 4 matrix of float64 values.

def row_major_trace(n):    return [(i * n + j) * 8 for i in range(n) for j in range(n)]def col_major_trace(n):    return [(i * n + j) * 8 for j in range(n) for i in range(n)]1print("row:", row_major_trace(4)[:6])print("col:", col_major_trace(4)[:6])

row: [0, 8, 16, 24, 32, 40]
col: [0, 32, 64, 96, 8, 40]

1. 1With i inside, the next read is a whole row (here 4 values = 32 bytes) further.

Row by row the addresses go up by 8 bytes, the stride of the first topic, and a line serves eight reads. Column by column they jump by 32 bytes here, and by the length of a whole row in a bigger matrix. We now count the misses of the two orders, in a small 4 KB cache, for matrices from 16 by 16 up to 73 by 73.

for n in (16, 32, 48, 73):    print(n, "x", n, "row:", count_misses(row_major_trace(n)), "col:", count_misses(col_major_trace(n)))

16 x 16 row: 32 col: 32
32 x 32 row: 128 col: 1024
48 x 48 row: 288 col: 2304
73 x 73 row: 667 col: 4573

A 16 by 16 matrix is 2 KB and fits in the cache, so after loading it both orders cost 32 misses, the same. From 32 by 32 the matrix no longer fits, and the column order misses on every single read: 1,024 misses against 128, a factor of eight, which is the eight values of a line. At 73 by 73 the factor is nearly seven. The shape of the loops did not change; the size of the data did.

Sometimes no loop order helps. The transpose B[j][i] = A[i][j] reads A and writes B, and if the loops go along the rows of A they go down the columns of B. Whichever loop is inside, one of the two matrices is walked with a big stride. The cure is loop tiling: process the matrix in small square tiles, so that the few lines a tile touches are used many times before they are evicted.

1. 1
   Choose a tile size tSmall enough that one tile of A and one of B fit in the cache.
2. 2
   Loop over the tilesTwo outer loops with step t: ii and jj give the corner of a tile.
3. 3
   Work inside the tileThe two inner loops visit only the t by t cells of this tile, as before.
4. 4
   Move to the next tileOnly now does the program leave the lines of the old tile behind.

The code below only lists the order in which the cells are visited, for a 4 by 4 matrix, without and with tiles of 2. With tiles, the program finishes the top-left 2 by 2 block before it touches any other cell.

def visit_order(n, tile=None):    step = n if tile is None else tile    order = []    for ii in range(0, n, step):1        for jj in range(0, n, step):            for i in range(ii, min(ii + step, n)):                for j in range(jj, min(jj + step, n)):                    order.append((i, j))    return orderprint("plain:", visit_order(4)[:8])print("tile 2:", visit_order(4, 2)[:8])

plain: [(0, 0), (0, 1), (0, 2), (0, 3), (1, 0), (1, 1), (1, 2), (1, 3)]
tile 2: [(0, 0), (0, 1), (1, 0), (1, 1), (0, 2), (0, 3), (1, 2), (1, 3)]

1. 1The two outer loops walk over the tiles; the two inner loops walk inside one tile.

Now the real transposition. We write the trace of the reads of A and the writes of B, with B stored right after A, and count the misses for a 73 by 73 matrix with no tile and with tiles of 4, 8, 16 and 32.

def transpose_trace(n, tile=None):    base_b = n * n * 8    step = n if tile is None else tile    trace = []    for ii in range(0, n, step):        for jj in range(0, n, step):            for i in range(ii, min(ii + step, n)):                for j in range(jj, min(jj + step, n)):                    trace.append((i * n + j) * 8)                    trace.append(base_b + (j * n + i) * 8)    return tracefor tile in (None, 4, 8, 16, 32):    print("tile", tile, "misses:", count_misses(transpose_trace(73, tile)))

tile None misses: 5996
tile 4 misses: 2603
tile 8 misses: 2036
tile 16 misses: 1960
tile 32 misses: 1714

Without tiling the transpose misses 5,996 times. With tiles of 4 it is already down to 2,603, and with tiles of 8 to 32 it stays between 2,036 and 1,714, about three times fewer. A tile that is too small wastes the edges of the lines it loads, so the gain grows with the tile until the tile stops fitting.

The lab follows the lecture: the two read orders, the transpose trace with an optional tile, the count of misses by order and by tile, the choice of the best tile and the study on the SPY matrix.

• Function · questionBuilds on
• row_major_trace(n)Which addresses does a row-by-row read touch?Plain arithmetic
• col_major_trace(n)And a column-by-column read?Plain arithmetic
• transpose_trace(n, tile=None)Which addresses does a transpose touch, tiled or not?Four nested loops
• compare_orders(n)How much worse is the column order?row_major_tracecol_major_trace
• misses_by_tile(n, tiles)How do tiles change the misses of the transpose?transpose_trace
• best_tile(n, tiles)Which tile is best?misses_by_tile
• study(bars)What does tiling gain on a matrix of SPY closes?compare_ordersbest_tile

Next topic: tiling hides the delay of a miss by reusing lines; prefetching hides it by asking for a line before it is needed. In Prefetching you add a stride detector to the cache and see which access patterns it can predict, and which, like a random walk through a table, it cannot.