Enable inverse() method for computing pseudo-inverse.

2 files changed, 26 insertions, 18 deletions

diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
index bee5bf47e..9bc768b7c 100644
--- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h
+++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
@@ -70,6 +70,7 @@ class CompleteOrthogonalDecomposition {
       MatrixType, typename internal::remove_all<
                       typename HCoeffsType::ConjugateReturnType>::type>
       HouseholderSequenceType;
+  typedef typename MatrixType::PlainObject PlainObject;
 
  private:
   typedef typename PermutationType::Index PermIndexType;
@@ -246,26 +247,15 @@ class CompleteOrthogonalDecomposition {
    */
   inline bool isInvertible() const { return m_cpqr.isInvertible(); }
 
-  /** \returns the inverse of the matrix of which *this is the complete
+  /** \returns the pseudo-inverse of the matrix of which *this is the complete
    * orthogonal decomposition.
-   *
-   * \note If this matrix is not invertible, the returned matrix has undefined
-   * coefficients. Use isInvertible() to first determine whether this matrix is
-   * invertible.
+   * \warning: Do not compute \c this->inverse()*rhs to solve a linear systems.
+   * It is more efficient and numerically stable to call \c this->solve(rhs).
    */
-
-  // TODO(rmlarsen): Add method for pseudo-inverse.
-  // inline const
-  // internal::solve_retval<CompleteOrthogonalDecomposition, typename
-  // MatrixType::IdentityReturnType>
-  // inverse() const
-  // {
-  //   eigen_assert(m_isInitialized && "CompleteOrthogonalDecomposition is not
-  //   initialized.");
-  //   return internal::solve_retval<CompleteOrthogonalDecomposition,typename
-  //   MatrixType::IdentityReturnType>
-  //            (*this, MatrixType::Identity(m_cpqr.rows(), m_cpqr.cols()));
-  // }
+  inline const Inverse<CompleteOrthogonalDecomposition> inverse() const
+  {
+    return Inverse<CompleteOrthogonalDecomposition>(*this);
+  }
 
   inline Index rows() const { return m_cpqr.rows(); }
   inline Index cols() const { return m_cpqr.cols(); }
@@ -513,6 +503,21 @@ void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl(
 }
 #endif
 
+namespace internal {
+
+template<typename DstXprType, typename MatrixType, typename Scalar>
+struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
+{
+  typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
+  typedef Inverse<CodType> SrcXprType;
+  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
+  {
+    dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows()));
+  }
+};
+
+} // end namespace internal
+
 /** \returns the matrix Q as a sequence of householder transformations */
 template <typename MatrixType>
 typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index 16f80d8b5..d8672ce33 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -57,6 +57,9 @@ void cod() {
   JacobiSVD<MatrixType> svd(matrix, ComputeThinU | ComputeThinV);
   MatrixType svd_solution = svd.solve(rhs);
   VERIFY_IS_APPROX(cod_solution, svd_solution);
+
+  MatrixType pinv = cod.inverse();
+  VERIFY_IS_APPROX(cod_solution, pinv * rhs);
 }
 
 template <typename MatrixType, int Cols2>