move PartialLU to the new API

Eigen/src/LU/LU.h

@@ -209,7 +209,7 @@ template<typename MatrixType> class LU
       *
       * \returns a solution.
       *
-      * \note_about_inexistant_solutions
+      * \note_about_checking_solutions
       *
       * \note_about_arbitrary_choice_of_solution
       * \note_about_using_kernel_to_study_multiple_solutions
@@ -374,15 +374,12 @@ template<typename MatrixType> class LU
     }
 
   protected:
-    bool m_isInitialized;
     MatrixType m_lu;
     IntColVectorType m_p;
     IntRowVectorType m_q;
-    int m_det_pq;
+    int m_det_pq, m_nonzero_pivots;
-    int m_nonzero_pivots;
+    RealScalar m_maxpivot, m_prescribedThreshold;
-    RealScalar m_maxpivot;
+    bool m_isInitialized, m_usePrescribedThreshold;
-    bool m_usePrescribedThreshold;
-    RealScalar m_prescribedThreshold;
 };
 
 template<typename MatrixType>

Eigen/src/LU/PartialLU.h

@@ -26,6 +26,8 @@
 #ifndef EIGEN_PARTIALLU_H
 #define EIGEN_PARTIALLU_H
 
+template<typename MatrixType, typename Rhs> struct ei_partiallu_solve_impl;
+
 /** \ingroup LU_Module
   *
   * \class PartialLU
@@ -112,26 +114,46 @@ template<typename MatrixType> class PartialLU
       return m_p;
     }
 
-    /** This method finds the solution x to the equation Ax=b, where A is the matrix of which
+    /** This method returns a solution x to the equation Ax=b, where A is the matrix of which
-      * *this is the LU decomposition. Since if this partial pivoting decomposition the matrix is assumed
+      * *this is the LU decomposition.
-      * to have full rank, such a solution is assumed to exist and to be unique.
-      *
-      * \warning Again, if your matrix may not have full rank, use class LU instead. See LU::solve().
       *
       * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
       *          the only requirement in order for the equation to make sense is that
       *          b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
-      * \param result a pointer to the vector or matrix in which to store the solution, if any exists.
+      *
-      *          Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
+      * \returns the solution.
-      *          If no solution exists, *result is left with undefined coefficients.
       *
       * Example: \include PartialLU_solve.cpp
       * Output: \verbinclude PartialLU_solve.out
       *
+      * Since this PartialLU 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 OtherDerived, typename ResultType>
+    template<typename Rhs>
-    void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+    inline const ei_partiallu_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());
+    }
+
+    /** \returns the inverse of the matrix of which *this is the LU decomposition.
+      *
+      * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
+      *          invertibility, use class LU instead.
+      *
+      * \sa MatrixBase::inverse(), LU::inverse()
+      */
+    inline const ei_partiallu_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> >
+               (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue());
+    }
 
     /** \returns the determinant of the matrix of which
       * *this is the LU decomposition. It has only linear complexity
@@ -148,34 +170,6 @@ template<typename MatrixType> class PartialLU
       */
     typename ei_traits<MatrixType>::Scalar determinant() const;
 
-    /** Computes the inverse of the matrix of which *this is the LU decomposition.
-      *
-      * \param result a pointer to the matrix into which to store the inverse. Resized if needed.
-      *
-      * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
-      *          invertibility, use class LU instead.
-      *
-      * \sa MatrixBase::computeInverse(), inverse()
-      */
-    inline void computeInverse(MatrixType *result) const
-    {
-      solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result);
-    }
-
-    /** \returns the inverse of the matrix of which *this is the LU decomposition.
-      *
-      * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
-      *          invertibility, use class LU instead.
-      *
-      * \sa computeInverse(), MatrixBase::inverse()
-      */
-    inline MatrixType inverse() const
-    {
-      MatrixType result;
-      computeInverse(&result);
-      return result;
-    }
-
   protected:
     MatrixType m_lu;
     IntColVectorType m_p;
@@ -411,36 +405,60 @@ typename ei_traits<MatrixType>::Scalar PartialLU<MatrixType>::determinant() cons
   return Scalar(m_det_p) * m_lu.diagonal().prod();
 }
 
-template<typename MatrixType>
+/***** Implementation of solve() *****************************************************/
-template<typename OtherDerived, typename ResultType>
+
-void PartialLU<MatrixType>::solve(
+template<typename MatrixType,typename Rhs>
-  const MatrixBase<OtherDerived>& b,
+struct ei_traits<ei_partiallu_solve_impl<MatrixType,Rhs> >
-  ResultType *result
-) const
 {
-  ei_assert(m_isInitialized && "PartialLU is not initialized.");
+  typedef Matrix<typename Rhs::Scalar,
+                 MatrixType::ColsAtCompileTime,
+                 Rhs::ColsAtCompileTime,
+                 Rhs::PlainMatrixType::Options,
+                 MatrixType::MaxColsAtCompileTime,
+                 Rhs::MaxColsAtCompileTime> ReturnMatrixType;
+};
 
-  /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
+template<typename MatrixType, typename Rhs>
-   * So we proceed as follows:
+struct ei_partiallu_solve_impl : public ReturnByValue<ei_partiallu_solve_impl<MatrixType, Rhs> >
-   * Step 1: compute c = Pb.
+{
-   * Step 2: replace c by the solution x to Lx = c.
+  typedef typename ei_cleantype<typename Rhs::Nested>::type RhsNested;
-   * Step 3: replace c by the solution x to Ux = c.
+  typedef PartialLU<MatrixType> LUType;
-   */
+  const LUType& m_lu;
+  const typename Rhs::Nested m_rhs;
 
-  const int size = m_lu.rows();
+  ei_partiallu_solve_impl(const LUType& lu, const Rhs& rhs)
-  ei_assert(b.rows() == size);
+    : m_lu(lu), m_rhs(rhs)
+  {}
 
-  result->resize(size, b.cols());
+  inline int rows() const { return m_lu.matrixLU().cols(); }
+  inline int cols() const { return m_rhs.cols(); }
 
-  // Step 1
+  template<typename Dest> void evalTo(Dest& dst) const
-  for(int i = 0; i < size; ++i) result->row(m_p.coeff(i)) = b.row(i);
+  {
+    /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
+    * So we proceed as follows:
+    * Step 1: compute c = Pb.
+    * Step 2: replace c by the solution x to Lx = c.
+    * Step 3: replace c by the solution x to Ux = c.
+    */
 
-  // Step 2
+    const int size = m_lu.matrixLU().rows();
-  m_lu.template triangularView<UnitLowerTriangular>().solveInPlace(*result);
+    ei_assert(m_rhs.rows() == size);
 
-  // Step 3
+    dst.resize(size, m_rhs.cols());
-  m_lu.template triangularView<UpperTriangular>().solveInPlace(*result);
+
-}
+    // Step 1
+    for(int i = 0; i < size; ++i) dst.row(m_lu.permutationP().coeff(i)) = m_rhs.row(i);
+
+    // Step 2
+    m_lu.matrixLU().template triangularView<UnitLowerTriangular>().solveInPlace(dst);
+
+    // Step 3
+    m_lu.matrixLU().template triangularView<UpperTriangular>().solveInPlace(dst);
+  }
+};
+
+/******** MatrixBase methods *******/
 
 /** \lu_module
   *

doc/Doxyfile.in

@@ -218,7 +218,7 @@ ALIASES                = "only_for_vectors=This is only for vectors (either row-
                          "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" \
                          "note_about_arbitrary_choice_of_solution=If there exists more than one solution, this method will arbitrarily choose one." \
                          "note_about_using_kernel_to_study_multiple_solutions=If you need a complete analysis of the space of solutions, take the one solution obtained by this method and add to it elements of the kernel, as determined by kernel()." \
-                         "note_about_inexistant_solutions=If only an approximate solution exists, then the result is only such an approximate solution. The case where no solution exists isn't treated separately: non-existent solutions are just \"very inaccurate solutions\". If you want to check whether a solution exists, just call this function to get a solution and then compute the error margin, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode As long as the input matrices have finite numeric coefficients (no \c inf or \c nan values), the result will also have finite numeric coefficients, there is no risk of getting \c nan values if no solution exists."
+                         "note_about_checking_solutions=This method just tries to find as good a solution as possible. If you want to check whether a solution "exists" or if it is accurate, just call this function to get a solution and then compute the error margin, or use MatrixBase::isApprox() directly, for instance like this: \code bool a_solution_exists = (A*result).isApprox(b, precision); \endcode The non-existence of a solution doesn't by itself mean that you'll get \c inf or \c nan values."
 
 # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C
 # sources only. Doxygen will then generate output that is more tailored for C.