Почему Eigen в 5 раз медленнее, чем ублас на следующем примере?

В версии Eigen я использую «истинные» матрицы и векторы фиксированного размера, лучший алгоритм (LDLT против LU в uBlas), он использует инструкции SIMD для внутреннего использования. Итак, почему это медленнее, чем uBlas на следующем примере?

Я уверен, что я делаю что-то не так — Эйген ДОЛЖЕН быть быстрее или, по крайней мере, сопоставимым.

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/progress.hpp>
#include <Eigen/Dense>
#include <iostream>

using namespace boost;
using namespace std;
const int n=9;
const int total=100000;

void test_ublas()
{
using namespace boost::numeric::ublas;
cout << "Boost.ublas ";
double r=1.0;
{
boost::progress_timer t;
for(int j=0;j!=total;++j)
{
//symmetric_matrix< double,lower,row_major,bounded_array<double,(1+n)*n/2> > A(n,n);
matrix<double,row_major,bounded_array<double,n*n> > A(n,n);
permutation_matrix< unsigned char,bounded_array<unsigned char,n> > P(n);
bounded_vector<double,n> v;
for(int i=0;i!=n;++i)
for(int k=0;k!=n;++k)
A(i,k)=0.0;
for(int i=0;i!=n;++i)
{
A(i,i)=1.0+i;
v[i]=i;
}
lu_factorize(A,P);
lu_substitute(A,P,v);
r+=inner_prod(v,v);
}
}
cout << r << endl;
}

void test_eigen()
{
using namespace Eigen;
cout << "Eigen ";
double r=1.0;
{
boost::progress_timer t;
for(int j=0;j!=total;++j)
{
Matrix<double,n,n> A;
Matrix<double,n,1> b;
for(int i=0;i!=n;++i)
{
A(i,i)=1.0+i;
b[i]=i;
}
Matrix<double,n,1> x=A.ldlt().solve(b);
r+=x.dot(x);
}
}
cout << r << endl;
}

int main()
{
test_ublas();
test_eigen();
}

Выход:

Boost.ublas 0.50 s

488184
Eigen 2.66 s

488184

(Visual Studio 2010 x64 Release)

РЕДАКТИРОВАТЬ:

За

const int n=4;
const int total=1000000;

Выход:

Boost.ublas 0.67 s

1.25695e+006
Eigen 0.40 s

5.4e+007

Я полагаю, такое поведение связано с тем, что версия uBlas вычисляет факторизацию на месте, в то время как версия Eigen создает COPY матрицы (LDLT), поэтому она хуже подходит для кеша.

Есть ли способ сделать вычисления на месте в Eigen? Или может быть есть другие способы улучшить это?

РЕДАКТИРОВАТЬ:

Следуя советам Fezvez и использую LLT вместо LDLT, я получаю:

Eigen 0.16 s

488184

Это хорошо, но все равно делает ненужное выделение стека матрицы:

sizeof(A.llt()) == 656

Я предпочитаю избегать этого — это должно быть еще быстрее.

РЕДАКТИРОВАТЬ:

Я удалил распределение, создав подклассы из LDLT (это внутренняя матрица защищена), и заполняя его напрямую. Теперь результат для ЛПНП:

Eigen 0.26 s
488209

Это работает, но это обходной путь — не реальное решение …

Подклассы из LLT также работают, но не дают такого большого эффекта.