c++ - SGEMM 结果不一致

Question

我正在使用英特尔 MKL 库中的 sgemm 函数在英特尔 CPU 上乘以大型矩阵。

我有一个单元测试，它获取一组数据并通过各种算法运行数据。已经证明，在使用这组数据的两次传递之间，如果不使用 sgemm(使用非优化算法而不是我公司某人编写的算法)，结果是完全相同的。

我们得到的结果与函数返回的矩阵中的最低有效数字不一致。然后，我们使用的算法类型可能会加剧此错误。

我通过切换到 dgemm 并使用 double 值而不是单精度值来避免效果的重要性。但是，我仍然对可能导致这种不一致的原因以及为什么乘以矩阵(使用我们自己未优化的算法)不会导致此问题感兴趣。

我目前的想法是，在乘以矩阵时，浮点乘法可能会乱序执行，并且因为这些浮点运算不是关联的，所以我们得到的值略有不同。

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 编译的=并行-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 代码。