对于矩阵乘法，本征+ MKL比Matlab慢

Question

我正在C ++程序中进行大量矩阵乘法运算，并且使用与英特尔MKL（2018.3.222）链接的Eigen（3.3.5）。我使用MKL的顺序版本，并且OpenMP被禁用。问题是它比Matlab慢。

一些示例代码：

#define NDEBUG
#define EIGEN_USE_MKL_ALL

#include <iostream>
#include <chrono>
#include <Core>

using namespace Eigen;
using namespace std;

int main(){
    MatrixXd jac = 100*MatrixXd::Random(10*1228, 2850);
    MatrixXd res = MatrixXd::Zero(2850, 2850);

    for (int i=0; i<10; i++){
        auto begin = chrono::high_resolution_clock::now();
        res.noalias() = jac.transpose()*jac;
        auto end = chrono::high_resolution_clock::now();

        cout<<"time: "<<chrono::duration_cast<chrono::milliseconds>(end-begin).count() <<endl;
    }

    return 0;
}

它报告平均大约8秒。在带有g ++ 6.4的Ubuntu 16.04上使用-O3编译且没有调试符号。

Matlab代码：

m=100*(-1+2*rand(10*1228, 2850));
res = zeros(2850, 2850);
tic; res=m'*m; toc

它报告约4秒，是原来的两倍。我在具有maxNumCompThreads（1）的同一系统上使用了Matlab R2017a。Matlab使用MKL 11.3。

如果没有MKL，仅使用Eigen，则大约需要18秒。如何使C ++的运行时间降低到与Matlab相同的值？谢谢。

稍后编辑：正如@Qubit所建议的那样，Matlab意识到我正在尝试将矩阵与其转置相乘，并进行了一些“隐藏”优化。当我在Matlab中将两个不同的矩阵相乘时，时间增加到了8秒。因此，现在的问题就变成了：我如何告诉Eigen这个矩阵产品是“特殊的”并且可以进一步优化？

后来的编辑2：我尝试这样做：

MatrixXd jac = 100*MatrixXd::Random(10*1228, 2850);
MatrixXd res = MatrixXd::Zero(2850, 2850);

auto begin = chrono::high_resolution_clock::now();
res.selfadjointView<Lower>().rankUpdate(jac.transpose(), 1);
res.triangularView<Upper>() = res.transpose();
auto end = chrono::high_resolution_clock::now();

MatrixXd oldSchool = jac.transpose()*jac;
if (oldSchool.isApprox(res)){
    cout<<"same result!"<<endl;
}
cout<<"time: "<<chrono::duration_cast<chrono::milliseconds>(end-begin).count() <<endl;

但是现在只需要9.4秒（这是传统产品不需要MKL的Eigen的一半时间）。禁用MKL对此时间没有时间影响，因此我相信'rankUpdate'方法不使用MKL吗？

最后编辑：我在本征头文件中发现了一个错误：

Core/products/GeneralMatrixMatrixTriangular_BLAS.h

在第55行。括号放错了位置。我改变了这个：

if ( lhs==rhs && ((UpLo&(Lower|Upper)==UpLo)) ) { \

对此：

if ( lhs==rhs && ((UpLo&(Lower|Upper))==UpLo) ) { \

现在，我的C ++版本和Matlab具有相同的执行速度（在我的系统上约为4秒）。

Answer 1

因为您已经解决了问题，所以要真正给出答案，但有一些评论：

1. 该问题Core/products/GeneralMatrixMatrixTriangular_BLAS.h已经在 devel 分支中修复，但事实证明它从未被移植到 3.3 分支。

2. 该问题现已在 3.3 分支中修复。该修复将成为 3.3.6 的一部分。

3. 单线程模式下内置 Eigen 和 MKL 之间的加速因子 x2 没有意义。确保通过编译启用 CPU 支持的所有-march=native功能-O3 -DNDEBUG。在我的 Haswell 2.6GHz 上，我得到 3.4 秒与 3 秒。