From 00c29c2caef8fb0c6b1d2ba5ecdf6780c0c766d4 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Thu, 26 May 2016 17:34:23 +0200 Subject: Store permutation's determinant as char. This also fixes some long to float conversion warnings --- Eigen/src/LU/FullPivLU.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'Eigen/src/LU/FullPivLU.h') diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 64b9eb7f1..c39f8e3d5 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -418,9 +418,10 @@ template class FullPivLU PermutationQType m_q; IntColVectorType m_rowsTranspositions; IntRowVectorType m_colsTranspositions; - Index m_det_pq, m_nonzero_pivots; + Index m_nonzero_pivots; RealScalar m_l1_norm; RealScalar m_maxpivot, m_prescribedThreshold; + char m_det_pq; bool m_isInitialized, m_usePrescribedThreshold; }; -- cgit v1.2.3 From 66e99ab6a1444d8e3d47211e4540837e6b982a3a Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 6 Jun 2016 15:11:41 +0200 Subject: Relax mixing-type constraints for binary coefficient-wise operators: - Replace internal::scalar_product_traits by Eigen::ScalarBinaryOpTraits - Remove the "functor_is_product_like" helper (was pretty ugly) - Currently, OP is not used, but it is available to the user for fine grained tuning - Currently, only the following operators have been generalized: *,/,+,-,=,*=,/=,+=,-= - TODO: generalize all other binray operators (comparisons,pow,etc.) - TODO: handle "scalar op array" operators (currently only * is handled) - TODO: move the handling of the "void" scalar type to ScalarBinaryOpTraits --- Eigen/src/Core/ArrayBase.h | 4 +- Eigen/src/Core/AssignEvaluator.h | 12 ++-- Eigen/src/Core/CwiseBinaryOp.h | 4 +- Eigen/src/Core/DiagonalMatrix.h | 6 +- Eigen/src/Core/Dot.h | 12 ++-- Eigen/src/Core/EigenBase.h | 4 +- Eigen/src/Core/MatrixBase.h | 4 +- Eigen/src/Core/NoAlias.h | 6 +- Eigen/src/Core/PlainObjectBase.h | 2 +- Eigen/src/Core/Product.h | 3 +- Eigen/src/Core/ProductEvaluators.h | 55 ++++++++++--------- Eigen/src/Core/Redux.h | 4 +- Eigen/src/Core/Ref.h | 2 +- Eigen/src/Core/SelfCwiseBinaryOp.h | 10 ++-- Eigen/src/Core/Solve.h | 13 +++-- Eigen/src/Core/TriangularMatrix.h | 20 +++---- Eigen/src/Core/VectorwiseOp.h | 4 +- Eigen/src/Core/functors/AssignmentFunctors.h | 54 +++++++++--------- Eigen/src/Core/functors/BinaryFunctors.h | 64 +++++++++------------- Eigen/src/Core/products/GeneralBlockPanelKernel.h | 4 +- Eigen/src/Core/products/GeneralMatrixMatrix.h | 4 +- .../Core/products/GeneralMatrixMatrixTriangular.h | 4 +- Eigen/src/Core/products/GeneralMatrixVector.h | 4 +- Eigen/src/Core/products/TriangularMatrixVector.h | 4 +- Eigen/src/Core/util/ForwardDeclarations.h | 5 +- Eigen/src/Core/util/Macros.h | 4 +- Eigen/src/Core/util/Meta.h | 64 +++++++++++++--------- Eigen/src/Core/util/XprHelper.h | 10 +--- Eigen/src/Geometry/AlignedBox.h | 6 +- Eigen/src/Geometry/Homogeneous.h | 12 ++-- Eigen/src/Householder/HouseholderSequence.h | 2 +- Eigen/src/IterativeLinearSolvers/SolveWithGuess.h | 4 +- Eigen/src/LU/FullPivLU.h | 6 +- Eigen/src/LU/InverseImpl.h | 6 +- Eigen/src/LU/PartialPivLU.h | 6 +- Eigen/src/QR/ColPivHouseholderQR.h | 4 +- Eigen/src/QR/CompleteOrthogonalDecomposition.h | 4 +- Eigen/src/QR/FullPivHouseholderQR.h | 4 +- Eigen/src/SparseCore/SparseAssign.h | 20 +++---- Eigen/src/SparseCore/SparseCwiseBinaryOp.h | 20 +++---- Eigen/src/SparseCore/SparseDenseProduct.h | 4 +- Eigen/src/SparseCore/SparseMatrix.h | 4 +- Eigen/src/SparseCore/SparseMatrixBase.h | 2 +- Eigen/src/SparseCore/SparseProduct.h | 12 ++-- Eigen/src/SparseCore/SparseSelfAdjointView.h | 10 ++-- Eigen/src/SparseQR/SparseQR.h | 8 +-- Eigen/src/plugins/ArrayCwiseBinaryOps.h | 4 +- Eigen/src/plugins/CommonCwiseBinaryOps.h | 18 +++++- Eigen/src/plugins/CommonCwiseUnaryOps.h | 6 +- Eigen/src/plugins/MatrixCwiseBinaryOps.h | 4 +- blas/PackedTriangularMatrixVector.h | 4 +- test/array.cpp | 2 +- test/array_for_matrix.cpp | 2 +- test/mixingtypes.cpp | 32 ++++++++--- test/vectorization_logic.cpp | 7 ++- unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h | 12 ++-- .../src/KroneckerProduct/KroneckerTensorProduct.h | 4 +- 57 files changed, 323 insertions(+), 292 deletions(-) (limited to 'Eigen/src/LU/FullPivLU.h') diff --git a/Eigen/src/Core/ArrayBase.h b/Eigen/src/Core/ArrayBase.h index 3d9c37bf6..62851a0c2 100644 --- a/Eigen/src/Core/ArrayBase.h +++ b/Eigen/src/Core/ArrayBase.h @@ -176,7 +176,7 @@ template EIGEN_STRONG_INLINE Derived & ArrayBase::operator-=(const ArrayBase &other) { - call_assignment(derived(), other.derived(), internal::sub_assign_op()); + call_assignment(derived(), other.derived(), internal::sub_assign_op()); return derived(); } @@ -189,7 +189,7 @@ template EIGEN_STRONG_INLINE Derived & ArrayBase::operator+=(const ArrayBase& other) { - call_assignment(derived(), other.derived(), internal::add_assign_op()); + call_assignment(derived(), other.derived(), internal::add_assign_op()); return derived(); } diff --git a/Eigen/src/Core/AssignEvaluator.h b/Eigen/src/Core/AssignEvaluator.h index 4b914ac0c..f966724cc 100644 --- a/Eigen/src/Core/AssignEvaluator.h +++ b/Eigen/src/Core/AssignEvaluator.h @@ -687,7 +687,7 @@ EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(const DstX template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_dense_assignment_loop(const DstXprType& dst, const SrcXprType& src) { - call_dense_assignment_loop(dst, src, internal::assign_op()); + call_dense_assignment_loop(dst, src, internal::assign_op()); } /*************************************************************************** @@ -722,13 +722,13 @@ template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_assignment(Dst& dst, const Src& src) { - call_assignment(dst, src, internal::assign_op()); + call_assignment(dst, src, internal::assign_op()); } template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_assignment(const Dst& dst, const Src& src) { - call_assignment(dst, src, internal::assign_op()); + call_assignment(dst, src, internal::assign_op()); } // Deal with "assume-aliasing" @@ -787,7 +787,7 @@ template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_assignment_no_alias(Dst& dst, const Src& src) { - call_assignment_no_alias(dst, src, internal::assign_op()); + call_assignment_no_alias(dst, src, internal::assign_op()); } template @@ -809,7 +809,7 @@ template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_assignment_no_alias_no_transpose(Dst& dst, const Src& src) { - call_assignment_no_alias_no_transpose(dst, src, internal::assign_op()); + call_assignment_no_alias_no_transpose(dst, src, internal::assign_op()); } // forward declaration @@ -838,7 +838,7 @@ template< typename DstXprType, typename SrcXprType, typename Functor, typename S struct Assignment { EIGEN_DEVICE_FUNC - static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static EIGEN_STRONG_INLINE void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) { eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); src.evalTo(dst); diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index 39820fd7d..aa3297354 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -160,7 +160,7 @@ template EIGEN_STRONG_INLINE Derived & MatrixBase::operator-=(const MatrixBase &other) { - call_assignment(derived(), other.derived(), internal::sub_assign_op()); + call_assignment(derived(), other.derived(), internal::sub_assign_op()); return derived(); } @@ -173,7 +173,7 @@ template EIGEN_STRONG_INLINE Derived & MatrixBase::operator+=(const MatrixBase& other) { - call_assignment(derived(), other.derived(), internal::add_assign_op()); + call_assignment(derived(), other.derived(), internal::add_assign_op()); return derived(); } diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h index 5a9e3abd4..aa619dd5c 100644 --- a/Eigen/src/Core/DiagonalMatrix.h +++ b/Eigen/src/Core/DiagonalMatrix.h @@ -320,16 +320,16 @@ template<> struct AssignmentKind { typedef Diagonal2De template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> struct Assignment { - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) { dst.setZero(); dst.diagonal() = src.diagonal(); } - static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &/*func*/) + static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &/*func*/) { dst.diagonal() += src.diagonal(); } - static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &/*func*/) + static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &/*func*/) { dst.diagonal() -= src.diagonal(); } }; diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index f3c869635..1d7f2262e 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -28,22 +28,24 @@ template struct dot_nocheck { - typedef typename scalar_product_traits::Scalar,typename traits::Scalar>::ReturnType ResScalar; + typedef scalar_conj_product_op::Scalar,typename traits::Scalar> conj_prod; + typedef typename conj_prod::result_type ResScalar; EIGEN_DEVICE_FUNC static inline ResScalar run(const MatrixBase& a, const MatrixBase& b) { - return a.template binaryExpr::Scalar,typename traits::Scalar> >(b).sum(); + return a.template binaryExpr(b).sum(); } }; template struct dot_nocheck { - typedef typename scalar_product_traits::Scalar,typename traits::Scalar>::ReturnType ResScalar; + typedef scalar_conj_product_op::Scalar,typename traits::Scalar> conj_prod; + typedef typename conj_prod::result_type ResScalar; EIGEN_DEVICE_FUNC static inline ResScalar run(const MatrixBase& a, const MatrixBase& b) { - return a.transpose().template binaryExpr::Scalar,typename traits::Scalar> >(b).sum(); + return a.transpose().template binaryExpr(b).sum(); } }; @@ -62,7 +64,7 @@ struct dot_nocheck template template EIGEN_DEVICE_FUNC -typename internal::scalar_product_traits::Scalar,typename internal::traits::Scalar>::ReturnType +typename ScalarBinaryOpTraits::Scalar,typename internal::traits::Scalar>::ReturnType MatrixBase::dot(const MatrixBase& other) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) diff --git a/Eigen/src/Core/EigenBase.h b/Eigen/src/Core/EigenBase.h index ba8e09674..f76995af9 100644 --- a/Eigen/src/Core/EigenBase.h +++ b/Eigen/src/Core/EigenBase.h @@ -138,7 +138,7 @@ template template Derived& DenseBase::operator+=(const EigenBase &other) { - call_assignment(derived(), other.derived(), internal::add_assign_op()); + call_assignment(derived(), other.derived(), internal::add_assign_op()); return derived(); } @@ -146,7 +146,7 @@ template template Derived& DenseBase::operator-=(const EigenBase &other) { - call_assignment(derived(), other.derived(), internal::sub_assign_op()); + call_assignment(derived(), other.derived(), internal::sub_assign_op()); return derived(); } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index b8b7f458f..f63505fef 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -193,7 +193,7 @@ template class MatrixBase template EIGEN_DEVICE_FUNC - typename internal::scalar_product_traits::Scalar,typename internal::traits::Scalar>::ReturnType + typename ScalarBinaryOpTraits::Scalar,typename internal::traits::Scalar>::ReturnType dot(const MatrixBase& other) const; EIGEN_DEVICE_FUNC RealScalar squaredNorm() const; @@ -381,7 +381,7 @@ template class MatrixBase #ifndef EIGEN_PARSED_BY_DOXYGEN /// \internal helper struct to form the return type of the cross product template struct cross_product_return_type { - typedef typename internal::scalar_product_traits::Scalar,typename internal::traits::Scalar>::ReturnType Scalar; + typedef typename ScalarBinaryOpTraits::Scalar,typename internal::traits::Scalar>::ReturnType Scalar; typedef Matrix type; }; #endif // EIGEN_PARSED_BY_DOXYGEN diff --git a/Eigen/src/Core/NoAlias.h b/Eigen/src/Core/NoAlias.h index ffb673cee..33908010b 100644 --- a/Eigen/src/Core/NoAlias.h +++ b/Eigen/src/Core/NoAlias.h @@ -39,7 +39,7 @@ class NoAlias EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ExpressionType& operator=(const StorageBase& other) { - call_assignment_no_alias(m_expression, other.derived(), internal::assign_op()); + call_assignment_no_alias(m_expression, other.derived(), internal::assign_op()); return m_expression; } @@ -47,7 +47,7 @@ class NoAlias EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ExpressionType& operator+=(const StorageBase& other) { - call_assignment_no_alias(m_expression, other.derived(), internal::add_assign_op()); + call_assignment_no_alias(m_expression, other.derived(), internal::add_assign_op()); return m_expression; } @@ -55,7 +55,7 @@ class NoAlias EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ExpressionType& operator-=(const StorageBase& other) { - call_assignment_no_alias(m_expression, other.derived(), internal::sub_assign_op()); + call_assignment_no_alias(m_expression, other.derived(), internal::sub_assign_op()); return m_expression; } diff --git a/Eigen/src/Core/PlainObjectBase.h b/Eigen/src/Core/PlainObjectBase.h index 570dbd53b..64f5eb052 100644 --- a/Eigen/src/Core/PlainObjectBase.h +++ b/Eigen/src/Core/PlainObjectBase.h @@ -718,7 +718,7 @@ class PlainObjectBase : public internal::dense_xpr_base::type //_resize_to_match(other); // the 'false' below means to enforce lazy evaluation. We don't use lazyAssign() because // it wouldn't allow to copy a row-vector into a column-vector. - internal::call_assignment_no_alias(this->derived(), other.derived(), internal::assign_op()); + internal::call_assignment_no_alias(this->derived(), other.derived(), internal::assign_op()); return this->derived(); } diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 8aa1de081..bad289761 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -18,11 +18,12 @@ namespace internal { // Determine the scalar of Product. This is normally the same as Lhs::Scalar times // Rhs::Scalar, but product with permutation matrices inherit the scalar of the other factor. +// TODO: this could be removed once ScalarBinaryOpTraits handles void. template::Shape, typename RhsShape = typename evaluator_traits::Shape > struct product_result_scalar { - typedef typename scalar_product_traits::ReturnType Scalar; + typedef typename ScalarBinaryOpTraits::ReturnType Scalar; }; template diff --git a/Eigen/src/Core/ProductEvaluators.h b/Eigen/src/Core/ProductEvaluators.h index cc7166062..7f041e5dd 100644 --- a/Eigen/src/Core/ProductEvaluators.h +++ b/Eigen/src/Core/ProductEvaluators.h @@ -124,12 +124,12 @@ protected: // Dense = Product template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar> -struct Assignment, internal::assign_op, Dense2Dense, +struct Assignment, internal::assign_op::Scalar>, Dense2Dense, typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct),Scalar>::type> { typedef Product SrcXprType; static EIGEN_STRONG_INLINE - void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { // FIXME shall we handle nested_eval here? generic_product_impl::evalTo(dst, src.lhs(), src.rhs()); @@ -138,12 +138,12 @@ struct Assignment, internal::assign_op -struct Assignment, internal::add_assign_op, Dense2Dense, +struct Assignment, internal::add_assign_op::Scalar>, Dense2Dense, typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct),Scalar>::type> { typedef Product SrcXprType; static EIGEN_STRONG_INLINE - void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &) + void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &) { // FIXME shall we handle nested_eval here? generic_product_impl::addTo(dst, src.lhs(), src.rhs()); @@ -152,12 +152,12 @@ struct Assignment, internal::add_assign_op< // Dense -= Product template< typename DstXprType, typename Lhs, typename Rhs, int Options, typename Scalar> -struct Assignment, internal::sub_assign_op, Dense2Dense, +struct Assignment, internal::sub_assign_op::Scalar>, Dense2Dense, typename enable_if<(Options==DefaultProduct || Options==AliasFreeProduct),Scalar>::type> { typedef Product SrcXprType; static EIGEN_STRONG_INLINE - void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &) + void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &) { // FIXME shall we handle nested_eval here? generic_product_impl::subTo(dst, src.lhs(), src.rhs()); @@ -187,37 +187,38 @@ struct Assignment" as well. template -struct evaluator_assume_aliasing, const OtherXpr, +struct evaluator_assume_aliasing::Scalar>, const OtherXpr, const Product >, DenseShape > { static const bool value = true; }; -template +template struct assignment_from_xpr_plus_product { - typedef CwiseBinaryOp, const OtherXpr, const ProductType> SrcXprType; + typedef CwiseBinaryOp, const OtherXpr, const ProductType> SrcXprType; + template static EIGEN_STRONG_INLINE - void run(DstXprType &dst, const SrcXprType &src, const Func1& func) + void run(DstXprType &dst, const SrcXprType &src, const InitialFunc& /*func*/) { - call_assignment_no_alias(dst, src.lhs(), func); + call_assignment_no_alias(dst, src.lhs(), Func1()); call_assignment_no_alias(dst, src.rhs(), Func2()); } }; -template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename Scalar> -struct Assignment, const OtherXpr, - const Product >, internal::assign_op, Dense2Dense> - : assignment_from_xpr_plus_product, Scalar, internal::assign_op, internal::add_assign_op > +template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar> +struct Assignment, const OtherXpr, + const Product >, internal::assign_op, Dense2Dense> + : assignment_from_xpr_plus_product, internal::assign_op, internal::add_assign_op > {}; -template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename Scalar> -struct Assignment, const OtherXpr, - const Product >, internal::add_assign_op, Dense2Dense> - : assignment_from_xpr_plus_product, Scalar, internal::add_assign_op, internal::add_assign_op > +template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar> +struct Assignment, const OtherXpr, + const Product >, internal::add_assign_op, Dense2Dense> + : assignment_from_xpr_plus_product, internal::add_assign_op, internal::add_assign_op > {}; -template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename Scalar> -struct Assignment, const OtherXpr, - const Product >, internal::sub_assign_op, Dense2Dense> - : assignment_from_xpr_plus_product, Scalar, internal::sub_assign_op, internal::sub_assign_op > +template< typename DstXprType, typename OtherXpr, typename Lhs, typename Rhs, typename DstScalar, typename SrcScalar, typename OtherScalar,typename ProdScalar> +struct Assignment, const OtherXpr, + const Product >, internal::sub_assign_op, Dense2Dense> + : assignment_from_xpr_plus_product, internal::sub_assign_op, internal::sub_assign_op > {}; //---------------------------------------- @@ -369,21 +370,21 @@ struct generic_product_impl { // Same as: dst.noalias() = lhs.lazyProduct(rhs); // but easier on the compiler side - call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::assign_op()); + call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::assign_op()); } template static EIGEN_STRONG_INLINE void addTo(Dst& dst, const Lhs& lhs, const Rhs& rhs) { // dst.noalias() += lhs.lazyProduct(rhs); - call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::add_assign_op()); + call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::add_assign_op()); } template static EIGEN_STRONG_INLINE void subTo(Dst& dst, const Lhs& lhs, const Rhs& rhs) { // dst.noalias() -= lhs.lazyProduct(rhs); - call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::sub_assign_op()); + call_assignment_no_alias(dst, lhs.lazyProduct(rhs), internal::sub_assign_op()); } // template @@ -735,7 +736,7 @@ template { - typedef typename scalar_product_traits::ReturnType Scalar; + typedef typename ScalarBinaryOpTraits::ReturnType Scalar; public: enum { CoeffReadCost = NumTraits::MulCost + evaluator::CoeffReadCost + evaluator::CoeffReadCost, diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h index 7984cd6e1..ec969d9b9 100644 --- a/Eigen/src/Core/Redux.h +++ b/Eigen/src/Core/Redux.h @@ -450,7 +450,7 @@ DenseBase::sum() const { if(SizeAtCompileTime==0 || (SizeAtCompileTime==Dynamic && size()==0)) return Scalar(0); - return derived().redux(Eigen::internal::scalar_sum_op()); + return derived().redux(Eigen::internal::scalar_sum_op()); } /** \returns the mean of all coefficients of *this @@ -465,7 +465,7 @@ DenseBase::mean() const #pragma warning push #pragma warning ( disable : 2259 ) #endif - return Scalar(derived().redux(Eigen::internal::scalar_sum_op())) / Scalar(this->size()); + return Scalar(derived().redux(Eigen::internal::scalar_sum_op())) / Scalar(this->size()); #ifdef __INTEL_COMPILER #pragma warning pop #endif diff --git a/Eigen/src/Core/Ref.h b/Eigen/src/Core/Ref.h index 6e94181f3..17065fdd5 100644 --- a/Eigen/src/Core/Ref.h +++ b/Eigen/src/Core/Ref.h @@ -262,7 +262,7 @@ template class Ref< template EIGEN_DEVICE_FUNC void construct(const Expression& expr, internal::false_type) { - internal::call_assignment_no_alias(m_object,expr,internal::assign_op()); + internal::call_assignment_no_alias(m_object,expr,internal::assign_op()); Base::construct(m_object); } diff --git a/Eigen/src/Core/SelfCwiseBinaryOp.h b/Eigen/src/Core/SelfCwiseBinaryOp.h index 78fff1549..719ed72a5 100644 --- a/Eigen/src/Core/SelfCwiseBinaryOp.h +++ b/Eigen/src/Core/SelfCwiseBinaryOp.h @@ -12,11 +12,13 @@ namespace Eigen { +// TODO generalize the scalar type of 'other' + template EIGEN_STRONG_INLINE Derived& DenseBase::operator*=(const Scalar& other) { typedef typename Derived::PlainObject PlainObject; - internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::mul_assign_op()); + internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::mul_assign_op()); return derived(); } @@ -24,7 +26,7 @@ template EIGEN_STRONG_INLINE Derived& ArrayBase::operator+=(const Scalar& other) { typedef typename Derived::PlainObject PlainObject; - internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::add_assign_op()); + internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::add_assign_op()); return derived(); } @@ -32,7 +34,7 @@ template EIGEN_STRONG_INLINE Derived& ArrayBase::operator-=(const Scalar& other) { typedef typename Derived::PlainObject PlainObject; - internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::sub_assign_op()); + internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::sub_assign_op()); return derived(); } @@ -40,7 +42,7 @@ template EIGEN_STRONG_INLINE Derived& DenseBase::operator/=(const Scalar& other) { typedef typename Derived::PlainObject PlainObject; - internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::div_assign_op()); + internal::call_assignment(this->derived(), PlainObject::Constant(rows(),cols(),other), internal::div_assign_op()); return derived(); } diff --git a/Eigen/src/Core/Solve.h b/Eigen/src/Core/Solve.h index ba2ee53b8..038ad5b11 100644 --- a/Eigen/src/Core/Solve.h +++ b/Eigen/src/Core/Solve.h @@ -134,10 +134,10 @@ protected: // 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> +struct Assignment, internal::assign_op, Dense2Dense, Scalar> { typedef Solve SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + 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); @@ -146,10 +146,10 @@ struct Assignment, internal::assign_op -struct Assignment,RhsType>, internal::assign_op, Dense2Dense, Scalar> +struct Assignment,RhsType>, internal::assign_op, Dense2Dense, Scalar> { typedef Solve,RhsType> SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { src.dec().nestedExpression().template _solve_impl_transposed(src.rhs(), dst); } @@ -157,10 +157,11 @@ struct Assignment,RhsType>, internal: // Specialization for "dst = dec.adjoint().solve(rhs)" template -struct Assignment, const Transpose >,RhsType>, internal::assign_op, Dense2Dense, Scalar> +struct Assignment, const Transpose >,RhsType>, + internal::assign_op, Dense2Dense, Scalar> { typedef Solve, const Transpose >,RhsType> SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { src.dec().nestedExpression().nestedExpression().template _solve_impl_transposed(src.rhs(), dst); } diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 5c5e5028e..8731e9127 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -367,14 +367,14 @@ template class TriangularViewImpl<_Mat template EIGEN_DEVICE_FUNC TriangularViewType& operator+=(const DenseBase& other) { - internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op()); + internal::call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op()); return derived(); } /** \sa MatrixBase::operator-=() */ template EIGEN_DEVICE_FUNC TriangularViewType& operator-=(const DenseBase& other) { - internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op()); + internal::call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op()); return derived(); } @@ -552,7 +552,7 @@ template inline TriangularView& TriangularViewImpl::operator=(const MatrixBase& other) { - internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op()); + internal::call_assignment_no_alias(derived(), other.derived(), internal::assign_op()); return derived(); } @@ -804,7 +804,7 @@ template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src) { - call_triangular_assignment_loop(dst, src, internal::assign_op()); + call_triangular_assignment_loop(dst, src, internal::assign_op()); } template<> struct AssignmentKind { typedef Triangular2Triangular Kind; }; @@ -933,10 +933,10 @@ namespace internal { // Triangular = Product template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> -struct Assignment, internal::assign_op, Dense2Triangular, Scalar> +struct Assignment, internal::assign_op::Scalar>, Dense2Triangular, Scalar> { typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { dst.setZero(); dst._assignProduct(src, 1); @@ -945,10 +945,10 @@ struct Assignment, internal::assign_ // Triangular += Product template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> -struct Assignment, internal::add_assign_op, Dense2Triangular, Scalar> +struct Assignment, internal::add_assign_op::Scalar>, Dense2Triangular, Scalar> { typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &) { dst._assignProduct(src, 1); } @@ -956,10 +956,10 @@ struct Assignment, internal::add_ass // Triangular -= Product template< typename DstXprType, typename Lhs, typename Rhs, typename Scalar> -struct Assignment, internal::sub_assign_op, Dense2Triangular, Scalar> +struct Assignment, internal::sub_assign_op::Scalar>, Dense2Triangular, Scalar> { typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &) { dst._assignProduct(src, -1); } diff --git a/Eigen/src/Core/VectorwiseOp.h b/Eigen/src/Core/VectorwiseOp.h index 193891189..00a4a8c39 100644 --- a/Eigen/src/Core/VectorwiseOp.h +++ b/Eigen/src/Core/VectorwiseOp.h @@ -540,7 +540,7 @@ template class VectorwiseOp /** Returns the expression of the sum of the vector \a other to each subvector of \c *this */ template EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC - CwiseBinaryOp, + CwiseBinaryOp, const ExpressionTypeNestedCleaned, const typename ExtendedType::Type> operator+(const DenseBase& other) const @@ -553,7 +553,7 @@ template class VectorwiseOp /** Returns the expression of the difference between each subvector of \c *this and the vector \a other */ template EIGEN_DEVICE_FUNC - CwiseBinaryOp, + CwiseBinaryOp, const ExpressionTypeNestedCleaned, const typename ExtendedType::Type> operator-(const DenseBase& other) const diff --git a/Eigen/src/Core/functors/AssignmentFunctors.h b/Eigen/src/Core/functors/AssignmentFunctors.h index 51fef50e8..9b373c783 100644 --- a/Eigen/src/Core/functors/AssignmentFunctors.h +++ b/Eigen/src/Core/functors/AssignmentFunctors.h @@ -18,20 +18,24 @@ namespace internal { * \brief Template functor for scalar/packet assignment * */ -template struct assign_op { +template struct assign_op { EIGEN_EMPTY_STRUCT_CTOR(assign_op) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void assignCoeff(Scalar& a, const Scalar& b) const { a = b; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void assignCoeff(DstScalar& a, const SrcScalar& b) const { a = b; } template - EIGEN_STRONG_INLINE void assignPacket(Scalar* a, const Packet& b) const - { internal::pstoret(a,b); } + EIGEN_STRONG_INLINE void assignPacket(DstScalar* a, const Packet& b) const + { internal::pstoret(a,b); } }; -template -struct functor_traits > { + +// Empty overload for void type (used by PermutationMatrix +template struct assign_op {}; + +template +struct functor_traits > { enum { - Cost = NumTraits::ReadCost, - PacketAccess = packet_traits::Vectorizable + Cost = NumTraits::ReadCost, + PacketAccess = is_same::value && packet_traits::Vectorizable && packet_traits::Vectorizable }; }; @@ -39,20 +43,20 @@ struct functor_traits > { * \brief Template functor for scalar/packet assignment with addition * */ -template struct add_assign_op { +template struct add_assign_op { EIGEN_EMPTY_STRUCT_CTOR(add_assign_op) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void assignCoeff(Scalar& a, const Scalar& b) const { a += b; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void assignCoeff(DstScalar& a, const SrcScalar& b) const { a += b; } template - EIGEN_STRONG_INLINE void assignPacket(Scalar* a, const Packet& b) const - { internal::pstoret(a,internal::padd(internal::ploadt(a),b)); } + EIGEN_STRONG_INLINE void assignPacket(DstScalar* a, const Packet& b) const + { internal::pstoret(a,internal::padd(internal::ploadt(a),b)); } }; -template -struct functor_traits > { +template +struct functor_traits > { enum { - Cost = NumTraits::ReadCost + NumTraits::AddCost, - PacketAccess = packet_traits::HasAdd + Cost = NumTraits::ReadCost + NumTraits::AddCost, + PacketAccess = is_same::value && packet_traits::HasAdd }; }; @@ -60,20 +64,20 @@ struct functor_traits > { * \brief Template functor for scalar/packet assignment with subtraction * */ -template struct sub_assign_op { +template struct sub_assign_op { EIGEN_EMPTY_STRUCT_CTOR(sub_assign_op) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void assignCoeff(Scalar& a, const Scalar& b) const { a -= b; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void assignCoeff(DstScalar& a, const SrcScalar& b) const { a -= b; } template - EIGEN_STRONG_INLINE void assignPacket(Scalar* a, const Packet& b) const - { internal::pstoret(a,internal::psub(internal::ploadt(a),b)); } + EIGEN_STRONG_INLINE void assignPacket(DstScalar* a, const Packet& b) const + { internal::pstoret(a,internal::psub(internal::ploadt(a),b)); } }; -template -struct functor_traits > { +template +struct functor_traits > { enum { - Cost = NumTraits::ReadCost + NumTraits::AddCost, - PacketAccess = packet_traits::HasSub + Cost = NumTraits::ReadCost + NumTraits::AddCost, + PacketAccess = is_same::value && packet_traits::HasSub }; }; @@ -98,7 +102,6 @@ struct functor_traits > { PacketAccess = is_same::value && packet_traits::HasMul }; }; -template struct functor_is_product_like > { enum { ret = 1 }; }; /** \internal * \brief Template functor for scalar/packet assignment with diviving @@ -120,7 +123,6 @@ struct functor_traits > { PacketAccess = is_same::value && packet_traits::HasDiv }; }; -template struct functor_is_product_like > { enum { ret = 1 }; }; /** \internal * \brief Template functor for scalar/packet assignment with swapping diff --git a/Eigen/src/Core/functors/BinaryFunctors.h b/Eigen/src/Core/functors/BinaryFunctors.h index 6eb5b91ce..98fcebae5 100644 --- a/Eigen/src/Core/functors/BinaryFunctors.h +++ b/Eigen/src/Core/functors/BinaryFunctors.h @@ -21,22 +21,23 @@ namespace internal { * * \sa class CwiseBinaryOp, MatrixBase::operator+, class VectorwiseOp, DenseBase::sum() */ -template struct scalar_sum_op { -// typedef Scalar result_type; +template struct scalar_sum_op { + typedef typename ScalarBinaryOpTraits::ReturnType result_type; EIGEN_EMPTY_STRUCT_CTOR(scalar_sum_op) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return a + b; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const result_type operator() (const LhsScalar& a, const RhsScalar& b) const { return a + b; } template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const { return internal::padd(a,b); } template - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar predux(const Packet& a) const + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const result_type predux(const Packet& a) const { return internal::predux(a); } }; -template -struct functor_traits > { +template +struct functor_traits > { enum { - Cost = NumTraits::AddCost, - PacketAccess = packet_traits::HasAdd + Cost = (NumTraits::AddCost+NumTraits::AddCost)/2, // rough estimate! + PacketAccess = is_same::value && packet_traits::HasAdd && packet_traits::HasAdd + // TODO vectorize mixed sum }; }; @@ -45,7 +46,7 @@ struct functor_traits > { * This is required to solve Bug 426. * \sa DenseBase::count(), DenseBase::any(), ArrayBase::cast(), MatrixBase::cast() */ -template<> struct scalar_sum_op : scalar_sum_op { +template<> struct scalar_sum_op : scalar_sum_op { EIGEN_DEPRECATED scalar_sum_op() {} }; @@ -57,11 +58,7 @@ template<> struct scalar_sum_op : scalar_sum_op { * \sa class CwiseBinaryOp, Cwise::operator*(), class VectorwiseOp, MatrixBase::redux() */ template struct scalar_product_op { - enum { - // TODO vectorize mixed product - Vectorizable = is_same::value && packet_traits::HasMul && packet_traits::HasMul - }; - typedef typename scalar_product_traits::ReturnType result_type; + typedef typename ScalarBinaryOpTraits::ReturnType result_type; EIGEN_EMPTY_STRUCT_CTOR(scalar_product_op) EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const result_type operator() (const LhsScalar& a, const RhsScalar& b) const { return a * b; } template @@ -75,7 +72,8 @@ template struct functor_traits > { enum { Cost = (NumTraits::MulCost + NumTraits::MulCost)/2, // rough estimate! - PacketAccess = scalar_product_op::Vectorizable + PacketAccess = is_same::value && packet_traits::HasMul && packet_traits::HasMul + // TODO vectorize mixed product }; }; @@ -90,7 +88,7 @@ template struct scalar_conj_product_op { Conj = NumTraits::IsComplex }; - typedef typename scalar_product_traits::ReturnType result_type; + typedef typename ScalarBinaryOpTraits::ReturnType result_type; EIGEN_EMPTY_STRUCT_CTOR(scalar_conj_product_op) EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const result_type operator() (const LhsScalar& a, const RhsScalar& b) const @@ -269,18 +267,19 @@ struct functor_traits > { * * \sa class CwiseBinaryOp, MatrixBase::operator- */ -template struct scalar_difference_op { +template struct scalar_difference_op { + typedef typename ScalarBinaryOpTraits::ReturnType result_type; EIGEN_EMPTY_STRUCT_CTOR(scalar_difference_op) - EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return a - b; } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const result_type operator() (const LhsScalar& a, const RhsScalar& b) const { return a - b; } template EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a, const Packet& b) const { return internal::psub(a,b); } }; -template -struct functor_traits > { +template +struct functor_traits > { enum { - Cost = NumTraits::AddCost, - PacketAccess = packet_traits::HasSub + Cost = (NumTraits::AddCost+NumTraits::AddCost)/2, + PacketAccess = is_same::value && packet_traits::HasSub && packet_traits::HasSub }; }; @@ -290,11 +289,7 @@ struct functor_traits > { * \sa class CwiseBinaryOp, Cwise::operator/() */ template struct scalar_quotient_op { - enum { - // TODO vectorize mixed product - Vectorizable = is_same::value && packet_traits::HasDiv && packet_traits::HasDiv - }; - typedef typename scalar_product_traits::ReturnType result_type; + typedef typename ScalarBinaryOpTraits::ReturnType result_type; EIGEN_EMPTY_STRUCT_CTOR(scalar_quotient_op) EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const result_type operator() (const LhsScalar& a, const RhsScalar& b) const { return a / b; } template @@ -305,7 +300,7 @@ template struct functor_traits > { typedef typename scalar_quotient_op::result_type result_type; enum { - PacketAccess = scalar_quotient_op::Vectorizable, + PacketAccess = is_same::value && packet_traits::HasDiv && packet_traits::HasDiv, Cost = NumTraits::template Div::Cost }; }; @@ -446,7 +441,7 @@ struct functor_traits > template struct scalar_multiple2_op { - typedef typename scalar_product_traits::ReturnType result_type; + typedef typename ScalarBinaryOpTraits::ReturnType result_type; EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE scalar_multiple2_op(const scalar_multiple2_op& other) : m_other(other.m_other) { } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE scalar_multiple2_op(const Scalar2& other) : m_other(other) { } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE result_type operator() (const Scalar1& a) const { return a * m_other; } @@ -481,7 +476,7 @@ struct functor_traits > template struct scalar_quotient2_op { - typedef typename scalar_product_traits::ReturnType result_type; + typedef typename ScalarBinaryOpTraits::ReturnType result_type; EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE scalar_quotient2_op(const scalar_quotient2_op& other) : m_other(other.m_other) { } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE scalar_quotient2_op(const Scalar2& other) : m_other(other) { } EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE result_type operator() (const Scalar1& a) const { return a / m_other; } @@ -491,15 +486,6 @@ template struct functor_traits > { enum { Cost = 2 * NumTraits::MulCost, PacketAccess = false }; }; -// In Eigen, any binary op (Product, CwiseBinaryOp) require the Lhs and Rhs to have the same scalar type, except for multiplication -// where the mixing of different types is handled by scalar_product_traits -// In particular, real * complex is allowed. -// FIXME move this to functor_traits adding a functor_default -template struct functor_is_product_like { enum { ret = 0 }; }; -template struct functor_is_product_like > { enum { ret = 1 }; }; -template struct functor_is_product_like > { enum { ret = 1 }; }; -template struct functor_is_product_like > { enum { ret = 1 }; }; - /** \internal * \brief Template functor to add a scalar to a fixed other one diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 253c03462..63a9fc462 100644 --- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -363,7 +363,7 @@ class gebp_traits public: typedef _LhsScalar LhsScalar; typedef _RhsScalar RhsScalar; - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; enum { ConjLhs = _ConjLhs, @@ -478,7 +478,7 @@ class gebp_traits, RealScalar, _ConjLhs, false> public: typedef std::complex LhsScalar; typedef RealScalar RhsScalar; - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; enum { ConjLhs = _ConjLhs, diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix.h b/Eigen/src/Core/products/GeneralMatrixMatrix.h index 7528fef24..b1465c3b5 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrix.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrix.h @@ -25,7 +25,7 @@ struct general_matrix_matrix_product Traits; - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; static EIGEN_STRONG_INLINE void run( Index rows, Index cols, Index depth, const LhsScalar* lhs, Index lhsStride, @@ -55,7 +55,7 @@ struct general_matrix_matrix_product Traits; -typedef typename scalar_product_traits::ReturnType ResScalar; +typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; static void run(Index rows, Index cols, Index depth, const LhsScalar* _lhs, Index lhsStride, const RhsScalar* _rhs, Index rhsStride, diff --git a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h index 80ba89465..29d6dc721 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrixTriangular.h @@ -40,7 +40,7 @@ template struct general_matrix_matrix_triangular_product { - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; static EIGEN_STRONG_INLINE void run(Index size, Index depth,const LhsScalar* lhs, Index lhsStride, const RhsScalar* rhs, Index rhsStride, ResScalar* res, Index resStride, const ResScalar& alpha, level3_blocking& blocking) @@ -57,7 +57,7 @@ template struct general_matrix_matrix_triangular_product { - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; static EIGEN_STRONG_INLINE void run(Index size, Index depth,const LhsScalar* _lhs, Index lhsStride, const RhsScalar* _rhs, Index rhsStride, ResScalar* _res, Index resStride, const ResScalar& alpha, level3_blocking& blocking) diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h index fc8886511..4a5cf3fb6 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/Eigen/src/Core/products/GeneralMatrixVector.h @@ -58,7 +58,7 @@ namespace internal { template struct general_matrix_vector_product { - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; enum { Vectorizable = packet_traits::Vectorizable && packet_traits::Vectorizable @@ -334,7 +334,7 @@ EIGEN_DONT_INLINE void general_matrix_vector_product struct general_matrix_vector_product { -typedef typename scalar_product_traits::ReturnType ResScalar; +typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; enum { Vectorizable = packet_traits::Vectorizable && packet_traits::Vectorizable diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index f79840aa7..c11a983c7 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h @@ -20,7 +20,7 @@ struct triangular_matrix_vector_product; template struct triangular_matrix_vector_product { - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; enum { IsLower = ((Mode&Lower)==Lower), HasUnitDiag = (Mode & UnitDiag)==UnitDiag, @@ -91,7 +91,7 @@ EIGEN_DONT_INLINE void triangular_matrix_vector_product struct triangular_matrix_vector_product { - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; enum { IsLower = ((Mode&Lower)==Lower), HasUnitDiag = (Mode & UnitDiag)==UnitDiag, diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 42e2e75b9..045e22658 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -131,6 +131,7 @@ template class ArrayWrapper; template class MatrixWrapper; template class SolverBase; template class InnerIterator; +template struct ScalarBinaryOpTraits; namespace internal { template struct kernel_retval_base; @@ -175,8 +176,8 @@ namespace internal { // with optional conjugation of the arguments. template struct conj_helper; -template struct scalar_sum_op; -template struct scalar_difference_op; +template struct scalar_sum_op; +template struct scalar_difference_op; template struct scalar_conj_product_op; template struct scalar_opposite_op; template struct scalar_conjugate_op; diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index c9a0b9893..35547fdda 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -885,9 +885,9 @@ namespace Eigen { } // the expression type of a cwise product -#define EIGEN_CWISE_PRODUCT_RETURN_TYPE(LHS,RHS) \ +#define EIGEN_CWISE_BINARY_RETURN_TYPE(LHS,RHS,OPNAME) \ CwiseBinaryOp< \ - internal::scalar_product_op< \ + EIGEN_CAT(EIGEN_CAT(internal::scalar_,OPNAME),_op)< \ typename internal::traits::Scalar, \ typename internal::traits::Scalar \ >, \ diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 7ecd59add..af661c313 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -375,33 +375,6 @@ template struct scalar_product_traits enum { Defined = 0 }; }; -template struct scalar_product_traits -{ - enum { - // Cost = NumTraits::MulCost, - Defined = 1 - }; - typedef T ReturnType; -}; - -template struct scalar_product_traits > -{ - enum { - // Cost = 2*NumTraits::MulCost, - Defined = 1 - }; - typedef std::complex ReturnType; -}; - -template struct scalar_product_traits, T> -{ - enum { - // Cost = 2*NumTraits::MulCost, - Defined = 1 - }; - typedef std::complex ReturnType; -}; - // FIXME quick workaround around current limitation of result_of // template // struct result_of(ArgType0,ArgType1)> { @@ -434,6 +407,43 @@ T div_ceil(const T &a, const T &b) } // end namespace numext + +/** \class ScalarBinaryOpTraits + * \ingroup Core_Module + * + * \brief Determines whether the given binary operation of two numeric types is allowed and what the scalar return type is. + * + * \sa CwiseBinaryOp + */ +template +struct ScalarBinaryOpTraits +#ifndef EIGEN_PARSED_BY_DOXYGEN + // for backward compatibility, use the hints given by the (deprecated) internal::scalar_product_traits class. + : internal::scalar_product_traits +#endif // EIGEN_PARSED_BY_DOXYGEN +{}; + +template +struct ScalarBinaryOpTraits +{ + enum { Defined = 1 }; + typedef T ReturnType; +}; + +template +struct ScalarBinaryOpTraits,BinaryOp> +{ + enum { Defined = 1 }; + typedef std::complex ReturnType; +}; + +template +struct ScalarBinaryOpTraits, T,BinaryOp> +{ + enum { Defined = 1 }; + typedef std::complex ReturnType; +}; + } // end namespace Eigen #endif // EIGEN_META_H diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 3605de6fd..4fd4a9b0d 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -649,17 +649,13 @@ std::string demangle_flags(int f) } // end namespace internal -// we require Lhs and Rhs to have the same scalar type. Currently there is no example of a binary functor -// that would take two operands of different types. If there were such an example, then this check should be -// moved to the BinaryOp functors, on a per-case basis. This would however require a change in the BinaryOp functors, as -// currently they take only one typename Scalar template parameter. +// We require Lhs and Rhs to have "compatible" scalar types. // It is tempting to always allow mixing different types but remember that this is often impossible in the vectorized paths. // So allowing mixing different types gives very unexpected errors when enabling vectorization, when the user tries to // add together a float matrix and a double matrix. +// Treat "void" as a special case. Needed for permutation products. TODO: this should be handled by ScalarBinaryOpTraits #define EIGEN_CHECK_BINARY_COMPATIBILIY(BINOP,LHS,RHS) \ - EIGEN_STATIC_ASSERT((internal::functor_is_product_like::ret \ - ? int(internal::scalar_product_traits::Defined) \ - : int(internal::is_same_or_void::value)), \ + EIGEN_STATIC_ASSERT(int(internal::is_same_or_void::value) || int(ScalarBinaryOpTraits::Defined), \ YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) } // end namespace Eigen diff --git a/Eigen/src/Geometry/AlignedBox.h b/Eigen/src/Geometry/AlignedBox.h index 03f1a11f8..aeb043a6c 100644 --- a/Eigen/src/Geometry/AlignedBox.h +++ b/Eigen/src/Geometry/AlignedBox.h @@ -112,7 +112,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) /** \returns the center of the box */ inline const CwiseUnaryOp, - const CwiseBinaryOp, const VectorType, const VectorType> > + const CwiseBinaryOp, const VectorType, const VectorType> > center() const { return (m_min+m_max)/2; } @@ -120,7 +120,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) * Note that this function does not get the same * result for integral or floating scalar types: see */ - inline const CwiseBinaryOp< internal::scalar_difference_op, const VectorType, const VectorType> sizes() const + inline const CwiseBinaryOp< internal::scalar_difference_op, const VectorType, const VectorType> sizes() const { return m_max - m_min; } /** \returns the volume of the bounding box */ @@ -131,7 +131,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) * if the length of the diagonal is needed: diagonal().norm() * will provide it. */ - inline CwiseBinaryOp< internal::scalar_difference_op, const VectorType, const VectorType> diagonal() const + inline CwiseBinaryOp< internal::scalar_difference_op, const VectorType, const VectorType> diagonal() const { return sizes(); } /** \returns the vertex of the bounding box at the corner defined by diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h index cd52b5470..1c35ca486 100644 --- a/Eigen/src/Geometry/Homogeneous.h +++ b/Eigen/src/Geometry/Homogeneous.h @@ -329,10 +329,10 @@ protected: // dense = homogeneous template< typename DstXprType, typename ArgType, typename Scalar> -struct Assignment, internal::assign_op, Dense2Dense, Scalar> +struct Assignment, internal::assign_op, Dense2Dense, Scalar> { typedef Homogeneous SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { dst.template topRows(src.nestedExpression().rows()) = src.nestedExpression(); dst.row(dst.rows()-1).setOnes(); @@ -341,10 +341,10 @@ struct Assignment, internal::assign_op // dense = homogeneous template< typename DstXprType, typename ArgType, typename Scalar> -struct Assignment, internal::assign_op, Dense2Dense, Scalar> +struct Assignment, internal::assign_op, Dense2Dense, Scalar> { typedef Homogeneous SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { dst.template leftCols(src.nestedExpression().cols()) = src.nestedExpression(); dst.col(dst.cols()-1).setOnes(); @@ -373,7 +373,7 @@ struct homogeneous_right_product_refactoring_helper typedef typename Rhs::ConstRowXpr ConstantColumn; typedef Replicate ConstantBlock; typedef Product LinearProduct; - typedef CwiseBinaryOp, const LinearProduct, const ConstantBlock> Xpr; + typedef CwiseBinaryOp, const LinearProduct, const ConstantBlock> Xpr; }; template @@ -414,7 +414,7 @@ struct homogeneous_left_product_refactoring_helper typedef typename Lhs::ConstColXpr ConstantColumn; typedef Replicate ConstantBlock; typedef Product LinearProduct; - typedef CwiseBinaryOp, const LinearProduct, const ConstantBlock> Xpr; + typedef CwiseBinaryOp, const LinearProduct, const ConstantBlock> Xpr; }; template diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h index a57f81764..3ce0a693d 100644 --- a/Eigen/src/Householder/HouseholderSequence.h +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -108,7 +108,7 @@ struct hseq_side_dependent_impl template struct matrix_type_times_scalar_type { - typedef typename scalar_product_traits::ReturnType + typedef typename ScalarBinaryOpTraits::ReturnType ResultScalar; typedef Matrix Type; diff --git a/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h b/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h index 35923be3d..7d67d3ce2 100644 --- a/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h +++ b/Eigen/src/IterativeLinearSolvers/SolveWithGuess.h @@ -91,10 +91,10 @@ protected: // Specialization for "dst = dec.solveWithGuess(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> +struct Assignment, internal::assign_op, Dense2Dense, Scalar> { typedef SolveWithGuess SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { // FIXME shall we resize dst here? dst = src.guess(); diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index c39f8e3d5..2d01b18c6 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -839,12 +839,12 @@ namespace internal { /***** Implementation of inverse() *****************************************************/ -template -struct Assignment >, internal::assign_op, Dense2Dense, Scalar> +template +struct Assignment >, internal::assign_op::Scalar>, Dense2Dense> { typedef FullPivLU LuType; typedef Inverse SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); } diff --git a/Eigen/src/LU/InverseImpl.h b/Eigen/src/LU/InverseImpl.h index e202a55cb..3134632e1 100644 --- a/Eigen/src/LU/InverseImpl.h +++ b/Eigen/src/LU/InverseImpl.h @@ -286,11 +286,11 @@ struct compute_inverse_and_det_with_check namespace internal { // Specialization for "dense = dense_xpr.inverse()" -template -struct Assignment, internal::assign_op, Dense2Dense, Scalar> +template +struct Assignment, internal::assign_op, Dense2Dense> { typedef Inverse SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + 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); diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index b68916287..ac2902261 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -525,12 +525,12 @@ MatrixType PartialPivLU::reconstructedMatrix() const namespace internal { /***** Implementation of inverse() *****************************************************/ -template -struct Assignment >, internal::assign_op, Dense2Dense, Scalar> +template +struct Assignment >, internal::assign_op::Scalar>, Dense2Dense> { typedef PartialPivLU LuType; typedef Inverse SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); } diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 7c559f952..525ee8c18 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -598,11 +598,11 @@ void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType & namespace internal { template -struct Assignment >, internal::assign_op, Dense2Dense, Scalar> +struct Assignment >, internal::assign_op, Dense2Dense, Scalar> { typedef ColPivHouseholderQR QrType; typedef Inverse SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); } diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h index 230d0d23c..52bcc2173 100644 --- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h +++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -510,11 +510,11 @@ void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl( namespace internal { template -struct Assignment >, internal::assign_op, Dense2Dense, Scalar> +struct Assignment >, internal::assign_op, Dense2Dense, Scalar> { typedef CompleteOrthogonalDecomposition CodType; typedef Inverse SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.rows())); } diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 32a10f3fe..4f55d52a5 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -560,11 +560,11 @@ void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType namespace internal { template -struct Assignment >, internal::assign_op, Dense2Dense, Scalar> +struct Assignment >, internal::assign_op, Dense2Dense, Scalar> { typedef FullPivHouseholderQR QrType; typedef Inverse SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); } diff --git a/Eigen/src/SparseCore/SparseAssign.h b/Eigen/src/SparseCore/SparseAssign.h index 4a8dd12e4..b284fa9e4 100644 --- a/Eigen/src/SparseCore/SparseAssign.h +++ b/Eigen/src/SparseCore/SparseAssign.h @@ -34,8 +34,8 @@ template inline Derived& SparseMatrixBase::operator=(const SparseMatrixBase& other) { // by default sparse evaluation do not alias, so we can safely bypass the generic call_assignment routine - internal::Assignment > - ::run(derived(), other.derived(), internal::assign_op()); + internal::Assignment > + ::run(derived(), other.derived(), internal::assign_op()); return derived(); } @@ -127,7 +127,7 @@ void assign_sparse_to_sparse(DstXprType &dst, const SrcXprType &src) template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> struct Assignment { - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) { assign_sparse_to_sparse(dst.derived(), src.derived()); } @@ -141,7 +141,7 @@ struct Assignment { eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); - if(internal::is_same >::value) + if(internal::is_same >::value) dst.setZero(); internal::evaluator srcEval(src); @@ -156,10 +156,10 @@ struct Assignment // Specialization for "dst = dec.solve(rhs)" // NOTE we need to specialize it for Sparse2Sparse to avoid ambiguous specialization error template -struct Assignment, internal::assign_op, Sparse2Sparse, Scalar> +struct Assignment, internal::assign_op, Sparse2Sparse, Scalar> { typedef Solve SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { src.dec()._solve_impl(src.rhs(), dst); } @@ -176,7 +176,7 @@ struct Assignment typedef Array ArrayXI; typedef Array ArrayXS; template - static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &/*func*/) { Index size = src.diagonal().size(); dst.makeCompressed(); @@ -187,15 +187,15 @@ struct Assignment } template - static void run(SparseMatrixBase &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static void run(SparseMatrixBase &dst, const SrcXprType &src, const internal::assign_op &/*func*/) { dst.diagonal() = src.diagonal(); } - static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &/*func*/) + static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &/*func*/) { dst.diagonal() += src.diagonal(); } - static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &/*func*/) + static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &/*func*/) { dst.diagonal() -= src.diagonal(); } }; } // end namespace internal diff --git a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h index d422f3cbe..dd21eb8c5 100644 --- a/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ b/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -579,7 +579,7 @@ template template Derived& SparseMatrixBase::operator+=(const DiagonalBase& other) { - call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op()); + call_assignment_no_alias(derived(), other.derived(), internal::add_assign_op()); return derived(); } @@ -587,7 +587,7 @@ template template Derived& SparseMatrixBase::operator-=(const DiagonalBase& other) { - call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op()); + call_assignment_no_alias(derived(), other.derived(), internal::sub_assign_op()); return derived(); } @@ -600,31 +600,31 @@ SparseMatrixBase::cwiseProduct(const MatrixBase &other) c } template -EIGEN_STRONG_INLINE const CwiseBinaryOp, const DenseDerived, const SparseDerived> +EIGEN_STRONG_INLINE const CwiseBinaryOp, const DenseDerived, const SparseDerived> operator+(const MatrixBase &a, const SparseMatrixBase &b) { - return CwiseBinaryOp, const DenseDerived, const SparseDerived>(a.derived(), b.derived()); + return CwiseBinaryOp, const DenseDerived, const SparseDerived>(a.derived(), b.derived()); } template -EIGEN_STRONG_INLINE const CwiseBinaryOp, const SparseDerived, const DenseDerived> +EIGEN_STRONG_INLINE const CwiseBinaryOp, const SparseDerived, const DenseDerived> operator+(const SparseMatrixBase &a, const MatrixBase &b) { - return CwiseBinaryOp, const SparseDerived, const DenseDerived>(a.derived(), b.derived()); + return CwiseBinaryOp, const SparseDerived, const DenseDerived>(a.derived(), b.derived()); } template -EIGEN_STRONG_INLINE const CwiseBinaryOp, const DenseDerived, const SparseDerived> +EIGEN_STRONG_INLINE const CwiseBinaryOp, const DenseDerived, const SparseDerived> operator-(const MatrixBase &a, const SparseMatrixBase &b) { - return CwiseBinaryOp, const DenseDerived, const SparseDerived>(a.derived(), b.derived()); + return CwiseBinaryOp, const DenseDerived, const SparseDerived>(a.derived(), b.derived()); } template -EIGEN_STRONG_INLINE const CwiseBinaryOp, const SparseDerived, const DenseDerived> +EIGEN_STRONG_INLINE const CwiseBinaryOp, const SparseDerived, const DenseDerived> operator-(const SparseMatrixBase &a, const MatrixBase &b) { - return CwiseBinaryOp, const SparseDerived, const DenseDerived>(a.derived(), b.derived()); + return CwiseBinaryOp, const SparseDerived, const DenseDerived>(a.derived(), b.derived()); } } // end namespace Eigen diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h index 476796dd7..0547db596 100644 --- a/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/Eigen/src/SparseCore/SparseDenseProduct.h @@ -74,7 +74,7 @@ struct sparse_time_dense_product_impl let's disable it for now as it is conflicting with generic scalar*matrix and matrix*scalar operators // template -// struct scalar_product_traits > +// struct ScalarBinaryOpTraits > // { // enum { // Defined = 1 @@ -97,7 +97,7 @@ struct sparse_time_dense_product_impl::ReturnType rhs_j(alpha * rhs.coeff(j,c)); + typename ScalarBinaryOpTraits::ReturnType rhs_j(alpha * rhs.coeff(j,c)); for(LhsInnerIterator it(lhsEval,j); it ;++it) res.coeffRef(it.index(),c) += it.value() * rhs_j; } diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index a78bd57c3..531fea399 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -440,7 +440,7 @@ class SparseMatrix template void setFromTriplets(const InputIterators& begin, const InputIterators& end, DupFunctor dup_func); - void sumupDuplicates() { collapseDuplicates(internal::scalar_sum_op()); } + void sumupDuplicates() { collapseDuplicates(internal::scalar_sum_op()); } template void collapseDuplicates(DupFunctor dup_func = DupFunctor()); @@ -979,7 +979,7 @@ template template void SparseMatrix::setFromTriplets(const InputIterators& begin, const InputIterators& end) { - internal::set_from_triplets >(begin, end, *this, internal::scalar_sum_op()); + internal::set_from_triplets >(begin, end, *this, internal::scalar_sum_op()); } /** The same as setFromTriplets but when duplicates are met the functor \a dup_func is applied: diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index 24df36884..45f64e7f2 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -256,7 +256,7 @@ template class SparseMatrixBase Derived& operator/=(const Scalar& other); template struct CwiseProductDenseReturnType { - typedef CwiseBinaryOp::Scalar, typename internal::traits::Scalar >::ReturnType>, diff --git a/Eigen/src/SparseCore/SparseProduct.h b/Eigen/src/SparseCore/SparseProduct.h index b23003bb1..84e69903b 100644 --- a/Eigen/src/SparseCore/SparseProduct.h +++ b/Eigen/src/SparseCore/SparseProduct.h @@ -99,10 +99,10 @@ struct generic_product_impl -struct Assignment, internal::assign_op, Sparse2Dense> +struct Assignment, internal::assign_op::Scalar>, Sparse2Dense> { typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &) { generic_product_impl::evalTo(dst,src.lhs(),src.rhs()); } @@ -110,10 +110,10 @@ struct Assignment, internal::assig // dense += sparse-product (can be sparse*sparse, sparse*perm, etc.) template< typename DstXprType, typename Lhs, typename Rhs> -struct Assignment, internal::add_assign_op, Sparse2Dense> +struct Assignment, internal::add_assign_op::Scalar>, Sparse2Dense> { typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::add_assign_op &) { generic_product_impl::addTo(dst,src.lhs(),src.rhs()); } @@ -121,10 +121,10 @@ struct Assignment, internal::add_a // dense -= sparse-product (can be sparse*sparse, sparse*perm, etc.) template< typename DstXprType, typename Lhs, typename Rhs> -struct Assignment, internal::sub_assign_op, Sparse2Dense> +struct Assignment, internal::sub_assign_op::Scalar>, Sparse2Dense> { typedef Product SrcXprType; - static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &) + static void run(DstXprType &dst, const SrcXprType &src, const internal::sub_assign_op &) { generic_product_impl::subTo(dst,src.lhs(),src.rhs()); } diff --git a/Eigen/src/SparseCore/SparseSelfAdjointView.h b/Eigen/src/SparseCore/SparseSelfAdjointView.h index b92bb17e2..4f0c84d88 100644 --- a/Eigen/src/SparseCore/SparseSelfAdjointView.h +++ b/Eigen/src/SparseCore/SparseSelfAdjointView.h @@ -223,13 +223,13 @@ struct Assignment - static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &/*func*/) { internal::permute_symm_to_fullsymm(src.matrix(), dst); } template - static void run(DynamicSparseMatrix& dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static void run(DynamicSparseMatrix& dst, const SrcXprType &src, const internal::assign_op &/*func*/) { // TODO directly evaluate into dst; SparseMatrix tmp(dst.rows(),dst.cols()); @@ -586,12 +586,12 @@ class SparseSymmetricPermutationProduct namespace internal { template -struct Assignment, internal::assign_op, Sparse2Sparse> +struct Assignment, internal::assign_op, Sparse2Sparse> { typedef SparseSymmetricPermutationProduct SrcXprType; typedef typename DstXprType::StorageIndex DstIndex; template - static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &) + static void run(SparseMatrix &dst, const SrcXprType &src, const internal::assign_op &) { // internal::permute_symm_to_fullsymm(m_matrix,_dest,m_perm.indices().data()); SparseMatrix tmp; @@ -600,7 +600,7 @@ struct Assignment } template - static void run(SparseSelfAdjointView& dst, const SrcXprType &src, const internal::assign_op &) + static void run(SparseSelfAdjointView& dst, const SrcXprType &src, const internal::assign_op &) { internal::permute_symm_to_symm(src.matrix(),dst.matrix(),src.perm().indices().data()); } diff --git a/Eigen/src/SparseQR/SparseQR.h b/Eigen/src/SparseQR/SparseQR.h index acd7f7e10..2d4498b03 100644 --- a/Eigen/src/SparseQR/SparseQR.h +++ b/Eigen/src/SparseQR/SparseQR.h @@ -705,12 +705,12 @@ struct evaluator_traits > }; template< typename DstXprType, typename SparseQRType> -struct Assignment, internal::assign_op, Sparse2Sparse> +struct Assignment, internal::assign_op, Sparse2Sparse> { typedef SparseQRMatrixQReturnType SrcXprType; typedef typename DstXprType::Scalar Scalar; typedef typename DstXprType::StorageIndex StorageIndex; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) { typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows()); idMat.setIdentity(); @@ -721,12 +721,12 @@ struct Assignment, internal: }; template< typename DstXprType, typename SparseQRType> -struct Assignment, internal::assign_op, Sparse2Dense> +struct Assignment, internal::assign_op, Sparse2Dense> { typedef SparseQRMatrixQReturnType SrcXprType; typedef typename DstXprType::Scalar Scalar; typedef typename DstXprType::StorageIndex StorageIndex; - static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) + static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op &/*func*/) { dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows()); } diff --git a/Eigen/src/plugins/ArrayCwiseBinaryOps.h b/Eigen/src/plugins/ArrayCwiseBinaryOps.h index c3f8c2575..c6ea12c90 100644 --- a/Eigen/src/plugins/ArrayCwiseBinaryOps.h +++ b/Eigen/src/plugins/ArrayCwiseBinaryOps.h @@ -4,10 +4,10 @@ */ template EIGEN_DEVICE_FUNC -EIGEN_STRONG_INLINE const EIGEN_CWISE_PRODUCT_RETURN_TYPE(Derived,OtherDerived) +EIGEN_STRONG_INLINE const EIGEN_CWISE_BINARY_RETURN_TYPE(Derived,OtherDerived,product) operator*(const EIGEN_CURRENT_STORAGE_BASE_CLASS &other) const { - return EIGEN_CWISE_PRODUCT_RETURN_TYPE(Derived,OtherDerived)(derived(), other.derived()); + return EIGEN_CWISE_BINARY_RETURN_TYPE(Derived,OtherDerived,product)(derived(), other.derived()); } /** \returns an expression of the coefficient wise quotient of \c *this and \a other diff --git a/Eigen/src/plugins/CommonCwiseBinaryOps.h b/Eigen/src/plugins/CommonCwiseBinaryOps.h index a8fa287c9..3c13b7f06 100644 --- a/Eigen/src/plugins/CommonCwiseBinaryOps.h +++ b/Eigen/src/plugins/CommonCwiseBinaryOps.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008-2009 Gael Guennebaud +// Copyright (C) 2008-2016 Gael Guennebaud // Copyright (C) 2006-2008 Benoit Jacob // // This Source Code Form is subject to the terms of the Mozilla @@ -16,7 +16,13 @@ * * \sa class CwiseBinaryOp, operator-=() */ -EIGEN_MAKE_CWISE_BINARY_OP(operator-,internal::scalar_difference_op) +template +EIGEN_DEVICE_FUNC +EIGEN_STRONG_INLINE const EIGEN_CWISE_BINARY_RETURN_TYPE(Derived,OtherDerived,difference) +operator-(const EIGEN_CURRENT_STORAGE_BASE_CLASS &other) const +{ + return EIGEN_CWISE_BINARY_RETURN_TYPE(Derived,OtherDerived,difference)(derived(), other.derived()); +} /** \returns an expression of the sum of \c *this and \a other * @@ -24,7 +30,13 @@ EIGEN_MAKE_CWISE_BINARY_OP(operator-,internal::scalar_difference_op) * * \sa class CwiseBinaryOp, operator+=() */ -EIGEN_MAKE_CWISE_BINARY_OP(operator+,internal::scalar_sum_op) +template +EIGEN_DEVICE_FUNC +EIGEN_STRONG_INLINE const EIGEN_CWISE_BINARY_RETURN_TYPE(Derived,OtherDerived,sum) +operator+(const EIGEN_CURRENT_STORAGE_BASE_CLASS &other) const +{ + return EIGEN_CWISE_BINARY_RETURN_TYPE(Derived,OtherDerived,sum)(derived(), other.derived()); +} /** \returns an expression of a custom coefficient-wise operator \a func of *this and \a other * diff --git a/Eigen/src/plugins/CommonCwiseUnaryOps.h b/Eigen/src/plugins/CommonCwiseUnaryOps.h index 67ec601b9..927167aff 100644 --- a/Eigen/src/plugins/CommonCwiseUnaryOps.h +++ b/Eigen/src/plugins/CommonCwiseUnaryOps.h @@ -73,7 +73,7 @@ operator/(const Scalar& scalar) const /** Overloaded for efficiently multipling with compatible scalar types */ template EIGEN_DEVICE_FUNC inline -typename internal::enable_if::Defined, +typename internal::enable_if::Defined, const CwiseUnaryOp, const Derived> >::type operator*(const T& scalar) const { @@ -91,7 +91,7 @@ operator*(const Scalar& scalar, const StorageBaseType& matrix) template EIGEN_DEVICE_FUNC inline friend -typename internal::enable_if::Defined, +typename internal::enable_if::Defined, const CwiseUnaryOp, const Derived> >::type operator*(const T& scalar, const StorageBaseType& matrix) { @@ -104,7 +104,7 @@ operator*(const T& scalar, const StorageBaseType& matrix) template EIGEN_DEVICE_FUNC inline -typename internal::enable_if::Defined, +typename internal::enable_if::Defined, const CwiseUnaryOp, const Derived> >::type operator/(const T& scalar) const { diff --git a/Eigen/src/plugins/MatrixCwiseBinaryOps.h b/Eigen/src/plugins/MatrixCwiseBinaryOps.h index 6dd2e1192..59581e618 100644 --- a/Eigen/src/plugins/MatrixCwiseBinaryOps.h +++ b/Eigen/src/plugins/MatrixCwiseBinaryOps.h @@ -19,10 +19,10 @@ */ template EIGEN_DEVICE_FUNC -EIGEN_STRONG_INLINE const EIGEN_CWISE_PRODUCT_RETURN_TYPE(Derived,OtherDerived) +EIGEN_STRONG_INLINE const EIGEN_CWISE_BINARY_RETURN_TYPE(Derived,OtherDerived,product) cwiseProduct(const EIGEN_CURRENT_STORAGE_BASE_CLASS &other) const { - return EIGEN_CWISE_PRODUCT_RETURN_TYPE(Derived,OtherDerived)(derived(), other.derived()); + return EIGEN_CWISE_BINARY_RETURN_TYPE(Derived,OtherDerived,product)(derived(), other.derived()); } /** \returns an expression of the coefficient-wise == operator of *this and \a other diff --git a/blas/PackedTriangularMatrixVector.h b/blas/PackedTriangularMatrixVector.h index e9886d56f..0039536a8 100644 --- a/blas/PackedTriangularMatrixVector.h +++ b/blas/PackedTriangularMatrixVector.h @@ -18,7 +18,7 @@ struct packed_triangular_matrix_vector_product; template struct packed_triangular_matrix_vector_product { - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; enum { IsLower = (Mode & Lower) ==Lower, HasUnitDiag = (Mode & UnitDiag)==UnitDiag, @@ -47,7 +47,7 @@ struct packed_triangular_matrix_vector_product struct packed_triangular_matrix_vector_product { - typedef typename scalar_product_traits::ReturnType ResScalar; + typedef typename ScalarBinaryOpTraits::ReturnType ResScalar; enum { IsLower = (Mode & Lower) ==Lower, HasUnitDiag = (Mode & UnitDiag)==UnitDiag, diff --git a/test/array.cpp b/test/array.cpp index 4cd4f262b..bd470d5f7 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -72,7 +72,7 @@ template void array(const ArrayType& m) VERIFY_IS_MUCH_SMALLER_THAN(abs(m1.rowwise().sum().sum() - m1.sum()), m1.abs().sum()); if (!internal::isMuchSmallerThan(abs(m1.sum() - (m1+m2).sum()), m1.abs().sum(), test_precision())) VERIFY_IS_NOT_APPROX(((m1+m2).rowwise().sum()).sum(), m1.sum()); - VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op())); + VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op())); // vector-wise ops m3 = m1; diff --git a/test/array_for_matrix.cpp b/test/array_for_matrix.cpp index 75e6a778f..97e03be83 100644 --- a/test/array_for_matrix.cpp +++ b/test/array_for_matrix.cpp @@ -45,7 +45,7 @@ template void array_for_matrix(const MatrixType& m) VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum().sum() - m1.sum(), m1.squaredNorm()); VERIFY_IS_MUCH_SMALLER_THAN(m1.colwise().sum() + m2.colwise().sum() - (m1+m2).colwise().sum(), (m1+m2).squaredNorm()); VERIFY_IS_MUCH_SMALLER_THAN(m1.rowwise().sum() - m2.rowwise().sum() - (m1-m2).rowwise().sum(), (m1-m2).squaredNorm()); - VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op())); + VERIFY_IS_APPROX(m1.colwise().sum(), m1.colwise().redux(internal::scalar_sum_op())); // vector-wise ops m3 = m1; diff --git a/test/mixingtypes.cpp b/test/mixingtypes.cpp index dbcf468ea..66d9b777a 100644 --- a/test/mixingtypes.cpp +++ b/test/mixingtypes.cpp @@ -42,6 +42,7 @@ template void mixingtypes(int size = SizeAtCompileType) Mat_f mf = Mat_f::Random(size,size); Mat_d md = mf.template cast(); + //Mat_d rd = md; Mat_cf mcf = Mat_cf::Random(size,size); Mat_cd mcd = mcf.template cast >(); Mat_cd rcd = mcd; @@ -56,16 +57,12 @@ template void mixingtypes(int size = SizeAtCompileType) mf+mf; - VERIFY_RAISES_ASSERT(mf+md); -#if !EIGEN_HAS_STD_RESULT_OF - // this one does not even compile with C++11 - VERIFY_RAISES_ASSERT(mf+mcf); -#endif + +// VERIFY_RAISES_ASSERT(mf+md); // does not even compile #ifdef EIGEN_DONT_VECTORIZE VERIFY_RAISES_ASSERT(vf=vd); VERIFY_RAISES_ASSERT(vf+=vd); - VERIFY_RAISES_ASSERT(mcd=md); #endif // check scalar products @@ -186,16 +183,35 @@ template void mixingtypes(int size = SizeAtCompileType) Mat_cd((scd * md.template cast().eval() * mcd).template triangularView())); - VERIFY_IS_APPROX( md.array() * mcd.array(), md.template cast().eval().array() * mcd.array() ); - VERIFY_IS_APPROX( mcd.array() * md.array(), mcd.array() * md.template cast().eval().array() ); + + VERIFY_IS_APPROX( md.array() * mcd.array(), md.template cast().eval().array() * mcd.array() ); + VERIFY_IS_APPROX( mcd.array() * md.array(), mcd.array() * md.template cast().eval().array() ); + + VERIFY_IS_APPROX( md.array() + mcd.array(), md.template cast().eval().array() + mcd.array() ); + VERIFY_IS_APPROX( mcd.array() + md.array(), mcd.array() + md.template cast().eval().array() ); + + VERIFY_IS_APPROX( md.array() - mcd.array(), md.template cast().eval().array() - mcd.array() ); + VERIFY_IS_APPROX( mcd.array() - md.array(), mcd.array() - md.template cast().eval().array() ); // VERIFY_IS_APPROX( md.array() / mcd.array(), md.template cast().eval().array() / mcd.array() ); VERIFY_IS_APPROX( mcd.array() / md.array(), mcd.array() / md.template cast().eval().array() ); + rcd = mcd; + VERIFY_IS_APPROX( rcd = md, md.template cast().eval() ); + rcd = mcd; + VERIFY_IS_APPROX( rcd += md, mcd + md.template cast().eval() ); + rcd = mcd; + VERIFY_IS_APPROX( rcd -= md, mcd - md.template cast().eval() ); rcd = mcd; VERIFY_IS_APPROX( rcd.array() *= md.array(), mcd.array() * md.template cast().eval().array() ); rcd = mcd; VERIFY_IS_APPROX( rcd.array() /= md.array(), mcd.array() / md.template cast().eval().array() ); + + rcd = mcd; + VERIFY_IS_APPROX( rcd += md + mcd*md, mcd + (md.template cast().eval()) + mcd*(md.template cast().eval())); + + rcd = mcd; + VERIFY_IS_APPROX( rcd += mcd + md*md, mcd + mcd + ((md*md).template cast().eval()) ); } void test_mixingtypes() diff --git a/test/vectorization_logic.cpp b/test/vectorization_logic.cpp index 24a7641ff..b7c2df64b 100644 --- a/test/vectorization_logic.cpp +++ b/test/vectorization_logic.cpp @@ -29,7 +29,7 @@ using internal::demangle_unrolling; template bool test_assign(const Dst&, const Src&, int traversal, int unrolling) { - typedef internal::copy_using_evaluator_traits,internal::evaluator, internal::assign_op > traits; + typedef internal::copy_using_evaluator_traits,internal::evaluator, internal::assign_op > traits; bool res = traits::Traversal==traversal; if(unrolling==InnerUnrolling+CompleteUnrolling) res = res && (int(traits::Unrolling)==InnerUnrolling || int(traits::Unrolling)==CompleteUnrolling); @@ -53,7 +53,7 @@ bool test_assign(const Dst&, const Src&, int traversal, int unrolling) template bool test_assign(int traversal, int unrolling) { - typedef internal::copy_using_evaluator_traits,internal::evaluator, internal::assign_op > traits; + typedef internal::copy_using_evaluator_traits,internal::evaluator, internal::assign_op > traits; bool res = traits::Traversal==traversal && traits::Unrolling==unrolling; if(!res) { @@ -73,7 +73,8 @@ bool test_assign(int traversal, int unrolling) template bool test_redux(const Xpr&, int traversal, int unrolling) { - typedef internal::redux_traits,internal::redux_evaluator > traits; + typedef typename Xpr::Scalar Scalar; + typedef internal::redux_traits,internal::redux_evaluator > traits; bool res = traits::Traversal==traversal && traits::Unrolling==unrolling; if(!res) diff --git a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h index 089042751..ba61288a3 100755 --- a/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h +++ b/unsupported/Eigen/src/AutoDiff/AutoDiffScalar.h @@ -501,36 +501,36 @@ struct make_coherent_impl -struct scalar_product_traits,A_Scalar> +struct ScalarBinaryOpTraits,A_Scalar> { enum { Defined = 1 }; typedef Matrix ReturnType; }; template -struct scalar_product_traits > +struct ScalarBinaryOpTraits > { enum { Defined = 1 }; typedef Matrix ReturnType; }; template -struct scalar_product_traits,typename DerType::Scalar> +struct ScalarBinaryOpTraits,typename DerType::Scalar> { enum { Defined = 1 }; typedef AutoDiffScalar ReturnType; }; template -struct scalar_product_traits > +struct ScalarBinaryOpTraits > { enum { Defined = 1 }; typedef AutoDiffScalar ReturnType; }; -} // end namespace internal - #define EIGEN_AUTODIFF_DECLARE_GLOBAL_UNARY(FUNC,CODE) \ template \ inline const Eigen::AutoDiffScalar::type>::Scalar>, const typename Eigen::internal::remove_all::type> > \ diff --git a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h index bf9727c21..582fa8512 100644 --- a/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h +++ b/unsupported/Eigen/src/KroneckerProduct/KroneckerTensorProduct.h @@ -203,7 +203,7 @@ struct traits > { typedef typename remove_all<_Lhs>::type Lhs; typedef typename remove_all<_Rhs>::type Rhs; - typedef typename scalar_product_traits::ReturnType Scalar; + typedef typename ScalarBinaryOpTraits::ReturnType Scalar; typedef typename promote_index_type::type StorageIndex; enum { @@ -222,7 +222,7 @@ struct traits > typedef MatrixXpr XprKind; typedef typename remove_all<_Lhs>::type Lhs; typedef typename remove_all<_Rhs>::type Rhs; - typedef typename scalar_product_traits::ReturnType Scalar; + typedef typename ScalarBinaryOpTraits::ReturnType Scalar; typedef typename cwise_promote_storage_type::StorageKind, typename traits::StorageKind, scalar_product_op >::ret StorageKind; typedef typename promote_index_type::type StorageIndex; -- cgit v1.2.3 From 32a41ee659686fe1fb76156f7a55287acf14d4bb Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 4 Jul 2016 15:13:35 +0200 Subject: bug #707: add inplace decomposition through Ref<> for Cholesky, LU and QR decompositions. --- Eigen/src/Cholesky/LDLT.h | 21 ++++- Eigen/src/Cholesky/LLT.h | 16 +++- Eigen/src/LU/FullPivLU.h | 44 +++++++--- Eigen/src/LU/PartialPivLU.h | 62 ++++++++++++-- Eigen/src/QR/ColPivHouseholderQR.h | 37 ++++++--- Eigen/src/QR/CompleteOrthogonalDecomposition.h | 47 +++++++---- Eigen/src/QR/FullPivHouseholderQR.h | 27 ++++-- Eigen/src/QR/HouseholderQR.h | 36 ++++++-- test/CMakeLists.txt | 1 + test/inplace_decomposition.cpp | 110 +++++++++++++++++++++++++ 10 files changed, 337 insertions(+), 64 deletions(-) create mode 100644 test/inplace_decomposition.cpp (limited to 'Eigen/src/LU/FullPivLU.h') diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 538aff956..a31b3d6aa 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -52,7 +52,6 @@ template class LDLT enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here! MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, UpLo = _UpLo @@ -61,7 +60,7 @@ template class LDLT typedef typename NumTraits::Real RealScalar; typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 typedef typename MatrixType::StorageIndex StorageIndex; - typedef Matrix TmpMatrixType; + typedef Matrix TmpMatrixType; typedef Transpositions TranspositionType; typedef PermutationMatrix PermutationType; @@ -97,6 +96,7 @@ template class LDLT /** \brief Constructor with decomposition * * This calculates the decomposition for the input \a matrix. + * * \sa LDLT(Index size) */ template @@ -110,6 +110,23 @@ template class LDLT compute(matrix.derived()); } + /** \brief Constructs a LDLT factorization from a given matrix + * + * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * + * \sa LDLT(const EigenBase&) + */ + template + explicit LDLT(EigenBase& matrix) + : m_matrix(matrix.derived()), + m_transpositions(matrix.rows()), + m_temporary(matrix.rows()), + m_sign(internal::ZeroSign), + m_isInitialized(false) + { + compute(matrix.derived()); + } + /** Clear any existing decomposition * \sa rankUpdate(w,sigma) */ diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 19578b216..ad163c749 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -54,7 +54,6 @@ template class LLT enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; typedef typename MatrixType::Scalar Scalar; @@ -95,6 +94,21 @@ template class LLT compute(matrix.derived()); } + /** \brief Constructs a LDLT factorization from a given matrix + * + * This overloaded constructor is provided for inplace solving when + * \c MatrixType is a Eigen::Ref. + * + * \sa LLT(const EigenBase&) + */ + template + explicit LLT(EigenBase& matrix) + : m_matrix(matrix.derived()), + m_isInitialized(false) + { + compute(matrix.derived()); + } + /** \returns a view of the upper triangular matrix U */ inline typename Traits::MatrixU matrixU() const { diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 2d01b18c6..113b8c7b8 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -97,6 +97,15 @@ template class FullPivLU template explicit FullPivLU(const EigenBase& matrix); + /** \brief Constructs a LU factorization from a given matrix + * + * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * + * \sa FullPivLU(const EigenBase&) + */ + template + explicit FullPivLU(EigenBase& matrix); + /** Computes the LU decomposition of the given matrix. * * \param matrix the matrix of which to compute the LU decomposition. @@ -105,7 +114,11 @@ template class FullPivLU * \returns a reference to *this */ template - FullPivLU& compute(const EigenBase& matrix); + FullPivLU& compute(const EigenBase& matrix) { + m_lu = matrix.derived(); + computeInPlace(); + return *this; + } /** \returns the LU decomposition matrix: the upper-triangular part is U, the * unit-lower-triangular part is L (at least for square matrices; in the non-square @@ -459,25 +472,28 @@ FullPivLU::FullPivLU(const EigenBase& matrix) template template -FullPivLU& FullPivLU::compute(const EigenBase& matrix) +FullPivLU::FullPivLU(EigenBase& matrix) + : m_lu(matrix.derived()), + m_p(matrix.rows()), + m_q(matrix.cols()), + m_rowsTranspositions(matrix.rows()), + m_colsTranspositions(matrix.cols()), + m_isInitialized(false), + m_usePrescribedThreshold(false) { - check_template_parameters(); - - // the permutations are stored as int indices, so just to be sure: - eigen_assert(matrix.rows()<=NumTraits::highest() && matrix.cols()<=NumTraits::highest()); - - m_lu = matrix.derived(); - m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff(); - computeInPlace(); - - m_isInitialized = true; - return *this; } template void FullPivLU::computeInPlace() { + check_template_parameters(); + + // the permutations are stored as int indices, so just to be sure: + eigen_assert(m_lu.rows()<=NumTraits::highest() && m_lu.cols()<=NumTraits::highest()); + + m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff(); + const Index size = m_lu.diagonalSize(); const Index rows = m_lu.rows(); const Index cols = m_lu.cols(); @@ -557,6 +573,8 @@ void FullPivLU::computeInPlace() m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k)); m_det_pq = (number_of_transpositions%2) ? -1 : 1; + + m_isInitialized = true; } template diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index ac2902261..c862d9692 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -26,6 +26,17 @@ template struct traits > }; }; +template +struct enable_if_ref; +// { +// typedef Derived type; +// }; + +template +struct enable_if_ref,Derived> { + typedef Derived type; +}; + } // end namespace internal /** \ingroup LU_Module @@ -102,8 +113,29 @@ template class PartialPivLU template explicit PartialPivLU(const EigenBase& matrix); + /** Constructor for inplace decomposition + * + * \param matrix the matrix of which to compute the LU decomposition. + * + * If \c MatrixType is an Eigen::Ref, then the storage of \a matrix will be shared + * between \a matrix and \c *this and the decomposition will take place in-place. + * The memory of \a matrix will be used througrough the lifetime of \c *this. In + * particular, further calls to \c this->compute(A) will still operate on the memory + * of \a matrix meaning. This also implies that the sizes of \c A must match the + * ones of \a matrix. + * + * \warning The matrix should have full rank (e.g. if it's square, it should be invertible). + * If you need to deal with non-full rank, use class FullPivLU instead. + */ + template + explicit PartialPivLU(EigenBase& matrix); + template - PartialPivLU& compute(const EigenBase& matrix); + PartialPivLU& compute(const EigenBase& matrix) { + m_lu = matrix.derived(); + compute(); + return *this; + } /** \returns the LU decomposition matrix: the upper-triangular part is U, the * unit-lower-triangular part is L (at least for square matrices; in the non-square @@ -251,6 +283,8 @@ template class PartialPivLU EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); } + void compute(); + MatrixType m_lu; PermutationType m_p; TranspositionType m_rowsTranspositions; @@ -284,7 +318,7 @@ PartialPivLU::PartialPivLU(Index size) template template PartialPivLU::PartialPivLU(const EigenBase& matrix) - : m_lu(matrix.rows(), matrix.rows()), + : m_lu(matrix.rows(),matrix.cols()), m_p(matrix.rows()), m_rowsTranspositions(matrix.rows()), m_l1_norm(0), @@ -294,6 +328,19 @@ PartialPivLU::PartialPivLU(const EigenBase& matrix) compute(matrix.derived()); } +template +template +PartialPivLU::PartialPivLU(EigenBase& matrix) + : m_lu(matrix.derived()), + m_p(matrix.rows()), + m_rowsTranspositions(matrix.rows()), + m_l1_norm(0), + m_det_p(0), + m_isInitialized(false) +{ + compute(); +} + namespace internal { /** \internal This is the blocked version of fullpivlu_unblocked() */ @@ -470,19 +517,17 @@ void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, t } // end namespace internal template -template -PartialPivLU& PartialPivLU::compute(const EigenBase& matrix) +void PartialPivLU::compute() { check_template_parameters(); // the row permutation is stored as int indices, so just to be sure: - eigen_assert(matrix.rows()::highest()); + eigen_assert(m_lu.rows()::highest()); - m_lu = matrix.derived(); m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff(); - eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); - const Index size = matrix.rows(); + eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); + const Index size = m_lu.rows(); m_rowsTranspositions.resize(size); @@ -493,7 +538,6 @@ PartialPivLU& PartialPivLU::compute(const EigenBase diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index e847bc434..db50b5675 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -51,7 +51,6 @@ template class ColPivHouseholderQR enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; @@ -59,7 +58,6 @@ template class ColPivHouseholderQR typedef typename MatrixType::RealScalar RealScalar; // FIXME should be int typedef typename MatrixType::StorageIndex StorageIndex; - typedef Matrix MatrixQType; typedef typename internal::plain_diag_type::type HCoeffsType; typedef PermutationMatrix PermutationType; typedef typename internal::plain_row_type::type IntRowVectorType; @@ -135,6 +133,27 @@ template class ColPivHouseholderQR compute(matrix.derived()); } + /** \brief Constructs a QR factorization from a given matrix + * + * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * + * \sa ColPivHouseholderQR(const EigenBase&) + */ + template + explicit ColPivHouseholderQR(EigenBase& matrix) + : m_qr(matrix.derived()), + m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), + m_colsPermutation(PermIndexType(matrix.cols())), + m_colsTranspositions(matrix.cols()), + m_temp(matrix.cols()), + m_colNormsUpdated(matrix.cols()), + m_colNormsDirect(matrix.cols()), + m_isInitialized(false), + m_usePrescribedThreshold(false) + { + computeInPlace(); + } + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which * *this is the QR decomposition, if any exists. * @@ -453,21 +472,19 @@ template template ColPivHouseholderQR& ColPivHouseholderQR::compute(const EigenBase& matrix) { - check_template_parameters(); - - // the column permutation is stored as int indices, so just to be sure: - eigen_assert(matrix.cols()<=NumTraits::highest()); - - m_qr = matrix; - + m_qr = matrix.derived(); computeInPlace(); - return *this; } template void ColPivHouseholderQR::computeInPlace() { + check_template_parameters(); + + // the column permutation is stored as int indices, so just to be sure: + eigen_assert(m_qr.cols()<=NumTraits::highest()); + using std::abs; Index rows = m_qr.rows(); diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h index 398e1aa77..967eb35dd 100644 --- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h +++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -48,16 +48,12 @@ class CompleteOrthogonalDecomposition { enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::StorageIndex StorageIndex; - typedef Matrix - MatrixQType; typedef typename internal::plain_diag_type::type HCoeffsType; typedef PermutationMatrix PermutationType; @@ -114,10 +110,27 @@ class CompleteOrthogonalDecomposition { explicit CompleteOrthogonalDecomposition(const EigenBase& matrix) : m_cpqr(matrix.rows(), matrix.cols()), m_zCoeffs((std::min)(matrix.rows(), matrix.cols())), - m_temp(matrix.cols()) { + m_temp(matrix.cols()) + { compute(matrix.derived()); } + /** \brief Constructs a complete orthogonal decomposition from a given matrix + * + * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * + * \sa CompleteOrthogonalDecomposition(const EigenBase&) + */ + template + explicit CompleteOrthogonalDecomposition(EigenBase& matrix) + : m_cpqr(matrix.derived()), + m_zCoeffs((std::min)(matrix.rows(), matrix.cols())), + m_temp(matrix.cols()) + { + computeInPlace(); + } + + /** This method computes the minimum-norm solution X to a least squares * problem \f[\mathrm{minimize} ||A X - B|| \f], where \b A is the matrix of * which \c *this is the complete orthogonal decomposition. @@ -165,7 +178,12 @@ class CompleteOrthogonalDecomposition { const MatrixType& matrixT() const { return m_cpqr.matrixQR(); } template - CompleteOrthogonalDecomposition& compute(const EigenBase& matrix); + CompleteOrthogonalDecomposition& compute(const EigenBase& matrix) { + // Compute the column pivoted QR factorization A P = Q R. + m_cpqr.compute(matrix); + computeInPlace(); + return *this; + } /** \returns a const reference to the column permutation matrix */ const PermutationType& colsPermutation() const { @@ -354,6 +372,8 @@ class CompleteOrthogonalDecomposition { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); } + void computeInPlace(); + /** Overwrites \b rhs with \f$ \mathbf{Z}^* * \mathbf{rhs} \f$. */ template @@ -384,20 +404,16 @@ CompleteOrthogonalDecomposition::logAbsDeterminant() const { * CompleteOrthogonalDecomposition(const MatrixType&) */ template -template -CompleteOrthogonalDecomposition& CompleteOrthogonalDecomposition< - MatrixType>::compute(const EigenBase& matrix) { +void CompleteOrthogonalDecomposition::computeInPlace() +{ check_template_parameters(); // the column permutation is stored as int indices, so just to be sure: - eigen_assert(matrix.cols() <= NumTraits::highest()); - - // Compute the column pivoted QR factorization A P = Q R. - m_cpqr.compute(matrix); + eigen_assert(m_cpqr.cols() <= NumTraits::highest()); const Index rank = m_cpqr.rank(); - const Index cols = matrix.cols(); - const Index rows = matrix.rows(); + const Index cols = m_cpqr.cols(); + const Index rows = m_cpqr.rows(); m_zCoeffs.resize((std::min)(rows, cols)); m_temp.resize(cols); @@ -443,7 +459,6 @@ CompleteOrthogonalDecomposition& CompleteOrthogonalDecomposition< } } } - return *this; } template diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index e21966056..5c2f57d04 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -60,7 +60,6 @@ template class FullPivHouseholderQR enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; @@ -135,6 +134,26 @@ template class FullPivHouseholderQR compute(matrix.derived()); } + /** \brief Constructs a QR factorization from a given matrix + * + * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * + * \sa FullPivHouseholderQR(const EigenBase&) + */ + template + explicit FullPivHouseholderQR(EigenBase& matrix) + : m_qr(matrix.derived()), + m_hCoeffs((std::min)(matrix.rows(), matrix.cols())), + m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())), + m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())), + m_cols_permutation(matrix.cols()), + m_temp(matrix.cols()), + m_isInitialized(false), + m_usePrescribedThreshold(false) + { + computeInPlace(); + } + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which * \c *this is the QR decomposition. * @@ -430,18 +449,16 @@ template template FullPivHouseholderQR& FullPivHouseholderQR::compute(const EigenBase& matrix) { - check_template_parameters(); - m_qr = matrix.derived(); - computeInPlace(); - return *this; } template void FullPivHouseholderQR::computeInPlace() { + check_template_parameters(); + using std::abs; Index rows = m_qr.rows(); Index cols = m_qr.cols(); diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 03bc8e6cd..f2a9cc080 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -47,7 +47,6 @@ template class HouseholderQR enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; @@ -102,6 +101,24 @@ template class HouseholderQR compute(matrix.derived()); } + + /** \brief Constructs a QR factorization from a given matrix + * + * This overloaded constructor is provided for inplace solving when + * \c MatrixType is a Eigen::Ref. + * + * \sa HouseholderQR(const EigenBase&) + */ + template + explicit HouseholderQR(EigenBase& matrix) + : m_qr(matrix.derived()), + m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), + m_temp(matrix.cols()), + m_isInitialized(false) + { + computeInPlace(); + } + /** This method finds a solution x to the equation Ax=b, where A is the matrix of which * *this is the QR decomposition, if any exists. * @@ -151,7 +168,11 @@ template class HouseholderQR } template - HouseholderQR& compute(const EigenBase& matrix); + HouseholderQR& compute(const EigenBase& matrix) { + m_qr = matrix.derived(); + computeInPlace(); + return *this; + } /** \returns the absolute value of the determinant of the matrix of which * *this is the QR decomposition. It has only linear complexity @@ -203,6 +224,8 @@ template class HouseholderQR { EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); } + + void computeInPlace(); MatrixType m_qr; HCoeffsType m_hCoeffs; @@ -354,16 +377,14 @@ void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) c * \sa class HouseholderQR, HouseholderQR(const MatrixType&) */ template -template -HouseholderQR& HouseholderQR::compute(const EigenBase& matrix) +void HouseholderQR::computeInPlace() { check_template_parameters(); - Index rows = matrix.rows(); - Index cols = matrix.cols(); + Index rows = m_qr.rows(); + Index cols = m_qr.cols(); Index size = (std::min)(rows,cols); - m_qr = matrix.derived(); m_hCoeffs.resize(size); m_temp.resize(cols); @@ -371,7 +392,6 @@ HouseholderQR& HouseholderQR::compute(const EigenBase::run(m_qr, m_hCoeffs, 48, m_temp.data()); m_isInitialized = true; - return *this; } #ifndef __CUDACC__ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 9d49f1e97..27a9ec4ac 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -258,6 +258,7 @@ ei_add_test(rvalue_types) ei_add_test(dense_storage) ei_add_test(ctorleak) ei_add_test(mpl2only) +ei_add_test(inplace_decomposition) add_executable(bug1213 bug1213.cpp bug1213_main.cpp) diff --git a/test/inplace_decomposition.cpp b/test/inplace_decomposition.cpp new file mode 100644 index 000000000..4900312f1 --- /dev/null +++ b/test/inplace_decomposition.cpp @@ -0,0 +1,110 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2016 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/. + +#include "main.h" +#include +#include +#include + +// This file test inplace decomposition through Ref<>, as supported by Cholesky, LU, and QR decompositions. + +template void inplace(bool square = false, bool SPD = false) +{ + typedef typename MatrixType::Scalar Scalar; + typedef Matrix RhsType; + typedef Matrix ResType; + + Index rows = MatrixType::RowsAtCompileTime==Dynamic ? internal::random(2,EIGEN_TEST_MAX_SIZE/2) : MatrixType::RowsAtCompileTime; + Index cols = MatrixType::ColsAtCompileTime==Dynamic ? (square?rows:internal::random(2,rows)) : MatrixType::ColsAtCompileTime; + + MatrixType A = MatrixType::Random(rows,cols); + RhsType b = RhsType::Random(rows); + ResType x(cols); + + if(SPD) + { + assert(square); + A.topRows(cols) = A.topRows(cols).adjoint() * A.topRows(cols); + A.diagonal().array() += 1e-3; + } + + MatrixType A0 = A; + MatrixType A1 = A; + + DecType dec(A); + + // Check that the content of A has been modified + VERIFY_IS_NOT_APPROX( A, A0 ); + + // Check that the decomposition is correct: + if(rows==cols) + { + VERIFY_IS_APPROX( A0 * (x = dec.solve(b)), b ); + } + else + { + VERIFY_IS_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); + } + + // Check that modifying A breaks the current dec: + A.setRandom(); + if(rows==cols) + { + VERIFY_IS_NOT_APPROX( A0 * (x = dec.solve(b)), b ); + } + else + { + VERIFY_IS_NOT_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); + } + + // Check that calling compute(A1) does not modify A1: + A = A0; + dec.compute(A1); + VERIFY_IS_EQUAL(A0,A1); + VERIFY_IS_NOT_APPROX( A, A0 ); + if(rows==cols) + { + VERIFY_IS_APPROX( A0 * (x = dec.solve(b)), b ); + } + else + { + VERIFY_IS_APPROX( A0.transpose() * A0 * (x = dec.solve(b)), A0.transpose() * b ); + } +} + + +void test_inplace_decomposition() +{ + EIGEN_UNUSED typedef Matrix Matrix43d; + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST_1(( inplace >, MatrixXd>(true,true) )); + CALL_SUBTEST_1(( inplace >, Matrix4d>(true,true) )); + + CALL_SUBTEST_2(( inplace >, MatrixXd>(true,true) )); + CALL_SUBTEST_2(( inplace >, Matrix4d>(true,true) )); + + CALL_SUBTEST_3(( inplace >, MatrixXd>(true,false) )); + CALL_SUBTEST_3(( inplace >, Matrix4d>(true,false) )); + + CALL_SUBTEST_4(( inplace >, MatrixXd>(true,false) )); + CALL_SUBTEST_4(( inplace >, Matrix4d>(true,false) )); + + CALL_SUBTEST_5(( inplace >, MatrixXd>(false,false) )); + CALL_SUBTEST_5(( inplace >, Matrix43d>(false,false) )); + + CALL_SUBTEST_6(( inplace >, MatrixXd>(false,false) )); + CALL_SUBTEST_6(( inplace >, Matrix43d>(false,false) )); + + CALL_SUBTEST_7(( inplace >, MatrixXd>(false,false) )); + CALL_SUBTEST_7(( inplace >, Matrix43d>(false,false) )); + + CALL_SUBTEST_8(( inplace >, MatrixXd>(false,false) )); + CALL_SUBTEST_8(( inplace >, Matrix43d>(false,false) )); + } +} -- cgit v1.2.3 From 7f7839c12f9bbd88bee1f3c9d7f5dbbe556045d4 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 4 Jul 2016 17:18:26 +0200 Subject: Add documentation and exemples for inplace decomposition. --- Eigen/src/Cholesky/LDLT.h | 4 +- Eigen/src/Cholesky/LLT.h | 4 +- Eigen/src/LU/FullPivLU.h | 4 +- Eigen/src/LU/PartialPivLU.h | 11 +-- Eigen/src/QR/ColPivHouseholderQR.h | 4 +- Eigen/src/QR/CompleteOrthogonalDecomposition.h | 4 +- Eigen/src/QR/FullPivHouseholderQR.h | 4 +- Eigen/src/QR/HouseholderQR.h | 4 +- doc/InplaceDecomposition.dox | 114 +++++++++++++++++++++++++ doc/Manual.dox | 2 + doc/examples/TutorialInplaceLU.cpp | 61 +++++++++++++ 11 files changed, 201 insertions(+), 15 deletions(-) create mode 100644 doc/InplaceDecomposition.dox create mode 100644 doc/examples/TutorialInplaceLU.cpp (limited to 'Eigen/src/LU/FullPivLU.h') diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index a31b3d6aa..69e176110 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -43,6 +43,8 @@ namespace internal { * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky * decomposition to determine whether a system of equations has a solution. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt(), class LLT */ template class LDLT @@ -112,7 +114,7 @@ template class LDLT /** \brief Constructs a LDLT factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa LDLT(const EigenBase&) */ diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index ad163c749..bd966656d 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -41,6 +41,8 @@ template struct LLT_Traits; * Example: \include LLT_example.cpp * Output: \verbinclude LLT_example.out * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT */ /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) @@ -96,7 +98,7 @@ template class LLT /** \brief Constructs a LDLT factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when * \c MatrixType is a Eigen::Ref. * * \sa LLT(const EigenBase&) diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 113b8c7b8..1632d3ac3 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -52,6 +52,8 @@ template struct traits > * \include class_FullPivLU.cpp * Output: \verbinclude class_FullPivLU.out * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse() */ template class FullPivLU @@ -99,7 +101,7 @@ template class FullPivLU /** \brief Constructs a LU factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa FullPivLU(const EigenBase&) */ diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index c862d9692..87ac6a281 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -68,6 +68,8 @@ struct enable_if_ref,Derived> { * * The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP(). * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU */ template class PartialPivLU @@ -113,17 +115,10 @@ template class PartialPivLU template explicit PartialPivLU(const EigenBase& matrix); - /** Constructor for inplace decomposition + /** Constructor for \link InplaceDecomposition inplace decomposition \endlink * * \param matrix the matrix of which to compute the LU decomposition. * - * If \c MatrixType is an Eigen::Ref, then the storage of \a matrix will be shared - * between \a matrix and \c *this and the decomposition will take place in-place. - * The memory of \a matrix will be used througrough the lifetime of \c *this. In - * particular, further calls to \c this->compute(A) will still operate on the memory - * of \a matrix meaning. This also implies that the sizes of \c A must match the - * ones of \a matrix. - * * \warning The matrix should have full rank (e.g. if it's square, it should be invertible). * If you need to deal with non-full rank, use class FullPivLU instead. */ diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index db50b5675..ee6989897 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -41,6 +41,8 @@ template struct traits > * This decomposition performs column pivoting in order to be rank-revealing and improve * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::colPivHouseholderQr() */ template class ColPivHouseholderQR @@ -135,7 +137,7 @@ template class ColPivHouseholderQR /** \brief Constructs a QR factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa ColPivHouseholderQR(const EigenBase&) */ diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h index 967eb35dd..69ffe4f42 100644 --- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h +++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -39,6 +39,8 @@ struct traits > * \b Q and \b Z are unitary matrices and \b T an upper triangular matrix of * size rank-by-rank. \b A may be rank deficient. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::completeOrthogonalDecomposition() */ template @@ -117,7 +119,7 @@ class CompleteOrthogonalDecomposition { /** \brief Constructs a complete orthogonal decomposition from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa CompleteOrthogonalDecomposition(const EigenBase&) */ diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 5c2f57d04..4c85a4693 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -50,6 +50,8 @@ struct traits > * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::fullPivHouseholderQr() */ template class FullPivHouseholderQR @@ -136,7 +138,7 @@ template class FullPivHouseholderQR /** \brief Constructs a QR factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when \c MatrixType is a Eigen::Ref. + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. * * \sa FullPivHouseholderQR(const EigenBase&) */ diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index f2a9cc080..2e64fa7fe 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -37,6 +37,8 @@ namespace Eigen { * This Householder QR decomposition is faster, but less numerically stable and less feature-full than * FullPivHouseholderQR or ColPivHouseholderQR. * + * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. + * * \sa MatrixBase::householderQr() */ template class HouseholderQR @@ -104,7 +106,7 @@ template class HouseholderQR /** \brief Constructs a QR factorization from a given matrix * - * This overloaded constructor is provided for inplace solving when + * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when * \c MatrixType is a Eigen::Ref. * * \sa HouseholderQR(const EigenBase&) diff --git a/doc/InplaceDecomposition.dox b/doc/InplaceDecomposition.dox new file mode 100644 index 000000000..fc89db4a1 --- /dev/null +++ b/doc/InplaceDecomposition.dox @@ -0,0 +1,114 @@ +namespace Eigen { + +/** \eigenManualPage InplaceDecomposition Inplace matrix decompositions + +Starting from %Eigen 3.3, the LU, Cholesky, and QR decompositions can operate \em inplace, that is, directly within the given input matrix. +This feature is especially useful when dealing with huae matrices, and or when the available memory is very limited (embedded systems). + +To this end, the respective decomposition class must be instantiated with a Ref<> matrix type, and the decomposition object must be constructed with the input matrix as argument. As an example, let us consider an inplace LU decomposition with partial pivoting. + +Let's start with the basic inclusions, and declaration of a 2x2 matrix \c A: + + + + + + + +
codeoutput
\snippet TutorialInplaceLU.cpp init + \snippet TutorialInplaceLU.out init +
+ +No surprise here! Then, let's declare our inplace LU object \c lu, and check the content of the matrix \c A: + + + + + + +
\snippet TutorialInplaceLU.cpp declaration + \snippet TutorialInplaceLU.out declaration +
+ +Here, the \c lu object computes and stores the \c L and \c U factors within the memory held by the matrix \c A. +The coefficients of \c A have thus been destroyed during the factorization, and replaced by the L and U factors as one can verify: + + + + + + +
\snippet TutorialInplaceLU.cpp matrixLU + \snippet TutorialInplaceLU.out matrixLU +
+ +Then, one can use the \c lu object as usual, for instance to solve the Ax=b problem: + + + + + +
\snippet TutorialInplaceLU.cpp solve + \snippet TutorialInplaceLU.out solve +
+ +Here, since the content of the original matrix \c A has been lost, we had to declared a new matrix \c A0 to verify the result. +Since the memory is shared between \c A and \c lu, modifying the matrix \c A will make \c lu invalid. +This can easily be verified by modifying the content of \c A and trying to solve the initial problem again: + + + + + + +
\snippet TutorialInplaceLU.cpp modifyA + \snippet TutorialInplaceLU.out modifyA +
+ +Note that there is no shared pointer under the hood, it is the \b responsibility \b of \b the \b user to keep the input matrix \c A in life as long as \c lu is living. + +If one wants to update the factorization with the modified A, one has to call the compute method as usual: + + + + + +
\snippet TutorialInplaceLU.cpp recompute + \snippet TutorialInplaceLU.out recompute +
+ +Note that calling compute does not change the memory which is referenced by the \c lu object. Therefore, if the compute method is called with another matrix \c A1 different than \c A, then the content of \c A1 won't be modified. This is still the content of \c A that will be used to store the L and U factors of the matrix \c A1. +This can easily be verified as follows: + + + + + +
\snippet TutorialInplaceLU.cpp recompute_bis0 + \snippet TutorialInplaceLU.out recompute_bis0 +
+The matrix \c A1 is unchanged, and one can thus solve A1*x=b, and directly check the residual without any copy of \c A1: + + + + + +
\snippet TutorialInplaceLU.cpp recompute_bis1 + \snippet TutorialInplaceLU.out recompute_bis1 +
+ + +Here is the list of matrix decomposition supported this inplace mechanism: + +- class LLT +- class LDLT +- class PartialPivLU +- class FullPivLU +- class HouseholderQR +- class ColPivHouseholderQR +- class FullPivHouseholderQR +- class CompleteOrthogonalDecomposition + +*/ + +} \ No newline at end of file diff --git a/doc/Manual.dox b/doc/Manual.dox index 8940b0070..5167cb9e5 100644 --- a/doc/Manual.dox +++ b/doc/Manual.dox @@ -106,6 +106,8 @@ namespace Eigen { \ingroup DenseLinearSolvers_chapter */ /** \addtogroup LeastSquares \ingroup DenseLinearSolvers_chapter */ +/** \addtogroup InplaceDecomposition + \ingroup DenseLinearSolvers_chapter */ /** \addtogroup DenseLinearSolvers_Reference \ingroup DenseLinearSolvers_chapter */ diff --git a/doc/examples/TutorialInplaceLU.cpp b/doc/examples/TutorialInplaceLU.cpp new file mode 100644 index 000000000..cb9c59b60 --- /dev/null +++ b/doc/examples/TutorialInplaceLU.cpp @@ -0,0 +1,61 @@ +#include +struct init { + init() { std::cout << "[" << "init" << "]" << std::endl; } +}; +init init_obj; +// [init] +#include +#include + +using namespace std; +using namespace Eigen; + +int main() +{ + MatrixXd A(2,2); + A << 2, -1, 1, 3; + cout << "Here is the input matrix A before decomposition:\n" << A << endl; +cout << "[init]" << endl; + +cout << "[declaration]" << endl; + PartialPivLU > lu(A); + cout << "Here is the input matrix A after decomposition:\n" << A << endl; +cout << "[declaration]" << endl; + +cout << "[matrixLU]" << endl; + cout << "Here is the matrix storing the L and U factors:\n" << lu.matrixLU() << endl; +cout << "[matrixLU]" << endl; + +cout << "[solve]" << endl; + MatrixXd A0(2,2); A0 << 2, -1, 1, 3; + VectorXd b(2); b << 1, 2; + VectorXd x = lu.solve(b); + cout << "Residual: " << (A0 * x - b).norm() << endl; +cout << "[solve]" << endl; + +cout << "[modifyA]" << endl; + A << 3, 4, -2, 1; + x = lu.solve(b); + cout << "Residual: " << (A0 * x - b).norm() << endl; +cout << "[modifyA]" << endl; + +cout << "[recompute]" << endl; + A0 = A; // save A + lu.compute(A); + x = lu.solve(b); + cout << "Residual: " << (A0 * x - b).norm() << endl; +cout << "[recompute]" << endl; + +cout << "[recompute_bis0]" << endl; + MatrixXd A1(2,2); + A1 << 5,-2,3,4; + lu.compute(A1); + cout << "Here is the input matrix A1 after decomposition:\n" << A1 << endl; +cout << "[recompute_bis0]" << endl; + +cout << "[recompute_bis1]" << endl; + x = lu.solve(b); + cout << "Residual: " << (A1 * x - b).norm() << endl; +cout << "[recompute_bis1]" << endl; + +} -- cgit v1.2.3 From 3cb914f33285f25de4b3e34ca33bbebc3e0a55d3 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 6 Sep 2016 09:55:50 +0200 Subject: bug #1266: remove CUDA guards on MatrixBase:: definitions. (those used to break old nvcc versions that we propably don't care anymore) --- Eigen/src/Cholesky/LDLT.h | 2 -- Eigen/src/Cholesky/LLT.h | 2 -- Eigen/src/LU/FullPivLU.h | 2 -- Eigen/src/LU/PartialPivLU.h | 4 ---- Eigen/src/QR/ColPivHouseholderQR.h | 2 -- Eigen/src/QR/CompleteOrthogonalDecomposition.h | 2 -- Eigen/src/QR/FullPivHouseholderQR.h | 2 -- Eigen/src/QR/HouseholderQR.h | 2 -- Eigen/src/SVD/JacobiSVD.h | 2 -- test/cuda_basic.cu | 3 ++- 10 files changed, 2 insertions(+), 21 deletions(-) (limited to 'Eigen/src/LU/FullPivLU.h') diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 795d19dce..fcee7b2e3 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -642,7 +642,6 @@ MatrixType LDLT::reconstructedMatrix() const return res; } -#ifndef __CUDACC__ /** \cholesky_module * \returns the Cholesky decomposition with full pivoting without square root of \c *this * \sa MatrixBase::ldlt() @@ -664,7 +663,6 @@ MatrixBase::ldlt() const { return LDLT(derived()); } -#endif // __CUDACC__ } // end namespace Eigen diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index bd966656d..ddf4875ab 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -507,7 +507,6 @@ MatrixType LLT::reconstructedMatrix() const return matrixL() * matrixL().adjoint().toDenseMatrix(); } -#ifndef __CUDACC__ /** \cholesky_module * \returns the LLT decomposition of \c *this * \sa SelfAdjointView::llt() @@ -529,7 +528,6 @@ SelfAdjointView::llt() const { return LLT(m_matrix); } -#endif // __CUDACC__ } // end namespace Eigen diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 1632d3ac3..fbb5a61bf 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -879,14 +879,12 @@ struct Assignment >, internal::assign_ * * \sa class FullPivLU */ -#ifndef __CUDACC__ template inline const FullPivLU::PlainObject> MatrixBase::fullPivLu() const { return FullPivLU(eval()); } -#endif } // end namespace Eigen diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 87ac6a281..6c35e88d3 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -584,14 +584,12 @@ struct Assignment >, internal::assi * * \sa class PartialPivLU */ -#ifndef __CUDACC__ template inline const PartialPivLU::PlainObject> MatrixBase::partialPivLu() const { return PartialPivLU(eval()); } -#endif /** \lu_module * @@ -601,14 +599,12 @@ MatrixBase::partialPivLu() const * * \sa class PartialPivLU */ -#ifndef __CUDACC__ template inline const PartialPivLU::PlainObject> MatrixBase::lu() const { return PartialPivLU(eval()); } -#endif } // end namespace Eigen diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 35e19b85b..9650781d6 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -637,7 +637,6 @@ typename ColPivHouseholderQR::HouseholderSequenceType ColPivHousehol return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); } -#ifndef __CUDACC__ /** \return the column-pivoting Householder QR decomposition of \c *this. * * \sa class ColPivHouseholderQR @@ -648,7 +647,6 @@ MatrixBase::colPivHouseholderQr() const { return ColPivHouseholderQR(eval()); } -#endif // __CUDACC__ } // end namespace Eigen diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h index f299d3c00..41e4ecdfd 100644 --- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h +++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h @@ -547,7 +547,6 @@ CompleteOrthogonalDecomposition::householderQ() const { return m_cpqr.householderQ(); } -#ifndef __CUDACC__ /** \return the complete orthogonal decomposition of \c *this. * * \sa class CompleteOrthogonalDecomposition @@ -557,7 +556,6 @@ const CompleteOrthogonalDecomposition::PlainObject> MatrixBase::completeOrthogonalDecomposition() const { return CompleteOrthogonalDecomposition(eval()); } -#endif // __CUDACC__ } // end namespace Eigen diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index a46d3f9f3..e0e15100d 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -660,7 +660,6 @@ inline typename FullPivHouseholderQR::MatrixQReturnType FullPivHouse return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions); } -#ifndef __CUDACC__ /** \return the full-pivoting Householder QR decomposition of \c *this. * * \sa class FullPivHouseholderQR @@ -671,7 +670,6 @@ MatrixBase::fullPivHouseholderQr() const { return FullPivHouseholderQR(eval()); } -#endif // __CUDACC__ } // end namespace Eigen diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 5e69cfee9..3513d995c 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -393,7 +393,6 @@ void HouseholderQR::computeInPlace() m_isInitialized = true; } -#ifndef __CUDACC__ /** \return the Householder QR decomposition of \c *this. * * \sa class HouseholderQR @@ -404,7 +403,6 @@ MatrixBase::householderQr() const { return HouseholderQR(eval()); } -#endif // __CUDACC__ } // end namespace Eigen diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 605c1a2a6..78dfd1d59 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -783,7 +783,6 @@ JacobiSVD::compute(const MatrixType& matrix, unsig return *this; } -#ifndef __CUDACC__ /** \svd_module * * \return the singular value decomposition of \c *this computed by two-sided @@ -797,7 +796,6 @@ MatrixBase::jacobiSvd(unsigned int computationOptions) const { return JacobiSVD(*this, computationOptions); } -#endif // __CUDACC__ } // end namespace Eigen diff --git a/test/cuda_basic.cu b/test/cuda_basic.cu index 0a515ce7b..cb2e4167a 100644 --- a/test/cuda_basic.cu +++ b/test/cuda_basic.cu @@ -26,7 +26,8 @@ #include "main.h" #include "cuda_common.h" -#include +// Check that dense modules can be properly parsed by nvcc +#include // struct Foo{ // EIGEN_DEVICE_FUNC -- cgit v1.2.3 From 73c8f2f697e4722bc8626a7412132fb13e75790e Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Tue, 13 Sep 2016 07:58:39 +0200 Subject: bug #1285: fix regression introduced in changeset 00c29c2caef8fb0c6b1d2ba5ecdf6780c0c766d4 --- Eigen/src/LU/FullPivLU.h | 2 +- Eigen/src/LU/PartialPivLU.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) (limited to 'Eigen/src/LU/FullPivLU.h') diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index fbb5a61bf..2b30fc146 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -436,7 +436,7 @@ template class FullPivLU Index m_nonzero_pivots; RealScalar m_l1_norm; RealScalar m_maxpivot, m_prescribedThreshold; - char m_det_pq; + signed char m_det_pq; bool m_isInitialized, m_usePrescribedThreshold; }; diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 6c35e88d3..d43961887 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -284,7 +284,7 @@ template class PartialPivLU PermutationType m_p; TranspositionType m_rowsTranspositions; RealScalar m_l1_norm; - char m_det_p; + signed char m_det_p; bool m_isInitialized; }; -- cgit v1.2.3 From c0d56a543e170d221e58b177b245c0e57ce1993a Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Wed, 14 Sep 2016 14:06:21 -0700 Subject: Added several missing EIGEN_DEVICE_FUNC qualifiers --- Eigen/src/Core/Inverse.h | 2 +- Eigen/src/Core/MathFunctions.h | 3 ++- Eigen/src/Core/MatrixBase.h | 2 +- Eigen/src/LU/FullPivLU.h | 8 ++++---- Eigen/src/LU/InverseImpl.h | 2 +- Eigen/src/LU/PartialPivLU.h | 4 ++-- 6 files changed, 11 insertions(+), 10 deletions(-) (limited to 'Eigen/src/LU/FullPivLU.h') diff --git a/Eigen/src/Core/Inverse.h b/Eigen/src/Core/Inverse.h index f3ec84990..f303aebf9 100644 --- a/Eigen/src/Core/Inverse.h +++ b/Eigen/src/Core/Inverse.h @@ -50,7 +50,7 @@ public: typedef typename internal::ref_selector::type Nested; typedef typename internal::remove_all::type NestedExpression; - explicit Inverse(const XprType &xpr) + explicit EIGEN_DEVICE_FUNC Inverse(const XprType &xpr) : m_xpr(xpr) {} diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index bf3044b96..fa322aca7 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -1312,11 +1312,12 @@ template struct scalar_fuzzy_default_impl { typedef typename NumTraits::Real RealScalar; - template + template EIGEN_DEVICE_FUNC static inline bool isMuchSmallerThan(const Scalar& x, const OtherScalar& y, const RealScalar& prec) { return numext::abs2(x) <= numext::abs2(y) * prec * prec; } + EIGEN_DEVICE_FUNC static inline bool isApprox(const Scalar& x, const Scalar& y, const RealScalar& prec) { return numext::abs2(x - y) <= numext::mini(numext::abs2(x), numext::abs2(y)) * prec * prec; diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 334a4d71e..976b130e3 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -98,7 +98,7 @@ template class MatrixBase /** \returns the size of the main diagonal, which is min(rows(),cols()). * \sa rows(), cols(), SizeAtCompileTime. */ EIGEN_DEVICE_FUNC - inline Index diagonalSize() const { return (std::min)(rows(),cols()); } + inline Index diagonalSize() const { return (numext::mini)(rows(),cols()); } typedef typename Base::PlainObject PlainObject; diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 2b30fc146..ebcd5c208 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -156,7 +156,7 @@ template class FullPivLU * * \sa permutationQ() */ - inline const PermutationPType& permutationP() const + EIGEN_DEVICE_FUNC inline const PermutationPType& permutationP() const { eigen_assert(m_isInitialized && "LU is not initialized."); return m_p; @@ -406,8 +406,8 @@ template class FullPivLU MatrixType reconstructedMatrix() const; - inline Index rows() const { return m_lu.rows(); } - inline Index cols() const { return m_lu.cols(); } + EIGEN_DEVICE_FUNC inline Index rows() const { return m_lu.rows(); } + EIGEN_DEVICE_FUNC inline Index cols() const { return m_lu.cols(); } #ifndef EIGEN_PARSED_BY_DOXYGEN template @@ -879,7 +879,7 @@ struct Assignment >, internal::assign_ * * \sa class FullPivLU */ -template +template EIGEN_DEVICE_FUNC inline const FullPivLU::PlainObject> MatrixBase::fullPivLu() const { diff --git a/Eigen/src/LU/InverseImpl.h b/Eigen/src/LU/InverseImpl.h index 3134632e1..147f9496c 100644 --- a/Eigen/src/LU/InverseImpl.h +++ b/Eigen/src/LU/InverseImpl.h @@ -327,7 +327,7 @@ struct Assignment, internal::assign_op +template EIGEN_DEVICE_FUNC inline const Inverse MatrixBase::inverse() const { EIGEN_STATIC_ASSERT(!NumTraits::IsInteger,THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES) diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index d43961887..13394fffa 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -584,7 +584,7 @@ struct Assignment >, internal::assi * * \sa class PartialPivLU */ -template +template EIGEN_DEVICE_FUNC inline const PartialPivLU::PlainObject> MatrixBase::partialPivLu() const { @@ -599,7 +599,7 @@ MatrixBase::partialPivLu() const * * \sa class PartialPivLU */ -template +template EIGEN_DEVICE_FUNC inline const PartialPivLU::PlainObject> MatrixBase::lu() const { -- cgit v1.2.3 From 59e9edfbf1e1a58d1e2401ec2f05b6fdd19fd87c Mon Sep 17 00:00:00 2001 From: Benoit Steiner Date: Mon, 19 Sep 2016 14:13:20 -0700 Subject: Removed EIGEN_DEVICE_FUNC qualifers for the lu(), fullPivLu(), partialPivLu(), and inverse() functions since they aren't ready to run on GPU --- Eigen/src/Core/MatrixBase.h | 4 ---- Eigen/src/LU/FullPivLU.h | 2 +- Eigen/src/LU/InverseImpl.h | 2 +- Eigen/src/LU/PartialPivLU.h | 4 ++-- 4 files changed, 4 insertions(+), 8 deletions(-) (limited to 'Eigen/src/LU/FullPivLU.h') diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 19a7483ba..d56df8249 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -330,15 +330,11 @@ template class MatrixBase /////////// LU module /////////// - EIGEN_DEVICE_FUNC inline const FullPivLU fullPivLu() const; - EIGEN_DEVICE_FUNC inline const PartialPivLU partialPivLu() const; - EIGEN_DEVICE_FUNC inline const PartialPivLU lu() const; - EIGEN_DEVICE_FUNC inline const Inverse inverse() const; template diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index ebcd5c208..03b6af706 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -879,7 +879,7 @@ struct Assignment >, internal::assign_ * * \sa class FullPivLU */ -template EIGEN_DEVICE_FUNC +template inline const FullPivLU::PlainObject> MatrixBase::fullPivLu() const { diff --git a/Eigen/src/LU/InverseImpl.h b/Eigen/src/LU/InverseImpl.h index 147f9496c..3134632e1 100644 --- a/Eigen/src/LU/InverseImpl.h +++ b/Eigen/src/LU/InverseImpl.h @@ -327,7 +327,7 @@ struct Assignment, internal::assign_op EIGEN_DEVICE_FUNC +template inline const Inverse MatrixBase::inverse() const { EIGEN_STATIC_ASSERT(!NumTraits::IsInteger,THIS_FUNCTION_IS_NOT_FOR_INTEGER_NUMERIC_TYPES) diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 13394fffa..d43961887 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -584,7 +584,7 @@ struct Assignment >, internal::assi * * \sa class PartialPivLU */ -template EIGEN_DEVICE_FUNC +template inline const PartialPivLU::PlainObject> MatrixBase::partialPivLu() const { @@ -599,7 +599,7 @@ MatrixBase::partialPivLu() const * * \sa class PartialPivLU */ -template EIGEN_DEVICE_FUNC +template inline const PartialPivLU::PlainObject> MatrixBase::lu() const { -- cgit v1.2.3