Почему 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 также работают, но не дают такого большого эффекта.

b=A.ldlt().solve(b);
r+=b.dot(b);

Boost.ublas 0.15 s
488184

Eigen 0.08 s
488184

Boost.ublas 0.47 s
Eigen 0.13 s

There were  100000  Eigen standard LDLT linear solvers applied in  12.658  seconds
There were  100000  Eigen standard LLT linear solvers applied in  4.652  seconds

There were  100000  Eigen in place LDLT linear solvers applied in  12.7  seconds
There were  100000  Eigen in place LLT linear solvers applied in  4.396  seconds