From ecd2c8f37b8023b56d00cec4aebec7d2f3157e3f Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 20 Feb 2014 14:18:24 +0100 Subject: [PATCH] Add general Inverse<> expression with evaluator --- Eigen/src/Core/MatrixBase.h | 5 + Eigen/src/Core/Solve.h | 6 +- Eigen/src/Core/util/ForwardDeclarations.h | 1 + Eigen/src/LU/Determinant.h | 4 + Eigen/src/LU/Inverse.h | 152 +++++++++++++++++++++- 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 class MatrixBase } #endif + #ifdef EIGEN_TEST_EVALUATORS + EIGEN_DEVICE_FUNC + const Inverse inverse() const; + #else EIGEN_DEVICE_FUNC const internal::inverse_impl inverse() const; + #endif template 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 class SolveImpl : public MatrixBase > 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 class CwiseBinaryOp; template class SelfCwiseBinaryOp; // TODO deprecated template class ProductBase; // TODO deprecated template class Solve; +template class Inverse; namespace internal { template 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 inline typename internal::traits::Scalar MatrixBase::determinant() const { eigen_assert(rows() == cols()); +#ifdef EIGEN_TEST_EVALUATORS + typedef typename internal::nested_eval::type Nested; +#else typedef typename internal::nested::type Nested; +#endif return internal::determinant_impl::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 static inline void run(const MatrixType& matrix, ResultType& result) { typedef typename MatrixType::Scalar Scalar; +#ifdef EIGEN_TEST_EVALUATORS + typename internal::evaluator::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 @@ -279,6 +284,7 @@ struct compute_inverse_and_det_with_check *** MatrixBase methods *** *************************/ +#ifndef EIGEN_TEST_EVALUATORS template struct traits > { @@ -313,9 +319,141 @@ struct inverse_impl : public ReturnByValue > compute_inverse::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 class InverseImpl; + +namespace internal { + +template +struct traits > + : traits +{ + typedef typename XprType::PlainObject PlainObject; + typedef traits 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 +class Inverse : public InverseImpl::StorageKind> +{ +public: + typedef typename XprType::Index Index; + typedef typename XprType::PlainObject PlainObject; + typedef typename internal::nested::type XprTypeNested; + typedef typename internal::remove_all::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 +class InverseImpl + : public MatrixBase > +{ + typedef Inverse Derived; + +public: + + typedef MatrixBase 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 +struct evaluator > + : public evaluator::PlainObject>::type +{ + typedef Inverse InverseType; + typedef typename InverseType::PlainObject PlainObject; + typedef typename evaluator::type Base; + + typedef evaluator type; + typedef evaluator nestedType; + + evaluator(const InverseType& inv_xpr) + : m_result(inv_xpr.rows(), inv_xpr.cols()) + { + ::new (static_cast(this)) Base(m_result); + + typedef typename internal::nested_eval::type ActualXprType; + typedef typename internal::remove_all::type ActualXprTypeCleanded; + + ActualXprType actual_xpr(inv_xpr.nestedExpression()); + + compute_inverse::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 +struct Assignment, internal::assign_op, Dense2Dense, Scalar> +{ + typedef Inverse SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + { + // 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::type ActualXprType; + typedef typename internal::remove_all::type ActualXprTypeCleanded; + + ActualXprType actual_xpr(src.nestedExpression()); + + compute_inverse::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 > * * \sa computeInverseAndDetWithCheck() */ +#ifdef EIGEN_TEST_EVALUATORS +template +inline const Inverse MatrixBase::inverse() const +{ + EIGEN_STATIC_ASSERT(!NumTraits::IsInteger,THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES) + eigen_assert(rows() == cols()); + return Inverse(derived()); +} +#else template inline const internal::inverse_impl MatrixBase::inverse() const { @@ -340,6 +487,7 @@ inline const internal::inverse_impl MatrixBase::inverse() cons eigen_assert(rows() == cols()); return internal::inverse_impl(derived()); } +#endif /** \lu_module *