From 8496f86f84a78e1d8b2424ad5048c5fdab37f58f Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Wed, 13 Nov 2019 21:16:53 +0100
Subject: [PATCH] Enable CompleteOrthogonalDecomposition::pseudoInverse with
 non-square fixed-size matrices.

---
 Eigen/src/Core/Inverse.h                       |  7 +++----
 Eigen/src/QR/CompleteOrthogonalDecomposition.h | 10 +++++++++-
 test/permutationmatrices.cpp                   |  6 +++++-
 test/qr_colpivoting.cpp                        |  6 +++++-
 4 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/Eigen/src/Core/Inverse.h b/Eigen/src/Core/Inverse.h
index b76f0439d..7352d8037 100644
--- a/Eigen/src/Core/Inverse.h
+++ b/Eigen/src/Core/Inverse.h
@@ -1,7 +1,7 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2014-2019 Gael Guennebaud <gael.guennebaud@inria.fr>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -44,7 +44,6 @@ class Inverse : public InverseImpl<XprType,typename internal::traits<XprType>::S
 {
 public:
   typedef typename XprType::StorageIndex StorageIndex;
-  typedef typename XprType::PlainObject                       PlainObject;
   typedef typename XprType::Scalar                            Scalar;
   typedef typename internal::ref_selector<XprType>::type      XprTypeNested;
   typedef typename internal::remove_all<XprTypeNested>::type  XprTypeNestedCleaned;
@@ -55,8 +54,8 @@ public:
     : m_xpr(xpr)
   {}
 
-  EIGEN_DEVICE_FUNC Index rows() const { return m_xpr.rows(); }
-  EIGEN_DEVICE_FUNC Index cols() const { return m_xpr.cols(); }
+  EIGEN_DEVICE_FUNC Index rows() const { return m_xpr.cols(); }
+  EIGEN_DEVICE_FUNC Index cols() const { return m_xpr.rows(); }
 
   EIGEN_DEVICE_FUNC const XprTypeNestedCleaned& nestedExpression() const { return m_xpr; }
 
diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
index 2fc3c871a..486d3373a 100644
--- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h
+++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
@@ -592,6 +592,13 @@ void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl_transposed(const
 
 namespace internal {
 
+template<typename MatrixType>
+struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > >
+  : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
+{
+  enum { Flags = 0 };
+};
+
 template<typename DstXprType, typename MatrixType>
 struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
 {
@@ -599,7 +606,8 @@ struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType
   typedef Inverse<CodType> SrcXprType;
   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
   {
-    dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows()));
+    typedef Matrix<typename CodType::Scalar, CodType::RowsAtCompileTime, CodType::RowsAtCompileTime, 0, CodType::MaxRowsAtCompileTime, CodType::MaxRowsAtCompileTime> IdentityMatrixType;
+    dst = src.nestedExpression().solve(IdentityMatrixType::Identity(src.cols(), src.cols()));
   }
 };
 
diff --git a/test/permutationmatrices.cpp b/test/permutationmatrices.cpp
index a26b4ad45..d4b68b2d4 100644
--- a/test/permutationmatrices.cpp
+++ b/test/permutationmatrices.cpp
@@ -54,7 +54,11 @@ template<typename MatrixType> void permutationmatrices(const MatrixType& m)
   m_permuted = m_original;
   VERIFY_EVALUATION_COUNT(m_permuted = lp * m_permuted * rp, 1);
   VERIFY_IS_APPROX(m_permuted, lm*m_original*rm);
-  
+
+  LeftPermutationType lpi;
+  lpi = lp.inverse();
+  VERIFY_IS_APPROX(lpi*m_permuted,lp.inverse()*m_permuted);
+
   VERIFY_IS_APPROX(lp.inverse()*m_permuted*rp.inverse(), m_original);
   VERIFY_IS_APPROX(lv.asPermutation().inverse()*m_permuted*rv.asPermutation().inverse(), m_original);
   VERIFY_IS_APPROX(MapLeftPerm(lv.data(),lv.size()).inverse()*m_permuted*MapRightPerm(rv.data(),rv.size()).inverse(), m_original);
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index a563b5470..06f16438f 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -70,10 +70,11 @@ void cod_fixedsize() {
     Cols = MatrixType::ColsAtCompileTime
   };
   typedef typename MatrixType::Scalar Scalar;
+  typedef CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols> > COD;
   int rank = internal::random<int>(1, (std::min)(int(Rows), int(Cols)) - 1);
   Matrix<Scalar, Rows, Cols> matrix;
   createRandomPIMatrixOfRank(rank, Rows, Cols, matrix);
-  CompleteOrthogonalDecomposition<Matrix<Scalar, Rows, Cols> > cod(matrix);
+  COD cod(matrix);
   VERIFY(rank == cod.rank());
   VERIFY(Cols - cod.rank() == cod.dimensionOfKernel());
   VERIFY(cod.isInjective() == (rank == Rows));
@@ -90,6 +91,9 @@ void cod_fixedsize() {
   JacobiSVD<MatrixType> svd(matrix, ComputeFullU | ComputeFullV);
   Matrix<Scalar, Cols, Cols2> svd_solution = svd.solve(rhs);
   VERIFY_IS_APPROX(cod_solution, svd_solution);
+
+  typename Inverse<COD>::PlainObject pinv = cod.pseudoInverse();
+  VERIFY_IS_APPROX(cod_solution, pinv * rhs);
 }
 
 template<typename MatrixType> void qr()