Foundations —
Where the Data Lives

Concept map

One idea, three pictures. The ladder is how far data travels, the line is how it travels, the roofline is whether the work per byte pays for the trip. Click a node to jump to it.

This is not an MPI lecture. It is the ground floor for OpenMP: from Mon Oct 19 the data still moves, silently. Weeks 9–13 are discovering what OpenMP moves for you.

The ladder: distance is time

Grab it from your hand, from the desk, from the shelf, from the basement, from the neighbour's basement, by mail. Each rung several times slower; top to bottom, about 100×.

The ladder, measured: chase.cpp

Same program. Only the fit changes. A random cyclic chain of pointers, one per 64-byte line, so every load depends on the previous one and the prefetcher cannot guess.

chase.cpp · average ns per load, min of 3 runs, g++ -O2, test VM
lscpu: 64 KiB L1d and 2 MiB L2 over 2 instances, 33 MiB L3
   4 –  32 KB    1.52 – 1.55 ns
  48 – 256 KB    4.28 – 4.66 ns
        1 MB    15.1 ns
   2 –   4 MB   25.4 – 26.5 ns
        8 MB   120 ns
  16 – 256 MB   157 – 189 ns

Touch every double. Then touch every 8th.

traverse.cpp, part 2. A 256 MB array of doubles. Loop 1: x[i] *= 1.000001 for every i. Loop 2: the same, for every 8th i. One-eighth of the work.

const size_t M = 32*1024*1024;                      // 32 Mi doubles = 256 MiB

for (int s : {1, 2, 4, 8, 16, 32, 64}) {
    for (size_t i = 0; i < M; i += s) {
        B[i] *= 1.000001;
    }
}

min of 5 runs, test VM
stride  1:  0.0371 s   touches 33,554,432 doubles
stride  2:  0.0306 s   touches 16,777,216
stride  4:  0.0285 s   touches  8,388,608
stride  8:  0.0274 s   touches  4,194,304
stride 16:  0.0245 s   touches  2,097,152
stride 32:  0.0142 s   touches  1,048,576
stride 64:  0.0060 s   touches    524,288

Sum an N × N matrix by rows. Then by columns.

traverse.cpp, part 1. Row-major doubles, A[i·N + j]. Same adds, same answer, same N² loads. Rows: next element 8 bytes away. Columns: next element 8·N bytes away.

for (int t = 0; t < 5; t++) {                       // min of 5
    s = 0;
    for (i) for (j) s += A[(size_t)i*N + j];       // row:    stride 1

    s = 0;
    for (j) for (i) s += A[(size_t)i*N + j];       // column: stride N doubles
}

$ ./traverse N           (min of 5, each size)
N =  256  (0.5 MB)  row 0.0001 s  column 0.0001 s   1.0x
N =  512  (  2 MB)  row 0.0003 s  column 0.0008 s   2.6x
N = 1024  (  8 MB)  row 0.0014 s  column 0.0112 s   7.9x
N = 2048  ( 34 MB)  row 0.0085 s  column 0.0656 s   7.7x
N = 4096  (134 MB)  row 0.0246 s  column 0.2536 s  10.3x

The 64-byte line is the unit

One fact, three consequences. 64 bytes on x86 and on Cardinal; Apple M-series reports 128.

Blocking's reuse is what the roofline will call arithmetic intensity, two screens from here.

NUMA: a Cardinal node is two sockets

First touch decides the socket. The OS places a page on the socket of the thread that first writes it, not the one that allocates it. Other socket: roughly 1.5–2× a local access. Correctness tests never see this.

The same picture in code

double* b = (double*)malloc(n * 8);  // address only, no pages

for (long i = 0; i < n; i++) {       // serial init:
    b[i] = 0.0;                      // left picture
}

// "omp parallel for" (Mon Oct 19): split the loop
// across threads, each takes a contiguous part
#pragma omp parallel for schedule(static)
for (long i = 0; i < n; i++) {
    b[i] = 0.0;                      // right picture
}

#pragma omp parallel for schedule(static)
for (long i = 0; i < n; i++) {       // same split:
    b[i] = f(b[i]);                  // every read local
}

The OpenMP lines are the first pragma you have seen; Monday Oct 19 explains them. The comment is the whole gloss for today.

The roofline: two speed limits. The slower one wins.

Flat roof: how fast the cores can compute. Slope: how fast memory can feed them. Across: FLOPs done per byte moved (arithmetic intensity). Which limit is yours?

AXPY, and reuse moves you right

The simplest loop in scientific computing: y[i] = a·x[i] + y[i]. Three numbers: FLOPs, bytes, fraction of peak.

time check, n = 10^8, model
memory:      2.4 GB / 0.5 TB/s        = 4.8 ms
arithmetic:  2 × 10^8 FLOP / 6.1 TF   = 0.033 ms
memory ≈ 145 × the arithmetic

Run it on Cardinal

The program behind the stride and row/column screens is traverse.cpp, the instructor's program; it lives in the deck's _deck_sources, not on a student branch. The listing is on this page. Run it on Cardinal for your own numbers, and say where it ran.

1. The program

   This is the part the deck shows. Around it: std::chrono::steady_clock, min of 5 runs per case, and the sums printed so the compiler cannot delete the loops. chase.cpp (the pointer chase) sits next to it in _deck_sources and is not listed here.

   // traverse.cpp, part 1: N x N row-major doubles in A
   for (int t = 0; t < 5; t++) {                       // min of 5
       s = 0;
       for (i) for (j) s += A[(size_t)i*N + j];       // row:    stride 1
   
       s = 0;
       for (j) for (i) s += A[(size_t)i*N + j];       // column: stride N doubles
   }
   
   // traverse.cpp, part 2: 256 MiB of doubles in B
   const size_t M = 32*1024*1024;                      // 32 Mi doubles = 256 MiB
   
   for (int s : {1, 2, 4, 8, 16, 32, 64}) {
       for (size_t i = 0; i < M; i += s) {
           B[i] *= 1.000001;
       }
   }

2. Log in, load the compiler, build, check the line size

   Modules do not persist between sessions; load it every time.

   ssh <your-username>@cardinal.osc.edu
   module load gcc/12.3.0
   g++ -O2 -o traverse traverse.cpp
   getconf LEVEL1_DCACHE_LINESIZE

   copy

3. Get cores of your own before you time anything

   The deck ran it on a login node, which is fine for a first look; a login node is shared, so for numbers you would write down use an allocation.

   salloc -N 1 -n 1 -c 8 -t 00:15:00 -A PCS0313
   module load gcc/12.3.0
   srun -c 8 ./traverse 4096

   copy

Experiments

getconf LEVEL1_DCACHE_LINESIZE
lscpu | grep -i cache

srun -c 8 ./traverse 4096
srun -c 8 ./traverse 256

for N in 256 512 1024 2048 4096; do
    srun -c 8 ./traverse $N
done

srun -c 8 ./traverse 4096        # part 2 runs after part 1

numactl -H