<style>
.reveal {
font-size: 24px;
}
.reveal section img {
background:none;
border:none;
box-shadow:none;
}
pre.graphviz {
box-shadow: none;
}
</style>
# Memory & Cache
*(AMH chapters 8-9)*
Performance Engineering, Lecture 4
Sergey Slotin
May 14, 2022
---
```nasm
; *c = *a + *b
mov eax, DWORD PTR [rsi]
add eax, DWORD PTR [rdi]
mov DWORD PTR [rdx], eax
```
`add` = one cycle, `mov` = depends
----
----
----
- *total size* $M$
- *block size* $B$
- *latency* (time to fetch one byte)
- *bandwidth* (> B-times-latency implies concurrent operations)
- *cost* (the chip itself, power, maintenance, etc.) <!-- .element: class="fragment" data-fragment-index="1" -->
- *volatility* (whether the data is preserved on shutdown) <!-- .element: class="fragment" data-fragment-index="2" -->
----
| Type | $M$ | $B$ | Latency | Bandwidth | $/GB/mo |
|:-----|:---------|-----|---------|-----------|:------------------|
| L1 | 10K | 64B | 2ns | 80G/s | - |
| L2 | 100K | 64B | 5ns | 40G/s | - |
| L3 | 1M/core | 64B | 20ns | 20G/s | - |
| RAM | GBs | 64B | 100ns | 10G/s | 1.5 |
| SSD | TBs | 4K | 0.1ms | 5G/s | 0.17 |
| HDD | TBs | - | 10ms | 1G/s | 0.04 |
| S3 | $\infty$ | - | 150ms | $\infty$ | 0.02 + i/o |
For comparison: one virtual CPU core is $5-15/mo
(Approximate 2022 figures, [Google Cloud pricing](https://cloud.google.com/products/calculator?skip_cache=true))
----
{%youtube 3owqvmMf6No %}
---
### External Memory Model
- Dataset of size $N$ stored in *external* memory
- We can read and write in blocks of $B$ elements
- We can store $M$ elements in *internal* memory (= up to $\left \lfloor \frac{M}{B} \right \rfloor$ blocks)
- We only care about I/O operations: any in-between computations are free
- We additionally assume $N \gg M \gg B$
----
### Array scan
$$
\underbrace{a_1, a_2, a_3,} _ {B_1}
\underbrace{a_4, a_5, a_6,} _ {B_2}
\ldots
\underbrace{a_{n-3}, a_{n-2}, a_{n-1}}_{B_{m-1}}
$$
----
$$
SCAN(N) \stackrel{\text{def}}{=} O\left(\left \lceil \frac{N}{B} \right \rceil \right) \; \text{IOPS}
$$
----
```cpp
FILE *input = fopen("input.bin", "rb");
const int M = 1024;
int buffer[M], sum = 0;
// while the file is not fully processed
while (true) {
// read up to M of 4-byte elements from the input stream
int n = fread(buffer, 4, M, input);
// ^ the number of elements that were actually read
// if we can't read any more elements, finish
if (n == 0)
break;
// sum elements in-memory
for (int i = 0; i < n; i++)
sum += buffer[i];
}
fclose(input);
printf("%d\n", sum);
```
(Buffering also happens when you just read a stream)
---
### External Sorting
----
```cpp
void merge(int *a, int *b, int *c, int n, int m) {
int i = 0, j = 0;
for (int k = 0; k < n + m; k++) {
if (i < n && (j == m || a[i] < b[j]))
c[k] = a[i++];
else
c[k] = b[j++];
}
}
```
$SCAN(N+M)$
----
K-way merging?
Still $SCAN(N)$ as long as we can fit $k$ blocks in memory
($M \geq B^{1+ε}$ for $\epsilon > 0$)
<!-- .element: class="fragment" data-fragment-index="1" -->
----
### Merge Sorting
The first $O(\log M)$ layers of mergesort are free, and there are only $O(\log_2 \frac{N}{M})$ non-zero-cost layers, each mergeable in $O(\frac{N}{B})$ IOPS in total
<!-- .element: class="fragment" data-fragment-index="1" -->
$$
O\left(\frac{N}{B} \log_2 \frac{N}{M}\right)
$$
<!-- .element: class="fragment" data-fragment-index="2" -->
----
### $k$-way Mergesort
What if we split and merge $k$ arrays instead of two?
$$
SORT(N) \stackrel{\text{def}}{=} O\left(\frac{N}{B} \log_{\frac{M}{B}} \frac{N}{M} \right)
$$
<!-- .element: class="fragment" data-fragment-index="1" -->
----
### Joining
**Problem.** Given two lists of tuples $(x_i, a_{x_i})$ and $(y_i, b_{y_i})$,
output a list $(k, a_{x_k}, b_{y_k})$ such that $x_k = y_k$
**Solution:** Sort two lists and use two pointers to merge them <!-- .element: class="fragment" data-fragment-index="1" -->
(Beneficial to keep the arrays sorted) <!-- .element: class="fragment" data-fragment-index="1" -->
---
- *Temporal locality*: same data is accessed within a small time period
- *Spatial locality*: neighboring elements are accessed simultaneously
----
### breadth-first vs depth-first recursion
(E. g. merge sorting)
- Breadth-first: read the entire array $O(\log N)$ times <!-- .element: class="fragment" data-fragment-index="1" -->
- Depth-first: read the entire array $O(\log \frac{N}{M})$ times (on each cache level) <!-- .element: class="fragment" data-fragment-index="2" -->
----
### Dynamic Programming
----
```python
@lru_cache
def f(n, w):
# check if we have no items to choose
if n == 0:
return 0
# check if we can't pick the last item (note zero-based indexing)
if c[n - 1] > w:
return f(n - 1, w)
# otherwise, we can either pick the last item or not
return max(f(n - 1, w), c[n - 1] + f(n - 1, w - c[n - 1]))
```
----
```cpp
int f[N + 1][W + 1] = {0}; // this zero-fills the array
for (int n = 1; n <= N; n++)
for (int w = 0; w <= W; w++)
f[n][w] = c[n - 1] > w ?
f[n - 1][w] :
max(f[n - 1][k], c[n - 1] + f[n - 1][w - c[n - 1]]);
```
----
```cpp
bool f[W + 1] = {0};
f[0] = 1;
for (int n = 0; n < N; n++)
for (int x = W - c[n]; x >= 0; x--)
f[x + c[n]] |= f[x];
```
----
```cpp
std::bitset<W + 1> b;
b[0] = 1;
for (int n = 0; n < N; n++)
b |= b << c[n];
```
----
### Sparse Table
$$
t[k][i] = \min \{ a_i, a_{i+1}, \ldots, a_{i+2^k-1} \}
$$
----
```cpp
int rmq(int l, int r) { // half-interval [l; r)
int t = __lg(r - l);
return min(mn[t][l], mn[t][r - (1 << t)]);
}
```
----
$$
t[k][i] = \min(t[k-1][i], t[k-1][i+2^{k-1}])
$$
There are $2×2=4$ ways to build it
----
```cpp
int mn[logn][maxn];
memcpy(mn[0], a, sizeof a);
for (int l = 0; l < logn - 1; l++)
for (int i = 0; i + (2 << l) <= n; i++)
mn[l + 1][i] = min(mn[l][i], mn[l][i + (1 << l)]);
```
---
---
## Bandwidth
Like throughput, but for data transmission
----
```cpp
int a[N];
for (int t = 0; t < K; t++)
for (int i = 0; i < N; i++)
a[i]++;
```
----
----
----
Read only:
```cpp
for (int i = 0; i < N; i++)
s += a[i];
```
Write only:
```cpp
for (int i = 0; i < N; i++)
a[i] = 0;
```
----
Why is read-only faster?
Write needs to update cache, doubling traffic between CPU and RAM <!-- .element: class="fragment" data-fragment-index="1" -->
----
```cpp
const __m256i zeros = _mm256_set1_epi32(0);
for (int i = 0; i + 7 < N; i += 8)
_mm256_store_si256((__m256i*) &a[i], zeros);
```
```cpp
const __m256i zeros = _mm256_set1_epi32(0);
for (int i = 0; i + 7 < N; i += 8)
_mm256_stream_si256((__m256i*) &a[i], zeros); // "stream" operations bypass cache
```
<!-- .element: class="fragment" data-fragment-index="1" -->
----
- the memory controller doesn't have to switch the between reading and writing <!-- .element: class="fragment" data-fragment-index="1" -->
- the instruction sequence is simpler, allowing for more pending operations <!-- .element: class="fragment" data-fragment-index="2" -->
- the memory controller can "fire and forget" non-temporal write requests <!-- .element: class="fragment" data-fragment-index="3" -->
You can speed up `memset` and `memcpy` if you don't need output array right away <!-- .element: class="fragment" data-fragment-index="4" -->
---
## Latency
Time to fetch one byte
----
```cpp
int p[N], q[N];
// generate a random permutation
iota(p, p + N, 0);
random_shuffle(p, p + N);
// make a guaranteed single-cycle permutation out of it
int k = p[N - 1];
for (int i = 0; i < N; i++)
k = q[k] = p[i];
for (int t = 0; t < K; t++)
for (int i = 0; i < N; i++)
k = q[k];
```
This anti-pattern is called *pointer chasing*
----
----
(Not *real* latency because of random cache line sharing, but close enough)
----
----
---
## Cache lines
The main memory is split into blocks of $B=64$ bytes
----
Tweakable parameter `D`:
```cpp
for (int i = 0; i < N; i += D)
a[i]++;
```
----
$D=1$ increments using two vector instructions, while $D=16$ needs just one
(bottlenecked by write)
----
$D=8$ also needs one vector write
----
```cpp
struct padded_int {
int val;
int padding[15];
};
padded_int q[N / 16];
// constructing a cycle from a random permutation
// ...
for (int i = 0; i < N / 16; i++)
k = q[k].val;
```
----
Each element is more likely to be kicked out of the cache
---
## Memory Sharing
RAM and L3 cache are shared between CPU cores
----
`lstopo` output
----
- `parallel` to run a command in parallel
- `taskset` to fix restrict them to specific cores (CPU affinity)
```bash
parallel taskset -c 0,1,2,3 ./run ::: {0..3}
```
----
----
```bash
parallel taskset -c 0,1 ./run ::: {0..1} # L3 cache sharing
parallel taskset -c 0,4 ./run ::: {0..1} # no L3 cache sharing
```
----
Workload affinity matters in NUMA systems
----
---
## Memory-Level Parallelism
Memory operations can overlap
(This is why pointer chasing is so slow)
----
----
A more direct experiment:
```cpp
const int M = N / D;
int p[M], q[D][M];
for (int d = 0; d < D; d++) {
iota(p, p + M, 0);
random_shuffle(p, p + M);
k[d] = p[M - 1];
for (int i = 0; i < M; i++)
k[d] = q[d][k[d]] = p[i];
}
for (int i = 0; i < M; i++)
for (int d = 0; d < D; d++)
k[d] = q[d][k[d]];
```
----
----
When # of "memory lanes" is fewer than # of registers:
```nasm
dec edx
movsx rdi, DWORD PTR q[0+rdi*4]
movsx rsi, DWORD PTR q[1048576+rsi*4]
movsx rcx, DWORD PTR q[2097152+rcx*4]
movsx rax, DWORD PTR q[3145728+rax*4]
jne .L9
```
Register spill otherwise:
```nasm
mov edx, DWORD PTR q[0+rdx*4]
mov DWORD PTR [rbp-128+rax*4], edx
```
---
## Prefetching
You can request the data you are likely to read ahead of time
----
If the access pattern is simple, the hardware will do the prefetching for you:
```cpp
int p[15], q[N];
iota(p, p + 15, 1);
for (int i = 0; i + 16 < N; i += 16) {
random_shuffle(p, p + 15);
int k = i;
for (int j = 0; j < 15; j++)
k = q[k] = i + p[j];
q[k] = i + 16;
}
```
Flat 3ns (as if it was in L1 cache) regardless of array size
----
- Hardware prefetcher activates when "natural" traffic is low
- It can only detect simple access patterns, e. g. linear iteration
- You can explicitly prefetch needed cache line with software: <!-- .element: class="fragment" data-fragment-index="1" -->
```cpp
__builtin_prefetch(&a[k]);
```
<!-- .element: class="fragment" data-fragment-index="1" -->
----
```cpp
const int n = find_prime(N); // largest prime not exceeding N
for (int i = 0; i < n; i++)
q[i] = (2 * i + 1) % n;
```
If $n$ is a prime, the period will be exactly $n$ (an LCG property)
----
We can now peek ahead:
```cpp
int k = 0;
for (int t = 0; t < K; t++) {
for (int i = 0; i < n; i++) {
__builtin_prefetch(&q[(2 * k + 1) % n]);
k = q[k];
}
}
```
----
----
$$
\begin{aligned}
f(x) &= 2 \cdot x + 1
\\ f^2(x) &= 4 \cdot x + 2 + 1
\\ f^3(x) &= 8 \cdot x + 4 + 2 + 1
\\ &\ldots
\\ f^k(x) &= 2^k \cdot x + (2^k - 1)
\end{aligned}
$$
----
Prefetch $D$-th element ahead:
```cpp
__builtin_prefetch(&q[((1 << D) * k + (1 << D) - 1) % n]);
```
(Due to MLP, this will concurrently execute $(D + 1)$ reads on average)
----
---
## Cache associativity
How are caches implemented in hardware?
----
Loop A:
```cpp
for (int i = 0; i < N; i += 256)
a[i]++;
```
Loop B:
```cpp
for (int i = 0; i < N; i += 257)
a[i]++;
```
Which one is faster? **B, by 10x!** <!-- .element: class="fragment" data-fragment-index="1" -->
----
----
Fully associative cache
----
Direct-mapped cache
----
Set-associative cache (2-way associative)
----
- *offset* — the index of the word within a 64B cache line ($\log_2 64 = 6$ bits)
- *index* — the index of the cache line set (the next $12$ bits as there are $2^{12}$ cache lines in the L3 cache)
- *tag* — the rest of the memory address, which is used to tell the memory blocks stored in the cache lines apart
When we iterate in steps of 256, all accessed elements map to the same cache line <!-- .element: class="fragment" data-fragment-index="1" -->
---
## Paging
The memory a process sees is not mapped directly to physical memory
----
```cpp
const int N = (1 << 13);
int a[D * N];
for (int i = 0; i < D * N; i += D)
a[i] += 1;
```
(Same as in the bandwidth experiment, but we adjust array size)
----
----
----
16G RAM / 4K page size = 4M pages
We can **increase page size** and/or introduce a separate **cache for the page table**
----
- The L1 TLB has 64 entries, and if the page size is 4K, then it can handle $64 \times 4K = 512K$ of active memory without going to the L2 TLB
- The L2 TLB has 2048 entries, and it can handle $2048 \times 4K = 8M$ of memory without going to the page table
$$
8K \times 256 \times 4B = 8M
$$
<!-- .element: class="fragment" data-fragment-index="1" -->
----
Page size is usually 4K but can be increased
```cpp
#include <sys/mman.h>
void *ptr = std::aligned_alloc(page_size, array_size);
madvise(ptr, array_size, MADV_HUGEPAGE);
```
----
----
---
### AoS vs SoA
Pointer chasing, but we have to read $D$ elements to compute the next pointer
----
```cpp
const int M = N / D; // # of memory accesses
int p[M], q[M][D]; // SoA is reversed (q[D][M])
iota(p, p + M, 0);
random_shuffle(p, p + M);
int k = p[M - 1];
for (int i = 0; i < M; i++)
q[k][0] = p[i];
for (int j = 1; j < D; j++)
q[i][0] ^= (q[j][i] = rand());
k = q[k][0];
}
for (int i = 0; i < M; i++) {
int x = 0;
for (int j = 0; j < D; j++)
x ^= q[k][j];
k = x;
}
```
----
----
----
```cpp
struct padded_int {
int val;
int padding[15];
};
const int M = N / D / 16;
padded_int q[M][D];
```
----
----
This must be a RAM-specific issue
----
{"metaMigratedAt":"2023-06-17T00:49:03.131Z","metaMigratedFrom":"YAML","title":"Memory & Cache","breaks":true,"slideOptions":"{\"theme\":\"white\"}","contributors":"[{\"id\":\"045b9308-fa89-4b5e-a7f0-d8e7d0b849fd\",\"add\":18680,\"del\":1581}]"}