*port the Cholesky module to the new solve() API

*improve documentation
2026-04-10 11:34:33 +08:00 · 2009-10-29 21:11:05 -04:00
parent e8dd552257
commit 6b48e932e9
8 changed files with 169 additions and 97 deletions
--- a/Eigen/src/Cholesky/LDLT.h
+++ b/Eigen/src/Cholesky/LDLT.h
@@ -27,6 +27,8 @@
 #ifndef EIGEN_LDLT_H
 #define EIGEN_LDLT_H

+template<typename MatrixType, typename Rhs> struct ei_ldlt_solve_impl;
+
 /** \ingroup cholesky_Module
  *
  * \class LDLT
@@ -43,8 +45,8 @@
  * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root
  * on D also stabilizes the computation.
  *
-  * Remember that Cholesky decompositions are not rank-revealing.  Also, do not use a Cholesky decomposition to determine
-  * whether a system of equations has a solution.
+  * Remember that Cholesky decompositions are not rank-revealing.  Also, do not use a Cholesky
+	* decomposition to determine whether a system of equations has a solution.
  *
  * \sa MatrixBase::ldlt(), class LLT
  */
@@ -117,14 +119,37 @@ template<typename MatrixType> class LDLT
      return m_sign == -1;
    }

-    template<typename RhsDerived, typename ResultType>
-    bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const;
-
+    /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A.
+      *
+      * \note_about_checking_solutions
+      *
+      * \sa solveInPlace(), MatrixBase::ldlt()
+      */
+    template<typename Rhs>
+    inline const ei_ldlt_solve_impl<MatrixType, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(m_isInitialized && "LDLT is not initialized.");
+      ei_assert(m_matrix.rows()==b.rows()
+                && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
+      return ei_ldlt_solve_impl<MatrixType, Rhs>(*this, b.derived());
+    }
+    
    template<typename Derived>
    bool solveInPlace(MatrixBase<Derived> &bAndX) const;

    LDLT& compute(const MatrixType& matrix);

+    /** \returns the LDLT decomposition matrix
+      *
+      * TODO: document the storage layout
+      */
+    inline const MatrixType& matrixLDLT() const
+    {
+      ei_assert(m_isInitialized && "LDLT is not initialized.");
+      return m_matrix;
+    }
+    
  protected:
    /** \internal
      * Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U.
@@ -134,7 +159,7 @@ template<typename MatrixType> class LDLT
      */
    MatrixType m_matrix;
    IntColVectorType m_p;
-    IntColVectorType m_transpositions;
+    IntColVectorType m_transpositions; // FIXME do we really need to store permanently the transpositions?
    int m_sign;
    bool m_isInitialized;
 };
@@ -238,27 +263,38 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a)
  return *this;
 }

-/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
-  * The result is stored in \a result
-  *
-  * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
-  *
-  * In other words, it computes \f$ b = A^{-1} b \f$ with
-  * \f$ P^T{L^{*}}^{-1} D^{-1} L^{-1} P b \f$ from right to left.
-  *
-  * \sa LDLT::solveInPlace(), MatrixBase::ldlt()
-  */
-template<typename MatrixType>
-template<typename RhsDerived, typename ResultType>
-bool LDLT<MatrixType>
-::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const
+template<typename MatrixType,typename Rhs>
+struct ei_traits<ei_ldlt_solve_impl<MatrixType,Rhs> >
 {
-  ei_assert(m_isInitialized && "LDLT is not initialized.");
-  const int size = m_matrix.rows();
-  ei_assert(size==b.rows() && "LDLT::solve(): invalid number of rows of the right hand side matrix b");
-  *result = b;
-  return solveInPlace(*result);
-}
+  typedef Matrix<typename Rhs::Scalar,
+                 MatrixType::ColsAtCompileTime,
+                 Rhs::ColsAtCompileTime,
+                 Rhs::PlainMatrixType::Options,
+                 MatrixType::MaxColsAtCompileTime,
+                 Rhs::MaxColsAtCompileTime> ReturnMatrixType;
+};
+
+template<typename MatrixType, typename Rhs>
+struct ei_ldlt_solve_impl : public ReturnByValue<ei_ldlt_solve_impl<MatrixType, Rhs> >
+{
+  typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
+  typedef LDLT<MatrixType> LDLTType;
+  const LDLTType& m_ldlt;
+  const typename Rhs::Nested m_rhs;
+
+  ei_ldlt_solve_impl(const LDLTType& ldlt, const Rhs& rhs)
+    : m_ldlt(ldlt), m_rhs(rhs)
+  {}
+
+  inline int rows() const { return m_ldlt.matrixLDLT().cols(); }
+  inline int cols() const { return m_rhs.cols(); }
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    dst = m_rhs;
+    m_ldlt.solveInPlace(dst);
+  }
+};

 /** This is the \em in-place version of solve().
  *
--- a/Eigen/src/Cholesky/LLT.h
+++ b/Eigen/src/Cholesky/LLT.h
@@ -26,6 +26,7 @@
 #define EIGEN_LLT_H

 template<typename MatrixType, int UpLo> struct LLT_Traits;
+template<typename MatrixType, int UpLo, typename Rhs> struct ei_llt_solve_impl;

 /** \ingroup cholesky_Module
  *
@@ -99,14 +100,41 @@ template<typename MatrixType, int _UpLo> class LLT
      return Traits::getL(m_matrix);
    }

-    template<typename RhsDerived, typename ResultType>
-    bool solve(const MatrixBase<RhsDerived> &b, ResultType *result) const;
-
+    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
+      *
+      * Since this LLT class assumes anyway that the matrix A is invertible, the solution
+      * theoretically exists and is unique regardless of b.
+      *
+      * Example: \include LLT_solve.cpp
+      * Output: \verbinclude LLT_solve.out
+      *
+      * \sa solveInPlace(), MatrixBase::llt()
+      */
+    template<typename Rhs>
+    inline const ei_llt_solve_impl<MatrixType, UpLo, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(m_isInitialized && "LLT is not initialized.");
+      ei_assert(m_matrix.rows()==b.rows()
+                && "LLT::solve(): invalid number of rows of the right hand side matrix b");
+      return ei_llt_solve_impl<MatrixType, UpLo, Rhs>(*this, b.derived());
+    }
+    
    template<typename Derived>
    bool solveInPlace(MatrixBase<Derived> &bAndX) const;

    LLT& compute(const MatrixType& matrix);

+    /** \returns the LLT decomposition matrix
+      *
+      * TODO: document the storage layout
+      */
+    inline const MatrixType& matrixLLT() const
+    {
+      ei_assert(m_isInitialized && "LLT is not initialized.");
+      return m_matrix;
+    }
+    
  protected:
    /** \internal
      * Used to compute and store L
@@ -229,28 +257,38 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a)
  return *this;
 }

-/** Computes the solution x of \f$ A x = b \f$ using the current decomposition of A.
-  * The result is stored in \a result
-  *
-  * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD.
-  *
-  * In other words, it computes \f$ b = A^{-1} b \f$ with
-  * \f$ {L^{*}}^{-1} L^{-1} b \f$ from right to left.
-  *
-  * Example: \include LLT_solve.cpp
-  * Output: \verbinclude LLT_solve.out
-  *
-  * \sa LLT::solveInPlace(), MatrixBase::llt()
-  */
-template<typename MatrixType, int _UpLo>
-template<typename RhsDerived, typename ResultType>
-bool LLT<MatrixType,_UpLo>::solve(const MatrixBase<RhsDerived> &b, ResultType *result) const
+template<typename MatrixType, int UpLo, typename Rhs>
+struct ei_traits<ei_llt_solve_impl<MatrixType,UpLo,Rhs> >
 {
-  ei_assert(m_isInitialized && "LLT is not initialized.");
-  const int size = m_matrix.rows();
-  ei_assert(size==b.rows() && "LLT::solve(): invalid number of rows of the right hand side matrix b");
-  return solveInPlace((*result) = b);
-}
+  typedef Matrix<typename Rhs::Scalar,
+                 MatrixType::ColsAtCompileTime,
+                 Rhs::ColsAtCompileTime,
+                 Rhs::PlainMatrixType::Options,
+                 MatrixType::MaxColsAtCompileTime,
+                 Rhs::MaxColsAtCompileTime> ReturnMatrixType;
+};
+
+template<typename MatrixType, int UpLo, typename Rhs>
+struct ei_llt_solve_impl : public ReturnByValue<ei_llt_solve_impl<MatrixType,UpLo,Rhs> >
+{
+  typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
+  typedef LLT<MatrixType,UpLo> LLTType;
+  const LLTType& m_llt;
+  const typename Rhs::Nested m_rhs;
+
+  ei_llt_solve_impl(const LLTType& llt, const Rhs& rhs)
+    : m_llt(llt), m_rhs(rhs)
+  {}
+
+  inline int rows() const { return m_llt.matrixLLT().cols(); }
+  inline int cols() const { return m_rhs.cols(); }
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    dst = m_rhs;
+    m_llt.solveInPlace(dst);
+  }
+};

 /** This is the \em in-place version of solve().
  *
--- a/Eigen/src/LU/FullPivLU.h
+++ b/Eigen/src/LU/FullPivLU.h
@@ -25,9 +25,9 @@
 #ifndef EIGEN_LU_H
 #define EIGEN_LU_H

-template<typename MatrixType, typename Rhs> struct ei_lu_solve_impl;
-template<typename MatrixType> struct ei_lu_kernel_impl;
-template<typename MatrixType> struct ei_lu_image_impl;
+template<typename MatrixType, typename Rhs> struct ei_fullpivlu_solve_impl;
+template<typename MatrixType> struct ei_fullpivlu_kernel_impl;
+template<typename MatrixType> struct ei_fullpivlu_image_impl;

 /** \ingroup LU_Module
  *
@@ -167,10 +167,10 @@ template<typename MatrixType> class FullPivLU
      *
      * \sa image()
      */
-    inline const ei_lu_kernel_impl<MatrixType> kernel() const
+    inline const ei_fullpivlu_kernel_impl<MatrixType> kernel() const
    {
      ei_assert(m_isInitialized && "LU is not initialized.");
-      return ei_lu_kernel_impl<MatrixType>(*this);
+      return ei_fullpivlu_kernel_impl<MatrixType>(*this);
    }

    /** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
@@ -193,11 +193,11 @@ template<typename MatrixType> class FullPivLU
      * \sa kernel()
      */
    template<typename OriginalMatrixType>
-    inline const ei_lu_image_impl<MatrixType>
+    inline const ei_fullpivlu_image_impl<MatrixType>
      image(const MatrixBase<OriginalMatrixType>& originalMatrix) const
    {
      ei_assert(m_isInitialized && "LU is not initialized.");
-      return ei_lu_image_impl<MatrixType>(*this, originalMatrix.derived());
+      return ei_fullpivlu_image_impl<MatrixType>(*this, originalMatrix.derived());
    }

    /** This method returns a solution x to the equation Ax=b, where A is the matrix of which
@@ -220,11 +220,11 @@ template<typename MatrixType> class FullPivLU
      * \sa TriangularView::solve(), kernel(), inverse()
      */
    template<typename Rhs>
-    inline const ei_lu_solve_impl<MatrixType, Rhs>
+    inline const ei_fullpivlu_solve_impl<MatrixType, Rhs>
    solve(const MatrixBase<Rhs>& b) const
    {
      ei_assert(m_isInitialized && "LU is not initialized.");
-      return ei_lu_solve_impl<MatrixType, Rhs>(*this, b.derived());
+      return ei_fullpivlu_solve_impl<MatrixType, Rhs>(*this, b.derived());
    }

    /** \returns the determinant of the matrix of which
@@ -365,11 +365,11 @@ template<typename MatrixType> class FullPivLU
      *
      * \sa MatrixBase::inverse()
      */
-    inline const ei_lu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const
+    inline const ei_fullpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const
    {
      ei_assert(m_isInitialized && "LU is not initialized.");
      ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
-      return ei_lu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> >
+      return ei_fullpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> >
               (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue());
    }

@@ -493,7 +493,7 @@ typename ei_traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() cons
 /********* Implementation of kernel() **************************************************/

 template<typename MatrixType>
-struct ei_traits<ei_lu_kernel_impl<MatrixType> >
+struct ei_traits<ei_fullpivlu_kernel_impl<MatrixType> >
 {
  typedef Matrix<
    typename MatrixType::Scalar,
@@ -509,7 +509,7 @@ struct ei_traits<ei_lu_kernel_impl<MatrixType> >
 };

 template<typename MatrixType>
-struct ei_lu_kernel_impl : public ReturnByValue<ei_lu_kernel_impl<MatrixType> >
+struct ei_fullpivlu_kernel_impl : public ReturnByValue<ei_fullpivlu_kernel_impl<MatrixType> >
 {
  typedef FullPivLU<MatrixType> LUType;
  typedef typename MatrixType::Scalar Scalar;
@@ -517,7 +517,7 @@ struct ei_lu_kernel_impl : public ReturnByValue<ei_lu_kernel_impl<MatrixType> >
  const LUType& m_lu;
  int m_rank, m_cols;
  
-  ei_lu_kernel_impl(const LUType& lu)
+  ei_fullpivlu_kernel_impl(const LUType& lu)
    : m_lu(lu),
      m_rank(lu.rank()),
      m_cols(m_rank==lu.matrixLU().cols() ? 1 : lu.matrixLU().cols() - m_rank){}
@@ -599,7 +599,7 @@ struct ei_lu_kernel_impl : public ReturnByValue<ei_lu_kernel_impl<MatrixType> >
 /***** Implementation of image() *****************************************************/

 template<typename MatrixType>
-struct ei_traits<ei_lu_image_impl<MatrixType> >
+struct ei_traits<ei_fullpivlu_image_impl<MatrixType> >
 {
  typedef Matrix<
    typename MatrixType::Scalar,
@@ -613,7 +613,7 @@ struct ei_traits<ei_lu_image_impl<MatrixType> >
 };

 template<typename MatrixType>
-struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> >
+struct ei_fullpivlu_image_impl : public ReturnByValue<ei_fullpivlu_image_impl<MatrixType> >
 {
  typedef FullPivLU<MatrixType> LUType;
  typedef typename MatrixType::RealScalar RealScalar;
@@ -621,7 +621,7 @@ struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> >
  int m_rank, m_cols;
  const MatrixType& m_originalMatrix;
  
-  ei_lu_image_impl(const LUType& lu, const MatrixType& originalMatrix)
+  ei_fullpivlu_image_impl(const LUType& lu, const MatrixType& originalMatrix)
    : m_lu(lu), m_rank(lu.rank()),
      m_cols(m_rank == 0 ? 1 : m_rank),
      m_originalMatrix(originalMatrix) {}
@@ -656,7 +656,7 @@ struct ei_lu_image_impl : public ReturnByValue<ei_lu_image_impl<MatrixType> >
 /***** Implementation of solve() *****************************************************/

 template<typename MatrixType,typename Rhs>
-struct ei_traits<ei_lu_solve_impl<MatrixType,Rhs> >
+struct ei_traits<ei_fullpivlu_solve_impl<MatrixType,Rhs> >
 {
  typedef Matrix<typename Rhs::Scalar,
                 MatrixType::ColsAtCompileTime,
@@ -667,14 +667,14 @@ struct ei_traits<ei_lu_solve_impl<MatrixType,Rhs> >
 };

 template<typename MatrixType, typename Rhs>
-struct ei_lu_solve_impl : public ReturnByValue<ei_lu_solve_impl<MatrixType, Rhs> >
+struct ei_fullpivlu_solve_impl : public ReturnByValue<ei_fullpivlu_solve_impl<MatrixType, Rhs> >
 {
  typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
  typedef FullPivLU<MatrixType> LUType;
  const LUType& m_lu;
  const typename Rhs::Nested m_rhs;
  
-  ei_lu_solve_impl(const LUType& lu, const Rhs& rhs)
+  ei_fullpivlu_solve_impl(const LUType& lu, const Rhs& rhs)
    : m_lu(lu), m_rhs(rhs)
  {}

--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -26,7 +26,7 @@
 #ifndef EIGEN_PARTIALLU_H
 #define EIGEN_PARTIALLU_H

-template<typename MatrixType, typename Rhs> struct ei_partiallu_solve_impl;
+template<typename MatrixType, typename Rhs> struct ei_partialpivlu_solve_impl;

 /** \ingroup LU_Module
  *
@@ -116,7 +116,7 @@ template<typename MatrixType> class PartialPivLU
      return m_p;
    }

-    /** This method returns a solution x to the equation Ax=b, where A is the matrix of which
+    /** This method returns the solution x to the equation Ax=b, where A is the matrix of which
      * *this is the LU decomposition.
      *
      * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
@@ -131,16 +131,14 @@ template<typename MatrixType> class PartialPivLU
      * Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution
      * theoretically exists and is unique regardless of b.
      *
-      * \note_about_checking_solutions
-      *
      * \sa TriangularView::solve(), inverse(), computeInverse()
      */
    template<typename Rhs>
-    inline const ei_partiallu_solve_impl<MatrixType, Rhs>
+    inline const ei_partialpivlu_solve_impl<MatrixType, Rhs>
    solve(const MatrixBase<Rhs>& b) const
    {
-      ei_assert(m_isInitialized && "LU is not initialized.");
-      return ei_partiallu_solve_impl<MatrixType, Rhs>(*this, b.derived());
+      ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
+      return ei_partialpivlu_solve_impl<MatrixType, Rhs>(*this, b.derived());
    }

    /** \returns the inverse of the matrix of which *this is the LU decomposition.
@@ -150,10 +148,10 @@ template<typename MatrixType> class PartialPivLU
      *
      * \sa MatrixBase::inverse(), LU::inverse()
      */
-    inline const ei_partiallu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const
+    inline const ei_partialpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const
    {
-      ei_assert(m_isInitialized && "LU is not initialized.");
-      return ei_partiallu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> >
+      ei_assert(m_isInitialized && "PartialPivLU is not initialized.");
+      return ei_partialpivlu_solve_impl<MatrixType,NestByValue<typename MatrixType::IdentityReturnType> >
               (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue());
    }

@@ -200,7 +198,7 @@ PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)



-/** This is the blocked version of ei_lu_unblocked() */
+/** This is the blocked version of ei_fullpivlu_unblocked() */
 template<typename Scalar, int StorageOrder>
 struct ei_partial_lu_impl
 {
@@ -410,7 +408,7 @@ typename ei_traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() c
 /***** Implementation of solve() *****************************************************/

 template<typename MatrixType,typename Rhs>
-struct ei_traits<ei_partiallu_solve_impl<MatrixType,Rhs> >
+struct ei_traits<ei_partialpivlu_solve_impl<MatrixType,Rhs> >
 {
  typedef Matrix<typename Rhs::Scalar,
                 MatrixType::ColsAtCompileTime,
@@ -421,14 +419,14 @@ struct ei_traits<ei_partiallu_solve_impl<MatrixType,Rhs> >
 };

 template<typename MatrixType, typename Rhs>
-struct ei_partiallu_solve_impl : public ReturnByValue<ei_partiallu_solve_impl<MatrixType, Rhs> >
+struct ei_partialpivlu_solve_impl : public ReturnByValue<ei_partialpivlu_solve_impl<MatrixType, Rhs> >
 {
  typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
  typedef PartialPivLU<MatrixType> LUType;
  const LUType& m_lu;
  const typename Rhs::Nested m_rhs;

-  ei_partiallu_solve_impl(const LUType& lu, const Rhs& rhs)
+  ei_partialpivlu_solve_impl(const LUType& lu, const Rhs& rhs)
    : m_lu(lu), m_rhs(rhs)
  {}