In 2008 the IEEE 754 standard was revised to include the fused multiply-add operation (FMA). The FMA operation computes $\text{rn}(X\times Y+Z)$ with only one rounding step. Without the FMA operation the result would have to be computed as $\text{rn}\left(\text{rn}\right(X\times Y)+Z)$ with two rounding steps, one for multiply and one for add. Because the FMA uses only a single rounding step the result is computed more accurately.

Let's consider an example to illustrate how the FMA operation works using decimal arithmetic first for clarity. Let's compute $x^2-1$ with four digits of precision after the decimal point, or a total of five digits of precision including the leading digit before the decimal point.

For $x=1.0008$ , the correct mathematical result is $x^2-1=1.60064\times {10}^{-4}$. The closest number using only four digits after the decimal point is $1.6006\times {10}^{-4}$. In this case $\text{rn}(x^2-1)=1.6006\times {10}^{-4}$ which corresponds to the fused multiply-add operation $\text{rn}(x\times x+(-1\left)\right)$. The alternative is to compute separate multiply and add steps. For the multiply, $x^2=1.00160064$, so $\text{rn}\left(x^2\right)=1.0016$. The final result is $\text{rn}\left(\text{rn}\right(x^2)-1)=1.6000\times {10}^{-4}$.

Rounding the multiply and add separately yields a result that is off by 0.00064. The corresponding FMA computation is wrong by only 0.00004, and its result is closest to the correct mathematical answer. The results are summarized below:

$\begin{array}{rll}x& =& 1.0008\\ x^2& =& 1.00160064\\ x^2-1& =& 1.60064\times {10}^{-4}& \text{true value}\\ \text{rn}(x^2-1)& =& 1.6006\times {10}^{-4}& \text{fused multiply-add}\\ \text{rn}\left(x^2\right)& =& 1.0016\times {10}^{-4}\\ \text{rn}\left(\text{rn}\right(x^2)-1)& =& 1.6000\times {10}^{-4}& \text{multiply, then add}\end{array}$

Below is another example, using binary single precision values:

$\begin{array}{rllll}A& =& & 2^0& \times 1.00000000000000000000001\\ B& =& -& 2^0& \times 1.00000000000000000000010\\ \text{rn}(A\times A+B)& =& & 2^{-46}& \times 1.00000000000000000000000\\ \text{rn}\left(\text{rn}\right(A\times A)+B)& =& & 0\end{array}$

In this particular case, computing $\text{rn}\left(\text{rn}\right(A\times A)+B)$ as an IEEE 754 multiply followed by an IEEE 754 add loses all bits of precision, and the computed result is 0. The alternative of computing the FMA $\text{rn}(A\times A+B)$ provides a result equal to the mathematical value. In general, the fused-multiply-add operation generates more accurate results than computing one multiply followed by one add. The choice of whether or not to use the fused operation depends on whether the platform provides the operation and also on how the code is compiled.

Figure 1 shows CUDA C code and output corresponding to inputs A and B and operations from the example above. The code is executed on two different hardware platforms: an x86-class CPU using SSE in single precision, and an NVIDIA GPU with compute capability 2.0. At the time this paper is written (Spring 2011) there are no commercially available x86 CPUs which offer hardware FMA. Because of this, the computed result in single precision in SSE would be 0. NVIDIA GPUs with compute capability 2.0 do offer hardware FMAs, so the result of executing this code will be the more accurate one by default. However, both results are correct according to the IEEE 754 standard. The code fragment was compiled without any special intrinsics or compiler options for either platform.

The fused multiply-add helps avoid loss of precision during subtractive cancellation. Subtractive cancellation occurs during the addition of quantities of similar magnitude with opposite signs. In this case many of the leading bits cancel, leaving fewer meaningful bits of precision in the result. The fused multiply-add computes a double-width product during the multiplication. Thus even if subtractive cancellation occurs during the addition there are still enough valid bits remaining in the product to get a precise result with no loss of precision.

union  {
    float  f;
    unsigned  int  i
} a,  b;
float  r;

a.i = 0x3F800001;
b.i = 0xBF800002;
r = a.f  * a.f  + b.f;

printf("a %.8g\n", a.f); 
printf("b %.8g\n", b.f); 
printf("r %.8g\n",   r);

a:  1.0000001 
b: -1.0000002 
r:  0

a:  1.0000001 
b: -1.0000002
r:  1.4210855e-14

We present three algorithms which differ in how the multiplications, additions, and possibly fused multiply-adds are organized. These algorithms are presented in Figure 2, Figure 3, and Figure 4. Each of the three algorithms is represented graphically. Individual operation are shown as a circle with arrows pointing from arguments to operations.

The simplest way to compute the dot product is using a short loop as shown in Figure 2. The multiplications and additions are done separately.

Figure 2. Serial Method to Compute Vectors Dot Product. The serial method uses a simple loop with separate multiplies and adds to compute the do t product of the vectors. The final result can be represented as ((((a₁ x b₁) + (a₂ x b₂)) + (a₃ x b₃)) + (a₄ x b₄)).

Figure 3. FMA Method to Compute Vector Dot Product. The FMA method uses a simple loop with fused multiply-adds to compute the dot product of the vectors. The final result can be represented as a₄ x b₄ = (a₃ x b₃ + (a₂ x b₂ + (a₁ x b₁ + 0))).

A simple improvement to the algorithm is to use the fused multiply-add to do the multiply and addition in one step to improve accuracy. Figure 3 shows this version.

Yet another way to compute the dot product is to use a divide-and-conquer strategy in which we first find the dot products of the first half and the second half of the vectors, then combine these results using addition. This is a recursive strategy; the base case is the dot product of vectors of length 1 which is a single multiply. Figure 4 graphically illustrates this approach. We call this algorithm the parallel algorithm because the two sub-problems can be computed in parallel as they have no dependencies. The algorithm does not require a parallel implementation, however; it can still be implemented with a single thread.

All three algorithms for computing a dot product use IEEE 754 arithmetic and can be implemented on any system that supports the IEEE standard. In fact, an implementation of the serial algorithm on multiple systems will give exactly the same result. So will implementations of the FMA or parallel algorithms. However, results computed by an implementation of the serial algorithm may differ from those computed by an implementation of the other two algorithms.

Figure 4. The Parallel Method to Reduce Individual Elements Products into a Final Sum. The parallel method uses a tree to reduce all the products of individual elements into a final sum. The final result can be represented as ((a₁ x b₁) + (a₂ x b₂)) + ((a₃ x b₃) + (a₄ x b₄)).

For example, consider the vectors:

• a = [1.907607, -.7862027, 1.148311, .9604002]
• b = [-.9355000, -.6915108, 1.724470, -.7097529]

whose elements are randomly chosen values between -1 and 2. The accuracy of each algorithm corresponding to these inputs is shown in Figure 5.

The main points to notice from the table are that each algorithm yields a different result, and they are all slightly different from the correct mathematical dot product. In this example the FMA version is the most accurate, and the parallel algorithm is more accurate than the serial algorithm. In our experience these results are typical; fused multiply-add significantly increases the accuracy of results, and parallel tree reductions for summation are usually much more accurate than serial summation.

Devices with compute capability 2.0 and above support both single and double precision IEEE 754 including fused multiply-add in both single and double precision. Operations such as square root and division will result in the floating point value closest to the correct mathematical result in both single and double precision, by default.

The IEEE 754 standard defines four rounding modes: round-to-nearest, round towards positive, round towards negative, and round towards zero. CUDA supports all four modes. By default, operations use round-to-nearest. Compiler intrinsics like the ones listed in the tables below can be used to select other rounding modes for individual operations.

In general, the fused multiply-add operation is faster and more accurate than performing separate multiply and add operations. However, on occasion you may wish to disable the merging of multiplies and adds into fused multiply-add instructions. To inhibit this optimization one can write the multiplies and additions using intrinsics with explicit rounding mode as shown in the previous tables. Operations written directly as intrinsics are guaranteed to remain independent and will not be merged into fused multiply-add instructions. With CUDA Fortran it is possible to disable FMA merging via a compiler flag.

Compiler flags relevant to IEEE 754 operations are -ftz={true|false}, -prec-div={true|false}, and -prec-sqrt={true|false}. These flags control single precision operations on devices of compute capability of 2.0 or later.

The default IEEE 754 mode means that single precision operations are correctly rounded and support denormals, as per the IEEE 754 standard. In the fast mode denormal numbers are flushed to zero, and the operations division and square root are not computed to the nearest floating point value. The flags have no effect on double precision or on devices of compute capability below 2.0.

NVIDIA GPUs differ from the x86 architecture in that rounding modes are encoded within each floating point instruction instead of dynamically using a floating point control word. Trap handlers for floating point exceptions are not supported. On the GPU there is no status flag to indicate when calculations have overflowed, underflowed, or have involved inexact arithmetic. Like SSE, the precision of each GPU operation is encoded in the instruction (for x87 the precision is controlled dynamically by the floating point control word).

So far we have only considered simple math operations such as addition, multiplication, division, and square root. These operations are simple enough that computing the best floating point result (e.g., the closest in round-to-nearest) is reasonable. For other mathematical operations computing the best floating point result is harder.

The problem is called the table maker's dilemma. To guarantee the correctly rounded result, it is not generally enough to compute the function to a fixed high accuracy. There might still be rare cases where the error in the high accuracy result affects the rounding step at the lower accuracy.

It is possible to solve the dilemma for particular functions by doing mathematical analysis and formal proofs [4], but most math libraries choose instead to give up the guarantee of correct rounding. Instead they provide implementations of math functions and document bounds on the relative error of the functions over the input range. For example, the double precision sin function in CUDA is guaranteed to be accurate to within 2 units in the last place (ulp) of the correctly rounded result. In other words, the difference between the computed result and the mathematical result is at most ±2 with respect to the least significant bit position of the fraction part of the floating point result.

For most inputs the sin function produces the correctly rounded result. precisions, libraries and hardware. Take for example the C code sequence shown in Figure 6. We compiled the code sequence on a 64-bit x86 platform using gcc version 4.4.3 (Ubuntu 4.3.3-4ubuntu5).

This shows that the result of computing cos(5992555.0) using a common library differs depending on whether the code is compiled in 32-bit mode or 64-bit mode.

The consequence is that different math libraries cannot be expected to compute exactly the same result for a given input. This applies to GPU programming as well. Functions compiled for the GPU will use the NVIDIA CUDA math library implementation while functions compiled for the CPU will use the host compiler math library implementation (e.g., glibc on Linux). Because these implementations are independent and neither is guaranteed to be correctly rounded, the results will often differ slightly.

volatile float x = 5992555.0;
printf("cos(%f): %.10g\n", x, cos(x));

gcc test.c -lm -m64
cos(5992555.000000): 3.320904615e-07
          
gcc test.c -lm -m32
cos(5992555.000000): 3.320904692e-07

One of the unfortunate realities of C compilers is that they are often poor at preserving IEEE 754 semantics of floating point operations [6]. This can be particularly confusing on platforms that support x87 and SSE operations. Just like CUDA operations, SSE operations are performed on single or double precision values, while x87 operations often use an additional internal 80-bit precision format. Sometimes the results of a computation using x87 can depend on whether an intermediate result was allocated to a register or stored to memory. Values stored to memory are rounded to the declared precision (e.g., single precision for float and double precision for double). Values kept in registers can remain in extended precision. Also, x87 instructions will often be used by default for 32-bit compiles but SSE instructions will be used by default for 64-bit compiles.

Because of these issues, guaranteeing a specific precision level on the CPU can sometimes be tricky. When comparing CPU results to results computed on the GPU, it is generally best to compare using SSE instructions. SSE instructions follow IEEE 754 for single and doubleprecision.

On 32-bit x86 targets without SSE it can be helpful to declare variables using volatile and force floating point values to be stored to memory (/Op in Visual Studio and -ffloat-store in gcc). This moves results from extended precision registers into memory, where the precision is precisely single or double precision. Alternately, the x87 control word can be updated to set the precision to 24 or 53 bits using the assembly instruction fldcw or a compiler option such as -mpc32 or-mpc64 in gcc.

As we have shown in Chapter 3, the final values computed using IEEE 754 arithmetic can depend on implementation choices such as whether to use fused multiply-add or whether additions are organized in series or parallel. These differences affect computation on the CPU and on the GPU.

One way such differences can arise is from differences between the number of concurrent threads involved in a computation. On the GPU, a common design pattern is to have all threads in a block coordinate to do a parallel reduction on data within the block, followed by a serial reduction of the results from each block. Changing the number of threads per block reorganizes the reduction; if the reduction is addition, then the change rearranges parentheses in the long string of additions.

Even if the same general strategy such as parallel reduction is used on the CPU and GPU, it is common to have widely different numbers of threads on the GPU compared to the CPU. For example, the GPU implementation might launch blocks with 128 threads per block, while the CPU implementation might use 4 threads in total.

The same inputs will give the same results for individual IEEE 754 operations to a given precision on the CPU and GPU. As we have explained, there are many reasons why the same sequence of operations may not be performed on the CPU and GPU. The GPU has fused multiply-add while the CPU does not. Parallelizing algorithms may rearrange operations, yielding different numeric results. The CPU may be computing results in a precision higher than expected. Finally, many common mathematical functions are not required by the IEEE 754 standard to be correctly rounded so should not be expected to yield identical results between implementations.

When porting numeric code from the CPU to the GPU of course it makes sense to use the x86 CPU results as a reference. But differences between the CPU result and GPU result must be interpreted carefully. Differences are not automatically evidence that the result computed by the GPU is wrong or that there is a problem on the GPU.

Computing results in a high precision and then comparing to results computed in a lower precision can be helpful to see if the lower precision is adequate for a particular application. However, rounding high precision results to a lower precision is not equivalent to performing the entire computation in lower precision. This can sometimes be a problem when using x87 and comparing results against the GPU. The results of the CPU may be computed to an unexpectedly high extended precision for some or all of the operations. The GPU result will be computed using single or double precision only.