2016-03-29 437 views
0

我想要Eigen3解決一個線性系統A * X = B與就地喬列斯基分解。我無法承受任何尺寸爲A的臨時物品,但我可以在此過程中自由銷燬AEigen LDLT Cholesky分解就地

不幸的是,

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()不必要地用於該系統的溶液

我無法找到與上述代碼類似的Eigen LDLT功能掛鉤方式,因爲代碼結構非常不同。

所以我的問題是:有沒有一種方法可以使用Eigen3來解決一個線性系統使用LDLT分解使用沒有更多的臨時空間比對角矩陣D

回答

0

一種選擇是將分配一個活體肝移植求解只有一次,並調用計算方法:

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

如果這也不是一個選項,你可以常量投用活體肝移植的對象存儲矩陣:

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

戲劇與A,然後:

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

這是醜陋的,但這種只要MatType是列專業應該工作。

+0

不幸的是,我真的不得不用我自己的記憶,所以這兩者都不起作用。我認爲我需要的是'internal :: ldlt_inplace :: unblocked()',但是它不像設置LLT那麼明顯。 – Stefan

+0

另外,如果我可以得到'unblocked()'設置,我們的矩陣'A'是行很重要的,但是因爲'A'是對稱的,所以我應該只能使用'A.transpose()',否? – Stefan

+0

如果矩陣是完整的,那麼,是的,你可以看到一個行主要的上三角部分作爲列主要的下三角部分。你只需要分配一個PermutationMatrix來將它傳遞給'internal :: ldlt_inplace :: unblocked()'。主要問題是您將不得不重新編寫解決步驟。 – ggael