2

Intel MKL ライブラリの sgemm 関数を使用して、Intel CPU で大きな行列を乗算しています。

一連のデータを取得し、さまざまなアルゴリズムでデータを実行する単体テストがあります。sgemm が使用されていない場合 (私の会社で誰かが書いた最適化されていないアルゴリズムが使用されている場合)、このデータ セットのパス間で結果が完全に同一であることが証明されています。

関数によって返された行列の最下位桁で一貫性のない結果が得られています。このエラーは、使用するアルゴリズムの種類によってさらに悪化する可能性があります。

dgemm に切り替え、単精度ではなく倍精度の値を使用することで、この効果の重要性を回避しました。ただし、この不一致の原因と、(最適化されていない独自のアルゴリズムを使用して) 行列を乗算してもこの問題が発生しない理由については、まだ興味があります。

私の現在の考えでは、行列を乗算する場合、浮動小数点乗算は順不同で実行される可能性があり、これらの浮動小数点演算は連想的ではないため、微妙に異なる値が得られます。

4

1 に答える 1

2

私はこれに興味を持ち、仮説を自分でテストするためのコードを少し書きましたが、SIMD は「標準モード」とは異なる結果をもたらすようです。

次のスニペットは、Mac OS X 10.8.3 で ICC 13.0.2 を使用してコンパイルされました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 が与えられた場合、コードは3 つの異なる方法を使用して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 コードを少し単純化するためにあります。

于 2013-06-25T00:46:53.877 に答える