Собственная ЛНПЛ разложения Холецкого на месте

Я пытаюсь заставить Eigen3 решить линейную систему A * X = B с на месте разложения Холецкого. Я не могу позволить себе иметь временные по размеру 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?

0

Решение

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

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

Если это тоже не вариант, вы можете преобразовать матрицу, хранящуюся в объекте ldlt:

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

играет с A, а потом:

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

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

0

Другие решения

Других решений пока нет …

По вопросам рекламы [email protected]