代码之家  ›  专栏  ›  技术社区  ›  Olivier

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

  •  1
  • Olivier  · 技术社区  · 8 年前

    最近,我开始使用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与:

    我得到以下结果(在Intel(R)Core(TM)i5-3210M CPU上,2.50GHz)

    • 对于双重类型:

    具有MKL的犰狳=>64.0s

    具有MKL的特征值=>72.2秒

    仅本征值=>68.7秒

    纯MKL=>64.9秒

    • 对于浮动类型:

    具有MKL的犰狳=>38.2s

    具有MKL的特征值=>61.1秒

    仅本征值=>42.6s

    纯MKL=>38.9秒

    注意:我为一个不会使用非常大的矩阵的项目运行这个测试,我不需要在这个级别上进行并行化,我最大的矩阵可能是25列2000行,而且我需要在更高的级别上并行,因此我希望避免任何可能会降低代码速度的嵌套并行。

    2 回复  |  直到 8 年前
        1
  •  2
  •   ggael    8 年前

    正如我在评论中所说,确保在基准测试时禁用turbo boost。

    作为补充说明和未来参考,您当前的特征码将调用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();
    

    对于这样小的矩阵,我看不出有什么不同。

        2
  •  0
  •   Olivier    8 年前

    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.9秒

    具有MKL的特征值=>79.8秒

    仅本征值=>71.1秒

    纯MKL=>81.1s

    • 对于浮动类型:

    具有MKL的犰狳=>47.2s

    仅本征值=>51.8秒

    纯MKL=>48.0s