使用 MKL 编译时 Eigen C++ 运行 变慢
Eigen C++ running slower when compiled with MKL
我最近开始使用 Eigen(版本 3.3.1),运行在 OLS 回归核心的简单矩阵运算上针对犰狳进行基准测试,即计算乘积的逆对于一个单独的矩阵,我注意到使用 MKL 库编译时 Eigen 比没有它时慢 运行ning 对于这种操作。我想知道我的编译说明是否有误。我还尝试实现直接调用 MKL BLAS 和 LAPACK 例程的此操作,并获得更快的结果,与犰狳一样快。我无法解释如此糟糕的性能,尤其是对于 float 类型。
我写了下面的代码来实现这个基准:
#define ARMA_DONT_USE_WRAPPER
#define ARMA_NO_DEBUG
#include <armadillo>
#define EIGEN_NO_DEBUG
#define EIGEN_NO_STATIC_ASSERT
#define EIGEN_USE_MKL_ALL
#include <Eigen/Dense>
template <typename T>
using Matrix = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
#ifdef USE_FLOAT
using T = float;
#else
using T = double;
#endif
int main()
{
arma::wall_clock timer;
int niter = 1000000;
int n = 1000;
int k = 20;
arma::Mat<T> Xa = arma::cumsum(arma::randn<arma::Mat<T>>(n, k));
Matrix<T> Xe = Matrix<T>::Map(Xa.memptr(), Xa.n_rows, Xa.n_cols);
// Armadillo compiled with MKL
timer.tic();
for (int i = 0; i < niter; ++i) {
arma::Mat<T> iX2a = (Xa.t() * Xa).i();
}
std::cout << "...Elapsed time: " << timer.toc() << "\n";
// Eigen compiled with MKL
timer.tic();
for (int i = 0; i < niter; ++i) {
Matrix<T> iX2e = (Xe.transpose() * Xe).inverse();
}
std::cout << "...Elapsed time: " << timer.toc() << "\n";*/
// Eigen Matrix with MKL routines
timer.tic();
for (int i = 0; i < niter; ++i) {
Matrix<T> iX2e = Matrix<T>::Zero(k, k);
// first stage => computing square matrix trans(X) * X
#ifdef USE_FLOAT
cblas_ssyrk(CblasColMajor, CblasLower, CblasTrans, k, n, 1.0, &Xe(0,0), n, 0.0, &iX2e(0,0), k);
#else
cblas_dsyrk(CblasColMajor, CblasLower, CblasTrans, k, n, 1.0, &Xe(0,0), n, 0.0, &iX2e(0,0), k);
#endif
// getting upper part
for (int i = 0; i < k; ++i)
for (int j = i + 1; j < k; ++j)
iX2e(i, j) = iX2e(j, i);
// second stage => inverting square matrix
// initializing pivots
int* ipiv = new int[k];
// factorizing matrix
#ifdef USE_FLOAT
LAPACKE_sgetrf(LAPACK_COL_MAJOR, k, k, &iX2e(0,0), k, ipiv);
#else
LAPACKE_dgetrf(LAPACK_COL_MAJOR, k, k, &iX2e(0,0), k, ipiv);
#endif
// computing the matrix inverse
#ifdef USE_FLOAT
LAPACKE_sgetri(LAPACK_COL_MAJOR, k, &iX2e(0,0), k, ipiv);
#else
LAPACKE_dgetri(LAPACK_COL_MAJOR, k, &iX2e(0,0), k, ipiv);
#endif
delete[] ipiv;
}
std::cout << "...Elapsed time: " << timer.toc() << "\n";
}
我编译这个名为 test.cpp 的文件:
g++ -std=c++14 -Wall -O3 -march=native -DUSE_FLOAT test.cpp -o 运行 -L${MKLROOT}/lib/intel64 - Wl,--不需要 -lmkl_gf_lp64 -lmkl_sequential -lmkl_core
我得到以下结果(在 Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz 上)
- 双精度类型:
带 MKL 的犰狳 => 64.0s
MKL 的特征值 => 72.2s
仅本征 => 68.7s
纯 MKL => 64.9 秒
- 对于浮动类型:
带 MKL 的犰狳 => 38.2s
MKL 的特征值 => 61.1s
仅本征 => 42.6s
纯 MKL => 38.9s
注意:我运行这个测试是针对一个不会使用非常大矩阵的项目,我不需要在这个级别进行并行化,我最大的矩阵可能是 2000 行 25 列,而且我会需要在更高级别并行,所以我想避免任何类型的嵌套并行性,这可能会降低我的代码速度。
正如我在评论中所说,确保在进行基准测试时禁用涡轮增压。
作为旁注和供将来参考,您当前的 Eigen 代码将调用 gemm 而不是 syrk。您可以明确要求后者:
Matrix<T> tmp = Matrix<T>::Zero(k, k);
tmp.selfadjointView<Eigen::Lower>().rankUpdate(Xe.transpose());
tmp.triangularView<Eigen::Upper>() = tmp.transpose().triangularView<Eigen::Lower>();
iX2e = tmp.inverse();
对于这么小的矩阵,我真的看不出有什么区别。
我只是想补充一点,以防某些人可能对此有疑问,ggael 给出的表达式必须按如下方式编写,以防它是模板的一部分 function/class,否则编译器将难以制作类型推导
Matrix<T> tmp = Matrix<T>::Zero(k, k);
tmp.template selfadjointView<Eigen::Lower>().rankUpdate(Xe.transpose());
tmp.template triangularView<Eigen::Upper>() = tmp.transpose().template triangularView<Eigen::Lower>();
Matrix<T> iX2e = tmp.inverse();
通过此修改并关闭涡轮增压,我得到以下结果:
- 双精度型:
带 MKL 的犰狳 => 79.9s
MKL 的特征值 => 79.8s
仅本征 => 71.1s
纯 MKL => 81.1s
- 对于浮动类型:
带 MKL 的犰狳 => 47.2s
MKL 的特征值 => 50.9s
仅本征 => 51.8s
纯 MKL => 48.0s
我最近开始使用 Eigen(版本 3.3.1),运行在 OLS 回归核心的简单矩阵运算上针对犰狳进行基准测试,即计算乘积的逆对于一个单独的矩阵,我注意到使用 MKL 库编译时 Eigen 比没有它时慢 运行ning 对于这种操作。我想知道我的编译说明是否有误。我还尝试实现直接调用 MKL BLAS 和 LAPACK 例程的此操作,并获得更快的结果,与犰狳一样快。我无法解释如此糟糕的性能,尤其是对于 float 类型。
我写了下面的代码来实现这个基准:
#define ARMA_DONT_USE_WRAPPER
#define ARMA_NO_DEBUG
#include <armadillo>
#define EIGEN_NO_DEBUG
#define EIGEN_NO_STATIC_ASSERT
#define EIGEN_USE_MKL_ALL
#include <Eigen/Dense>
template <typename T>
using Matrix = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
#ifdef USE_FLOAT
using T = float;
#else
using T = double;
#endif
int main()
{
arma::wall_clock timer;
int niter = 1000000;
int n = 1000;
int k = 20;
arma::Mat<T> Xa = arma::cumsum(arma::randn<arma::Mat<T>>(n, k));
Matrix<T> Xe = Matrix<T>::Map(Xa.memptr(), Xa.n_rows, Xa.n_cols);
// Armadillo compiled with MKL
timer.tic();
for (int i = 0; i < niter; ++i) {
arma::Mat<T> iX2a = (Xa.t() * Xa).i();
}
std::cout << "...Elapsed time: " << timer.toc() << "\n";
// Eigen compiled with MKL
timer.tic();
for (int i = 0; i < niter; ++i) {
Matrix<T> iX2e = (Xe.transpose() * Xe).inverse();
}
std::cout << "...Elapsed time: " << timer.toc() << "\n";*/
// Eigen Matrix with MKL routines
timer.tic();
for (int i = 0; i < niter; ++i) {
Matrix<T> iX2e = Matrix<T>::Zero(k, k);
// first stage => computing square matrix trans(X) * X
#ifdef USE_FLOAT
cblas_ssyrk(CblasColMajor, CblasLower, CblasTrans, k, n, 1.0, &Xe(0,0), n, 0.0, &iX2e(0,0), k);
#else
cblas_dsyrk(CblasColMajor, CblasLower, CblasTrans, k, n, 1.0, &Xe(0,0), n, 0.0, &iX2e(0,0), k);
#endif
// getting upper part
for (int i = 0; i < k; ++i)
for (int j = i + 1; j < k; ++j)
iX2e(i, j) = iX2e(j, i);
// second stage => inverting square matrix
// initializing pivots
int* ipiv = new int[k];
// factorizing matrix
#ifdef USE_FLOAT
LAPACKE_sgetrf(LAPACK_COL_MAJOR, k, k, &iX2e(0,0), k, ipiv);
#else
LAPACKE_dgetrf(LAPACK_COL_MAJOR, k, k, &iX2e(0,0), k, ipiv);
#endif
// computing the matrix inverse
#ifdef USE_FLOAT
LAPACKE_sgetri(LAPACK_COL_MAJOR, k, &iX2e(0,0), k, ipiv);
#else
LAPACKE_dgetri(LAPACK_COL_MAJOR, k, &iX2e(0,0), k, ipiv);
#endif
delete[] ipiv;
}
std::cout << "...Elapsed time: " << timer.toc() << "\n";
}
我编译这个名为 test.cpp 的文件:
g++ -std=c++14 -Wall -O3 -march=native -DUSE_FLOAT test.cpp -o 运行 -L${MKLROOT}/lib/intel64 - Wl,--不需要 -lmkl_gf_lp64 -lmkl_sequential -lmkl_core
我得到以下结果(在 Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz 上)
- 双精度类型:
带 MKL 的犰狳 => 64.0s
MKL 的特征值 => 72.2s
仅本征 => 68.7s
纯 MKL => 64.9 秒
- 对于浮动类型:
带 MKL 的犰狳 => 38.2s
MKL 的特征值 => 61.1s
仅本征 => 42.6s
纯 MKL => 38.9s
注意:我运行这个测试是针对一个不会使用非常大矩阵的项目,我不需要在这个级别进行并行化,我最大的矩阵可能是 2000 行 25 列,而且我会需要在更高级别并行,所以我想避免任何类型的嵌套并行性,这可能会降低我的代码速度。
正如我在评论中所说,确保在进行基准测试时禁用涡轮增压。
作为旁注和供将来参考,您当前的 Eigen 代码将调用 gemm 而不是 syrk。您可以明确要求后者:
Matrix<T> tmp = Matrix<T>::Zero(k, k);
tmp.selfadjointView<Eigen::Lower>().rankUpdate(Xe.transpose());
tmp.triangularView<Eigen::Upper>() = tmp.transpose().triangularView<Eigen::Lower>();
iX2e = tmp.inverse();
对于这么小的矩阵,我真的看不出有什么区别。
我只是想补充一点,以防某些人可能对此有疑问,ggael 给出的表达式必须按如下方式编写,以防它是模板的一部分 function/class,否则编译器将难以制作类型推导
Matrix<T> tmp = Matrix<T>::Zero(k, k);
tmp.template selfadjointView<Eigen::Lower>().rankUpdate(Xe.transpose());
tmp.template triangularView<Eigen::Upper>() = tmp.transpose().template triangularView<Eigen::Lower>();
Matrix<T> iX2e = tmp.inverse();
通过此修改并关闭涡轮增压,我得到以下结果:
- 双精度型:
带 MKL 的犰狳 => 79.9s
MKL 的特征值 => 79.8s
仅本征 => 71.1s
纯 MKL => 81.1s
- 对于浮动类型:
带 MKL 的犰狳 => 47.2s
MKL 的特征值 => 50.9s
仅本征 => 51.8s
纯 MKL => 48.0s