Несоответствие результатов SGEMM
Я использую функцию sgemm в библиотеке Intel MKL для умножения больших матриц на процессорах Intel.
У меня есть модульный тест, который принимает набор данных и запускает данные с помощью различных алгоритмов. Было доказано, что между проходами с этим набором данных, если sgemm не используется (вместо этого используется неоптимизированный алгоритм, который кто-то написал в моей компании), результаты полностью идентичны.
Мы получаем противоречивые результаты с наименьшей значащей цифрой в матрицах, возвращаемых функцией. Эта ошибка может быть затем усугублена типом алгоритмов, которые мы используем.
Я избежал значимости эффекта, переключившись на dgemm и используя значения двойной точности, а не одиночные. Однако я все еще интересуюсь тем, что может вызвать это несоответствие и почему умножение матрицы (с нашим собственным неоптимизированным алгоритмом) не вызывает этой проблемы.
Мои нынешние мысли заключаются в том, что при умножении матриц умножение с плавающей запятой может выполняться не по порядку, и, поскольку эти операции с плавающей запятой не ассоциативны, мы получаем слегка разные значения.
Решение
Я заинтересовался этим и сам написал немного кода для проверки гипотезы, и кажется, что 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, код вычисляет AAT используя три разных метода, а затем сравнивайте результаты, распечатывая компонент, который отличается наибольшим количеством. На моей машине вывод выглядит следующим образом:
$ ./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.
Другие решения
Других решений пока нет …