Locality benchmark

/*

 * Locality benchmark: row-major vs column-major traversal of a matrix.

 *

 * Companion to "When CPUs Stopped Scaling: Why Hardware Got Complicated"

 * https://hed.am/briefs/when-cpus-stopped-scaling/

 *

 *   cc -O2 -Wall -Wextra -o bench locality-benchmark.c && ./bench

 *

 * sum_array_rows and sum_array_cols are exactly the two functions from the

 * brief. The rest of the file exists to time them accurately, and to answer

 * a question the brief's two-version comparison does not: whether the

 * column walk is slow because it wastes 31 of every 32 bytes it fetches, or

 * because it fetches roughly the same bytes as the row walk but badly. A

 * third function, sum_array_cols_scrambled, answers it. Several of the

 * file's choices are deliberate, and each has a silent failure mode: get it

 * wrong and the program still runs and still prints a plausible number.

 *

 * Why 645 MiB. The array has to be larger than the largest cache, or the

 * column walk is served from cache and there is no difference left to

 * measure. It also has to be small enough that the operating system never

 * pages or compresses it, since that would time the memory manager rather

 * than the cache. On an M3 the difference stops growing at roughly 80 MiB and

 * is flat above that, so 645 MiB sits well clear of the cache without

 * approaching the size at which paging begins.

 *

 * Why 13000 columns and not 16384. A cache selects where to store an address

 * from bits in the middle of that address, so addresses separated by a large

 * power of two map to the same location and displace one another. A row here

 * is 13000 * 4 = 52,000 bytes, which is not a power of two, so the addresses

 * a column walk visits are spread across the cache. At 16384 columns the rows

 * would be 65,536 bytes apart and every element of a column would map to the

 * same few locations. That would measure cache conflict rather than locality.

 *

 * Why the array is written before timing begins. malloc returns pages that

 * have not yet been given physical memory, and what a read of one does next

 * is platform-specific. Linux points every untouched page of a private

 * anonymous mapping at one shared page of zeroes, so the array would occupy a

 * single page, that page would stay in cache, and both loops would report

 * almost the same time. macOS commits a real page per fault instead, which

 * does not collapse the array but does charge the first pass a fault for

 * every page in it. Writing each element first removes both effects.

 *

 * Best of 24 rather than a mean: processors raise their clock frequency under

 * sustained load, so the earliest runs are slower than the steady-state rate.

 *

 * The ns/elem column is six wide because a slower machine than the one in the

 * published run reports a two-digit rate, which a narrower field pushes out of

 * its column and out of alignment with the rows above it.

 *

 * Why the scrambled pass, and why stride 997. A 128-byte line holds 32

 * consecutive columns of one row, so sum_array_cols may not pay for a fresh

 * fetch on every element -- it may pay once per line and reuse the rest,

 * since one such sweep touches about 1.6 MiB and the L2 is 16 MiB. The

 * scrambled pass visits the same N columns in the order j, j+997, j+2*997

 * mod N, ..., which touches the same bytes but decorrelates consecutive

 * visits from cache-line adjacency. 997 is prime and does not divide 13000,

 * which is what guarantees every column is visited exactly once.

 */

#define _POSIX_C_SOURCE 199309L

#include <stdio.h>

#include <stdlib.h>

#include <time.h>

#define M 13000

#define N 13000

#define REPS 24

#define STRIDE 997

// Version A: This is a great way to sum a matrix.

int sum_array_rows(int a[M][N])

{

    int i, j, sum = 0;

    for (i = 0; i < M; i++)

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

            sum += a[i][j];

    return sum;

}

// Version B: This is a terrible way to sum a matrix.

int sum_array_cols(int a[M][N])

{

    int i, j, sum = 0;

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

        for (i = 0; i < M; i++)

            sum += a[i][j];

    return sum;

}

// Version C: Same columns as B, visited in a different order.

int sum_array_cols_scrambled(int a[M][N])

{

    int i, sum = 0, j = 0;

    for (int k = 0; k < N; k++) {

        for (i = 0; i < M; i++)

            sum += a[i][j];

        j = (j + STRIDE) % N;

    }

    return sum;

}

/* Assigned but never read, and volatile for that reason: without it the

 * compiler can prove the sums are unused and remove every loop. */

static volatile int sink;

static double now(void)

{

    struct timespec ts;

    clock_gettime(CLOCK_MONOTONIC, &ts);

    return (double)ts.tv_sec + (double)ts.tv_nsec / 1e9;

}

int main(void)

{

    /* 13000 * 13000 fits in an int. The cast matters only if the matrix is

     * enlarged: beyond 46,340 columns the product exceeds INT_MAX, and

     * computing it as an int would allocate less memory than the loops read. */

    const size_t elems = (size_t)M * N;

    int(*a)[N] = malloc(elems * sizeof(int));

    if (!a) {

        perror("malloc");

        return 1;

    }

    for (size_t k = 0; k < elems; k++)

        ((int *)a)[k] = (int)(k & 1);

    double rows = 1e30, cols = 1e30, scrambled = 1e30;

    for (int r = 0; r < REPS; r++) {

        double t = now();

        sink = sum_array_rows(a);

        double d = now() - t;

        if (d < rows)

            rows = d;

        t = now();

        sink = sum_array_cols(a);

        d = now() - t;

        if (d < cols)

            cols = d;

        t = now();

        sink = sum_array_cols_scrambled(a);

        d = now() - t;

        if (d < scrambled)

            scrambled = d;

    }

    printf("%d x %d ints, %.0f MiB, best of %d\n\n"

           "  row-major               %7.3f s   %6.3f ns/elem\n"

           "  column-major            %7.3f s   %6.3f ns/elem\n"

           "  column-major, scrambled %7.3f s   %6.3f ns/elem\n\n"

           "  column-major / row-major       %6.2fx slower\n"

           "  scrambled / column-major       %6.2fx slower\n",

           M, N, (double)(elems * sizeof(int)) / 1048576.0, REPS, rows,

           rows * 1e9 / (double)elems, cols, cols * 1e9 / (double)elems,

           scrambled, scrambled * 1e9 / (double)elems, cols / rows,

           scrambled / cols);

    free(a);

    return 0;

}