SIMD today: packed extensions, GPUs, tensor cores and NPUs

Packed SIMD: the vector length written into the opcode

Start with the loop of the exam's strip-mining quiz, for (i=0; i<260; i++) Y[i] = X[i] + Y[i]; (lesson 21 takes its wrong options apart), and compile it not for VMIPS but for the CPU in your laptop.

From the mid-1990s the general-purpose CPUs (x86, ARM, PowerPC) added SIMD instruction sets (p. 414): Intel MMX (1996), then SSE, AVX and AVX-512; ARM NEON. A register of 128, 256 or 512 bits is treated as a short row of numbers packed side by side, and one instruction applies the same operation to all of them: 4 doubles in a 256-bit AVX register, 8 in a 512-bit AVX-512 register. The notes call these "SIMD ISA extensions", Hennessy & Patterson "multimedia SIMD extensions"; this lesson says packed SIMD.

For AVX, 260 = 65 × 4: the compiler emits a loop of 65 trips, one 4-wide add each. For AVX-512, 260 = 32 × 8 + 4: 32 trips of the 8-wide loop, and then four elements that no 8-wide add can take. They need a second loop, the scalar remainder loop:

// AVX-512: 8 doubles in every add
// main loop, 32 trips:
for (i = 0; i + 8 <= n; i += 8)
  Y[i:i+8] = X[i:i+8] + Y[i:i+8];
// remainder loop, 4 scalar trips:
for (; i < n; i++)
  Y[i] = X[i] + Y[i];

On a vector machine the same loop is the strip-mined loop of lesson 31. With MVL = 64, 260 = 4 × 64 + 4: five strips, each the same four vector instructions with VLR set to the strip's length, and no second loop:

// vector ISA: the length is in VLR
// strips: VLR = 64, 64, 64, 64, 4
for (i = 0; i < n; i += VLR) {
  VLR = min(n - i, MVL);
  // L.V, L.V, ADDVV.D, S.V:
  Y[i:i+VLR] = X[i:i+VLR]
             + Y[i:i+VLR];
}

The difference fits in one sentence, and Table 30 of the notes compares the two machines row by row: in packed SIMD the vector length is part of the opcode, in a vector ISA it is the content of a register.

Three rows — vector length, instruction count, loop handling — are one fact seen three times. The register width fixes the length, so one instruction covers only what fits in the register, and whatever does not fit at the end goes to the remainder loop. For a loop of \(n\) elements on a width of \(W\):

A vector instruction covers VLR elements whatever MVL is, and the last strip simply sets a shorter VLR: \(\lceil n/\mathrm{MVL} \rceil\) instructions per operation, no second loop.

The width in the opcode also freezes the binary. A program compiled for AVX keeps issuing 4-wide adds on a core with 512-bit hardware: half the lanes stay idle until someone recompiles it. A program compiled for AVX-512 does not run at all on a core without AVX-512. A vector binary runs on any number of lanes — the transparent scaling of lesson 31 — and gets longer strips when MVL grows, provided its loop takes MVL from the hardware instead of writing 64 into the code.

The two rows left, memory access and masking, are about what packed SIMD lacked at birth. It loaded only contiguous, aligned blocks. Strided access, gather/scatter and per-element masks — all in vector ISAs since the Cray era — arrived generations later, if at all: gather with AVX2 (strided loads only as gathers), scatter and the k-mask registers with AVX-512. The notes misdate gather:

What packed SIMD buys, and what it costs

Why did the CPU makers take this road anyway? The notes list three gains (p. 415): compatibility with scalar CPUs — "no need for separate vector hardware, just extra datapaths in the same ALU pipeline"; widespread adoption, so compilers and libraries target it; and a good fit for multimedia and data-parallel kernels: audio, video, ML inner loops. Hennessy & Patterson add the hardware side: short vectors add little state to save on a context switch, they are fed by the caches instead of the banked memory a vector unit needs, and a 32- or 64-byte access touches at most two pages, where one 64-element vector load can touch many.

The notes' list of losses (p. 415) mirrors the other column of Table 30: no performance portability to wider registers, compilers that must generate remainder loops, clumsy strided and sparse access. The fourth bullet carries the wrong label:

Packed SIMD writes the vector length into the opcode; a vector ISA keeps it in a register (VLR). Remainder loops and recompilation for every new width both follow from that one choice.

A GPU is a vector processor that juggles threads

Here is DAXPY, Y = a·X + Y, written for an NVIDIA GPU in CUDA. There is no loop and no vector instruction: you write the body for one element and launch as many copies of it as there are elements.

__global__
void daxpy(int n, double a,
           double *x, double *y)
{
  // which element am I?
  int i = blockIdx.x * blockDim.x
          + threadIdx.x;
  if (i < n)   // the tail idles
    y[i] = a * x[i] + y[i];
}
// host: n threads, 256 per block
daxpy<<<(n + 255) / 256, 256>>>
     (n, a, d_x, d_y);

Each copy is a CUDA thread with its own index i and its own registers. The launch line creates a grid of blocks of 256 threads, and the hardware cuts every block into groups of 32 consecutive threads, the warps. A warp is the real unit of execution: one instruction is fetched for it, and its 32 threads execute that instruction in lockstep on the SM's SIMD lanes, each on its own registers. In the words of lessons 30–31, a warp instruction is a vector instruction with VL = 32: on 16 lanes it takes 32 / 16 = 2 cycles. NVIDIA calls the style SIMT: you program threads, the machine runs vectors.

The hardware that runs the warps is NVIDIA's streaming multiprocessor (SM); Hennessy & Patterson call it a multithreaded SIMD processor. Figure 30 of the notes (p. 392) draws it next to a vector processor, and the widget turns one into the other.

The tail and the if: masks instead of VLR

A GPU has no vector-length register. The kernel's if (i < n) does its job: in the last block, the threads with i ≥ n are switched off by the hardware's per-thread mask and their lanes idle — the masked tail of the first widget. The same mechanism handles any if. When the threads of a warp disagree, the warp runs both paths, each under the mask of its own threads: the vector mask register of lesson 31, set by the hardware. This is branch divergence, and it costs what masking costs. A warp split 8 against 24 across an if-else issues both bodies with part of its lanes idle; a warp whose 32 threads take 32 different paths runs at 1/32 of its speed.

Memory: coalescing finds the stride at run time

VMIPS knows the access pattern from the instruction: L.V means unit stride. A GPU receives 32 independent addresses, one per thread. The address coalescing unit compares them and merges them into as few wide transactions as it can: 32 threads reading 32 consecutive doubles cost a couple of transactions, 32 threads walking down a column of a row-major matrix cost 32. It is the bad-stride story of lesson 31 again. The cure is a traversal order or a data layout in which consecutive threads touch consecutive words; padding the rows, the cure for bank conflicts, does not merge a column's addresses.

Latency: switch warps instead of waiting

A load that misses goes to DRAM and takes hundreds of cycles (the memory wall of lesson 22). A vector machine pays that latency once per 64-element vector and streams behind it. A GPU has short vectors, 32 elements, and no out-of-order engine. Instead, every SM keeps dozens of warps resident, and its warp scheduler picks, every cycle, one whose operands are ready — the fine-grained multithreading of lesson 29. The switch costs nothing because the registers of every resident warp are already in the register file: that is what the 1024 registers per lane of Figure 30 are for.

How many warps are enough? Each warp offers a run of \(r\) instructions before its next stall, and the other warps must fill the \(t_{wait}\) cycles of its wait:

In Fig. 1, 1 + 6/2 = 4 warps are exactly enough. With runs of 10 instructions and a 400-cycle DRAM access you need 41. That is why an SM is built to hold dozens of warps at once, and why a single GPU thread is slow while the GPU is fast.

A streaming multiprocessor is a multithreaded SIMD processor: a warp is one vector instruction stream of length 32 with per-lane registers; masks replace the VLR, coalescing replaces the strided load, and switching warps replaces long vectors as the way to hide memory latency.

Host and device: two memories, explicit copies

The kernel is half of a CUDA program. The other half runs on the CPU, the host, and its main job is moving data, because the GPU, the device, has its own DRAM (GDDR or HBM) on its own board, linked to the host by PCIe or NVLink. DAXPY on the host side:

size_t bytes = n * sizeof(double);
double *d_x, *d_y;
// 1. reserve device memory
cudaMalloc(&d_x, bytes);
cudaMalloc(&d_y, bytes);
// 2. copy the inputs in
cudaMemcpy(d_x, x, bytes,
           cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, bytes,
           cudaMemcpyHostToDevice);
// 3. run the kernel on the device
daxpy<<<(n + 255) / 256, 256>>>
     (n, a, d_x, d_y);
// 4. copy the result out
cudaMemcpy(y, d_y, bytes,
           cudaMemcpyDeviceToHost);

The pointers d_x and d_y hold addresses in the device's memory. In the basic model the host reaches that memory only through cudaMemcpy, and the kernel cannot see x and y in host memory. That is the quiz of the June 2025 paper (Quizzes, Question 3, p. 11): "In the GP-GPU programming model, the CPU (host) and GPU (device) have separate memory address spaces." True. Later CUDA versions added unified memory, which migrates pages between the two memories automatically, but it is a convenience laid on top: the basic model, and the hardware, keep two memories, and every byte crosses the link.

The split has a price: the copies cost time, and the GPU's own bandwidth serves only data that is already on the device, so a kernel must do enough work per byte copied to pay for the trip.

In the GP-GPU programming model the host and the device have separate memory address spaces: the program allocates device memory, copies the data in, launches the kernel and copies the results back.

Tensor cores: make the operand a matrix

Multiply two 4 × 4 matrices and add a third: D = A × B + C. Each entry of D is a row of A times a column of B plus one entry of C — four multiply-adds, FMAs. Sixteen entries make 64 FMAs.

On ordinary lanes those 64 FMAs are 64 lane operations, each reading its own a and b from the register file: 128 operand reads. A tensor core, a functional unit NVIDIA added to every SM with the Volta generation (2017), takes the three 4 × 4 matrices and returns D as one operation: 64 FMAs per clock, A and B in FP16, C and D in FP16 or FP32. It reads the 32 values of A and B once and uses each of them 4 times inside the unit. A Tesla V100 has 8 tensor cores per SM, 640 in all, for up to 125 tensor TFLOPS (NVIDIA's figure). Software sees a warp-level instruction: the 32 threads of a warp hold 16 × 16 tiles of the matrices in their registers and issue the multiply together.

Why does it pay? For three reasons you have met before. One operation of the unit does 64 FMAs: the vector idea, in two dimensions. Each value read is used 4 times, so the register-file traffic per FMA drops by four. And FP16 halves the bytes per operand. The accumulation stays in FP32 because a product of two FP16 numbers loses little, while a sum of thousands of them would lose precision or overflow: accumulate wide, store narrow.

The cost is that a tensor core does nothing for work that is not a matrix product. DAXPY gains nothing, and a model from Hennessy & Patterson shows why.

The roofline: when is a faster unit useless?

The roofline model (H&P, ch. 4) caps the flop rate a kernel can reach by two ceilings: the peak of the arithmetic units, and the memory bandwidth times the kernel's arithmetic intensity — the flops it performs per byte it moves from DRAM:

DAXPY does 2 flops per element, a multiply and an add, and moves 24 bytes: X[i] and Y[i] loaded, Y[i] stored, 8 bytes each. A product of two \(n \times n\) matrices does \(2n^3\) flops on \(3n^2\) doubles, read and written once each:

DAXPY stays at 1/12 flop per byte whatever \(n\) is; the matrix product grows with \(n\), to 85 for \(n = 1024\) (four times more in FP16, the tensor cores' format). The slope meets a roof at its ridge point. A kernel left of the ridge, like DAXPY, is memory-bound and cannot use a higher roof; only a kernel to its right, like the large matrix product, can.

A tensor core turns the operand from a vector into a matrix tile: 64 FMAs per operation, each input used 4 times, in reduced precision. It raises the compute roof, so only kernels of high arithmetic intensity — matrix products — gain.

NPUs: the whole chip becomes the matrix unit

A neural processing unit (NPU) is a domain-specific accelerator for neural networks: a large array of low-precision multiply-accumulate units (MACs) fed from on-chip memory, with none of the general-purpose machinery around it. The reference design is Google's first TPU (Jouppi et al., ISCA 2017): a 256 × 256 array of 8-bit MACs at 700 MHz, 28 MiB of on-chip memory managed by software, and about a dozen CISC-style instructions (Read_Weights, MatrixMultiply, Activate …) sent by the host over PCIe. Each MAC counts as two operations, so the peak is 92 TOPS, tera-operations per second:

The array is systolic: data moves between neighbouring cells in lockstep, pumped like a heartbeat, so every value fetched from memory is used again and again on its way across. The TPU's dataflow is weight-stationary: the weights are loaded into the cells and stay there; the activations enter from the left, each row one cycle after the row above, and move right; the partial sums move down and leave at the bottom into the accumulators. The widget runs it on a 3 × 3 array.

The systolic array is to a matrix what the vector pipeline of lesson 30 is to a vector: the latency of the wavefront is paid once, then results stream out, one per column per cycle.

Hennessy & Patterson's five guidelines

Hennessy & Patterson (6th ed., ch. 7) distil such designs into five guidelines for domain-specific architectures. The TPU follows each of them:

1. Dedicated memories to minimise data movement — a 24 MiB unified buffer for activations and 4 MiB of accumulators, both managed by software; no caches.
2. Spend the resources saved on advanced microarchitecture on more arithmetic units or bigger memories — 65,536 MACs, and control logic under 2 % of the die.
3. The easiest form of parallelism that matches the domain — SIMD in two dimensions, a systolic array; no MIMD.
4. Reduce data size and type to the simplest the domain needs — 8-bit integers. Quantising Google's Inception network shrank it from 91 MB to 23 MB.
5. A domain-specific programming language — models written in TensorFlow, compiled for the TPU.

The NPUs of phones and laptops (Apple's Neural Engine, Qualcomm's Hexagon, the NPU blocks of recent Intel and AMD chips) apply the same recipe at a few watts: 8-bit or FP16 MAC arrays and large local SRAM.

Back to the puzzle

The TPU of the hook beat the CPUs and GPUs of its day without caches, branch prediction, out-of-order execution or multithreading. That machinery exists for programs that hide their future: branches whose direction is unknown until they execute, loads whose latency depends on a cache, independent instructions that the hardware must discover inside a window of a few hundred (lesson 8). A neural-network layer announces everything in advance — a matrix product of known size, a fixed dataflow, no data-dependent branches. Hardware that rediscovers what the program already states wastes silicon and energy; the TPU spends them on multipliers and on-chip memory instead. It is the choice Cray made in lesson 30, taken as far as it goes.

Every machine of this lesson is a step on the same road. Packed SIMD writes the parallelism into the opcode, a GPU into thousands of threads, a tensor core into a matrix tile, and an NPU into the whole chip.

An NPU is the vector idea pushed to two dimensions and to one domain: a systolic array that reuses every operand hundreds of times, fed from software-managed on-chip memory, in low precision, with none of the hardware that discovers parallelism at run time.