Les implémentations BLAS garantissent-elles exactement le même résultat?

Étant donné deux implémentations BLAS différentes, pouvons-nous nous attendre à ce qu'elles effectuent exactement les mêmes calculs en virgule flottante et renvoient les mêmes résultats? Ou peut-il arriver, par exemple, que l'on calcule un produit scalaire comme et un comme donnant ainsi éventuellement un résultat différent en virgule flottante IEEE arithmétique?

Non, ce n'est pas garanti. Si vous utilisez un NETLIB BLAS sans aucune optimisation, il est généralement vrai que les résultats sont les mêmes. Mais pour toute utilisation pratique de BLAS et LAPACK, on ​​utilise un BLAS parallèle hautement optimisé. La parallélisation provoque, même si elle ne fonctionne qu'en parallèle à l'intérieur des registres vectoriels d'un CPU, que l'ordre de façon dont les termes simples sont évalués change et que l'ordre de la sommation change également. Maintenant, il résulte de la propriété associative manquante dans la norme IEEE que les résultats ne sont pas les mêmes. Donc, exactement ce que vous avez mentionné peut arriver.

Dans le NETLIB BLAS, le produit scalaire n'est qu'une boucle for déroulée par un facteur 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

et il appartient au compilateur si chaque multiplication est ajoutée à DTEMP immédiatement ou si les 5 composants sont additionnés en premier et ensuite ajoutés à DTEMP. Dans OpenBLAS, cela dépend de l'architecture d'un noyau plus compliqué:

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

qui divise le produit scalaire en petits produits scalaires de longueur 4 et les résume.

En utilisant les autres implémentations BLAS typiques comme ATLAS, MKL, ESSL, ... ce problème reste le même car chaque implémentation BLAS utilise différentes optimisations pour obtenir du code rapide. Mais pour autant que je sache, il faut un exemple artificiel pour provoquer des résultats vraiment défectueux.

S'il est nécessaire que la bibliothèque BLAS renvoie pour les mêmes résultats (au niveau du bit la même), il faut utiliser une bibliothèque BLAS reproductible telle que:

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

La réponse courte

Si les deux implémentations BLAS sont écrites pour exécuter les opérations dans le même ordre exact, et que les bibliothèques ont été compilées en utilisant les mêmes drapeaux de compilateur et avec le même compilateur, alors elles vous donneront le même résultat. L'arithmétique en virgule flottante n'est pas aléatoire, donc deux implémentations identiques donneront des résultats identiques.

Cependant, il existe une variété de choses qui peuvent briser ce comportement pour des raisons de performances ...

La réponse plus longue

L'IEEE spécifie également l' ordre dans lequel ces opérations sont effectuées, en plus du comportement de chaque opération. Cependant, si vous compilez votre implémentation BLAS avec des options comme "-ffast-math", le compilateur peut effectuer des transformations qui seraient vraies en arithmétique exacte mais pas "correctes" en virgule flottante IEEE. L'exemple canonique est la non-associativité de l'addition en virgule flottante, comme vous l'avez souligné. Avec les paramètres d'optimisation les plus agressifs, l'associativité sera supposée, et le processeur fera autant de cela en parallèle que possible en réorganisant les opérations.

$a+b*c$

En général, non. Laissant de côté l'associativité, le choix des drapeaux du compilateur (par exemple, les instructions SIMD étant activées, l'utilisation de la fusion multipliée , etc.) ou le matériel (par exemple, si une précision étendue est utilisée) peut produire des résultats différents.

Il y a quelques efforts pour obtenir des implémentations BLAS reproductibles. Voir ReproBLAS et ExBLAS pour plus d'informations.