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

SimplicialCholesky*: s/LLt/LLT and s/LDLt/LDLT for consistency with dense names

Browse Source
This commit is contained in:
Gael Guennebaud 2012-02-27 14:28:07 +01:00
parent ece30e9e6f
commit 5effdba2c6
1 changed files with 61 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

122
Eigen/src/SparseCholesky/SimplicialCholesky.h
View File

@ -64,8 +64,8 @@ LDL License:
#define EIGEN_SIMPLICIAL_CHOLESKY_H #define EIGEN_SIMPLICIAL_CHOLESKY_H
enum SimplicialCholeskyMode { enum SimplicialCholeskyMode {
SimplicialCholeskyLLt, SimplicialCholeskyLLT,
SimplicialCholeskyLDLt SimplicialCholeskyLDLT
}; };
/** \ingroup SparseCholesky_Module /** \ingroup SparseCholesky_Module
@ -142,7 +142,7 @@ class SimplicialCholeskyBase
inline const internal::solve_retval<SimplicialCholeskyBase, Rhs> inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
solve(const MatrixBase<Rhs>& b) const solve(const MatrixBase<Rhs>& b) const
{ {
eigen_assert(m_isInitialized && "Simplicial LLt or LDLt is not initialized."); eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
eigen_assert(rows()==b.rows() eigen_assert(rows()==b.rows()
&& "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b"); && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
@ -156,7 +156,7 @@ class SimplicialCholeskyBase
inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs> inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
solve(const SparseMatrixBase<Rhs>& b) const solve(const SparseMatrixBase<Rhs>& b) const
{ {
eigen_assert(m_isInitialized && "Simplicial LLt or LDLt is not initialized."); eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
eigen_assert(rows()==b.rows() eigen_assert(rows()==b.rows()
&& "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b"); && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
@ -256,10 +256,10 @@ class SimplicialCholeskyBase
protected: protected:
template<bool DoLDLt> template<bool DoLDLT>
void factorize(const MatrixType& a); void factorize(const MatrixType& a);
void analyzePattern(const MatrixType& a, bool doLDLt); void analyzePattern(const MatrixType& a, bool doLDLT);
/** keeps off-diagonal entries; drops diagonal entries */ /** keeps off-diagonal entries; drops diagonal entries */
struct keep_diag { struct keep_diag {
@ -275,7 +275,7 @@ class SimplicialCholeskyBase
bool m_analysisIsOk; bool m_analysisIsOk;
CholMatrixType m_matrix; CholMatrixType m_matrix;
VectorType m_diag; // the diagonal coefficients (LDLt mode) VectorType m_diag; // the diagonal coefficients (LDLT mode)
VectorXi m_parent; // elimination tree VectorXi m_parent; // elimination tree
VectorXi m_nonZerosPerCol; VectorXi m_nonZerosPerCol;
PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation
@ -285,13 +285,13 @@ class SimplicialCholeskyBase
RealScalar m_shiftScale; RealScalar m_shiftScale;
}; };
template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLt; template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLt; template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky; template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
namespace internal { namespace internal {
template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLt<_MatrixType,_UpLo> > template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> >
{ {
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
@ -304,7 +304,7 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLt<_MatrixTyp
static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
}; };
template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLt<_MatrixType,_UpLo> > template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
{ {
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
@ -326,8 +326,8 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_Matr
} }
/** \ingroup SparseCholesky_Module /** \ingroup SparseCholesky_Module
* \class SimplicialLLt * \class SimplicialLLT
* \brief A direct sparse LLt Cholesky factorizations * \brief A direct sparse LLT Cholesky factorizations
* *
* This class provides a LL^T Cholesky factorizations of sparse matrices that are * This class provides a LL^T Cholesky factorizations of sparse matrices that are
* selfadjoint and positive definite. The factorization allows for solving A.X = B where * selfadjoint and positive definite. The factorization allows for solving A.X = B where
@ -337,39 +337,39 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_Matr
* \tparam _UpLo the triangular part that will be used for the computations. It can be Lower * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower. * or Upper. Default is Lower.
* *
* \sa class SimplicialLDLt * \sa class SimplicialLDLT
*/ */
template<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo>
class SimplicialLLt : public SimplicialCholeskyBase<SimplicialLLt<_MatrixType,_UpLo> > class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
typedef SimplicialCholeskyBase<SimplicialLLt> Base; typedef SimplicialCholeskyBase<SimplicialLLT> Base;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
typedef Matrix<Scalar,Dynamic,1> VectorType; typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef internal::traits<SimplicialLLt> Traits; typedef internal::traits<SimplicialLLT> Traits;
typedef typename Traits::MatrixL MatrixL; typedef typename Traits::MatrixL MatrixL;
typedef typename Traits::MatrixU MatrixU; typedef typename Traits::MatrixU MatrixU;
public: public:
/** Default constructor */ /** Default constructor */
SimplicialLLt() : Base() {} SimplicialLLT() : Base() {}
/** Constructs and performs the LLt factorization of \a matrix */ /** Constructs and performs the LLT factorization of \a matrix */
SimplicialLLt(const MatrixType& matrix) SimplicialLLT(const MatrixType& matrix)
: Base(matrix) {} : Base(matrix) {}
/** \returns an expression of the factor L */ /** \returns an expression of the factor L */
inline const MatrixL matrixL() const { inline const MatrixL matrixL() const {
eigen_assert(Base::m_factorizationIsOk && "Simplicial LLt not factorized"); eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
return Traits::getL(Base::m_matrix); return Traits::getL(Base::m_matrix);
} }
/** \returns an expression of the factor U (= L^*) */ /** \returns an expression of the factor U (= L^*) */
inline const MatrixU matrixU() const { inline const MatrixU matrixU() const {
eigen_assert(Base::m_factorizationIsOk && "Simplicial LLt not factorized"); eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
return Traits::getU(Base::m_matrix); return Traits::getU(Base::m_matrix);
} }
@ -404,8 +404,8 @@ public:
}; };
/** \ingroup SparseCholesky_Module /** \ingroup SparseCholesky_Module
* \class SimplicialLDLt * \class SimplicialLDLT
* \brief A direct sparse LDLt Cholesky factorizations without square root. * \brief A direct sparse LDLT Cholesky factorizations without square root.
* *
* This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are
* selfadjoint and positive definite. The factorization allows for solving A.X = B where * selfadjoint and positive definite. The factorization allows for solving A.X = B where
@ -415,45 +415,45 @@ public:
* \tparam _UpLo the triangular part that will be used for the computations. It can be Lower * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower. * or Upper. Default is Lower.
* *
* \sa class SimplicialLLt * \sa class SimplicialLLT
*/ */
template<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo>
class SimplicialLDLt : public SimplicialCholeskyBase<SimplicialLDLt<_MatrixType,_UpLo> > class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
typedef SimplicialCholeskyBase<SimplicialLDLt> Base; typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
typedef Matrix<Scalar,Dynamic,1> VectorType; typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef internal::traits<SimplicialLDLt> Traits; typedef internal::traits<SimplicialLDLT> Traits;
typedef typename Traits::MatrixL MatrixL; typedef typename Traits::MatrixL MatrixL;
typedef typename Traits::MatrixU MatrixU; typedef typename Traits::MatrixU MatrixU;
public: public:
/** Default constructor */ /** Default constructor */
SimplicialLDLt() : Base() {} SimplicialLDLT() : Base() {}
/** Constructs and performs the LLt factorization of \a matrix */ /** Constructs and performs the LLT factorization of \a matrix */
SimplicialLDLt(const MatrixType& matrix) SimplicialLDLT(const MatrixType& matrix)
: Base(matrix) {} : Base(matrix) {}
/** \returns a vector expression of the diagonal D */ /** \returns a vector expression of the diagonal D */
inline const VectorType vectorD() const { inline const VectorType vectorD() const {
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLt not factorized"); eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
return Base::m_diag; return Base::m_diag;
} }
/** \returns an expression of the factor L */ /** \returns an expression of the factor L */
inline const MatrixL matrixL() const { inline const MatrixL matrixL() const {
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLt not factorized"); eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
return Traits::getL(Base::m_matrix); return Traits::getL(Base::m_matrix);
} }
/** \returns an expression of the factor U (= L^*) */ /** \returns an expression of the factor U (= L^*) */
inline const MatrixU matrixU() const { inline const MatrixU matrixU() const {
eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLt not factorized"); eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
return Traits::getU(Base::m_matrix); return Traits::getU(Base::m_matrix);
} }
@ -486,11 +486,11 @@ public:
} }
}; };
/** \deprecated use SimplicialLDLt or class SimplicialLLt /** \deprecated use SimplicialLDLT or class SimplicialLLT
* \ingroup SparseCholesky_Module * \ingroup SparseCholesky_Module
* \class SimplicialCholesky * \class SimplicialCholesky
* *
* \sa class SimplicialLDLt, class SimplicialLLt * \sa class SimplicialLDLT, class SimplicialLLT
*/ */
template<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo>
class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> > class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
@ -505,13 +505,13 @@ public:
typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
typedef Matrix<Scalar,Dynamic,1> VectorType; typedef Matrix<Scalar,Dynamic,1> VectorType;
typedef internal::traits<SimplicialCholesky> Traits; typedef internal::traits<SimplicialCholesky> Traits;
typedef internal::traits<SimplicialLDLt<MatrixType,UpLo> > LDLtTraits; typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
typedef internal::traits<SimplicialLLt<MatrixType,UpLo> > LLtTraits; typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
public: public:
SimplicialCholesky() : Base(), m_LDLt(true) {} SimplicialCholesky() : Base(), m_LDLT(true) {}
SimplicialCholesky(const MatrixType& matrix) SimplicialCholesky(const MatrixType& matrix)
: Base(), m_LDLt(true) : Base(), m_LDLT(true)
{ {
Base::compute(matrix); Base::compute(matrix);
} }
@ -520,11 +520,11 @@ public:
{ {
switch(mode) switch(mode)
{ {
case SimplicialCholeskyLLt: case SimplicialCholeskyLLT:
m_LDLt = false; m_LDLT = false;
break; break;
case SimplicialCholeskyLDLt: case SimplicialCholeskyLDLT:
m_LDLt = true; m_LDLT = true;
break; break;
default: default:
break; break;
@ -550,7 +550,7 @@ public:
*/ */
void analyzePattern(const MatrixType& a) void analyzePattern(const MatrixType& a)
{ {
Base::analyzePattern(a, m_LDLt); Base::analyzePattern(a, m_LDLT);
} }
/** Performs a numeric decomposition of \a matrix /** Performs a numeric decomposition of \a matrix
@ -561,7 +561,7 @@ public:
*/ */
void factorize(const MatrixType& a) void factorize(const MatrixType& a)
{ {
if(m_LDLt) if(m_LDLT)
Base::template factorize<true>(a); Base::template factorize<true>(a);
else else
Base::template factorize<false>(a); Base::template factorize<false>(a);
@ -584,10 +584,10 @@ public:
if(Base::m_matrix.nonZeros()>0) // otherwise L==I if(Base::m_matrix.nonZeros()>0) // otherwise L==I
{ {
if(m_LDLt) if(m_LDLT)
LDLtTraits::getL(Base::m_matrix).solveInPlace(dest); LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
else else
LLtTraits::getL(Base::m_matrix).solveInPlace(dest); LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
} }
if(Base::m_diag.size()>0) if(Base::m_diag.size()>0)
@ -595,10 +595,10 @@ public:
if (Base::m_matrix.nonZeros()>0) // otherwise I==I if (Base::m_matrix.nonZeros()>0) // otherwise I==I
{ {
if(m_LDLt) if(m_LDLT)
LDLtTraits::getU(Base::m_matrix).solveInPlace(dest); LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
else else
LLtTraits::getU(Base::m_matrix).solveInPlace(dest); LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
} }
if(Base::m_P.size()>0) if(Base::m_P.size()>0)
@ -607,7 +607,7 @@ public:
Scalar determinant() const Scalar determinant() const
{ {
if(m_LDLt) if(m_LDLT)
{ {
return Base::m_diag.prod(); return Base::m_diag.prod();
} }
@ -619,11 +619,11 @@ public:
} }
protected: protected:
bool m_LDLt; bool m_LDLT;
}; };
template<typename Derived> template<typename Derived>
void SimplicialCholeskyBase<Derived>::analyzePattern(const MatrixType& a, bool doLDLt) void SimplicialCholeskyBase<Derived>::analyzePattern(const MatrixType& a, bool doLDLT)
{ {
eigen_assert(a.rows()==a.cols()); eigen_assert(a.rows()==a.cols());
const Index size = a.rows(); const Index size = a.rows();
@ -678,7 +678,7 @@ void SimplicialCholeskyBase<Derived>::analyzePattern(const MatrixType& a, bool d
Index* Lp = m_matrix.outerIndexPtr(); Index* Lp = m_matrix.outerIndexPtr();
Lp[0] = 0; Lp[0] = 0;
for(Index k = 0; k < size; ++k) for(Index k = 0; k < size; ++k)
Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLt ? 0 : 1); Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1);
m_matrix.resizeNonZeros(Lp[size]); m_matrix.resizeNonZeros(Lp[size]);
@ -690,7 +690,7 @@ void SimplicialCholeskyBase<Derived>::analyzePattern(const MatrixType& a, bool d
template<typename Derived> template<typename Derived>
template<bool DoLDLt> template<bool DoLDLT>
void SimplicialCholeskyBase<Derived>::factorize(const MatrixType& a) void SimplicialCholeskyBase<Derived>::factorize(const MatrixType& a)
{ {
eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
@ -711,7 +711,7 @@ void SimplicialCholeskyBase<Derived>::factorize(const MatrixType& a)
ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv); ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_Pinv);
bool ok = true; bool ok = true;
m_diag.resize(DoLDLt ? size : 0); m_diag.resize(DoLDLT ? size : 0);
for(Index k = 0; k < size; ++k) for(Index k = 0; k < size; ++k)
{ {
@ -749,21 +749,21 @@ void SimplicialCholeskyBase<Derived>::factorize(const MatrixType& a)
/* the nonzero entry L(k,i) */ /* the nonzero entry L(k,i) */
Scalar l_ki; Scalar l_ki;
if(DoLDLt) if(DoLDLT)
l_ki = yi / m_diag[i]; l_ki = yi / m_diag[i];
else else
yi = l_ki = yi / Lx[Lp[i]]; yi = l_ki = yi / Lx[Lp[i]];
Index p2 = Lp[i] + m_nonZerosPerCol[i]; Index p2 = Lp[i] + m_nonZerosPerCol[i];
Index p; Index p;
for(p = Lp[i] + (DoLDLt ? 0 : 1); p < p2; ++p) for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p)
y[Li[p]] -= internal::conj(Lx[p]) * yi; y[Li[p]] -= internal::conj(Lx[p]) * yi;
d -= internal::real(l_ki * internal::conj(yi)); d -= internal::real(l_ki * internal::conj(yi));
Li[p] = k; /* store L(k,i) in column form of L */ Li[p] = k; /* store L(k,i) in column form of L */
Lx[p] = l_ki; Lx[p] = l_ki;
++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */
} }
if(DoLDLt) if(DoLDLT)
{ {
m_diag[k] = d; m_diag[k] = d;
if(d == RealScalar(0)) if(d == RealScalar(0))