Add general Inverse<> expression with evaluator

5 files changed, 163 insertions, 5 deletions

diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index 5a3675a20..72aa75da8 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -412,8 +412,13 @@ template<typename Derived> class MatrixBase
     }
     #endif
 
+    #ifdef EIGEN_TEST_EVALUATORS
+    EIGEN_DEVICE_FUNC
+    const Inverse<Derived> inverse() const;
+    #else
     EIGEN_DEVICE_FUNC
     const internal::inverse_impl<Derived> inverse() const;
+    #endif
     template<typename ResultType>
     void computeInverseAndDetWithCheck(
       ResultType& inverse,
diff --git a/Eigen/src/Core/Solve.h b/Eigen/src/Core/Solve.h
index 2da2aff56..cab0e916e 100644
--- a/Eigen/src/Core/Solve.h
+++ b/Eigen/src/Core/Solve.h
@@ -7,8 +7,8 @@
 // Public License v. 2.0. If a copy of the MPL was not distributed
 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-#ifndef EIGEN_INVERSE_H
-#define EIGEN_INVERSE_H
+#ifndef EIGEN_SOLVE_H
+#define EIGEN_SOLVE_H
 
 namespace Eigen {
 
@@ -81,7 +81,7 @@ protected:
 };
 
 
-// Specilaization of the Solve expression for dense results
+// Specialization of the Solve expression for dense results
 template<typename Decomposition, typename RhsType>
 class SolveImpl<Decomposition,RhsType,Dense>
   : public MatrixBase<Solve<Decomposition,RhsType> >
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index 8cf13abd9..ecb811910 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -96,6 +96,7 @@ template<typename BinaryOp,  typename Lhs, typename Rhs>  class CwiseBinaryOp;
 template<typename BinOp,     typename Lhs, typename Rhs>  class SelfCwiseBinaryOp;      // TODO deprecated
 template<typename Derived,   typename Lhs, typename Rhs>  class ProductBase;            // TODO deprecated
 template<typename Decomposition, typename Rhstype>        class Solve;
+template<typename XprType>                                class Inverse;
 
 namespace internal {
   template<typename Lhs, typename Rhs> struct product_tag;
diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h
index bb8e78a8a..9726bd96a 100644
--- a/Eigen/src/LU/Determinant.h
+++ b/Eigen/src/LU/Determinant.h
@@ -92,7 +92,11 @@ template<typename Derived>
 inline typename internal::traits<Derived>::Scalar MatrixBase<Derived>::determinant() const
 {
   eigen_assert(rows() == cols());
+#ifdef EIGEN_TEST_EVALUATORS
+  typedef typename internal::nested_eval<Derived,Base::RowsAtCompileTime>::type Nested;
+#else
   typedef typename internal::nested<Derived,Base::RowsAtCompileTime>::type Nested;
+#endif
   return internal::determinant_impl<typename internal::remove_all<Nested>::type>::run(derived());
 }
 
diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h
index 8d1364e0a..593ce6f8f 100644
--- a/Eigen/src/LU/Inverse.h
+++ b/Eigen/src/LU/Inverse.h
@@ -42,7 +42,12 @@ struct compute_inverse<MatrixType, ResultType, 1>
   static inline void run(const MatrixType& matrix, ResultType& result)
   {
     typedef typename MatrixType::Scalar Scalar;
+#ifdef EIGEN_TEST_EVALUATORS
+    typename internal::evaluator<MatrixType>::type matrixEval(matrix);
+    result.coeffRef(0,0) = Scalar(1) / matrixEval.coeff(0,0);
+#else
     result.coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0);
+#endif
   }
 };
 
@@ -75,10 +80,10 @@ inline void compute_inverse_size2_helper(
     const MatrixType& matrix, const typename ResultType::Scalar& invdet,
     ResultType& result)
 {
-  result.coeffRef(0,0) = matrix.coeff(1,1) * invdet;
+  result.coeffRef(0,0) =  matrix.coeff(1,1) * invdet;
   result.coeffRef(1,0) = -matrix.coeff(1,0) * invdet;
   result.coeffRef(0,1) = -matrix.coeff(0,1) * invdet;
-  result.coeffRef(1,1) = matrix.coeff(0,0) * invdet;
+  result.coeffRef(1,1) =  matrix.coeff(0,0) * invdet;
 }
 
 template<typename MatrixType, typename ResultType>
@@ -279,6 +284,7 @@ struct compute_inverse_and_det_with_check<MatrixType, ResultType, 4>
 *** MatrixBase methods ***
 *************************/
 
+#ifndef EIGEN_TEST_EVALUATORS
 template<typename MatrixType>
 struct traits<inverse_impl<MatrixType> >
 {
@@ -313,9 +319,141 @@ struct inverse_impl : public ReturnByValue<inverse_impl<MatrixType> >
     compute_inverse<MatrixTypeNestedCleaned, Dest>::run(m_matrix, dst);
   }
 };
+#endif
+} // end namespace internal
+
+#ifdef EIGEN_TEST_EVALUATORS
+
+// TODO move the general declaration in Core, and rename this file DenseInverseImpl.h, or something like this...
+
+template<typename XprType,typename StorageKind> class InverseImpl;
+
+namespace internal {
+
+template<typename XprType>
+struct traits<Inverse<XprType> >
+  : traits<typename XprType::PlainObject>
+{
+  typedef typename XprType::PlainObject PlainObject;
+  typedef traits<PlainObject> BaseTraits;
+  enum {
+    Flags = BaseTraits::Flags & RowMajorBit,
+    CoeffReadCost = Dynamic
+  };
+};
+
+} // end namespace internal
+
+/** \class Inverse
+  * \ingroup LU_Module
+  *
+  * \brief Expression of the inverse of another expression
+  *
+  * \tparam XprType the type of the expression we are taking the inverse
+  *
+  * This class represents an expression of A.inverse()
+  * and most of the time this is the only way it is used.
+  *
+  */
+template<typename XprType>
+class Inverse : public InverseImpl<XprType,typename internal::traits<XprType>::StorageKind>
+{
+public:
+  typedef typename XprType::Index Index;
+  typedef typename XprType::PlainObject                       PlainObject;
+  typedef typename internal::nested<XprType>::type            XprTypeNested;
+  typedef typename internal::remove_all<XprTypeNested>::type  XprTypeNestedCleaned;
+  
+  Inverse(const XprType &xpr)
+    : 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 const XprTypeNestedCleaned& nestedExpression() const { return m_xpr; }
+
+protected:
+  XprTypeNested &m_xpr;
+};
 
+// Specialization of the Inverse expression for dense expressions
+template<typename XprType>
+class InverseImpl<XprType,Dense>
+  : public MatrixBase<Inverse<XprType> >
+{
+  typedef Inverse<XprType> Derived;
+  
+public:
+  
+  typedef MatrixBase<Derived> Base;
+  EIGEN_DENSE_PUBLIC_INTERFACE(Derived)
+
+private:
+  
+  Scalar coeff(Index row, Index col) const;
+  Scalar coeff(Index i) const;
+};
+
+namespace internal {
+
+// Evaluator of Inverse -> eval into a temporary
+template<typename XprType>
+struct evaluator<Inverse<XprType> >
+  : public evaluator<typename Inverse<XprType>::PlainObject>::type
+{
+  typedef Inverse<XprType> InverseType;
+  typedef typename InverseType::PlainObject PlainObject;
+  typedef typename evaluator<PlainObject>::type Base;
+  
+  typedef evaluator type;
+  typedef evaluator nestedType;
+
+  evaluator(const InverseType& inv_xpr)
+    : m_result(inv_xpr.rows(), inv_xpr.cols())
+  {
+    ::new (static_cast<Base*>(this)) Base(m_result);
+    
+    typedef typename internal::nested_eval<XprType,XprType::ColsAtCompileTime>::type  ActualXprType;
+    typedef typename internal::remove_all<ActualXprType>::type                        ActualXprTypeCleanded;
+    
+    ActualXprType actual_xpr(inv_xpr.nestedExpression());
+    
+    compute_inverse<ActualXprTypeCleanded, PlainObject>::run(actual_xpr, m_result);
+  }
+  
+protected:
+  PlainObject m_result;
+};
+
+// Specialization for "dst = xpr.inverse()"
+// NOTE we need to specialize it for Dense2Dense to avoid ambiguous specialization error and a Sparse2Sparse specialization must exist somewhere
+template<typename DstXprType, typename XprType, typename Scalar>
+struct Assignment<DstXprType, Inverse<XprType>, internal::assign_op<Scalar>, Dense2Dense, Scalar>
+{
+  typedef Inverse<XprType> SrcXprType;
+  static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar> &)
+  {
+    // FIXME shall we resize dst here?
+    const int Size = EIGEN_PLAIN_ENUM_MIN(XprType::ColsAtCompileTime,DstXprType::ColsAtCompileTime);
+    EIGEN_ONLY_USED_FOR_DEBUG(Size);
+    eigen_assert(( (Size<=1) || (Size>4) || (extract_data(src.nestedExpression())!=extract_data(dst)))
+              && "Aliasing problem detected in inverse(), you need to do inverse().eval() here.");
+
+    typedef typename internal::nested_eval<XprType,XprType::ColsAtCompileTime>::type  ActualXprType;
+    typedef typename internal::remove_all<ActualXprType>::type                        ActualXprTypeCleanded;
+    
+    ActualXprType actual_xpr(src.nestedExpression());
+    
+    compute_inverse<ActualXprTypeCleanded, DstXprType>::run(actual_xpr, dst);
+  }
+};
+
+  
 } // end namespace internal
 
+#endif
+
 /** \lu_module
   *
   * \returns the matrix inverse of this matrix.
@@ -333,6 +471,15 @@ struct inverse_impl : public ReturnByValue<inverse_impl<MatrixType> >
   *
   * \sa computeInverseAndDetWithCheck()
   */
+#ifdef EIGEN_TEST_EVALUATORS
+template<typename Derived>
+inline const Inverse<Derived> MatrixBase<Derived>::inverse() const
+{
+  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsInteger,THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES)
+  eigen_assert(rows() == cols());
+  return Inverse<Derived>(derived());
+}
+#else
 template<typename Derived>
 inline const internal::inverse_impl<Derived> MatrixBase<Derived>::inverse() const
 {
@@ -340,6 +487,7 @@ inline const internal::inverse_impl<Derived> MatrixBase<Derived>::inverse() cons
   eigen_assert(rows() == cols());
   return internal::inverse_impl<Derived>(derived());
 }
+#endif
 
 /** \lu_module
   *