diff --git a/Eigen/src/Core/NoAlias.h b/Eigen/src/Core/NoAlias.h index 6ac525336..3c9c951f0 100644 --- a/Eigen/src/Core/NoAlias.h +++ b/Eigen/src/Core/NoAlias.h @@ -40,6 +40,8 @@ class NoAlias EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ExpressionType& operator=(const StorageBase& other) { + // TODO either call resize here or call "call_assignment" through m_expression.lazyAssign() ?? + m_expression.resizeLike(other.derived()); call_assignment(*this, other.derived(), internal::assign_op()); return m_expression; } diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h index 5dd480cad..c3a9f0db4 100644 --- a/Eigen/src/Core/ProductEvaluators.h +++ b/Eigen/src/Core/ProductEvaluators.h @@ -17,7 +17,7 @@ namespace Eigen { namespace internal { -// Like more general binary expressions, products need they own evaluator: +// Like more general binary expressions, products need their own evaluator: template< typename T, int ProductTag = internal::product_tag::ret, typename LhsShape = typename evaluator_traits::Shape, @@ -39,11 +39,14 @@ struct evaluator > evaluator(const XprType& xpr) : Base(xpr) {} }; -// Helper class to perform a dense product with the destination at hand. +// Helper class to perform a matrix product with the destination at hand. // Depending on the sizes of the factors, there are different evaluation strategies // as controlled by internal::product_type. -template::value> -struct dense_product_impl; +template< typename Lhs, typename Rhs, + typename LhsShape = typename evaluator_traits::Shape, + typename RhsShape = typename evaluator_traits::Shape, + int ProductType = internal::product_type::value> +struct generic_product_impl; template struct evaluator_traits > @@ -52,7 +55,7 @@ struct evaluator_traits > enum { AssumeAliasing = 1 }; }; -// The evaluator for default dense products creates a temporary and call dense_product_impl +// The evaluator for default dense products creates a temporary and call generic_product_impl template struct product_evaluator, ProductTag, DenseShape, DenseShape, typename Lhs::Scalar, typename Rhs::Scalar> : public evaluator::PlainObject>::type @@ -65,7 +68,7 @@ struct product_evaluator, ProductTag, DenseSha : m_result(xpr.rows(), xpr.cols()) { ::new (static_cast(this)) Base(m_result); - dense_product_impl::evalTo(m_result, xpr.lhs(), xpr.rhs()); + generic_product_impl::evalTo(m_result, xpr.lhs(), xpr.rhs()); } protected: @@ -79,7 +82,7 @@ struct Assignment, internal::assign_ typedef Product SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { - dense_product_impl::evalTo(dst, src.lhs(), src.rhs()); + generic_product_impl::evalTo(dst, src.lhs(), src.rhs()); } }; @@ -90,7 +93,7 @@ struct Assignment, internal::add_ass typedef Product SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &) { - dense_product_impl::addTo(dst, src.lhs(), src.rhs()); + generic_product_impl::addTo(dst, src.lhs(), src.rhs()); } }; @@ -101,12 +104,12 @@ struct Assignment, internal::sub_ass typedef Product SrcXprType; static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &) { - dense_product_impl::subTo(dst, src.lhs(), src.rhs()); + generic_product_impl::subTo(dst, src.lhs(), src.rhs()); } }; template -struct dense_product_impl +struct generic_product_impl { template static inline void evalTo(Dst& dst, const Lhs& lhs, const Rhs& rhs) @@ -128,7 +131,7 @@ struct dense_product_impl template -struct dense_product_impl +struct generic_product_impl { typedef typename Product::Scalar Scalar; @@ -165,7 +168,7 @@ struct dense_product_impl // This base class provides default implementations for evalTo, addTo, subTo, in terms of scaleAndAddTo template -struct dense_product_impl_base +struct generic_product_impl_base { typedef typename Product::Scalar Scalar; @@ -188,7 +191,8 @@ struct dense_product_impl_base }; template -struct dense_product_impl : dense_product_impl_base > +struct generic_product_impl + : generic_product_impl_base > { typedef typename Product::Scalar Scalar; enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight }; @@ -205,20 +209,21 @@ struct dense_product_impl : dense_product_impl_base -struct dense_product_impl : dense_product_impl_base > +struct generic_product_impl + : generic_product_impl_base > { typedef typename Product::Scalar Scalar; -// template -// static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) -// { -// // TODO bypass GeneralProduct class -// GeneralProduct(lhs,rhs).scaleAndAddTo(dst, alpha); -// } + template + static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) + { + // TODO bypass GeneralProduct class + GeneralProduct(lhs,rhs).scaleAndAddTo(dst, alpha); + } }; template -struct dense_product_impl +struct generic_product_impl { typedef typename Product::Scalar Scalar; @@ -249,8 +254,10 @@ struct dense_product_impl // { dst += alpha * lazyprod(lhs,rhs); } }; +// This specialization enforces the use of a coefficient-based evaluation strategy template -struct dense_product_impl : dense_product_impl {}; +struct generic_product_impl + : generic_product_impl {}; // Case 2: Evaluate coeff by coeff // @@ -547,6 +554,81 @@ struct etor_product_packet_impl } }; + +/*************************************************************************** +* Triangular products +***************************************************************************/ + +template +struct generic_product_impl + : generic_product_impl_base > +{ + typedef typename Product::Scalar Scalar; + + template + static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) + { + // TODO bypass TriangularProduct class + TriangularProduct(lhs.nestedExpression(),rhs).scaleAndAddTo(dst, alpha); + } +}; + +template +struct product_evaluator, ProductTag, TriangularShape, DenseShape, typename Lhs::Scalar, typename Rhs::Scalar> + : public evaluator::PlainObject>::type +{ + typedef Product XprType; + typedef typename XprType::PlainObject PlainObject; + typedef typename evaluator::type Base; + + product_evaluator(const XprType& xpr) + : m_result(xpr.rows(), xpr.cols()) + { + ::new (static_cast(this)) Base(m_result); + generic_product_impl::evalTo(m_result, xpr.lhs(), xpr.rhs()); + } + +protected: + PlainObject m_result; +}; + + +template +struct generic_product_impl +: generic_product_impl_base > +{ + typedef typename Product::Scalar Scalar; + + template + static void scaleAndAddTo(Dest& dst, const Lhs& lhs, const Rhs& rhs, const Scalar& alpha) + { + // TODO bypass TriangularProduct class + TriangularProduct(lhs,rhs.nestedExpression()).scaleAndAddTo(dst, alpha); + } +}; + +template +struct product_evaluator, ProductTag, DenseShape, TriangularShape, typename Lhs::Scalar, typename Rhs::Scalar> + : public evaluator::PlainObject>::type +{ + typedef Product XprType; + typedef typename XprType::PlainObject PlainObject; + typedef typename evaluator::type Base; + + product_evaluator(const XprType& xpr) + : m_result(xpr.rows(), xpr.cols()) + { + ::new (static_cast(this)) Base(m_result); + generic_product_impl::evalTo(m_result, xpr.lhs(), xpr.rhs()); + } + +protected: + PlainObject m_result; +}; + + + + } // end namespace internal } // end namespace Eigen diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 602f0c5d1..28191694d 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -49,6 +49,7 @@ template class TriangularBase : public EigenBase typedef typename internal::traits::Index Index; typedef typename internal::traits::DenseMatrixType DenseMatrixType; typedef DenseMatrixType DenseType; + typedef Derived const& Nested; EIGEN_DEVICE_FUNC inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); } @@ -204,6 +205,7 @@ template class TriangularView enum { Mode = _Mode, + Flags = internal::traits::Flags, TransposeMode = (Mode & Upper ? Lower : 0) | (Mode & Lower ? Upper : 0) | (Mode & (UnitDiag)) @@ -326,6 +328,27 @@ template class TriangularView return m_matrix.transpose(); } +#ifdef EIGEN_TEST_EVALUATORS + + /** Efficient triangular matrix times vector/matrix product */ + template + EIGEN_DEVICE_FUNC + const Product + operator*(const MatrixBase& rhs) const + { + return Product(*this, rhs.derived()); + } + + /** Efficient vector/matrix times triangular matrix product */ + template friend + EIGEN_DEVICE_FUNC + const Product + operator*(const MatrixBase& lhs, const TriangularView& rhs) + { + return Product(lhs.derived(),rhs); + } + +#else // EIGEN_TEST_EVALUATORS /** Efficient triangular matrix times vector/matrix product */ template EIGEN_DEVICE_FUNC @@ -347,6 +370,7 @@ template class TriangularView (lhs.derived(),rhs.m_matrix); } +#endif #ifdef EIGEN2_SUPPORT template diff --git a/Eigen/src/Core/products/CoeffBasedProduct.h b/Eigen/src/Core/products/CoeffBasedProduct.h index c987a30b0..5f8182aa9 100644 --- a/Eigen/src/Core/products/CoeffBasedProduct.h +++ b/Eigen/src/Core/products/CoeffBasedProduct.h @@ -49,8 +49,8 @@ struct traits > enum { LhsCoeffReadCost = _LhsNested::CoeffReadCost, RhsCoeffReadCost = _RhsNested::CoeffReadCost, - LhsFlags = _LhsNested::Flags, - RhsFlags = _RhsNested::Flags, + LhsFlags = traits<_LhsNested>::Flags, + RhsFlags = traits<_RhsNested>::Flags, RowsAtCompileTime = _LhsNested::RowsAtCompileTime, ColsAtCompileTime = _RhsNested::ColsAtCompileTime, diff --git a/test/evaluators.cpp b/test/evaluators.cpp index c07146260..d4b737348 100644 --- a/test/evaluators.cpp +++ b/test/evaluators.cpp @@ -451,5 +451,10 @@ void test_evaluators() C.triangularView().swap(D.triangularView()); swap_using_evaluator(B.triangularView(), A.triangularView()); VERIFY(B.isApprox(C) && "swap_using_evaluator(B.triangularView(), A.triangularView())"); + + + VERIFY_IS_APPROX_EVALUATOR2(B, prod(A.triangularView(),A), MatrixXd(A.triangularView()*A)); + + B.col(0).noalias() = prod( (2.1 * A.adjoint()).triangularView() , (A.row(0)).adjoint() ); } }