Непоследовательность результата SGEMM

Question

Я использую функцию sgemm в библиотеке Intel MKL для умножения больших матриц на процессоры Intel.

У меня есть единичный тест, который берет набор данных и управляет данными с помощью различных алгоритмов. Было доказано, что между проходами с этим набором данных, если 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; 
} 

Обратите внимание, что «случайная» матрица была сгенерирована с использованием стандартного семени, чтобы результат был полностью воспроизводимым. В принципе, учитывая матрицу 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 от того vectorised, но петля в matmul_simd явно вынужден быть vectorised #pragma simd. Матрица транспонирования - это просто, чтобы немного упростить код SIMD.