求问C++的Eigen矩阵运算库有没有提供两个矩阵对应元素相乘的方法

Posted 2023-05-01

参考技术A #includeusing namespace std;#define M 6#define N 3void mulMatri(int x[M][N],int y[N][M],int z[M][M],int m,int n);int main() int i,j; int x[M][N],y[N][M],z[M][M]; for(i=0;i>x[i][j]; for(i=0;i>y[i][j]; mulMatri( x,y,z,M,N); cout<追问

如果我没有看错的话，您的意思应该是用循环实现吧。这不是我想要的。不过还是谢谢您的回答。我已经找到方法了：MatrixXd类型相乘是矩阵积，转化成Array类型相乘就是对应元素相乘。

本回答被提问者和网友采纳 参考技术B 直接用 m.cwiseProduct(n) ，不用转为array 参考技术C 转化为array再乘, a.array()*b.array() ,很优雅的一个操作。矩阵和array可以互相转化，那样乘要求行和列相同。

使用 MKL 编译时，Eigen C++ 运行速度较慢

【中文标题】使用 MKL 编译时，Eigen C++ 运行速度较慢

【英文标题】：Eigen C++ running slower when compiled with MKL

【发布时间】：2016-12-11 18:02:56

【问题描述】：

我最近开始使用 Eigen（版本 3.3.1），在 OLS 回归核心的简单矩阵运算上针对 Armadillo 运行基准测试，即自行计算矩阵乘积的逆，我注意到使用 MKL 库编译时，Eigen 的运行速度比没有它时运行得慢。我想知道我的编译说明是否错误。我还尝试直接调用 MKL BLAS 和 LAPACK 例程来实现此操作，并获得了更快的结果，与犰狳一样快。我无法解释如此糟糕的性能，尤其是对于浮点类型。

我写了下面的代码来实现这个基准：

#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 run -L$MKLROOT/lib/intel64 -Wl,--no-as-needed -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.9s

对于浮点类型：

带 MKL 的犰狳 => 38.2 秒

MKL 的特征 => 61.1s

仅本征 => 42.6s

纯 MKL => 38.9s

注意：我为一个不会使用非常大矩阵的项目运行此测试，我不需要在这个级别进行并行化，我最大的矩阵可能是 25 列的 2000 行，而且我需要并行更高级别，所以我想避免任何可能减慢代码速度的嵌套并行性。

【问题讨论】：

问题是什么？

问题是为什么它在使用 MKL 编译时运行得这么慢，而我可以更快地直接调用 MKL 例程，我做错了什么吗？

很难说，我建议在分析器下运行它（英特尔的 VTune 非常适合，并且可以在您的项目中稍后显示并行度评估）。

我无法重现如此大的差异。确保在基准测试时禁用涡轮增压，否则你的数字毫无意义，因为你不能保证恒定的频率。

【参考方案1】：

正如我在评论中所说，请确保在进行基准测试时禁用涡轮增压。

作为旁注和供将来参考，您当前的 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();

不过，对于这么小的矩阵，我看不出有太大的差异。

【讨论】：

确实我的错，我忘了禁用涡轮增压，但是即使这样做我仍然有很大的差异，只有在我的代码中替换你写的我得到了更均匀的结果，没有外部的 Eigen库甚至最快的双倍和稍慢的浮点数。我没有意识到我在调用 gemm 而不是 syrk 因为大多数其他库似乎可以检测到同一对象何时在矩阵产品中使用两次。非常感谢我注意到我似乎比我想象的要了解更多关于 Eigen 的知识！

【参考方案2】：

我只是想补充一点，以防有些人可能对此有疑问，ggael 给出的表达式必须写成如下，以防它是模板函数/类的一部分，否则编译器将难以进行类型推断

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