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?

  • 1
    $\begingroup$ There are some complaints about BLAS quality in this thread, search for CBLAS in the page. That would suggest that not only do they not give the same result, but not all of them are accurate enough for any task ... $\endgroup$ – Szabolcs Feb 9 '17 at 20:11

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)

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:

| cite | improve this answer | |

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.

| cite | improve this answer | |
  • 3
    $\begingroup$ "Floating point arithmetic is not random". It's sad that this has to be explicitly stated, but it seems that too many people think it is... $\endgroup$ – pipe Feb 9 '17 at 16:23
  • $\begingroup$ Well, unpredictable and random look pretty similar, and a lot of intro programming classes say "never compare floats for equality," which gives the impression that the exact value can't be relied upon in the same way as integers. $\endgroup$ – Tyler Olsen Feb 9 '17 at 16:42
  • $\begingroup$ @TylerOlsen This isn't relevant to the question, and this isn't why those classes say such things, but IIRC, there used to be a class of compiler bugs where equality couldn't be relied on. E.g., if (x == 0) assert(x == 0) might sometimes fail, which from a certain point of view is as good as random. $\endgroup$ – Kirill Feb 9 '17 at 17:45
  • $\begingroup$ @Kirill Sorry, I was simply trying to make the point that many people never really learn how floating point works. As for the historical point, that's kind of terrifying, but it must have been resolved before I got into the game. $\endgroup$ – Tyler Olsen Feb 9 '17 at 17:48
  • $\begingroup$ @TylerOlsen Oops, my example is wrong. It should be if (x != 0) assert(x != 0), because of the extended-precision arithmetic. $\endgroup$ – Kirill Feb 9 '17 at 18:10

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.

| cite | improve this answer | |
  • 1
    $\begingroup$ See also the Conditional Numerical Reproducibility (CNR) feature in recent versions of Intel's MKL BLAS. Realize that getting this level of reproducibility will typically slow down your computations and may slow them down a lot! $\endgroup$ – Brian Borchers Feb 10 '17 at 21:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.