Floating-point performance on GPUs: precision, FLOPs, and numerical error

Prerequisites

This article assumes you have read Inside a modern NVIDIA GPU, specifically the sections on CUDA cores and tensor cores, and warp divergence for the SIMT execution model. You should know that a warp issues one instruction to 32 lanes at once and that GPUs support multiple floating-point precisions in hardware; this article explains what those precisions actually are, what they cost you numerically, and how to turn “cores times clock speed” into a peak FLOP/s number that means something.

Why floating-point is not one thing

Every number in this article’s code samples so far has been declared float or float32 without much comment. That convenience hides a real decision: GPUs support at least five distinct floating-point formats in hardware (FP64, FP32, TF32, FP16, and BF16), and the choice between them changes throughput, memory footprint, and the numerical correctness of your results simultaneously. Getting this decision wrong either wastes performance you left on the table by using more precision than you need, or silently corrupts results by using less precision than your computation can tolerate.

IEEE-754: sign, exponent, significand

A floating-point number is stored as three fields: a sign bit, an exponent (stored with a bias so it can represent both very large and very small magnitudes), and a significand (also called the mantissa, the fractional part of the number after an implicit leading 1). The value is reconstructed as:

value = (-1)^sign * 1.significand * 2^(exponent - bias)

graph LR
  S["Sign
1 bit"] --> E["Exponent
(biased)"]
  E --> M["Significand / mantissa
(fractional bits)"]
  style S fill:#e8744f,color:#fff
  style E fill:#4a90d9,color:#fff
  style M fill:#00CC96,color:#fff

The number of exponent bits determines the range (how large or small a number can be represented at all); the number of significand bits determines the precision (how many significant decimal digits the format can distinguish between).

Bit layout across formats

TF32 is not a full IEEE-754 type: it is NVIDIA’s tensor-core input format that stores a value in a 32-bit container but only reads 10 significand bits, giving it FP32’s exponent range (so it does not lose dynamic range) with roughly FP16’s precision. BF16 takes the opposite trade from FP16: it keeps FP32’s 8-bit exponent (so it has FP32’s range) but truncates the significand to only 7 bits, sacrificing precision to keep the same range, which makes it much more forgiving of the very large and very small magnitudes that show up during deep learning training.

Subnormals, infinity, and NaN

Every one of these formats reserves specific bit patterns for special values:

• Subnormals (denormals): when the exponent field is all zeros, the implicit leading 1 is dropped and the value is interpreted as (-1)^sign * 0.significand * 2^(1-bias). This lets the format represent numbers smaller than the smallest normal value, at reduced precision, instead of abruptly rounding to zero. Subnormal arithmetic is also frequently much slower in hardware than normal-range arithmetic, because it can fall back to a slower microcoded or software path.
• Infinity: exponent field all ones, significand all zeros. Produced by overflow (a value too large to represent) or by an operation like 1.0f / 0.0f.
• NaN (Not a Number): exponent field all ones, significand nonzero. Produced by indeterminate operations like 0.0f / 0.0f or sqrtf(-1.0f). NaN is “contagious”: almost any arithmetic operation involving a NaN produces another NaN, which is a useful debugging signal, since a single NaN in a large GPU buffer usually spreads visibly through downstream computation rather than vanishing silently.

Rounding and why 0.1 is not exact

Binary floating-point cannot represent every decimal fraction exactly, for the same reason decimal notation cannot represent 1/3 exactly: some fractions that terminate in one base do not terminate in another. 0.1 in decimal is an infinitely repeating binary fraction, so it gets rounded to the nearest representable FP32 or FP64 value, which is not exactly 0.1. The default IEEE-754 rounding mode is round-to-nearest, ties-to-even: round to the closest representable value, and if a value is exactly between two representable values, round to whichever of the two has an even last bit. This specific tie-breaking rule exists to avoid a statistical bias that would accumulate if ties always rounded the same direction.

>>> 0.1 + 0.2
0.30000000000000004

This is not a bug and not specific to Python; it is a direct consequence of the binary representation, and the exact same effect occurs identically in CUDA C++ float/double arithmetic. The practical lesson is to never compare floating-point results with ==; always compare with a tolerance (fabs(a - b) < epsilon, or the relative-tolerance form used later in this article).

Non-associativity: why order changes the answer

Real-number addition is associative: (a + b) + c == a + (b + c) always. Floating-point addition is not, because each addition rounds its result to the nearest representable value, and rounding a partial sum before adding the next term can lose different information depending on the order of operations.

a = 1e16, b = 1.0, c = -1e16

(a + b) + c:
  a + b  rounds to 1e16 (1.0 is too small to change a 1e16-magnitude float)
  1e16 + c = 1e16 - 1e16 = 0.0

a + (b + c):
  b + c  = 1.0 + (-1e16) rounds to -1e16
  a + (-1e16) = 1e16 - 1e16 = 0.0

But with different magnitudes, e.g. a = 1e8, b = 1.0, c = 1.0, d = -1e8:
  ((a + b) + c) + d  can differ from  a + (b + (c + d))
  depending on which additions lose the small terms to rounding first.

This matters enormously on a GPU because a parallel reduction (summing an array with a tree of concurrent partial sums across thousands of threads) adds the same set of numbers in a different order than a sequential CPU loop does. Both results are “correct” in the sense that both are valid IEEE-754 rounding outcomes, but they are not bit-identical, and a test that asserts exact equality between a GPU reduction and a CPU reference sum will fail even when both implementations are bug-free. Parallel prefix sum and reduction covers the algorithmic side of this; here, the point is that the numerical difference is expected, not a defect to chase.

FMA and the two-FLOP convention

Almost every floating-point-heavy kernel is dominated by the pattern result = a * b + c. CUDA cores (and tensor cores) implement this as a single hardware instruction, FMA (fused multiply-add), which computes the product and the addition with only one final rounding step, rather than rounding once after the multiply and again after the add. This makes FMA both faster (one instruction instead of two) and more accurate (one rounding error instead of two) than computing the multiply and add separately.

Because FMA does the work of a multiply and an add in one instruction, the convention in GPU computing is to count each FMA as two FLOPs (one multiply, one add) even though it is one hardware instruction. This convention is why peak FLOP/s figures for GPUs look larger than a naive “one instruction per cycle” count would suggest: the hardware genuinely does twice the arithmetic work per FMA instruction that it does per plain add or multiply instruction.

// The compiler will typically fuse this into a single FMA instruction:
float result = a * b + c;

// Explicit intrinsic form, if you want to guarantee FMA (and its single
// rounding step) rather than rely on compiler fusion:
float result = fmaf(a, b, c);

Computing peak FLOP/s

Peak (theoretical, not achieved) FLOP/s for CUDA-core arithmetic follows directly from the FMA convention:

Peak FLOP/s = (CUDA cores) * (clock speed in Hz) * (FLOPs per core per cycle)

Since each CUDA core can retire one FMA per cycle in the steady state, and one FMA counts as 2 FLOPs, the FLOPs-per-core-per-cycle term is 2 for FP32 CUDA-core throughput. As a worked example with illustrative, hypothetical numbers (always use cudaGetDeviceProperties and your GPU’s own datasheet for real figures, never a number copied from a tutorial):

Hypothetical device: 6912 FP32 CUDA cores, 1.4 GHz boost clock

Peak FP32 FLOP/s = 6912 cores * 1.4e9 Hz * 2 FLOPs/cycle
                  = 6912 * 1.4e9 * 2
                  = 1.93536e13 FLOP/s
                  = 19.3536 TFLOP/s

Tensor cores compute peak throughput differently, because a single tensor-core instruction performs many more FMA-equivalent operations per cycle by processing a small matrix multiply-accumulate in one shot rather than one scalar FMA at a time. This is why a GPU’s tensor-core peak FLOP/s (often quoted separately, and typically several times higher than its CUDA-core FP32 peak) cannot be derived from the CUDA-core core count at all; it depends on the tensor core’s specific operand shape and precision, which differs by architecture generation and is published per-precision on the datasheet (for example, separate peak numbers for FP16, BF16, TF32, and INT8 tensor-core throughput).

FLOP/s vs. bytes: arithmetic intensity and the roofline boundary

A kernel’s achievable throughput is bounded by whichever resource it exhausts first: compute (FLOP/s) or memory bandwidth (bytes/s). The ratio between how much arithmetic a kernel does and how much data it moves is its arithmetic intensity:

Arithmetic intensity (FLOPs/byte) = Total FLOPs performed / Total bytes moved (from/to DRAM)

The roofline model plots achievable performance (FLOP/s, y-axis) against arithmetic intensity (FLOPs/byte, x-axis). Two ceilings bound every kernel: a flat line at the device’s peak FLOP/s, and a diagonal line whose slope is the device’s peak memory bandwidth. The point where the diagonal memory-bandwidth line meets the flat compute-peak line is the roofline boundary (also called the machine balance point):

Roofline boundary (FLOPs/byte) = Peak FLOP/s / Peak memory bandwidth (bytes/s)

A kernel with arithmetic intensity below this boundary is memory-bound: no matter how fast the arithmetic units are, the kernel cannot outrun the rate at which data arrives from DRAM. A kernel above the boundary is compute-bound: the arithmetic units, not the memory system, are the limiting resource. Vector addition (one FMA-equivalent op per three array accesses) sits far to the left of this boundary on almost every GPU; a well-tiled dense matrix multiply sits far to the right. Debugging and profiling CUDA programs and performance case studies cover how Nsight Compute plots your actual kernel on this chart and how to interpret where it lands; this article’s job is to make sure the two axes of that chart (a FLOP/s number and an arithmetic-intensity number) are things you know how to compute yourself.

Tensor cores and accumulation precision

Tensor cores generally accept reduced-precision inputs but accumulate the running sum in higher precision, specifically to control the numerical error that reduced-precision multiplies would otherwise compound across many accumulation steps:

This is the mechanism that makes mixed-precision training viable at all: the individual multiplies are cheap and fast in reduced precision, but the accumulator (which sums potentially thousands of terms in a large matrix multiply’s inner dimension) does not compound reduced-precision rounding error thousands of times over, because it is carrying full FP32 precision throughout the accumulation.

Where reduced precision becomes unsafe

Reduced precision is not free of risk even with FP32 accumulation. Two failure modes are common enough to name explicitly:

• Overflow in the input format itself. FP16’s exponent range tops out around 65,504. A value that would be perfectly representable in FP32 (say, an unscaled gradient or activation in the millions) simply becomes Inf in FP16 before it ever reaches the accumulator. This is why mixed-precision training uses loss scaling: multiply the loss (and therefore the gradients) by a scale factor before the backward pass so that small gradient values do not underflow to zero in FP16, then divide back out before the optimizer step.
• Silent precision loss in reductions with many terms, even when every individual value is comfortably in range. Summing millions of FP16 or BF16 values sequentially, without a higher-precision accumulator, means each addition’s rounding error compounds; the running sum’s relative error grows with the number of terms rather than staying bounded, and can become large enough to change a result meaningfully at scale.

The rule of thumb: reduced precision is safe for the storage and multiply side of an operation when you have controlled for its dynamic range (via scaling or by choosing BF16 for its FP32-equivalent range), and it is safe for accumulation only when the hardware or your algorithm explicitly promotes the running sum to a higher-precision accumulator.

Reproducibility, tolerances, and controlling error

Because parallel floating-point results are not bit-identical to serial ones (non-associativity, above) and because different GPU architectures, library versions, or even different kernel launch configurations can choose different reduction orders internally, exact bit-for-bit reproducibility across runs is not something CUDA gives you by default. Practical strategies:

• Compare with tolerances, not equality. Use a combined absolute/relative tolerance check, the same pattern numpy.allclose and torch.allclose implement: abs(a - b) <= atol + rtol * abs(b). Choose rtol/atol based on the precision you actually used (FP16 needs a much looser tolerance than FP64) and the number of accumulation steps in your computation (more terms summed generally means more accumulated rounding error to tolerate).
• Kahan summation (compensated summation) tracks a running compensation term that captures the low-order bits lost to rounding on each addition and feeds them back into the next addition, dramatically reducing the growth of accumulated error in a long sequential sum, at the cost of roughly four times the arithmetic per addition.
• Pairwise (tree) summation sums an array by recursively summing halves and combining the two partial sums, rather than accumulating strictly left to right. This bounds the growth of rounding error to logarithmic in the number of terms rather than linear, and is naturally what a well-designed parallel reduction already does, which is one reason parallel reductions are not automatically less accurate than a naive sequential CPU sum, only differently rounded.
• Scaling inputs to keep values within a numerically well-behaved range for the precision in use (the same idea as loss scaling above, generalized beyond training) avoids both overflow and the precision cliff near the smallest representable magnitudes.
• Library-level determinism controls. Frameworks and libraries built on CUDA expose explicit switches when you need run-to-run reproducibility badly enough to accept a performance cost: for example, PyTorch’s torch.use_deterministic_algorithms(True) and the CUBLAS_WORKSPACE_CONFIG environment variable it depends on for deterministic cuBLAS behavior, or cuDNN’s deterministic algorithm selection flags. These typically force a fixed, often slower, algorithm variant instead of letting the library pick whichever variant is fastest for the current input shape (which is exactly the source of run-to-run nondeterminism in the first place).

In practice

Default to FP32 for anything where you have not specifically verified reduced precision is safe; it is the format every CUDA feature, library, and debugging tool supports without caveats. Move to FP16/BF16/TF32 deliberately, for the specific operations (usually large matrix multiplies inside a training or inference loop) where the throughput gain is large and you can either rely on hardware FP32 accumulation or add scaling. Never assert bit-exact equality between a GPU result and a CPU reference; assert a tolerance appropriate to the precision and the number of accumulation steps involved. When a result looks suspiciously far off, check for Inf/NaN first (both propagate visibly), then check for reduced-precision overflow or underflow before assuming your kernel logic itself is wrong.

What comes next

You now know what a floating-point number actually is at the bit level, why FMA counts as two FLOPs, how to turn hardware specs into a peak throughput number, and where the roofline boundary between memory-bound and compute-bound comes from. The next article, CUDA memory hierarchy, covers the other half of that roofline equation in depth: where data actually lives (registers, shared memory, caches, global memory) and how each memory space’s bandwidth and latency shape the arithmetic-intensity threshold you just learned to compute.