Я пытаюсь заставить 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
может быть только положительным полуопределенным, в этом случае требуется разложение LDLTLLT
разложение вычисляет sqrt()
излишне для решения системыЯ не мог найти способ подключиться к функциональности LDLT в Eigen, как в коде выше, так как код структурирован очень по-разному.
Итак, мой вопрос: есть ли способ использовать Eigen3 для решения линейной системы, используя разложения LDLT, используя не более царапин пространство, чем для диагональной матрицы D
?
Одним из вариантов является выделение решателя 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
это столбец мажор.
Других решений пока нет …