Are BLAS implementations guaranteed to give the exact same result?

Question

Given two different BLAS implementations, can we expect that they make the exact same floating point computations and return the same results? Or can it happen, for instance, that one computes a scalar product as $$ ((x_1y_1 + x_2y_2) + x_3y_3) + x_4y_4 $$ and one as $$ (x_1y_1 + x_2y_2) + (x_3y_3 + x_4y_4), $$ so possibly giving different result in IEEE floating point arithmetic?

Answer 1

No, that is not guaranteed. If you are using a NETLIB BLAS without any optimizations, it it mostly true that the results are the same. But for any practical usage of BLAS and LAPACK one uses a highly optimized an parallel BLAS. The parallelization causes, even if it only works in parallel inside the vector registers of a CPU, that the order how the single terms are evaluated changes and the order of the summation changes too. Now it follows form the missing associative property in the IEEE standard that the results are not the same. So exactly the thing you mentioned can happen.

In the NETLIB BLAS the scalar product is only a for loop unrolled by a factor 5:

DO I = MP1,N,5
          DTEMP = DTEMP + DX(I)*DY(I) + DX(I+1)*DY(I+1) +
     $            DX(I+2)*DY(I+2) + DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
END DO

and it is up to the compiler if each multiplication is added to DTEMP immediately or if all 5 components are summed up first and than added to DTEMP. In OpenBLAS it is depending on the architecture a more complicated kernel:

 __asm__  __volatile__
    (
    "vxorpd     %%ymm4, %%ymm4, %%ymm4               \n\t"
    "vxorpd     %%ymm5, %%ymm5, %%ymm5               \n\t"
    "vxorpd     %%ymm6, %%ymm6, %%ymm6               \n\t"
    "vxorpd     %%ymm7, %%ymm7, %%ymm7               \n\t"

    ".align 16                           \n\t"
    "1:                          \n\t"
        "vmovups                  (%2,%0,8), %%ymm12         \n\t"  // 2 * x
        "vmovups                32(%2,%0,8), %%ymm13         \n\t"  // 2 * x
        "vmovups                64(%2,%0,8), %%ymm14         \n\t"  // 2 * x
        "vmovups                96(%2,%0,8), %%ymm15         \n\t"  // 2 * x

    "vmulpd      (%3,%0,8), %%ymm12, %%ymm12 \n\t"  // 2 * y
    "vmulpd    32(%3,%0,8), %%ymm13, %%ymm13 \n\t"  // 2 * y
    "vmulpd    64(%3,%0,8), %%ymm14, %%ymm14 \n\t"  // 2 * y
    "vmulpd    96(%3,%0,8), %%ymm15, %%ymm15 \n\t"  // 2 * y

    "vaddpd      %%ymm4 , %%ymm12, %%ymm4 \n\t"  // 2 * y
    "vaddpd      %%ymm5 , %%ymm13, %%ymm5 \n\t"  // 2 * y
    "vaddpd      %%ymm6 , %%ymm14, %%ymm6 \n\t"  // 2 * y
    "vaddpd      %%ymm7 , %%ymm15, %%ymm7 \n\t"  // 2 * y

    "addq       $16 , %0	  	     \n\t"
	"subq	        $16 , %1            \n\t"      
    "jnz        1b                   \n\t"
...

which splits the scalar product in small scalar products of length 4 and sum them up.

Using the other typical BLAS implementations like ATLAS, MKL, ESSL,... this problem stays the same because each BLAS implementation uses different optimizations to get fast code. But as far as I know one need an artificial example to cause really faulty results.

If it is necessary that the BLAS library returns for the same results ( bit-wise the same) one have to use a reproducible BLAS library such as:

• ReproBLAS http://bebop.cs.berkeley.edu/reproblas/
• exBLAS https://exblas.lip6.fr/

Answer 2

The Short Answer

If the two BLAS implementations are written to carry out the operations in the exact same order, and the libraries were compiled using the same compiler flags and with the same compiler, then they'll give you the same result. Floating point arithmetic is not random, so two identical implementations will give identical results.

However, there are a variety of things that can break this behavior for the sake of performance...

The Longer Answer

IEEE also specifies the order in which these operations are carried out, in addition to how each operation should behave. However, if you compile your BLAS implementation with options like "-ffast-math", the compiler can perform transformations that would be true in exact arithmetic but not "correct" in IEEE floating point. The canonical example is the non-associativity of floating point addition, as you pointed out. With the more aggressive optimization settings, associativity will be assumed, and the processor will do as much of that in parallel as possible by re-ordering the operations.

The other standard-breaking behavior comes via the use of FMA (fused multiply-add) instructions. These are prominent in operations like matrix multiplication, and they have the potential to double the throughput of your routine. However, they perform the operation $a+b*c$ in a single operation, and it only incurs a single floating point rounding step. This deviates from the IEEE standard, which requires that this operation have two rounding steps. This makes the FMA result actually more accurate than the IEEE one, but it is technically standard-breaking behavior.

Answer 3

In general, no. Leaving associativity aside, the choice of compiler flags (for example, SIMD instructions being enabled, usage of fused multiply add, etc.) or the hardware (e.g., whether extended precision is being used) may produce different results.

There are some efforts to get reproducible BLAS implementations. See ReproBLAS and ExBLAS for more information.