为什么带有 MKL 的 Eigen C++ 不为这种大型矩阵乘法使用多线程？

Question

我正在做一些计算，包括大量（每次执行约 40000）具有复杂双元素（从随机分布中选择）的 4x4 矩阵的 QR 分解。我开始使用英特尔 MKL 函数直接编写代码。但经过一些研究，似乎使用 Eigen 会更简单，并且更容易维护代码。（部分原因是我发现很难在英特尔 MKL 中使用二维数组，并且需要注意内存管理）。

在转向 Eigen 之前，我开始进行一些性能检查。我采用了一个代码（来自早期关于 SO 的类似问题），用于将 10000x100000 矩阵与另一个 100000x1000 矩阵（选择大尺寸以产生并行化效果）相乘。我运行它在一个 36 核节点上。当我检查统计数据时，没有英特尔 MKL 指令的 Eigen（但使用 -O3 -fopenmp 编译）使用了所有内核并在 ~7 秒内完成了任务。

另一方面，

#define EIGEN_USE_MKL_ALL
#define EIGEN_VECTORIZE_SSE4_2

该代码耗时 28 秒且仅使用一个内核。这是我的编译说明

g++ -m64 -std=c++17 -fPIC -c -I. -I/apps/IntelParallelStudio/mkl/include -O2 -DNDEBUG -Wall -Wno-unused-variable -O3 -fopenmp -I /home/bart/work/eigen-3.4.0 -o eigen_test.o eigen_test.cc
g++ -m64 -std=c++17 -fPIC -I.  -I/apps/IntelParallelStudio/mkl/include -O2 -DNDEBUG -Wall -Wno-unused-variable -O3 -fopenmp -I /home/bart/work/eigen-3.4.0 eigen_test.o -o eigen_test  -L/apps/IntelParallelStudio/linux/lib/intel64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_rt -lmkl_core -liomp5 -lpthread

代码在这里，

//#define EIGEN_USE_MKL_ALL // Determine if use MKL
//#define EIGEN_VECTORIZE_SSE4_2

#include <iostream>
#include <Eigen/Dense>
using namespace Eigen;
int main()
{

  int n_a_rows = 10000;
  int n_a_cols = 10000;
  int n_b_rows = n_a_cols;
  int n_b_cols = 1000;

  MatrixXi a(n_a_rows, n_a_cols);

  for (int i = 0; i < n_a_rows; ++ i)
      for (int j = 0; j < n_a_cols; ++ j)
        a (i, j) = n_a_cols * i + j;

  MatrixXi b (n_b_rows, n_b_cols);
  for (int i = 0; i < n_b_rows; ++ i)
      for (int j = 0; j < n_b_cols; ++ j)
        b (i, j) = n_b_cols * i + j;

  MatrixXi d (n_a_rows, n_b_cols);

  clock_t begin = clock();

  d = a * b;

  clock_t end = clock();
  double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
  std::cout << "Time taken : " << elapsed_secs << std::endl;

}

在上一个与此主题相关的问题中，发现速度差异是涡轮增压（&差异不是很大）。我知道对于小矩阵，Eigen 可能比 MKL 更有效。但是我不明白为什么 Eigen+MKL 拒绝使用多核，即使我在编译过程中通过了 -liomp5。

提前致谢。 (CentOS 7 与 GCC 7.4.0 eigen 3.4.0)

Answer 1

请在执行程序之前设置以下 shell 变量。

export MKL_NUM_THREADS="$(nproc)"
export OMP_NUM_THREADS="$(nproc)"

还有构建命令（第一行运行和$define EIGEN_USE_MKL_ALL被注释掉了）

. /opt/intel/oneapi/setvars.sh
$CXX -I /usr/local/include/eigen3 eigen.cpp -o eigen_test -lblas
$CXX -I /usr/local/include/eigen3 eigen.cpp -o eigen_test -lmkl_rt

与 CXX 一起工作正常，如 clang++、g++ 和 icpx。如上所示设置环境很重要。在那种情况下 -lmkl_rt 就足够了。对代码稍作调整即可获得挂钟的净收益：

#define EIGEN_USE_BLAS
#define EIGEN_USE_MKL_ALL

#include <iostream>
#include <chrono>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std::chrono;
int main()
{

  int n_a_rows = 10000;
  int n_a_cols = 10000;
  int n_b_rows = n_a_cols;
  int n_b_cols = 1000;

  MatrixXd a(n_a_rows, n_a_cols);

  for (int i = 0; i < n_a_rows; ++ i)
      for (int j = 0; j < n_a_cols; ++ j)
        a (i, j) = n_a_cols * i + j;

  MatrixXd b (n_b_rows, n_b_cols);
  for (int i = 0; i < n_b_rows; ++ i)
      for (int j = 0; j < n_b_cols; ++ j)
        b (i, j) = n_b_cols * i + j;

  MatrixXd d (n_a_rows, n_b_cols);

  using wall_clock_t = std::chrono::high_resolution_clock;
  auto const start = wall_clock_t::now();
  clock_t begin = clock();

  d = a * b;

  clock_t end = clock();
  auto const wall = std::chrono::duration<double>(wall_clock_t::now() - start);
  
  double elapsed_secs = double(end - begin) / CLOCKS_PER_SEC;
  std::cout << "CPU time : " << elapsed_secs << std::endl;
  std::cout << "Wall time : " << wall.count() << std::endl;
  std::cout << "Speed up : " << elapsed_secs/wall.count() << std::endl;

}

我的 8 核 i7-4790K @4GHz 上的运行时间显示完美的并行化：

板载 blas：

CPU time : 12.5134
Wall time : 1.69036
Speed up : 7.40277

使用 MKL：

> ./eigen_test
CPU time : 11.4391
Wall time : 1.52542
Speed up : 7.49898

为什么带有 MKL 的 Eigen C++ 不为这种大型矩阵乘法使用多线程？

Why Eigen C++ with MKL doesn't use multi-threading for this large matrix multiplication?

c++

parallel-processing

linear-algebra

eigen3

intel-mkl