MOAB
4.9.3pre
|
#include <LDLT.h>
Static Public Member Functions | |
template<typename MatrixType , typename TranspositionType , typename Workspace > | |
static bool | unblocked (MatrixType &mat, TranspositionType &transpositions, Workspace &temp, SignMatrix &sign) |
template<typename MatrixType , typename WDerived > | |
static bool | updateInPlace (MatrixType &mat, MatrixBase< WDerived > &w, const typename MatrixType::RealScalar &sigma=1) |
template<typename MatrixType , typename TranspositionType , typename Workspace , typename WType > | |
static bool | update (MatrixType &mat, const TranspositionType &transpositions, Workspace &tmp, const WType &w, const typename MatrixType::RealScalar &sigma=1) |
static bool Eigen::internal::ldlt_inplace< Lower >::unblocked | ( | MatrixType & | mat, |
TranspositionType & | transpositions, | ||
Workspace & | temp, | ||
SignMatrix & | sign | ||
) | [inline, static] |
Definition at line 257 of file LDLT.h.
{ using std::abs; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename TranspositionType::StorageIndex IndexType; eigen_assert(mat.rows()==mat.cols()); const Index size = mat.rows(); if (size <= 1) { transpositions.setIdentity(); if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef; else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef; else sign = ZeroSign; return true; } for (Index k = 0; k < size; ++k) { // Find largest diagonal element Index index_of_biggest_in_corner; mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); index_of_biggest_in_corner += k; transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner); if(k != index_of_biggest_in_corner) { // apply the transposition while taking care to consider only // the lower triangular part Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k)); mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s)); std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner)); for(Index i=k+1;i<index_of_biggest_in_corner;++i) { Scalar tmp = mat.coeffRef(i,k); mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i)); mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp); } if(NumTraits<Scalar>::IsComplex) mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k)); } // partition the matrix: // A00 | - | - // lu = A10 | A11 | - // A20 | A21 | A22 Index rs = size - k - 1; Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1); Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k); Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k); if(k>0) { temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint(); mat.coeffRef(k,k) -= (A10 * temp.head(k)).value(); if(rs>0) A21.noalias() -= A20 * temp.head(k); } // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot // was smaller than the cutoff value. However, since LDLT is not rank-revealing // we should only make sure that we do not introduce INF or NaN values. // Remark that LAPACK also uses 0 as the cutoff value. RealScalar realAkk = numext::real(mat.coeffRef(k,k)); if((rs>0) && (abs(realAkk) > RealScalar(0))) A21 /= realAkk; if (sign == PositiveSemiDef) { if (realAkk < 0) sign = Indefinite; } else if (sign == NegativeSemiDef) { if (realAkk > 0) sign = Indefinite; } else if (sign == ZeroSign) { if (realAkk > 0) sign = PositiveSemiDef; else if (realAkk < 0) sign = NegativeSemiDef; } } return true; }
static bool Eigen::internal::ldlt_inplace< Lower >::update | ( | MatrixType & | mat, |
const TranspositionType & | transpositions, | ||
Workspace & | tmp, | ||
const WType & | w, | ||
const typename MatrixType::RealScalar & | sigma = 1 |
||
) | [inline, static] |
Definition at line 385 of file LDLT.h.
{ // Apply the permutation to the input w tmp = transpositions * w; return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma); }
static bool Eigen::internal::ldlt_inplace< Lower >::updateInPlace | ( | MatrixType & | mat, |
MatrixBase< WDerived > & | w, | ||
const typename MatrixType::RealScalar & | sigma = 1 |
||
) | [inline, static] |
Definition at line 347 of file LDLT.h.
{ using numext::isfinite; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; const Index size = mat.rows(); eigen_assert(mat.cols() == size && w.size()==size); RealScalar alpha = 1; // Apply the update for (Index j = 0; j < size; j++) { // Check for termination due to an original decomposition of low-rank if (!(isfinite)(alpha)) break; // Update the diagonal terms RealScalar dj = numext::real(mat.coeff(j,j)); Scalar wj = w.coeff(j); RealScalar swj2 = sigma*numext::abs2(wj); RealScalar gamma = dj*alpha + swj2; mat.coeffRef(j,j) += swj2/alpha; alpha += swj2/dj; // Update the terms of L Index rs = size-j-1; w.tail(rs) -= wj * mat.col(j).tail(rs); if(gamma != 0) mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs); } return true; }