Eigen LDLT Cholesky разложение на месте

Question

Я пытаюсь получить Eigen3 для решения линейной системы A * X = B с декомпозицией на месте Cholesky. Я не могу позволить себе иметь временные файлы размером A, нажатыми на стек, но я могу уничтожить A в процессе.

К сожалению,

A.llt().solveInPlace(B); 

может быть и речи, так как A.llt() неявно выталкивает временную матрицу размером A на стеке. Для LLT случае, я мог бы получить доступ к необходимой функциональности, как так:

// solve A * X = B in-place for positive-definite A 
template <typename AType, typename BType> 
void AllInPlaceSolve(AType& A, BType& B) 
{ 
    typedef Eigen::internal::LLT_Traits<AType, Eigen::Upper> TraitsType; 
    TraitsType::inplace_decomposition(A); 
    TraitsType::getL(A).solveInPlace(B); 
    TraitsType::getU(A).solveInPlace(B); 
} 

Это прекрасно работает, но я обеспокоен тем, что:

• Мои матрицы A может быть неотрицательно только, в котором случае разложение LDLT требуется
• разложение LLT вычисляет sqrt() без необходимости для решения системы

Я не мог найти способ подключить функциональность LDLT от Eigen аналогично приведенному выше коду, поскольку код структурирован очень по-разному.

Итак, мой вопрос: есть ли способ использовать Eigen3 для решения линейной системы с использованием декомпозиций LDLT, используя не более царапающее пространство, чем для диагональной матрицы D?

Answer 1

Один из вариантов состоит в выделении LDLT решатель только один раз, и вызвать метод вычислений:

LDLT<MatType> ldlt(size); 
// ... 
ldlt.compute(A); 
x = ldlt.solve(b); 

Если это тоже не вариант, вы можете Const отливать матрицу, сохраненную LDLT объекта:

LDLT<MatType> ldlt(MatType::Identity(size,size)); 
MatType& A = const_cast<MatType&>(ldlt.matrixLDLT()); 

играет с A, а затем:

ldlt.compute(A); 
x = ldlt.solve(b); 

Это некрасиво, но это должен работать до тех пор, пока MatType является основным столбцом.