As implementações do BLAS são garantidas para fornecer exatamente o mesmo resultado?

Dadas duas implementações BLAS diferentes, podemos esperar que elas façam exatamente os mesmos cálculos de ponto flutuante e retornem os mesmos resultados? Ou pode acontecer, por exemplo, que alguém calcule um produto escalar como e um como ( x 1 y 1 + x 2 y 2 ) + ( x 3 y 3 + x

Não, isso não é garantido. Se você estiver usando um NETLIB BLAS sem nenhuma otimização, é verdade que os resultados são os mesmos. Porém, para qualquer uso prático do BLAS e LAPACK, utiliza-se um BLAS paralelo altamente otimizado. A paralelização faz com que, mesmo que funcione paralelamente nos registros vetoriais de uma CPU, a ordem de avaliação dos termos únicos seja alterada e a ordem da soma também. Agora segue a propriedade associativa ausente no padrão IEEE que os resultados não são os mesmos. Então, exatamente o que você mencionou pode acontecer.

No NETLIB BLAS, o produto escalar é apenas um loop for desenrolado por um fator 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

e cabe ao compilador se cada multiplicação for adicionada ao DTEMP imediatamente ou se todos os 5 componentes forem somados primeiro e depois adicionados ao DTEMP. No OpenBLAS, dependendo da arquitetura, um kernel mais 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 o produto escalar em pequenos produtos escalares de comprimento 4 e os resume.

Usando outras implementações típicas do BLAS, como ATLAS, MKL, ESSL, ... esse problema permanece o mesmo, porque cada implementação do BLAS usa otimizações diferentes para obter código rápido. Mas, tanto quanto eu sei, é preciso um exemplo artificial para causar resultados realmente defeituosos.

Se for necessário que a biblioteca BLAS retorne para os mesmos resultados (bit a bit o mesmo), é necessário usar uma biblioteca BLAS reproduzível, como:

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

A resposta curta

Se as duas implementações do BLAS forem escritas para executar as operações exatamente na mesma ordem e as bibliotecas foram compiladas usando os mesmos sinalizadores do compilador e com o mesmo compilador, eles fornecerão o mesmo resultado. A aritmética de ponto flutuante não é aleatória, portanto, duas implementações idênticas fornecerão resultados idênticos.

No entanto, há uma variedade de coisas que podem quebrar esse comportamento em prol do desempenho ...

A resposta mais longa

O IEEE também especifica a ordem em que essas operações são realizadas, além de como cada operação deve se comportar. No entanto, se você compilar sua implementação do BLAS com opções como "-ffast-math", o compilador poderá executar transformações que seriam verdadeiras na aritmética exata, mas não "corretas" no ponto flutuante IEEE. O exemplo canônico é a não associatividade da adição de ponto flutuante, como você apontou. Com as configurações de otimização mais agressivas, a associatividade será assumida e o processador fará o máximo possível em paralelo ao reordenar as operações.

$a+b*c$

Em geral, não. Deixando de lado a associatividade, a escolha de sinalizadores do compilador (por exemplo, instruções SIMD ativadas, uso de adição de multiplicação por fusão etc.) ou o hardware (por exemplo, se a precisão estendida está sendo usada) pode produzir resultados diferentes.

Existem alguns esforços para obter implementações reproduzíveis do BLAS. Consulte ReproBLAS e ExBLAS para obter mais informações.