100 DAYS OF GPU Course preview
Day 71Module 8
in-technical-review

Mixed precision: FP16, BF16, TF32, FP8, FP4

Store 300,000 in a 16-bit float. In __half you get inf: the format tops out at 65,504. Store the same number in __nv_bfloat16 and you get 299,008, off by a third of a percent.

Now store 1.0007. __half keeps 1.0009765625, four correct digits; bfloat16 answers 1.0.

The same sixteen bits give each format a different range and precision. Neither is "half the accuracy of FP32".

Accuracy is not one number.

By the end of this page you can read a float format by its range and precision, convert day 16's matmul to FP16 storage with FP32 accumulation, and put a bound on its error that a program checks rather than relying on a broad claim.

Range and precision use different bits

A binary float has one sign bit, exponent bits, and mantissa bits. More exponent bits increase the range before inf. More mantissa bits reduce the gap between representable values.

Mixed precision uses different formats in different parts of one computation. The table shows how each format divides its bits.

Format bits exp mantissa largest finite gap at 1.0 tensor-core mma needs
FP32 32 8 23 3.4e38 1.19e-7 none (CUDA cores)
FP16 16 5 10 65,504 9.77e-4 sm_75
BF16 16 8 7 3.39e38 7.81e-3 sm_80
TF32 32 see below see below FP32's range at most 9.77e-4 sm_80
FP8 E4M3 8 4 3 448 0.125 sm_89
FP8 E5M2 8 5 2 57,344 0.25 sm_89
FP4 E2M1 4 2 1 6.0 0.5 sm_120a

The bit layouts are from the PTX ISA, section 5.2.3 "Alternate Floating-Point Data Formats" (https://docs.nvidia.com/cuda/parallel-thread-execution/index.html , checked 2026-08-29). The largest-finite and gap columns follow from those layouts by arithmetic: the gap at 1.0 is 2^-mantissa-bits, E4M3 uses its top exponent code on values rather than infinity so it reaches 448 with no inf at all, and E2M1 has neither infinity nor NaN.

The mma minimum capabilities are the Target ISA Notes of the same document's mma instruction, checked 2026-09-01: ".f16 floating point type mma operation with .m16n8k8 shape requires sm_75 or higher", ".bf16 alternate floating point type mma operation with .m16n8k8 and .m16n8k16 shapes requires sm_80 or higher", ".tf32 alternate floating point type mma operation with .m16n8k4 and .m16n8k8 shapes requires sm_80 or higher", ".e4m3 and .e5m2 alternate floating point type mma operation requires sm_89 or higher", and ".e3m2, .e2m3 and .e2m1 alternate floating point type mma operation requires sm_120a and is supported on sm_120f from PTX ISA version 8.8". Datacenter Blackwell reaches FP4 through a different instruction family entirely, which day 78 covers.

For the block-scaled forms that make 4-bit weights usable, see MX block scaling; day 98 runs them.

The ISA leaves TF32's layout open. The document says it is "a special 32-bit floating point format supported by the matrix multiply-and-accumulate instructions, with the same range as .f32 and reduced precision (>=10 bits)", and that "the internal layout of tf32 format is implementation defined".

Most explainers draw TF32 as 1 sign, 8 exponent and 10 mantissa bits. The ISA does not promise that, so this page and the widget below do not draw it.

Compare the FP16 and BF16 rows with the opening values. FP16 uses three more mantissa bits, so it resolves 1.0007 and overflows above 65,504. BF16 keeps FP32's eight exponent bits, so it stores 300,000 but cannot tell 1.0007 from 1.0.

The fp16-vs-bf16-range frame below shows that comparison. Toggle bits to see how each change affects the range and precision.

Why the accumulator decides how error grows

A common assumption is that a format's error is a fixed percentage, so FP16 storage means FP16-sized answers everywhere. That is true for a single rounding and false for a running sum, and a matmul is a running sum: every output element is a dot product over K terms.

An FP16 matmul has two sources of error. Storage rounding happens once per input: each value moves by at most one part in 2^11 when it becomes __half, so each product carries a relative error near eps = 2^-10. Those per-term errors do not share a sign, so a standard independent-error model estimates growth at sqrt(K) times eps rather than K times eps.

The accumulator type determines the other error. With FP32 accumulation, each add contributes rounding at about 2^-23, three orders of magnitude below the storage term.

Accumulate in FP16 and every partial sum is rounded back to 10 mantissa bits, at the scale of the running sum rather than of one term, so this term grows with K and can exceed the storage error. FP32 accumulation makes the growing term negligible and leaves the sqrt(K) term.

Day 66 turned that model into the gate this course uses: rtol = max(table value, 4 x eps x sqrt(K)), with the factor 4 as slack for a different-but-valid summation order, and atol = rtol x max|reference|. For FP16 at K = 256, 1024 and 2048 the rtol works out to 0.0625, 0.125 and 0.177. Those are ceilings for any valid order of the same arithmetic, not estimates of where a run lands; the Results section predicts that the measured error sits far under them.

Formats depend on hardware as well as their encodings. The mma minimum capabilities in the table mean a compute capability 7.5 card executes FP16 on its tensor cores and none of the others.

This page's program reads the capability at run time and prints a yes or no per format, so the transcript states what the card could not have been doing. On a compute capability 7.5 card, the report is one yes and four nos.

Three pipelines, one geometry

Full program in matmul_fp16.cu, under code/day71-mixed-precision/. It runs day 16's tiled kernel three ways, FP32 throughout, __half storage with a float accumulator, and __half storage with a __half accumulator, at 256, 1024 and 2048.

The reference is double precision over the original FP32 inputs. Not over the half-rounded copies. A reference built from rounded inputs would forgive the storage error, and the storage error is half of what this day bounds.

The inputs are chosen so FP16 has to round. Day 16 filled its matrices with small integers so a mismatch had to be an indexing bug. Those integers are also exact in __half: reuse them here and the half paths report an error of exactly zero, a test that cannot fail and therefore proves nothing. This program scales by 0.013 and 0.017, values no binary format represents, so every conversion causes the rounding being measured.

All three kernels share tiles, guards and launch shape, so a table row differs from the row above it in number format and nothing else. The half-storage kernel converts each value to float in a register at the moment of use:

__global__ void matmulF16StoreF32Acc(const __half* __restrict__ a,
                                     const __half* __restrict__ b,
                                     float* __restrict__ c, size_t m, size_t n,
                                     size_t k) {
    __shared__ __half tileA[kTileDim][kTileDim];
    __shared__ __half tileB[kTileDim][kTileDim];

    const size_t col =
        blockIdx.x * static_cast<size_t>(blockDim.x) + threadIdx.x;
    const size_t row =
        blockIdx.y * static_cast<size_t>(blockDim.y) + threadIdx.y;
    const size_t tiles = (k + kTileDim - 1) / kTileDim;

    float acc = 0.0f;
    for (size_t tileIdx = 0; tileIdx < tiles; ++tileIdx) {
        const size_t aCol = tileIdx * kTileDim + threadIdx.x;
        const size_t bRow = tileIdx * kTileDim + threadIdx.y;

        tileA[threadIdx.y][threadIdx.x] =
            (row < m && aCol < k) ? a[row * k + aCol] : __float2half(0.0f);
        tileB[threadIdx.y][threadIdx.x] =
            (bRow < k && col < n) ? b[bRow * n + col] : __float2half(0.0f);
        __syncthreads();

        for (int p = 0; p < kTileDim; ++p) {
            acc += __half2float(tileA[threadIdx.y][p]) *
                   __half2float(tileB[p][threadIdx.x]);
        }
        __syncthreads();
    }

    if (row < m && col < n) {
        c[row * n + col] = acc;
    }
}

The gate each kernel must pass prints its own arithmetic, so a transcript reader can recompute the bound:

static double kScaledRtol(double tableRtol, int mantissaBits, size_t kDim) {
    const double eps = std::ldexp(1.0, -mantissaBits);
    const double grown = 4.0 * eps * std::sqrt(static_cast<double>(kDim));
    return grown > tableRtol ? grown : tableRtol;
}

Timing is CUDA events, three warm-ups and the mean of ten runs per kernel, copies excluded, exactly the day 9 discipline. The __half tiles in shared memory cost half the bytes of the float pair, and the global loads cost half the sectors, which is the only mechanism behind any speedup here. Nothing in this program touches a tensor core, so the FP16 win here is bandwidth, not arithmetic.

The FP16-accumulate kernel is included for measurement. The program reports its max-error ratio against FP32 accumulation, but does not make that ordering a correctness gate: cancellation can reverse the ordering for one input even when both kernels satisfy their error bounds.

Results

Re-verified on the same Tesla T4 with driver 580.173.02 and CUDA 13.0 (V13.0.88). All format, error and correctness results reproduced exactly.

Timing changed at the two larger sizes: the FP32-accumulate half path improved at 1024 while the FP16-accumulate path slowed, and the 2048 FP32 baseline slowed. Both CUDA 12.6 and CUDA 13.0 transcripts are listed in front matter; the table below is the CUDA 13.0 run. Timings are means of ten runs after three warm-ups, with copies excluded.

case kernel max abs err gate frac time (ms) vs f32
256 matmulTiledF32 4.582e-06 0.0348 0.0832 1.00
256 matmulF16StoreF32Acc 2.595e-03 0.0036 0.0752 1.11
256 matmulF16StoreF16Acc 2.950e-02 0.0368 0.0745 1.12
1024 matmulTiledF32 1.740e-05 0.1253 5.7951 1.00
1024 matmulF16StoreF32Acc 3.642e-03 0.0026 4.1275 1.40
1024 matmulF16StoreF16Acc 1.867e-01 0.1657 3.8691 1.50
2048 matmulTiledF32 7.610e-05 0.1952 31.8951 1.00
2048 matmulF16StoreF32Acc 2.016e-03 0.0004 23.4591 1.36
2048 matmulF16StoreF16Acc 5.078e-01 0.1467 21.8186 1.46

The transcript settles four predictions.

  1. Refuted: the baseline did not reproduce day 44 within a few percent. The 2048 FP32 row was 31.8951 ms, farther from day 44's 25.350 ms than the CUDA 12.6 run's 27.8516 ms. The environment stamp changed, so the warning attached to the bet applies, but the numeric prediction still missed.
  2. Partly held, refuted overall. At 256 the half-storage, FP32-accumulate path was 1.11x, inside the predicted under-1.2x cache band. At 2048 it reached only 1.36x, below the predicted 1.4x to 2.0x DRAM band. The memory-bound explanation did not predict the magnitude.
  3. Partly held, refuted overall. FP32-accumulate max error was not monotonic: 2.595e-03, 3.642e-03 and 2.016e-03. At 2048 it was about 1.19e-04 of the maximum reference magnitude, below the predicted order 1e-3 relative. The other half held sharply: FP16 accumulation error rose from 2.950e-02 to 1.867e-01 to 5.078e-01, and its max-error ratio against FP32 accumulation rose from 11.4 to 51.3 to 252.
  4. Held. The format report printed FP16 yes and BF16, TF32, FP8 and FP4 no. All three cases passed every correctness gate and the process exited

These milliseconds belong to one T4. Across cards, compare the ratio between the two error columns and whether speed changes with the working set.

Run it yourself

Any card at sm_75 or above builds and runs the program:

nvcc -std=c++17 -O3 -arch=sm_75 -o matmul_fp16 matmul_fp16.cu

There is no Compiler Explorer embed because the double-precision CPU reference at 2048 is 8.6 billion multiply-adds, which does not fit a 20 second execution cap. On an sm_80+ card the same source builds unchanged and the format report's answers change, because the kernels themselves use nothing above sm_75.

Exercise

Do the conversion yourself instead of reading this one. Start from day 16's tiled_matmul.cu, convert it to __half storage with a float accumulator, gate it against a double reference with the K-scaled tolerance, and time it against the FP32 original at 2048.

Time: 30 to 45 minutes. Submit: your converted kernel, its max abs error and gate fraction next to the FP32 kernel's, the two times at 2048, and one sentence naming where the speedup came from.

Check: the shipped program is the reference implementation; your kernel must pass the same gate it prints, |got-ref| <= atol + rtol*|ref| with rtol = max(1e-2, 4 x 2^-10 x sqrt(K)), and any inf or NaN anywhere is an automatic fail with the index printed. The checker cannot detect one weak test: if you keep day 16's integer makeInputs, your error will be exactly zero and your gate proves nothing. Change the inputs so the format has to round.

Hint 1

Only the storage changes. Every declaration that holds a matrix element becomes __half; every piece of arithmetic stays float. Where is the one place a stored value must turn back into a float, and what is the cheapest place in the memory hierarchy for that to happen?

Hint 2

Two epsilons matter: 2^-10 for the stored inputs, 2^-23 for the accumulator. Write the error budget as one term per epsilon and ask how each grows with K. Which term comes from the accumulator, and what happens to it if the accumulator is __half too?

Solution

The diff is mechanical: const __half* parameters, __half tiles, __float2half(0.0f) for the guard fills, and acc += __half2float(...) * __half2float(...) in the inner loop, with acc still float. The full version is matmul_fp16.cu's matmulF16StoreF32Acc.

Every global load and every shared-memory read moves half the data. Since this kernel was bandwidth-limited, that lower byte count causes the speedup. The error stays at the storage term, sqrt(K) x 2^-10-ish, because the float accumulator's own rounding is a thousand times finer.

Storage precision is a bandwidth decision. Accumulator precision is an accuracy decision. Mixed precision lets you choose them separately.

Pitfalls

Your FP16 kernel reports zero error and you conclude FP16 adds no error. Your test inputs are exactly representable in half: integers, or multiples of a power of two, stay bit-exact through conversion and the matmul of exact halves in a float accumulator can be exact. The gate passed because it could not fail. Feed values that force rounding, the way this page's 0.013 and 0.017 scales do.

The error was small at 512 and large at 8192. A per-element tolerance that ignores K will do that: the storage term grows like sqrt(K) and an FP16 accumulator's term grows faster. Scale the gate with K, and if the growth exceeds sqrt(K), suspect the accumulator, not the storage. Day 68 covers why the answer also moves when the summation order does.

Outputs are inf and only at large sizes. 65,504 is a modest number. A dot product whose partial sums cross it in a half accumulator overflows even when every input and the true answer are small. FP32 accumulation raises the limit; rescaling inputs only changes where overflow occurs.

This program's inputs stay below 1 so that its half-accumulate row can show error growth without overflowing. Real data may not allow that choice.

Your "FP32" cuBLAS baseline uses TF32 without an explicit choice. On an sm_80+ card, letting the compute type change to CUBLAS_COMPUTE_32F_FAST_TF32 puts cuBLAS on tensor cores with about 10 mantissa bits while your kernel does real FP32. Comparisons and error budgets both become invalid. Set the compute type explicitly, the way day 44's harness pins CUBLAS_COMPUTE_32F (https://docs.nvidia.com/cuda/cublas/index.html#cublasmath-t , checked 2026-08-29).

You benchmarked FP16 and reported a tensor-core speedup. This page's kernels never touch a tensor core; any win is bandwidth. The instruction that would use them is day 72's subject, and day 47 already showed how easily a speedup gets attributed to the wrong mechanism when the binary is not checked.

Go deeper

Next

Day 72 uses these formats through the WMMA API, which runs matrix operations on tensor cores. Day 44's hand-written FP32 kernel reached 67.7 percent of cuBLAS on the measured card. The next lesson tests whether a tensor-core path beats that CUDA-core path.