aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen/src/Cholesky
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen/src/Cholesky')
-rw-r--r--Eigen/src/Cholesky/LDLT.h72
-rw-r--r--Eigen/src/Cholesky/LLT.h60
2 files changed, 74 insertions, 58 deletions
diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h
index 15a92f8fe..15925ba5d 100644
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -27,7 +27,9 @@
#ifndef EIGEN_LDLT_H
#define EIGEN_LDLT_H
+namespace internal {
template<typename MatrixType, int UpLo> struct LDLT_Traits;
+}
/** \ingroup cholesky_Module
*
@@ -74,7 +76,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
- typedef LDLT_Traits<MatrixType,UpLo> Traits;
+ typedef internal::LDLT_Traits<MatrixType,UpLo> Traits;
/** \brief Default Constructor.
*
@@ -108,14 +110,14 @@ template<typename _MatrixType, int _UpLo> class LDLT
/** \returns a view of the upper triangular matrix U */
inline typename Traits::MatrixU matrixU() const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
return Traits::getU(m_matrix);
}
/** \returns a view of the lower triangular matrix L */
inline typename Traits::MatrixL matrixL() const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
return Traits::getL(m_matrix);
}
@@ -123,28 +125,28 @@ template<typename _MatrixType, int _UpLo> class LDLT
*/
inline const TranspositionType& transpositionsP() const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_transpositions;
}
/** \returns the coefficients of the diagonal matrix D */
inline Diagonal<MatrixType,0> vectorD(void) const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix.diagonal();
}
/** \returns true if the matrix is positive (semidefinite) */
inline bool isPositive(void) const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == 1;
}
/** \returns true if the matrix is negative (semidefinite) */
inline bool isNegative(void) const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_sign == -1;
}
@@ -155,13 +157,13 @@ template<typename _MatrixType, int _UpLo> class LDLT
* \sa solveInPlace(), MatrixBase::ldlt()
*/
template<typename Rhs>
- inline const ei_solve_retval<LDLT, Rhs>
+ inline const internal::solve_retval<LDLT, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
- ei_assert(m_matrix.rows()==b.rows()
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_matrix.rows()==b.rows()
&& "LDLT::solve(): invalid number of rows of the right hand side matrix b");
- return ei_solve_retval<LDLT, Rhs>(*this, b.derived());
+ return internal::solve_retval<LDLT, Rhs>(*this, b.derived());
}
template<typename Derived>
@@ -175,7 +177,7 @@ template<typename _MatrixType, int _UpLo> class LDLT
*/
inline const MatrixType& matrixLDLT() const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
return m_matrix;
}
@@ -199,9 +201,11 @@ template<typename _MatrixType, int _UpLo> class LDLT
bool m_isInitialized;
};
-template<int UpLo> struct ei_ldlt_inplace;
+namespace internal {
+
+template<int UpLo> struct ldlt_inplace;
-template<> struct ei_ldlt_inplace<Lower>
+template<> struct ldlt_inplace<Lower>
{
template<typename MatrixType, typename TranspositionType, typename Workspace>
static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
@@ -209,14 +213,14 @@ template<> struct ei_ldlt_inplace<Lower>
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
- ei_assert(mat.rows()==mat.cols());
+ eigen_assert(mat.rows()==mat.cols());
const Index size = mat.rows();
if (size <= 1)
{
transpositions.setIdentity();
if(sign)
- *sign = ei_real(mat.coeff(0,0))>0 ? 1:-1;
+ *sign = real(mat.coeff(0,0))>0 ? 1:-1;
return true;
}
@@ -234,10 +238,10 @@ template<> struct ei_ldlt_inplace<Lower>
// The biggest overall is the point of reference to which further diagonals
// are compared; if any diagonal is negligible compared
// to the largest overall, the algorithm bails.
- cutoff = ei_abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
+ cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner);
if(sign)
- *sign = ei_real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
+ *sign = real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1;
}
// Finish early if the matrix is not full rank.
@@ -259,11 +263,11 @@ template<> struct ei_ldlt_inplace<Lower>
for(int i=k+1;i<index_of_biggest_in_corner;++i)
{
Scalar tmp = mat.coeffRef(i,k);
- mat.coeffRef(i,k) = ei_conj(mat.coeffRef(index_of_biggest_in_corner,i));
- mat.coeffRef(index_of_biggest_in_corner,i) = ei_conj(tmp);
+ mat.coeffRef(i,k) = conj(mat.coeffRef(index_of_biggest_in_corner,i));
+ mat.coeffRef(index_of_biggest_in_corner,i) = conj(tmp);
}
if(NumTraits<Scalar>::IsComplex)
- mat.coeffRef(index_of_biggest_in_corner,k) = ei_conj(mat.coeff(index_of_biggest_in_corner,k));
+ mat.coeffRef(index_of_biggest_in_corner,k) = conj(mat.coeff(index_of_biggest_in_corner,k));
}
// partition the matrix:
@@ -282,7 +286,7 @@ template<> struct ei_ldlt_inplace<Lower>
if(rs>0)
A21.noalias() -= A20 * temp.head(k);
}
- if((rs>0) && (ei_abs(mat.coeffRef(k,k)) > cutoff))
+ if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff))
A21 /= mat.coeffRef(k,k);
}
@@ -290,13 +294,13 @@ template<> struct ei_ldlt_inplace<Lower>
}
};
-template<> struct ei_ldlt_inplace<Upper>
+template<> struct ldlt_inplace<Upper>
{
template<typename MatrixType, typename TranspositionType, typename Workspace>
static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0)
{
Transpose<MatrixType> matt(mat);
- return ei_ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
+ return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign);
}
};
@@ -316,12 +320,14 @@ template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper>
inline static MatrixU getU(const MatrixType& m) { return m; }
};
+} // end namespace internal
+
/** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix
*/
template<typename MatrixType, int _UpLo>
LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
{
- ei_assert(a.rows()==a.cols());
+ eigen_assert(a.rows()==a.cols());
const Index size = a.rows();
m_matrix = a;
@@ -330,22 +336,23 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
m_isInitialized = false;
m_temporary.resize(size);
- ei_ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
+ internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign);
m_isInitialized = true;
return *this;
}
+namespace internal {
template<typename _MatrixType, int _UpLo, typename Rhs>
-struct ei_solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
- : ei_solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
+struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
+ : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs>
{
typedef LDLT<_MatrixType,_UpLo> LDLTType;
EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs)
template<typename Dest> void evalTo(Dest& dst) const
{
- ei_assert(rhs().rows() == dec().matrixLDLT().rows());
+ eigen_assert(rhs().rows() == dec().matrixLDLT().rows());
// dst = P b
dst = dec().transpositionsP() * rhs();
@@ -362,6 +369,7 @@ struct ei_solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
dst = dec().transpositionsP().transpose() * dst;
}
};
+}
/** \internal use x = ldlt_object.solve(x);
*
@@ -380,9 +388,9 @@ template<typename MatrixType,int _UpLo>
template<typename Derived>
bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
const Index size = m_matrix.rows();
- ei_assert(size == bAndX.rows());
+ eigen_assert(size == bAndX.rows());
bAndX = this->solve(bAndX);
@@ -395,7 +403,7 @@ bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
template<typename MatrixType, int _UpLo>
MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
{
- ei_assert(m_isInitialized && "LDLT is not initialized.");
+ eigen_assert(m_isInitialized && "LDLT is not initialized.");
const Index size = m_matrix.rows();
MatrixType res(size,size);
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h
index df135fce3..80a248546 100644
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -25,7 +25,9 @@
#ifndef EIGEN_LLT_H
#define EIGEN_LLT_H
+namespace internal{
template<typename MatrixType, int UpLo> struct LLT_Traits;
+}
/** \ingroup cholesky_Module
*
@@ -68,12 +70,12 @@ template<typename _MatrixType, int _UpLo> class LLT
typedef typename MatrixType::Index Index;
enum {
- PacketSize = ei_packet_traits<Scalar>::size,
+ PacketSize = internal::packet_traits<Scalar>::size,
AlignmentMask = int(PacketSize)-1,
UpLo = _UpLo
};
- typedef LLT_Traits<MatrixType,UpLo> Traits;
+ typedef internal::LLT_Traits<MatrixType,UpLo> Traits;
/**
* \brief Default Constructor.
@@ -102,14 +104,14 @@ template<typename _MatrixType, int _UpLo> class LLT
/** \returns a view of the upper triangular matrix U */
inline typename Traits::MatrixU matrixU() const
{
- ei_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
return Traits::getU(m_matrix);
}
/** \returns a view of the lower triangular matrix L */
inline typename Traits::MatrixL matrixL() const
{
- ei_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
return Traits::getL(m_matrix);
}
@@ -124,13 +126,13 @@ template<typename _MatrixType, int _UpLo> class LLT
* \sa solveInPlace(), MatrixBase::llt()
*/
template<typename Rhs>
- inline const ei_solve_retval<LLT, Rhs>
+ inline const internal::solve_retval<LLT, Rhs>
solve(const MatrixBase<Rhs>& b) const
{
- ei_assert(m_isInitialized && "LLT is not initialized.");
- ei_assert(m_matrix.rows()==b.rows()
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(m_matrix.rows()==b.rows()
&& "LLT::solve(): invalid number of rows of the right hand side matrix b");
- return ei_solve_retval<LLT, Rhs>(*this, b.derived());
+ return internal::solve_retval<LLT, Rhs>(*this, b.derived());
}
template<typename Derived>
@@ -144,7 +146,7 @@ template<typename _MatrixType, int _UpLo> class LLT
*/
inline const MatrixType& matrixLLT() const
{
- ei_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
return m_matrix;
}
@@ -158,7 +160,7 @@ template<typename _MatrixType, int _UpLo> class LLT
*/
ComputationInfo info() const
{
- ei_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
return m_info;
}
@@ -175,9 +177,11 @@ template<typename _MatrixType, int _UpLo> class LLT
ComputationInfo m_info;
};
-template<int UpLo> struct ei_llt_inplace;
+namespace internal {
+
+template<int UpLo> struct llt_inplace;
-template<> struct ei_llt_inplace<Lower>
+template<> struct llt_inplace<Lower>
{
template<typename MatrixType>
static bool unblocked(MatrixType& mat)
@@ -185,7 +189,7 @@ template<> struct ei_llt_inplace<Lower>
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef typename MatrixType::Index Index;
- ei_assert(mat.rows()==mat.cols());
+ eigen_assert(mat.rows()==mat.cols());
const Index size = mat.rows();
for(Index k = 0; k < size; ++k)
{
@@ -195,11 +199,11 @@ template<> struct ei_llt_inplace<Lower>
Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k);
Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k);
- RealScalar x = ei_real(mat.coeff(k,k));
+ RealScalar x = real(mat.coeff(k,k));
if (k>0) x -= mat.row(k).head(k).squaredNorm();
if (x<=RealScalar(0))
return false;
- mat.coeffRef(k,k) = x = ei_sqrt(x);
+ mat.coeffRef(k,k) = x = sqrt(x);
if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
if (rs>0) A21 *= RealScalar(1)/x;
}
@@ -210,7 +214,7 @@ template<> struct ei_llt_inplace<Lower>
static bool blocked(MatrixType& m)
{
typedef typename MatrixType::Index Index;
- ei_assert(m.rows()==m.cols());
+ eigen_assert(m.rows()==m.cols());
Index size = m.rows();
if(size<32)
return unblocked(m);
@@ -239,19 +243,19 @@ template<> struct ei_llt_inplace<Lower>
}
};
-template<> struct ei_llt_inplace<Upper>
+template<> struct llt_inplace<Upper>
{
template<typename MatrixType>
static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat)
{
Transpose<MatrixType> matt(mat);
- return ei_llt_inplace<Lower>::unblocked(matt);
+ return llt_inplace<Lower>::unblocked(matt);
}
template<typename MatrixType>
static EIGEN_STRONG_INLINE bool blocked(MatrixType& mat)
{
Transpose<MatrixType> matt(mat);
- return ei_llt_inplace<Lower>::blocked(matt);
+ return llt_inplace<Lower>::blocked(matt);
}
};
@@ -262,7 +266,7 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Lower>
inline static MatrixL getL(const MatrixType& m) { return m; }
inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); }
static bool inplace_decomposition(MatrixType& m)
- { return ei_llt_inplace<Lower>::blocked(m); }
+ { return llt_inplace<Lower>::blocked(m); }
};
template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
@@ -272,9 +276,11 @@ template<typename MatrixType> struct LLT_Traits<MatrixType,Upper>
inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); }
inline static MatrixU getU(const MatrixType& m) { return m; }
static bool inplace_decomposition(MatrixType& m)
- { return ei_llt_inplace<Upper>::blocked(m); }
+ { return llt_inplace<Upper>::blocked(m); }
};
+} // end namespace internal
+
/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix
*
*
@@ -295,9 +301,10 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
return *this;
}
+namespace internal {
template<typename _MatrixType, int UpLo, typename Rhs>
-struct ei_solve_retval<LLT<_MatrixType, UpLo>, Rhs>
- : ei_solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
+struct solve_retval<LLT<_MatrixType, UpLo>, Rhs>
+ : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs>
{
typedef LLT<_MatrixType,UpLo> LLTType;
EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs)
@@ -308,6 +315,7 @@ struct ei_solve_retval<LLT<_MatrixType, UpLo>, Rhs>
dec().solveInPlace(dst);
}
};
+}
/** \internal use x = llt_object.solve(x);
*
@@ -326,8 +334,8 @@ template<typename MatrixType, int _UpLo>
template<typename Derived>
void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
{
- ei_assert(m_isInitialized && "LLT is not initialized.");
- ei_assert(m_matrix.rows()==bAndX.rows());
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(m_matrix.rows()==bAndX.rows());
matrixL().solveInPlace(bAndX);
matrixU().solveInPlace(bAndX);
}
@@ -338,7 +346,7 @@ void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
template<typename MatrixType, int _UpLo>
MatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const
{
- ei_assert(m_isInitialized && "LLT is not initialized.");
+ eigen_assert(m_isInitialized && "LLT is not initialized.");
return matrixL() * matrixL().adjoint().toDenseMatrix();
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-17 17:04:03 +0000