加速本征 C++ 转置？

Question

我知道这个 'eigen speed-up' 问题经常出现，但是在阅读了很多问题并尝试了几个标志之后，与执行转置的传统方式相比，我无法更好地使用 c++ eigen。实际上使用阻塞更有效。以下是代码

#include <cstdio>
#include <ctime>
#include <cstdlib>
#include <iostream>
#include <Eigen/Dense>

#define min( a, b ) ( ((a) < (b)) ? (a) : (b) )

int main(){
    const int n = 10000;
    const int csize = 32;
    float **a, **b;
    clock_t cputime1, cputime2;
    int i,j,k,ii,jj,kk;
  
    // Allocating memory for array/matrix
    a = new float * [n];
    for (i=0; i<n; i++){
        a[i] = new float [n];
    }
    b = new float * [n];
    for (i=0; i<n; i++){
        b[i] = new float[n];
    }
    // eigen matrices
    Eigen::MatrixXf M1 = Eigen::MatrixXf::Constant(n, n, 0.0);
    Eigen::MatrixXf M2 = Eigen::MatrixXf::Constant(n, n, 0.0);
    
    // Filling matrices with zeros
    for(i=0; i<n; ++i)
        for (j=0; j<n; ++j)
            a[i][j] = 0;
    for(i=0; i<n; ++i)
        for (j=0; j<n; ++j)
            b[i][j] = 0;

    // Direct (inefficient) transposition
    cputime1 = clock();
    for (i=0; i<n; ++i)
        for (j=0; j<n; ++j)
            a[i][j] = b[j][i];
    cputime2 = clock() - cputime1;
    std::printf("Time for transposition: %f\n", ((double)cputime2)/CLOCKS_PER_SEC);

    // Transposition using cache-blocking
    cputime1 = clock();
    for (ii=0; ii<n; ii+=csize)
        for (jj=0; jj<n; jj+=csize)
            for (i=ii; i<min(n,ii+csize-1); ++i)
                for (j=jj; j<min(n,jj+csize-1); ++j)
                    a[i][j] = b[j][i];
    cputime2 = clock() - cputime1;
    std::printf("Time for transposition: %f\n", ((double)cputime2)/CLOCKS_PER_SEC);

    // eigen
    cputime1 = clock();
    M1.noalias() = M2.transpose();
    cputime2 = clock() - cputime1;
    std::printf("Time for transposition with eigen: %f\n", ((double)cputime2)/CLOCKS_PER_SEC);

    // use data
    std::cout << a[n/2][n/2] << std::endl;
    std::cout << b[n/2][n/2] << std::endl;
    std::cout << M1(n/2,n/2) << std::endl;

    return 0;
}

而我使用的编译命令是

g++ -fno-math-errno -ffast-math -march=native -fopenmp -O2 -msse2 -DNDEBUG  blocking_and_eigen.cpp

有结果

Time for transposition: 1.926674
Time for transposition: 0.280653
Time for transposition with eigen: 2.018217

我正在使用 eigen 3.4.0 和 g++ 11.2.0。

您对提高本征性能有什么建议吗？提前致谢

Answer 1

正如 INS 在评论中所建议的那样，矩阵的实际复制导致性能下降，我稍微修改了您的示例以使用一些数字而不是全零（以避免任何类型的优化）：

for(i=0; i<n; ++i) {
    for (j=0; j<n; ++j) {
        a[i][j] = i+j;
        M1(i,j) = i+j;
      }
}
for(i=0; i<n; ++i) {
    for (j=0; j<n; ++j) {
        b[i][j] = i+j;
        M1(i,j) = i+j;
    }
}

此外，我修改了最终的打印语句，对结果进行了全面检查（如果不正确，将针对 M2 执行检查）：

    for (i=0; i<n; ++i)
    for (j=0; j<n; ++j)
      if (a[i][j] != M1(i,j))
        std::cout << "Diff here! " << std::endl;

然后我尝试了几个测试：

预分配和赋值

Eigen::MatrixXf M2 = Eigen::MatrixXf::Constant(n, n, 0.0); ...这里有一些代码... M2 = M1.transpose();
复制构造函数

特征值::MatrixXf M2(M1.transpose());
到位

M1.transpose就地();
使用 auto 和 c++17 复制构造

自动 M2{ M1.transpose() };

这是最令人费解的，性能非常出色，我认为这个故事有两个部分，如果我为案例 2 和案例 4 打印 M2 的 typeid 名称，它们是不同的，并且名称被破坏但是它给我们一个线索：

N5本征6矩阵IfLin1ELin1ELi0ELin1ELin1EEE N5特征9转置INS_6矩阵IfLin1ELin1ELi0ELin1ELin1EEEEE

auto 关键字解析为特定于转置矩阵的不同类型。故事的第二部分是 M1 之后没有修改的事实，所以要么编译器移动它，要么很可能 EigenTransposeMatrix (https://eigen.tuxfamily.org/dox/classEigen_1_1Transpose.html) 只保留原始矩阵的引用并且它不复制它。

结果

Test	Direct (s)	Cache block (s)	eigen (s)
1	2.633	0.312	1.861
2	2.599	0.262	1.968
3	2.602	0.262	0.216
4	2.552	0.280	0.000002

加速本征 C++ 转置？

Speed-up eigen c++ transpose?

c++

performance

transpose

linear-algebra

eigen