Eigen vs Matlab:并行化矩阵乘法

Question

我想比较矩阵乘法中Matlab的速度与Intel(R)Core(TM)i7-4770 CPU @ 3.40GHz上的Eigen 3的速度.代码包括Eigen:

#include <iostream>
#include "Eigen/Dense"
#include <chrono>
#include <omp.h>


using namespace std;
using namespace Eigen;

const int dim=100;

int main()
{
    std::chrono::time_point<std::chrono::system_clock> start, end;

    int n;
    n = Eigen::nbThreads();
    cout<<n<<"\n";

    Matrix<double, Dynamic, Dynamic> m1(dim,dim);
    Matrix<double, Dynamic, Dynamic> m2(dim,dim);
    Matrix<double, Dynamic, Dynamic> m_res(dim,dim);

    start = std::chrono::system_clock::now();

    for (int i = 0 ; i <100000; ++i) {
        m1.setRandom(dim,dim);
        m2.setRandom(dim,dim);
        m_res=m1*m2;

    }

    end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end-start;

    std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";

    return 0;
}

它是用编译g++ -O3 -std=c++11 -fopenmp和执行的OMP_NUM_THREADS=8 ./prog.在Matlab我正在使用

function mat_test(N,dim)
%
% N:    how many tests
% dim:  dimension of the matrices

tic
parfor i=1:N
     A = rand(dim);
     B = rand(dim);
     C = A*B;
end
toc

结果是:Matlab为9s ,Eigen为36s.我在Eigen案中做错了什么？我可以排除矩阵的动态分配.此外,仅使用3个线程而不是8个线程.

编辑:

也许我没有说得足够清楚:任务是乘以10000倍dim = 100的双值矩阵,每次随机填充,不仅一​​次.使用Eigen尽可能快地完成.如果Eigen不能应付Matlab,你会建议什么选择？

Answer 1

下面是一个更好的代码版本,可以合理地使用Eigen.总结一下:

• 移动setRandom()基准测试循环外部.setRandom()调用rand()相当慢的系统功能.
• 使用.noalias()以避免临时(仅使得当右手侧是一种产品的意义上)的创建
• 将OMP_NUM_THREADS设置为真正的核心数,而不是超线程数.(在你的情况下为4)
• 你的CPU支持只有Eigen的devel分支支持的AVX和FMA(将变为3.3),所以使用devel分支并使用-mavx和-mfma编译器选项启用它们(与仅SSE相比,约为x3.5加速)

代码:

#include <iostream>
#include "Eigen/Dense"
#include <chrono>

using namespace std;
using namespace Eigen;

const int dim=100;

int main()
{
    std::chrono::time_point<std::chrono::system_clock> start, end;

    int n;
    n = Eigen::nbThreads();
    cout << n << "\n";

    Matrix<double, Dynamic, Dynamic> m1(dim,dim);
    Matrix<double, Dynamic, Dynamic> m2(dim,dim);
    Matrix<double, Dynamic, Dynamic> m_res(dim,dim);

    start = std::chrono::system_clock::now();

    m1.setRandom();
    m2.setRandom();

    for (int i = 0 ; i <100000; ++i) {
      m_res.noalias() = m1 * m2;
    }

    end = std::chrono::system_clock::now();
    std::chrono::duration<double> elapsed_seconds = end-start;

    std::cout << "elapsed time: " << elapsed_seconds.count() << "s\n";

    return 0;
}