Eigen: медленный доступ к столбцам матрицы 4

Question

Я использую Eigen для операций, подобных обновлению Cholesky, что подразумевает много AXPY (сумма плюс умножение на скаляр) на столбцах матрицы фиксированного размера, как правило, Matrix4d. Вкратце, это в 3 раза дороже, чтобы получить доступ к столбцам матрицы 4, чем в вектор 4.

Как правило, код ниже:

for(int i=0;i<4;++i) L.col(0) += x*y[i]; 

в 3 раза менее эффективен, чем ниже кода :

for(int i=0;i<4;++i) l4 += x*y[i]; 

, где L обычно представляет собой матрицу размером 4, X, Y и L4 являются векторами размера 4.

Кроме того, время, проведенное в первой строке кода не зависит от т atrix (или RowMajor от ColMajor).

На Intel i7 (2.5 ГГц) для векторных операций требуется около 0,007us, а для матричных операций - 0,02% (тайминги выполняются повторением 100000 раз по той же операции). Мое приложение потребует тысячи таких операций в таймингах, которые, надеюсь, намного ниже миллисекунды.

Вопрос: Я делаю что-то неправильно при доступе к столбцам моей матрицы 4x4? Есть ли чем-то, чтобы сделать первую строку кода более эффективной?

Полный код, используемый для таймингов ниже:

#include <iostream> 
#include <Eigen/Core> 
#include <vector> 
#include <sys/time.h> 

typedef Eigen::Matrix<double,4,1,Eigen::ColMajor> Vector4; 
//typedef Eigen::Matrix<double,4,4,Eigen::RowMajor,4,4> Matrix4; 
typedef Eigen::Matrix<double,4,4,Eigen::ColMajor,4,4> Matrix4; 

inline double operator- ( const struct timeval & t1,const struct timeval & t0) 
{ 
    /* TODO: double check the double conversion from long (on 64x). */ 
    return double(t1.tv_sec - t0.tv_sec)+1e-6*double(t1.tv_usec - t0.tv_usec); 
} 

void sumCols(Matrix4 & L, 
       Vector4 & x4, 
       Vector4 & y) 
{ 
    for(int i=0;i<4;++i) 
    { 
     L.col(0) += x4*y[i]; 
    } 
} 

void sumVec(Vector4 & L, 
      Vector4 & x4, 
      Vector4 & y) 
{ 
    for(int i=0;i<4;++i) 
    { 
     //L.tail(4-i) += x4.tail(4-i)*y[i]; 
     L   += x4   *y[i]; 
    } 
} 

int main() 
{ 
    using namespace Eigen; 

    const int NBT = 1000000; 

    struct timeval t0,t1; 

    std::vector<  Vector4> x4s(NBT); 
    std::vector<  Vector4> y4s(NBT); 
    std::vector<  Vector4> z4s(NBT); 
    std::vector<  Matrix4> L4s(NBT); 

    for(int i=0;i<NBT;++i) 
    { 
    x4s[i] = Vector4::Random(); 
    y4s[i] = Vector4::Random(); 
    L4s[i] = Matrix4::Random(); 
    } 

    int sample = int(z4s[55][2]/10*NBT); 
    std::cout << "*** SAMPLE = " << sample << std::endl; 

    gettimeofday(&t0,NULL); 
    for(int i=0;i<NBT;++i) 
    { 
     sumCols(L4s[i], x4s[i], y4s[i]); 
    } 
    gettimeofday(&t1,NULL); 
    std::cout << (t1-t0) << std::endl; 
    std::cout << "\t\t\t\t\t\t\tForce check" << L4s[sample](1,0) << std::endl; 

    gettimeofday(&t0,NULL); 
    for(int i=0;i<NBT;++i) 
    { 
     sumVec(z4s[i], x4s[i], y4s[i]); 
    } 
    gettimeofday(&t1,NULL); 
    std::cout << (t1-t0) << std::endl; 
    std::cout << "\t\t\t\t\t\t\tForce check" << z4s[sample][2] << std::endl; 

    return -1; 
} 

Answer 1

Как я уже сказал в комментариях, генерируемая сборка точно так же, как для функций.

Проблема заключается в том, что ваш бенчмарк смещен в том смысле, что L4s в 4 раза больше, чем z4s, и таким образом вы получаете больше промахов для кеша в матричном случае, чем в случае с вектором.