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?
Are BLAS implementations guaranteed to give the exact same result?
floating point;blas;ieee
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 DOand 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/