¿Se garantiza que las implementaciones de BLAS darán exactamente el mismo resultado?

Dadas dos implementaciones de BLAS diferentes, ¿podemos esperar que hagan exactamente los mismos cálculos de coma flotante y devuelvan los mismos resultados? ¿O puede suceder, por ejemplo, que uno calcule un producto escalar como y uno como por lo que posiblemente da un resultado diferente en el punto flotante IEEE ¿aritmética? ( x 1 y 1 + x 2 y 2 ) + ( x 3 y 3 + x 4 y 4 )

No, eso no está garantizado. Si está utilizando un NETLIB BLAS sin optimizaciones, es cierto que los resultados son los mismos. Pero para cualquier uso práctico de BLAS y LAPACK se usa un BLAS paralelo altamente optimizado. La paralelización provoca, incluso si solo funciona en paralelo dentro de los registros vectoriales de una CPU, que el orden en que se evalúan los términos individuales cambia y el orden de la suma también cambia. Ahora se deduce de la propiedad asociativa que falta en el estándar IEEE que los resultados no son los mismos. Entonces, exactamente lo que mencionaste puede suceder.

En NETLIB BLAS, el producto escalar es solo un bucle for desenrollado por un 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

y depende del compilador si cada multiplicación se agrega a DTEMP inmediatamente o si los 5 componentes se suman primero y luego se agregan a DTEMP. En OpenBLAS depende de la arquitectura un núcleo más complicado:

 __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"
...

que divide el producto escalar en pequeños productos escalares de longitud 4 y los resume.

Usando las otras implementaciones típicas de BLAS como ATLAS, MKL, ESSL, ... este problema permanece igual porque cada implementación de BLAS usa diferentes optimizaciones para obtener un código rápido. Pero hasta donde yo sé, uno necesita un ejemplo artificial para causar resultados realmente defectuosos.

Si es necesario que la biblioteca BLAS regrese para obtener los mismos resultados (en lo que respecta a los bits), se debe usar una biblioteca BLAS reproducible como:

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

La respuesta corta

Si las dos implementaciones de BLAS se escriben para llevar a cabo las operaciones en el mismo orden exacto, y las bibliotecas se compilan usando los mismos indicadores de compilación y con el mismo compilador, entonces le darán el mismo resultado. La aritmética de coma flotante no es aleatoria, por lo que dos implementaciones idénticas darán resultados idénticos.

Sin embargo, hay una variedad de cosas que pueden romper este comportamiento en aras del rendimiento ...

La respuesta más larga

IEEE también especifica el orden en que se llevan a cabo estas operaciones, además de cómo debe comportarse cada operación. Sin embargo, si compila su implementación BLAS con opciones como "-ffast-math", el compilador puede realizar transformaciones que serían verdaderas en aritmética exacta pero no "correctas" en coma flotante IEEE. El ejemplo canónico es la no asociatividad de la adición de punto flotante, como usted señaló. Con las configuraciones de optimización más agresivas, se asumirá la asociatividad, y el procesador hará tanto de eso en paralelo como sea posible al reordenar las operaciones.

$a+b*c$

En general, no. Dejando a un lado la asociatividad, la elección de los indicadores del compilador (por ejemplo, las instrucciones SIMD habilitadas, el uso de suma múltiple con fusibles , etc.) o el hardware (por ejemplo, si se está utilizando una precisión extendida ) puede producir resultados diferentes.

Hay algunos esfuerzos para obtener implementaciones BLAS reproducibles. Ver ReproBLAS y ExBLAS para más información.