GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-18 06:35:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

feature 319: fix LDLT::rankUpdate for complex/upper, simply the algortihm, update copyrights

Browse Source
This commit is contained in:
Gael Guennebaud 2011-12-09 23:08:38 +01:00
parent 2d7c3eea53
commit 38277e8a9b
1 changed files with 21 additions and 34 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
Eigen/src/Cholesky/LDLT.h
View File

@ -1,9 +1,10 @@
// This file is part of Eigen, a lightweight C++ template library // This file is part of Eigen, a lightweight C++ template library
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2009 Keir Mierle <mierle@gmail.com> // Copyright (C) 2009 Keir Mierle <mierle@gmail.com>
// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com>
// //
// Eigen is free software; you can redistribute it and/or // Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public // modify it under the terms of the GNU Lesser General Public
@ -115,7 +116,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
/** Clear any existing decomposition /** Clear any existing decomposition
* \sa rankUpdate(w,sigma) * \sa rankUpdate(w,sigma)
*/ */
void clear() void setZero()
{ {
m_isInitialized = false; m_isInitialized = false;
} }
@ -143,14 +144,14 @@ template<typename _MatrixType, int _UpLo> class LDLT
} }
/** \returns the coefficients of the diagonal matrix D */ /** \returns the coefficients of the diagonal matrix D */
inline Diagonal<const MatrixType> vectorD(void) const inline Diagonal<const MatrixType> vectorD() const
{ {
eigen_assert(m_isInitialized && "LDLT is not initialized."); eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix.diagonal(); return m_matrix.diagonal();
} }
/** \returns true if the matrix is positive (semidefinite) */ /** \returns true if the matrix is positive (semidefinite) */
inline bool isPositive(void) const inline bool isPositive() const
{ {
eigen_assert(m_isInitialized && "LDLT is not initialized."); eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == 1; return m_sign == 1;
@ -348,47 +349,33 @@ template<> struct ldlt_inplace<Lower>
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
typedef typename MatrixType::ColXpr ColXpr;
typedef typename internal::remove_all<ColXpr>::type ColXprCleaned;
typedef typename ColXprCleaned::SegmentReturnType ColXprSegment;
typedef typename MatrixType::Scalar Scalar;
// typedef Matrix<Scalar,Dynamic,1> TempVectorType;
typedef typename WDerived::SegmentReturnType TempVecSegment;
const Index size = mat.rows(); const Index size = mat.rows();
eigen_assert(mat.cols() == size && w.size()==size); eigen_assert(mat.cols() == size && w.size()==size);
// Prepare the update RealScalar alpha = 1;
RealScalar alpha,alphabar,temp,dtemp,gammatmp;
Scalar wtemp,gamma;
alpha = 1;
// Apply the update // Apply the update
for (Index j = 0; j < size; j++) for (Index j = 0; j < size; j++)
{ {
// Check for termination due to an original decomposition of low-rank // Check for termination due to an original decomposition of low-rank
if (!std::isfinite(alpha)) if (!std::isfinite(alpha))
break; break;
// Update the diagonal terms // Update the diagonal terms
dtemp = real(mat.diagonal().coeff(j)); RealScalar dj = real(mat.coeff(j,j));
wtemp = w.coeff(j); Scalar wj = w.coeff(j);
temp = sigma*real(wtemp*conj(wtemp)); RealScalar swj2 = sigma*abs2(wj);
alphabar = alpha + temp/dtemp; RealScalar gamma = dj*alpha + swj2;
gammatmp = dtemp*alpha + temp;
if (gammatmp != 0) mat.coeffRef(j,j) += swj2/alpha;
gamma = conj(wtemp)/gammatmp; // FIXME: guessing on conj here alpha += swj2/dj;
else
gamma = 0;
dtemp += temp/alpha;
alpha = alphabar;
mat.diagonal().coeffRef(j) = dtemp;
// Update the terms of L // Update the terms of L
w.tail(size-j-1) -= wtemp*mat.col(j).tail(size-j-1); Index rs = size-j-1;
mat.col(j).tail(size-j-1) += (sigma*gamma)*w.tail(size-j-1); w.tail(rs) -= wj * mat.col(j).tail(rs);
if(gamma != 0)
mat.col(j).tail(rs) += (sigma*conj(wj)/gamma)*w.tail(rs);
} }
return true; return true;
} }
@ -416,7 +403,7 @@ template<> struct ldlt_inplace<Upper>
static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, typename MatrixType::RealScalar sigma=1) static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, typename MatrixType::RealScalar sigma=1)
{ {
Transpose<MatrixType> matt(mat); Transpose<MatrixType> matt(mat);
return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w, sigma); return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma);
} }
}; };
@ -461,7 +448,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
/** Update the LDLT decomposition: given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T. /** Update the LDLT decomposition: given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T.
* \param w a vector to be incorporated into the decomposition. * \param w a vector to be incorporated into the decomposition.
* \param sigma a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1. * \param sigma a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1.
* \sa clear() * \sa setZero()
*/ */
template<typename MatrixType, int _UpLo> template<typename MatrixType, int _UpLo>
template<typename Derived> template<typename Derived>