From ccc41128fb936563651a3c5b25bfc80f303710bf Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Wed, 19 Feb 2014 11:33:29 +0100 Subject: [PATCH] Add a Solve expression for uniform treatment of solve() methods. --- Eigen/Core | 3 + Eigen/src/Core/CoreEvaluators.h | 4 +- Eigen/src/Core/Solve.h | 145 ++++++++++++++++++++++ Eigen/src/Core/TriangularMatrix.h | 17 +++ Eigen/src/Core/util/ForwardDeclarations.h | 3 +- 5 files changed, 169 insertions(+), 3 deletions(-) create mode 100644 Eigen/src/Core/Solve.h diff --git a/Eigen/Core b/Eigen/Core index c1d168ec4..e00cd7a16 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -361,6 +361,9 @@ using std::ptrdiff_t; #include "src/Core/Flagged.h" #include "src/Core/ProductBase.h" #include "src/Core/GeneralProduct.h" +#ifdef EIGEN_ENABLE_EVALUATORS +#include "src/Core/Solve.h" +#endif #include "src/Core/TriangularMatrix.h" #include "src/Core/SelfAdjointView.h" #include "src/Core/products/GeneralBlockPanelKernel.h" diff --git a/Eigen/src/Core/CoreEvaluators.h b/Eigen/src/Core/CoreEvaluators.h index d02bac0a8..726a6854a 100644 --- a/Eigen/src/Core/CoreEvaluators.h +++ b/Eigen/src/Core/CoreEvaluators.h @@ -495,14 +495,14 @@ struct evaluator > PacketScalar packet(Index row, Index col) const { return m_functor.packetOp(m_lhsImpl.template packet(row, col), - m_rhsImpl.template packet(row, col)); + m_rhsImpl.template packet(row, col)); } template PacketScalar packet(Index index) const { return m_functor.packetOp(m_lhsImpl.template packet(index), - m_rhsImpl.template packet(index)); + m_rhsImpl.template packet(index)); } protected: diff --git a/Eigen/src/Core/Solve.h b/Eigen/src/Core/Solve.h new file mode 100644 index 000000000..2da2aff56 --- /dev/null +++ b/Eigen/src/Core/Solve.h @@ -0,0 +1,145 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2014 Gael Guennebaud +// +// 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 +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_INVERSE_H +#define EIGEN_INVERSE_H + +namespace Eigen { + +template class SolveImpl; + +/** \class Solve + * \ingroup Core_Module + * + * \brief Pseudo expression representing a solving operation + * + * \tparam Decomposition the type of the matrix or decomposion object + * \tparam Rhstype the type of the right-hand side + * + * This class represents an expression of A.solve(B) + * and most of the time this is the only way it is used. + * + */ +namespace internal { + +// this solve_traits class permits to determine the evaluation type with respect to storage kind (Dense vs Sparse) +template struct solve_traits; + +template +struct solve_traits +{ + typedef typename Decomposition::MatrixType MatrixType; + typedef Matrix PlainObject; +}; + +template +struct traits > + : traits::StorageKind>::PlainObject> +{ + typedef typename solve_traits::StorageKind>::PlainObject PlainObject; + typedef traits BaseTraits; + enum { + Flags = BaseTraits::Flags & RowMajorBit, + CoeffReadCost = Dynamic + }; +}; + +} + + +template +class Solve : public SolveImpl::StorageKind> +{ +public: + typedef typename RhsType::Index Index; + typedef typename internal::traits::PlainObject PlainObject; + + Solve(const Decomposition &dec, const RhsType &rhs) + : m_dec(dec), m_rhs(rhs) + {} + + EIGEN_DEVICE_FUNC Index rows() const { return m_dec.rows(); } + EIGEN_DEVICE_FUNC Index cols() const { return m_rhs.cols(); } + + EIGEN_DEVICE_FUNC const Decomposition& dec() const { return m_dec; } + EIGEN_DEVICE_FUNC const RhsType& rhs() const { return m_rhs; } + +protected: + const Decomposition &m_dec; + const RhsType &m_rhs; +}; + + +// Specilaization of the Solve expression for dense results +template +class SolveImpl + : public MatrixBase > +{ + typedef Solve 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 Solve -> eval into a temporary +template +struct evaluator > + : public evaluator::PlainObject>::type +{ + typedef Solve SolveType; + typedef typename SolveType::PlainObject PlainObject; + typedef typename evaluator::type Base; + + typedef evaluator type; + typedef evaluator nestedType; + + evaluator(const SolveType& solve) + : m_result(solve.rows(), solve.cols()) + { + ::new (static_cast(this)) Base(m_result); + solve.dec()._solve_impl(solve.rhs(), m_result); + } + +protected: + PlainObject m_result; +}; + +// Specialization for "dst = dec.solve(rhs)" +// 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 Solve SrcXprType; + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + { + // FIXME shall we resize dst here? + src.dec()._solve_impl(src.rhs(), dst); + } +}; + +} // end namepsace internal + +} // end namespace Eigen + +#endif // EIGEN_SOLVE_H diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index ffa2faeaa..76ba76d5f 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -437,11 +437,19 @@ template class TriangularView EIGEN_DEVICE_FUNC void solveInPlace(const MatrixBase& other) const; +#ifdef EIGEN_TEST_EVALUATORS + template + EIGEN_DEVICE_FUNC + inline const Solve + solve(const MatrixBase& other) const + { return Solve(*this, other.derived()); } +#else // EIGEN_TEST_EVALUATORS template EIGEN_DEVICE_FUNC inline const internal::triangular_solve_retval solve(const MatrixBase& other) const { return solve(other); } +#endif // EIGEN_TEST_EVALUATORS template EIGEN_DEVICE_FUNC @@ -547,6 +555,15 @@ template class TriangularView #endif // EIGEN_TEST_EVALUATORS #ifdef EIGEN_TEST_EVALUATORS + + template + EIGEN_DEVICE_FUNC + EIGEN_STRONG_INLINE void _solve_impl(const RhsType &rhs, DstType &dst) const { + if(!(internal::is_same::value && internal::extract_data(dst) == internal::extract_data(rhs))) + dst = rhs; + this->template solveInPlace(dst); + } + template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TriangularView& _assignProduct(const ProductType& prod, const Scalar& alpha); diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index db0e2b033..8cf13abd9 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -94,7 +94,8 @@ template class CwiseUnaryOp; template class CwiseUnaryView; template class CwiseBinaryOp; template class SelfCwiseBinaryOp; // TODO deprecated -template class ProductBase; +template class ProductBase; // TODO deprecated +template class Solve; namespace internal { template struct product_tag;