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

17

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 ) ,

((X1y1+X2y2)+X3y3)+X4 4y4 4
(X1y1+X2y2)+(X3y3+X4 4y4 4),
Federico Poloni
fuente
1
Hay algunas quejas sobre la calidad de BLAS en este hilo , busque CBLAS en la página. Eso sugeriría que no solo no dan el mismo resultado, sino que no todos son lo suficientemente precisos para cualquier tarea ...
Szabolcs

Respuestas:

15

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:

MK aka Grisu
fuente
8

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.

un+siC

Tyler Olsen
fuente
3
"La aritmética de coma flotante no es aleatoria" . Es triste que esto tiene que ser declarado explícitamente, pero parece que mucha gente piensa que es ...
tubería
Bueno, impredecible y aleatorio se ven bastante similares, y muchas clases de programación introductoria dicen "nunca compares flotantes por igualdad", lo que da la impresión de que no se puede confiar en el valor exacto de la misma manera que los enteros.
Tyler Olsen
@TylerOlsen Esto no es relevante para la pregunta, y no es por eso que esas clases dicen esas cosas, pero IIRC, solía haber una clase de errores de compilación en los que no se podía confiar en la igualdad. Por ejemplo, a veces if (x == 0) assert(x == 0) puede fallar, lo que desde cierto punto de vista es tan bueno como aleatorio.
Kirill
@Kirill Lo siento, simplemente estaba tratando de señalar que muchas personas nunca aprenden realmente cómo funciona el punto flotante. En cuanto al punto histórico, eso es un poco aterrador, pero debe haber sido resuelto antes de entrar en el juego.
Tyler Olsen
@TylerOlsen Vaya, mi ejemplo está mal. Debería serlo if (x != 0) assert(x != 0), debido a la aritmética de precisión extendida.
Kirill
4

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.

Juan M. Bello-Rivas
fuente
1
Consulte también la función de Reproducibilidad Numérica Condicional (CNR) en versiones recientes de MKL BLAS de Intel. ¡Tenga en cuenta que obtener este nivel de reproducibilidad generalmente ralentizará sus cálculos y puede ralentizarlos mucho!
Brian Borchers