SGEMM結果不一致

Question

我正在使用英特爾MKL庫中的sgemm函數來乘以英特爾CPU上的大型矩陣。

我有一個單元測試，它需要一組數據並通過各種算法運行數據。已經證明，如果不使用sgemm（使用非優化算法，而不是有人在我公司寫的），結果是完全相同的。

我們得到的結果與函數返回的矩陣中的最小有效數字不一致。這個錯誤可以通過我們使用的算法類型加以複合。

我切換到DGEMM，並使用雙精度值，而不是單一避免的影響意義。不過，我仍然對可能導致這種不一致的原因感興趣，以及爲什麼乘以矩陣（使用我們自己的未優化算法）不會導致此問題。

我現在的想法是矩陣相乘時，浮點乘法可以亂序執行，因爲這些浮點運算都沒有關聯，我們得到微妙的不同值。

Answer 1

我了興趣，這一點，並編寫了一段代碼來測試假設自己，它似乎是比「標準模式」 SIMD給出了不同的結果。

以下片段與ICC 13.0.2在Mac OS X 10.8.3使用icpc -std=c++11 -O3 -ip -xAVX -fp-model source -fp-model precise -mkl=parallel -openmp編譯。

#include <cmath> 
#include <cstring> 
#include <iostream> 
#include <random> 
#include <utility> 

#include <immintrin.h> 
#include <mkl.h> 

template <typename type, size_t rows, size_t cols> 
class matrix 
{ 
private: 
    void *_data; 

public: 
    matrix() : 
     _data (_mm_malloc(sizeof(type) * rows * cols, 64)) 
    { 
     if (_data == nullptr) throw std::bad_alloc(); 
     else memset(_data, 0, sizeof(type) * rows * cols); 
    } 

    matrix(matrix<type, rows, cols> const& other) : 
     _data (_mm_malloc(sizeof(type) * rows * cols, 64)) 
    { 
     if (_data == nullptr) throw std::bad_alloc(); 
     else memcpy(_data, other._data, sizeof(type) * rows * cols); 
    } 

    ~matrix() 
    { 
     if (_data != nullptr) _mm_free(_data); 
    } 

    typedef type array_type[cols]; 
    array_type& operator[](size_t i) 
    { 
     return static_cast<array_type*>(_data)[i]; 
    } 

    typedef type const_array_type[cols]; 
    const_array_type& operator[](size_t i) const 
    { 
     return static_cast<const_array_type*>(_data)[i]; 
    } 
}; 

template <typename type, size_t m, size_t n> 
type max_diff(matrix<type, m, n> const& a, matrix<type, m, n> const& b) 
{ 
    type value = static_cast<type>(0); 
    for (size_t i = 0; i < m; ++i) 
    { 
     #pragma novector 
     for (size_t j = 0; j < n; ++j) 
     { 
      const type diff = a[i][j] - b[i][j]; 
      if (std::abs(diff) > value) value = std::abs(diff); 
     } 
    } 
    return value; 
} 

template <typename type, size_t m, size_t n, size_t k> 
matrix<type, m, n> matmul_loop(matrix<type, m, k> const& a, matrix<type, n, k> const& b) 
{ 
    matrix<type, m, n> out; 

    #pragma omp parallel for 
    for (size_t i = 0; i < m; ++i) 
    { 
     for (size_t j = 0; j < n; ++j) 
     { 
      for (size_t l = 0; l < k; ++l) 
      { 
       out[i][j] += a[i][l] * b[j][l]; 
      } 
     } 
    } 

    return out; 
} 

template <typename type, size_t m, size_t n, size_t k> 
matrix<type, m, n> matmul_simd(matrix<type, m, k> const& a, matrix<type, n, k> const& b) 
{ 
    matrix<type, m, n> out; 
    type *temp = static_cast<type*>(_mm_malloc(sizeof(type) * k, 64)); 

    #pragma omp parallel for 
    for (size_t i = 0; i < m; ++i) 
    { 
     for (size_t j = 0; j < n; ++j) 
     { 
      type temp = 0.; 

      #pragma vector aligned 
      #pragma ivdep 
      #pragma simd vectorlengthfor(type) 
      for (size_t l = 0; l < k; ++l) 
      { 
       temp += a[i][l] * b[j][l]; 
      } 

      out[i][j] = temp; 
     } 
    } 

    return out; 
} 

template <size_t m, size_t n, size_t k> 
matrix<float, m, n> matmul_sgemm(matrix<float, m, k> const& a, matrix<float, n, k> const& b) 
{ 
    matrix<float, m, n> out; 
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, 1., &a[0][0], m, &b[0][0], n, 0., &out[0][0], m); 
    return out; 
} 

int main() 
{ 
    std::mt19937_64 generator; 
    std::uniform_real_distribution<float> rand_dist(-1000.0,1000.0); 

    const size_t size = 4096; 

    matrix<float, size, size> mat; 
    for (size_t i = 0; i < size; ++i) 
    { 
     for (size_t j = 0; j < size; ++j) 
     { 
      mat[i][j] = rand_dist(generator); 
     } 
    } 

    matrix<float, size, size> result_loop = matmul_loop(mat, mat); 
    matrix<float, size, size> result_simd = matmul_simd(mat, mat); 
    matrix<float, size, size> result_sgemm = matmul_sgemm(mat, mat); 

    std::cout << "SIMD differs from LOOP by a maximum of " << max_diff(result_loop, result_simd) << std::endl; 
    std::cout << "SGEMM differs from LOOP by a maximum of " << max_diff(result_loop, result_sgemm) << std::endl; 
    std::cout << "SGEMM differs from SIMD by a maximum of " << max_diff(result_simd, result_sgemm) << std::endl; 

    return 0; 
} 

請注意，「隨機」矩陣是使用標準種子生成的，因此結果應該是完全可重現的。基本上，給定一個4096x4096矩陣A，代碼使用三種不同的方法計算AA ^T，然後比較結果，打印出相差最大的成分。在我的機器，輸出如下：

$ ./matmul 
SIMD differs from LOOP by a maximum of 6016 
SGEMM differs from LOOP by a maximum of 6016 
SGEMM differs from SIMD by a maximum of 512 

的編譯器標誌-fp-model source -fp-model precise防止matmul_loop被矢量化，但在matmul_simd循環顯然被迫矢量化#pragma simd。矩陣轉置只是爲了簡化SIMD代碼。