Apakah implementasi BLAS dijamin untuk memberikan hasil yang sama persis?

Dengan dua implementasi BLAS yang berbeda, dapatkah kita berharap mereka membuat perhitungan floating point yang sama persis dan mengembalikan hasil yang sama? Atau bisa terjadi, misalnya, bahwa salah satu Toedjoe menghitung produk skalar sebagai dan satu sebagai ( x 1 y 1 + x 2 y 2 ) + ( x 3 y 3 + x

Tidak, itu tidak dijamin. Jika Anda menggunakan NETLIB BLAS tanpa optimasi apa pun, sebagian besar hasilnya benar-benar sama. Tetapi untuk penggunaan praktis BLAS dan LAPACK, orang menggunakan BLAS paralel yang sangat dioptimalkan. Paralelisasi menyebabkan, bahkan jika itu hanya bekerja secara paralel di dalam register vektor CPU, bahwa urutan bagaimana istilah tunggal dievaluasi berubah dan urutan penjumlahan juga berubah. Sekarang mengikuti bentuk properti asosiatif yang hilang dalam standar IEEE bahwa hasilnya tidak sama. Jadi persis hal yang Anda sebutkan bisa terjadi.

Dalam NETLIB BLAS produk skalar hanya untuk loop yang tidak digerakkan oleh faktor 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

dan terserah kompiler jika setiap perkalian ditambahkan ke DTEMP segera atau jika semua 5 komponen dirangkum terlebih dahulu dan kemudian ditambahkan ke DTEMP. Di OpenBLAS tergantung pada arsitektur kernel yang lebih rumit:

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

yang membagi produk skalar menjadi produk skalar kecil dengan panjang 4 dan menjumlahkannya.

Menggunakan implementasi BLAS khas lainnya seperti ATLAS, MKL, ESSL, ... masalah ini tetap sama karena setiap implementasi BLAS menggunakan optimisasi yang berbeda untuk mendapatkan kode cepat. Tapi sejauh yang saya tahu orang perlu contoh buatan untuk menyebabkan hasil yang benar-benar salah.

Jika perlu bahwa perpustakaan BLAS kembali untuk hasil yang sama (agak bijaksana sama) kita harus menggunakan perpustakaan BLAS yang dapat direproduksi seperti:

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

Jawaban Singkat

Jika dua implementasi BLAS ditulis untuk menjalankan operasi dalam urutan yang sama persis, dan perpustakaan dikompilasi menggunakan flag compiler yang sama dan dengan compiler yang sama, maka mereka akan memberi Anda hasil yang sama. Aritmatika floating point tidak acak, sehingga dua implementasi yang identik akan memberikan hasil yang identik.

Namun, ada berbagai hal yang dapat mematahkan perilaku ini demi kinerja ...

Jawaban yang Lebih Panjang

IEEE juga menentukan urutan di mana operasi ini dilakukan, di samping bagaimana masing-masing operasi harus berperilaku. Namun, jika Anda mengkompilasi implementasi BLAS Anda dengan opsi-opsi seperti "-ast-matematika", kompiler dapat melakukan transformasi yang akan benar dalam aritmatika yang tepat tetapi tidak "benar" di IEEE floating point. Contoh kanonik adalah non-associativity dari penambahan floating point, seperti yang Anda tunjukkan. Dengan pengaturan optimisasi yang lebih agresif, asosiasi akan diasumsikan, dan prosesor akan melakukan sebanyak mungkin secara paralel dengan memesan ulang operasi.

$a+b*c$

Secara umum, tidak. Mengesampingkan asosiatif, pilihan flag compiler (misalnya, instruksi SIMD diaktifkan, penggunaan tambah banyak pengganda , dll.) Atau perangkat keras (mis., Apakah presisi yang diperluas digunakan) dapat menghasilkan hasil yang berbeda.

Ada beberapa upaya untuk mendapatkan implementasi BLAS yang dapat direproduksi. Lihat ReproBLAS dan ExBLAS untuk informasi lebih lanjut.