Czy wdrożenia BLAS gwarantują dokładnie ten sam wynik?

Biorąc pod uwagę dwie różne implementacje BLAS, czy możemy oczekiwać, że wykonają dokładnie takie same obliczenia zmiennoprzecinkowe i zwrócą te same wyniki? Lub może się zdarzyć, na przykład, że oblicza się produkt skalarny jako i jeden jako co może dać inny wynik w zmiennoprzecinkowym IEEE arytmetyka?

((x_{1} y_{1} + x_{2)} y_{2)}) + x_{3)} y_{3)}) + x_{4} y_{4}

$((x_1y_1 + x_2y_2) + x_3y_3) + x_4y_4$

(x_{1} y_{1} + x_{2)} y_{2)}) + (x_{3)} y_{3)} + x_{4} y_{4}),

$(x_1y_1 + x_2y_2) + (x_3y_3 + x_4y_4),$

floating-point blas

— Federico Poloni
źródło

W tym wątku pojawiają się skargi dotyczące jakości BLAS , poszukaj CBLAS na stronie. To sugerowałoby, że nie tylko nie dają tego samego rezultatu, ale nie wszystkie są wystarczająco dokładne do każdego zadania ...

— Szabolcs

Odpowiedzi:

Nie, to nie jest gwarantowane. Jeśli używasz NETLIB BLAS bez żadnych optymalizacji, to w większości przypadków prawdą jest, że wyniki są takie same. Ale do każdego praktycznego zastosowania BLAS i LAPACK stosuje się wysoce zoptymalizowany równoległy BLAS. Równoległość powoduje, nawet jeśli działa ona równolegle tylko w rejestrach wektorowych CPU, że zmienia się kolejność, w jakiej oceniane są poszczególne warunki, a także zmienia się kolejność sumowania. Teraz z brakującej właściwości asocjacyjnej w standardzie IEEE wynika, że wyniki nie są takie same. Może się zdarzyć dokładnie to, o czym wspomniałeś.

W NETLIB BLAS iloczyn skalarny jest tylko pętlą for rozwiniętą 5-krotnie:

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

i to od kompilatora zależy, czy każde zwielokrotnienie zostanie dodane do DTEMP natychmiast, czy wszystkie 5 składników zostanie najpierw zsumowane, a następnie dodane do DTEMP. W OpenBLAS jest zależne od architektury bardziej skomplikowane jądro:

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

który dzieli iloczyn skalarny na małe iloczyny skalarne o długości 4 i sumuje je.

Używając innych typowych implementacji BLAS, takich jak ATLAS, MKL, ESSL, ... ten problem pozostaje ten sam, ponieważ każda implementacja BLAS korzysta z różnych optymalizacji w celu uzyskania szybkiego kodu. Ale o ile mi wiadomo, potrzebny jest sztuczny przykład, aby spowodować naprawdę błędne wyniki.

Jeśli konieczne jest, aby biblioteka BLAS zwróciła te same wyniki (pod względem bitów to samo), należy użyć odtwarzalnej biblioteki BLAS, takiej jak:

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

— MK alias Grisu
źródło

Krótka odpowiedź

Jeśli dwie implementacje BLAS zostały napisane w celu przeprowadzenia operacji w dokładnie tej samej kolejności, a biblioteki zostały skompilowane przy użyciu tych samych flag kompilatora i przy użyciu tego samego kompilatora, to dają ten sam wynik. Arytmetyka zmiennoprzecinkowa nie jest przypadkowa, więc dwie identyczne implementacje dadzą identyczne wyniki.

Istnieje jednak wiele rzeczy, które mogą złamać to zachowanie ze względu na wydajność ...

Dłuższa odpowiedź

IEEE określa także kolejność wykonywania tych operacji, oprócz tego, jak powinna się zachowywać każda operacja. Jeśli jednak skompilujesz swoją implementację BLAS z opcjami takimi jak „-ffast-matematyka”, kompilator może wykonać transformacje, które byłyby prawdziwe w dokładnej arytmetyce, ale nie byłyby „poprawne” w zmiennoprzecinkowym IEEE. Jak zauważyłeś, kanonicznym przykładem jest brak asocjatywności dodawania zmiennoprzecinkowego. Przy bardziej agresywnych ustawieniach optymalizacji założona zostanie asocjatywność, a procesor wykona tyle czynności równolegle, jak to możliwe, poprzez ponowne uporządkowanie operacji.

$a+b*c$

— Tyler Olsen
źródło

„Arytmetyka zmiennoprzecinkowa nie jest przypadkowa” . To smutne, że trzeba to jednoznacznie stwierdzić, ale wydaje się, że zbyt wielu ludzi uważa, że to ...

— rura

Cóż, nieprzewidywalne i losowe wyglądają bardzo podobnie, a wiele klas programowania wstępnego mówi „nigdy nie porównuj liczb zmiennoprzecinkowych dla równości”, co sprawia wrażenie, że nie można polegać na takiej samej wartości jak liczby całkowite.

— Tyler Olsen,

@TylerOlsen Nie ma to znaczenia dla pytania i nie dlatego te klasy mówią takie rzeczy, ale IIRC istniała klasa błędów kompilatora, w których nie można było polegać na równości. Na przykład if (x == 0) assert(x == 0) może czasami zawieść, co z pewnego punktu widzenia jest tak dobre, jak losowe.

— Kirill,

@Kirill Niestety, po prostu starałem się podkreślić, że wiele osób tak naprawdę nigdy nie uczy się, jak działa zmiennoprzecinkowy. Jeśli chodzi o punkt historyczny, to trochę przerażające, ale musiało zostać rozwiązane, zanim wszedłem do gry.

— Tyler Olsen,

@TylerOlsen Ups, mój przykład jest zły. Powinno tak być z if (x != 0) assert(x != 0)powodu arytmetyki o rozszerzonej precyzji.

— Kirill,

Ogólnie nie. Pomijając skojarzenia, wybór flag kompilatora (na przykład włączenie instrukcji SIMD, użycie stopionego dodawania wielokrotnego itd.) Lub sprzętu (np. Czy używana jest zwiększona precyzja ) może dawać różne wyniki.

Istnieją pewne wysiłki, aby uzyskać odtwarzalne implementacje BLAS. Aby uzyskać więcej informacji, zobacz ReproBLAS i ExBLAS .

— Juan M. Bello-Rivas
źródło

Zobacz także funkcję warunkowej odtwarzalności numerycznej (CNR) w najnowszych wersjach Intel MKL BLAS. Uświadom sobie, że uzyskanie tego poziomu odtwarzalności zwykle spowalnia obliczenia i może je znacznie spowolnić!

— Brian Borchers,