/*
 * 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. Three of its
 * 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.
 */

#define _POSIX_C_SOURCE 199309L

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define M 13000
#define N 13000
#define REPS 24

// 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;
}

/* Assigned but never read, and volatile for that reason: without it the
 * compiler can prove the sums are unused and remove both loops. */
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;

    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;
    }

    printf("%d x %d ints, %.0f MiB, best of %d\n\n"
           "  row-major     %7.3f s   %5.3f ns/elem\n"
           "  column-major  %7.3f s   %5.3f ns/elem\n"
           "  ratio        %7.2fx slower\n",
           M, N, (double)(elems * sizeof(int)) / 1048576.0, REPS, rows,
           rows * 1e9 / (double)elems, cols, cols * 1e9 / (double)elems,
           cols / rows);

    free(a);
    return 0;
}