From 6cc9dc17f2b016883260c00ff05c0a819e24bd24 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Tue, 26 Jan 2010 18:45:23 -0500
Subject: [PATCH] In LU / Inverse, decouple the output type from the input
 type. This has long been done in the default branch

---
 Eigen/src/Core/MatrixBase.h |  3 +-
 Eigen/src/LU/Inverse.h      | 55 +++++++++++++++++++------------------
 Eigen/src/LU/LU.h           |  3 +-
 3 files changed, 32 insertions(+), 29 deletions(-)

diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index fd1b74735..de1b242a3 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -583,7 +583,8 @@ template<typename Derived> class MatrixBase
 
     const LU<PlainMatrixType> lu() const;
     const PlainMatrixType inverse() const;
-    void computeInverse(PlainMatrixType *result) const;
+    template<typename ResultType>
+    void computeInverse(MatrixBase<ResultType> *result) const;
     Scalar determinant() const;
 
 /////////// Cholesky module ///////////
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index 98e7fc5d0..9e1c8cfa6 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -54,10 +54,10 @@ bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixTy
   return true;
 }
 
-template<typename MatrixType>
-void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result)
+template<typename Derived, typename OtherDerived>
+void ei_compute_inverse_in_size3_case(const Derived& matrix, OtherDerived* result)
 {
-  typedef typename MatrixType::Scalar Scalar;
+  typedef typename Derived::Scalar Scalar;
   const Scalar det_minor00 = matrix.minor(0,0).determinant();
   const Scalar det_minor10 = matrix.minor(1,0).determinant();
   const Scalar det_minor20 = matrix.minor(2,0).determinant();
@@ -75,10 +75,10 @@ void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* resu
   result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet;
 }
 
-template<typename MatrixType, typename Scalar = typename MatrixType::Scalar>
+template<typename Derived, typename OtherDerived, typename Scalar = typename Derived::Scalar>
 struct ei_compute_inverse_in_size4_case
 {
-  static void run(const MatrixType& matrix, MatrixType& result)
+  static void run(const Derived& matrix, OtherDerived& result)
   {
     result.coeffRef(0,0) = matrix.minor(0,0).determinant();
     result.coeffRef(1,0) = -matrix.minor(0,1).determinant();
@@ -109,10 +109,10 @@ struct ei_compute_inverse_in_size4_case
 // here that if Intel makes this document publically available, with source code
 // and detailed explanations, it's because they want their CPUs to be fed with
 // good code, and therefore they presumably don't mind us using it in Eigen.
-template<typename MatrixType>
-struct ei_compute_inverse_in_size4_case<MatrixType, float>
+template<typename Derived, typename OtherDerived>
+struct ei_compute_inverse_in_size4_case<Derived, OtherDerived, float>
 {
-  static void run(const MatrixType& matrix, MatrixType& result)
+  static void run(const Derived& matrix, OtherDerived& result)
   {
     // Variables (Streaming SIMD Extensions registers) which will contain cofactors and, later, the
     // lines of the inverted matrix.
@@ -229,50 +229,50 @@ struct ei_compute_inverse_in_size4_case<MatrixType, float>
 *** Part 2 : selector and MatrixBase methods ***
 ***********************************************/
 
-template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime>
+template<typename Derived, typename OtherDerived, int Size = Derived::RowsAtCompileTime>
 struct ei_compute_inverse
 {
-  static inline void run(const MatrixType& matrix, MatrixType* result)
+  static inline void run(const Derived& matrix, OtherDerived* result)
   {
-    LU<MatrixType> lu(matrix);
+    LU<Derived> lu(matrix);
     lu.computeInverse(result);
   }
 };
 
-template<typename MatrixType>
-struct ei_compute_inverse<MatrixType, 1>
+template<typename Derived, typename OtherDerived>
+struct ei_compute_inverse<Derived, OtherDerived, 1>
 {
-  static inline void run(const MatrixType& matrix, MatrixType* result)
+  static inline void run(const Derived& matrix, OtherDerived* result)
   {
-    typedef typename MatrixType::Scalar Scalar;
+    typedef typename Derived::Scalar Scalar;
     result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
   }
 };
 
-template<typename MatrixType>
-struct ei_compute_inverse<MatrixType, 2>
+template<typename Derived, typename OtherDerived>
+struct ei_compute_inverse<Derived, OtherDerived, 2>
 {
-  static inline void run(const MatrixType& matrix, MatrixType* result)
+  static inline void run(const Derived& matrix, OtherDerived* result)
   {
     ei_compute_inverse_in_size2_case(matrix, result);
   }
 };
 
-template<typename MatrixType>
-struct ei_compute_inverse<MatrixType, 3>
+template<typename Derived, typename OtherDerived>
+struct ei_compute_inverse<Derived, OtherDerived, 3>
 {
-  static inline void run(const MatrixType& matrix, MatrixType* result)
+  static inline void run(const Derived& matrix, OtherDerived* result)
   {
     ei_compute_inverse_in_size3_case(matrix, result);
   }
 };
 
-template<typename MatrixType>
-struct ei_compute_inverse<MatrixType, 4>
+template<typename Derived, typename OtherDerived>
+struct ei_compute_inverse<Derived, OtherDerived, 4>
 {
-  static inline void run(const MatrixType& matrix, MatrixType* result)
+  static inline void run(const Derived& matrix, OtherDerived* result)
   {
-    ei_compute_inverse_in_size4_case<MatrixType>::run(matrix, *result);
+    ei_compute_inverse_in_size4_case<Derived, OtherDerived>::run(matrix, *result);
   }
 };
 
@@ -290,11 +290,12 @@ struct ei_compute_inverse<MatrixType, 4>
   * \sa inverse()
   */
 template<typename Derived>
-inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const
+template<typename OtherDerived>
+inline void MatrixBase<Derived>::computeInverse(MatrixBase<OtherDerived> *result) const
 {
   ei_assert(rows() == cols());
   EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT)
-  ei_compute_inverse<PlainMatrixType>::run(eval(), result);
+  ei_compute_inverse<PlainMatrixType, OtherDerived>::run(eval(), static_cast<OtherDerived*>(result));
 }
 
 /** \lu_module
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h
index 176e76a91..36ee886d4 100644
--- a/Eigen/src/LU/LU.h
+++ b/Eigen/src/LU/LU.h
@@ -297,7 +297,8 @@ template<typename MatrixType> class LU
       *
       * \sa MatrixBase::computeInverse(), inverse()
       */
-    inline void computeInverse(MatrixType *result) const
+    template<typename ResultType>
+    inline void computeInverse(ResultType *result) const
     {
       solve(MatrixType::Identity(m_lu.rows(), m_lu.cols()), result);
     }