特征:访问矩阵4列的速度很慢

我使用Eigen进行类似于Cholesky更新的操作，这意味着在固定大小的矩阵(通常是Matrix4d)的列上有大量的AXPY(和加上标量乘法)。简而言之，访问矩阵4的列的开销是访问向量4的3倍。

通常，下面的代码:

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

的效率是下面代码的三分之一:

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

其中L通常是大小为4的矩阵，x, y和14是大小为4的向量。

此外，在第一行代码中花费的时间并不取决于矩阵存储组织(RowMajor或ColMajor)。

在Intel i7 (2.5GHz)上，矢量运算大约需要0.007us，矩阵运算大约需要0.02us(计时是通过重复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 << "tttttttForce 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 << "tttttttForce check" <<  z4s[sample][2] << std::endl;
  return -1;
}