diff options
| -rw-r--r-- | Eigen/src/Core/AssignEvaluator.h | 13 | ||||
| -rw-r--r-- | Eigen/src/Core/SelfAdjointView.h | 5 | ||||
| -rw-r--r-- | Eigen/src/Core/Swap.h | 1 | ||||
| -rw-r--r-- | Eigen/src/Core/TriangularMatrix.h | 301 |
4 files changed, 167 insertions, 153 deletions
diff --git a/Eigen/src/Core/AssignEvaluator.h b/Eigen/src/Core/AssignEvaluator.h index 6d3739407..86cd1cd68 100644 --- a/Eigen/src/Core/AssignEvaluator.h +++ b/Eigen/src/Core/AssignEvaluator.h @@ -140,6 +140,7 @@ public: template<typename Kernel, int Index, int Stop> struct copy_using_evaluator_DefaultTraversal_CompleteUnrolling { + // FIXME: this is not very clean, perhaps this information should be provided by the kernel? typedef typename Kernel::DstEvaluatorType DstEvaluatorType; typedef typename DstEvaluatorType::XprType DstXprType; @@ -204,9 +205,10 @@ struct copy_using_evaluator_LinearTraversal_CompleteUnrolling<Kernel, Stop, Stop template<typename Kernel, int Index, int Stop> struct copy_using_evaluator_innervec_CompleteUnrolling { + // FIXME: this is not very clean, perhaps this information should be provided by the kernel? typedef typename Kernel::DstEvaluatorType DstEvaluatorType; typedef typename DstEvaluatorType::XprType DstXprType; - + enum { outer = Index / DstXprType::InnerSizeAtCompileTime, inner = Index % DstXprType::InnerSizeAtCompileTime, @@ -233,8 +235,7 @@ struct copy_using_evaluator_innervec_InnerUnrolling static EIGEN_STRONG_INLINE void run(Kernel &kernel, int outer) { kernel.template assignPacketByOuterInner<Aligned, Aligned>(outer, Index); - typedef typename Kernel::DstEvaluatorType::XprType DstXprType; - enum { NextIndex = Index + packet_traits<typename DstXprType::Scalar>::size }; + enum { NextIndex = Index + packet_traits<typename Kernel::Scalar>::size }; copy_using_evaluator_innervec_InnerUnrolling<Kernel, NextIndex, Stop>::run(kernel, outer); } }; @@ -542,12 +543,6 @@ public: m_functor.assignCoeff(m_dst.coeffRef(row,col), m_src.coeff(row,col)); } - /// This overload by-passes the source expression, i.e., dst(row,col) ?= value - void assignCoeff(Index row, Index col, const Scalar &value) - { - m_functor.assignCoeff(m_dst.coeffRef(row,col), value); - } - /// \sa assignCoeff(Index,Index) void assignCoeff(Index index) { diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index a5c6d5ceb..f34de8b39 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -244,6 +244,8 @@ template<typename _MatrixType, unsigned int UpLo> class SelfAdjointView namespace internal { +#ifndef EIGEN_TEST_EVALUATORS + template<typename Derived1, typename Derived2, int UnrollCount, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount, ClearOpposite> { @@ -336,6 +338,8 @@ struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dyn } }; +#endif // EIGEN_TEST_EVALUATORS + #ifdef EIGEN_ENABLE_EVALUATORS // TODO currently a selfadjoint expression has the form SelfAdjointView<.,.> // in the future selfadjoint-ness should be defined by the expression traits @@ -348,6 +352,7 @@ struct evaluator_traits<SelfAdjointView<MatrixType,Mode> > static const int AssumeAliasing = 0; }; + #endif // EIGEN_ENABLE_EVALUATORS } // end namespace internal diff --git a/Eigen/src/Core/Swap.h b/Eigen/src/Core/Swap.h index 117f667f6..e7d525572 100644 --- a/Eigen/src/Core/Swap.h +++ b/Eigen/src/Core/Swap.h @@ -148,6 +148,7 @@ template<typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT> class generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, swap_assign_op<typename DstEvaluatorTypeT::Scalar>, Specialized> : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, swap_assign_op<typename DstEvaluatorTypeT::Scalar>, BuiltIn> { +protected: typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, swap_assign_op<typename DstEvaluatorTypeT::Scalar>, BuiltIn> Base; typedef typename DstEvaluatorTypeT::PacketScalar PacketScalar; using Base::m_dst; diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 28191694d..c658e2290 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -298,8 +298,8 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView template<typename OtherDerived> EIGEN_DEVICE_FUNC - void lazyAssign(const MatrixBase<OtherDerived>& other); - + void lazyAssign(const MatrixBase<OtherDerived>& other); + /** \sa MatrixBase::conjugate() */ EIGEN_DEVICE_FUNC inline TriangularView<MatrixConjugateReturnType,Mode> conjugate() @@ -469,6 +469,8 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView return m_matrix.diagonal().prod(); } +#ifndef EIGEN_TEST_EVALUATORS + // TODO simplify the following: template<typename ProductDerived, typename Lhs, typename Rhs> EIGEN_DEVICE_FUNC @@ -514,7 +516,9 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView { return assignProduct(other,-other.alpha()); } - + +#endif // EIGEN_TEST_EVALUATORS + protected: template<typename ProductDerived, typename Lhs, typename Rhs> @@ -530,6 +534,8 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView namespace internal { +#ifndef EIGEN_TEST_EVALUATORS + template<typename Derived1, typename Derived2, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_selector { @@ -694,8 +700,52 @@ struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, Cl } }; +#endif // EIGEN_TEST_EVALUATORS + } // end namespace internal +#ifdef EIGEN_TEST_EVALUATORS + +// FIXME should we keep that possibility +template<typename MatrixType, unsigned int Mode> +template<typename OtherDerived> +inline TriangularView<MatrixType, Mode>& +TriangularView<MatrixType, Mode>::operator=(const MatrixBase<OtherDerived>& other) +{ + internal::call_assignment(*this, other.template triangularView<Mode>(), internal::assign_op<Scalar>()); + return *this; +} + +// FIXME should we keep that possibility +template<typename MatrixType, unsigned int Mode> +template<typename OtherDerived> +void TriangularView<MatrixType, Mode>::lazyAssign(const MatrixBase<OtherDerived>& other) +{ + internal::call_assignment(this->noalias(), other.template triangularView<Mode>()); +} + + + +template<typename MatrixType, unsigned int Mode> +template<typename OtherDerived> +inline TriangularView<MatrixType, Mode>& +TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& other) +{ + eigen_assert(Mode == int(OtherDerived::Mode)); + internal::call_assignment(*this, other.derived()); + return *this; +} + +template<typename MatrixType, unsigned int Mode> +template<typename OtherDerived> +void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDerived>& other) +{ + eigen_assert(Mode == int(OtherDerived::Mode)); + internal::call_assignment(this->noalias(), other.derived()); +} + +#else + // FIXME should we keep that possibility template<typename MatrixType, unsigned int Mode> template<typename OtherDerived> @@ -770,6 +820,8 @@ void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDeri >::run(m_matrix.const_cast_derived(), other.derived().nestedExpression()); } +#endif // EIGEN_TEST_EVALUATORS + /*************************************************************************** * Implementation of TriangularBase methods ***************************************************************************/ @@ -790,6 +842,8 @@ void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const evalToLazy(other.derived()); } +#ifndef EIGEN_TEST_EVALUATORS + /** Assigns a triangular or selfadjoint matrix to a dense matrix. * If the matrix is triangular, the opposite part is set to zero. */ template<typename Derived> @@ -811,6 +865,8 @@ void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const >::run(other.derived(), derived().nestedExpression()); } +#endif // EIGEN_TEST_EVALUATORS + /*************************************************************************** * Implementation of TriangularView methods ***************************************************************************/ @@ -962,15 +1018,15 @@ struct evaluator_traits<TriangularView<MatrixType,Mode> > template<typename MatrixType, unsigned int Mode, typename Kind> struct evaluator<TriangularView<MatrixType,Mode>, Kind, typename MatrixType::Scalar> - : evaluator<MatrixType> + : evaluator<typename internal::remove_all<MatrixType>::type> { typedef TriangularView<MatrixType,Mode> XprType; - typedef evaluator<MatrixType> Base; + typedef evaluator<typename internal::remove_all<MatrixType>::type> Base; typedef evaluator type; evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {} }; -// Additional assignement kinds: +// Additional assignment kinds: struct Triangular2Triangular {}; struct Triangular2Dense {}; struct Dense2Triangular {}; @@ -978,6 +1034,59 @@ struct Dense2Triangular {}; template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop; + +// Specialization of the dense assignment kernel for triangular matrices. +// The main additions are: +// - An overload of assignCoeff(i, j, val) that by-passes the source expression. Needed to set the opposite triangular part to 0. +// - When calling assignCoeff(i,j...), we guarantee that i!=j +// - New assignDiagonalCoeff(i), and assignDiagonalCoeff(i, val) for assigning coefficients on the diagonal. +template<int Mode, int ClearOpposite, typename DstEvaluatorTypeT, typename SrcEvaluatorTypeT, typename Functor, int Version = Specialized> +class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> +{ +protected: + typedef generic_dense_assignment_kernel<DstEvaluatorTypeT, SrcEvaluatorTypeT, Functor, Version> Base; + typedef typename Base::DstXprType DstXprType; + typedef typename Base::SrcXprType SrcXprType; + using Base::m_dst; + using Base::m_src; + using Base::m_functor; +public: + + typedef typename Base::DstEvaluatorType DstEvaluatorType; + typedef typename Base::SrcEvaluatorType SrcEvaluatorType; + typedef typename Base::Scalar Scalar; + typedef typename Base::Index Index; + typedef typename Base::AssignmentTraits AssignmentTraits; + + + triangular_dense_assignment_kernel(DstEvaluatorType &dst, const SrcEvaluatorType &src, const Functor &func, DstXprType& dstExpr) + : Base(dst, src, func, dstExpr) + {} + +#ifdef EIGEN_INTERNAL_DEBUGGING + void assignCoeff(Index row, Index col) + { + eigen_internal_assert(row!=col); + Base::assignCoeff(row,col); + } +#else + using Base::assignCoeff; +#endif + + void assignDiagonalCoeff(Index id) + { + if((Mode&UnitDiag) && ClearOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(1)); + else if((Mode&ZeroDiag) && ClearOpposite) m_functor.assignCoeff(m_dst.coeffRef(id,id), Scalar(0)); + else if((Mode&(UnitDiag|ZeroDiag))==0) Base::assignCoeff(id,id); + } + + void assignOppositeCoeff(Index row, Index col) + { + eigen_internal_assert(row!=col); + if(ClearOpposite) + m_functor.assignCoeff(m_dst.coeffRef(row,col), Scalar(0)); + } +}; template<int Mode, bool ClearOpposite, typename DstXprType, typename SrcXprType, typename Functor> void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func) @@ -994,13 +1103,13 @@ void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& sr DstEvaluatorType dstEvaluator(dst); SrcEvaluatorType srcEvaluator(src); - typedef generic_dense_assignment_kernel<DstEvaluatorType,SrcEvaluatorType,Functor> Kernel; + typedef triangular_dense_assignment_kernel<Mode,ClearOpposite,DstEvaluatorType,SrcEvaluatorType,Functor> Kernel; Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived()); enum { unroll = DstXprType::SizeAtCompileTime != Dynamic - && internal::traits<SrcXprType>::CoeffReadCost != Dynamic - && DstXprType::SizeAtCompileTime * internal::traits<SrcXprType>::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT +// && internal::traits<SrcXprType>::CoeffReadCost != Dynamic +// && DstXprType::SizeAtCompileTime * internal::traits<SrcXprType>::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT }; triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, ClearOpposite>::run(kernel); @@ -1050,9 +1159,13 @@ struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular, Scalar> template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop { + // FIXME: this is not very clean, perhaps this information should be provided by the kernel? + typedef typename Kernel::DstEvaluatorType DstEvaluatorType; + typedef typename DstEvaluatorType::XprType DstXprType; + enum { - col = (UnrollCount-1) / Kernel::DstEvaluatorType::RowsAtCompileTime, - row = (UnrollCount-1) % Kernel::DstEvaluatorType::RowsAtCompileTime + col = (UnrollCount-1) / DstXprType::RowsAtCompileTime, + row = (UnrollCount-1) % DstXprType::RowsAtCompileTime }; typedef typename Kernel::Scalar Scalar; @@ -1061,24 +1174,13 @@ struct triangular_assignment_loop static inline void run(Kernel &kernel) { triangular_assignment_loop<Kernel, Mode, UnrollCount-1, ClearOpposite>::run(kernel); - - // TODO should be a static assert - eigen_assert( Mode == Upper || Mode == Lower - || Mode == StrictlyUpper || Mode == StrictlyLower - || Mode == UnitUpper || Mode == UnitLower); - if((Mode == Upper && row <= col) - || (Mode == Lower && row >= col) - || (Mode == StrictlyUpper && row < col) - || (Mode == StrictlyLower && row > col) - || (Mode == UnitUpper && row < col) - || (Mode == UnitLower && row > col)) - kernel.assignCoeff(row, col); - else if(ClearOpposite) - { - if((Mode&UnitDiag) && row==col) kernel.assignCoeff(row, col, Scalar(1)); - else kernel.assignCoeff(row, col, Scalar(0)); - } + if(row==col) + kernel.assignDiagonalCoeff(row); + else if( ((Mode&Lower) && row>col) || ((Mode&Upper) && row<col) ) + kernel.assignCoeff(row,col); + else + kernel.assignOppositeCoeff(row,col); } }; @@ -1090,120 +1192,14 @@ struct triangular_assignment_loop<Kernel, Mode, 0, ClearOpposite> static inline void run(Kernel &) {} }; -// TODO: merge the following 6 variants into a single one (or max two), -// and perhaps write them using sub-expressions -// TODO: expreiment with a recursive assignment procedure splitting the current -// triangular part into one rectangular and two triangular parts. -template<typename Kernel, bool ClearOpposite> -struct triangular_assignment_loop<Kernel, Upper, Dynamic, ClearOpposite> -{ - typedef typename Kernel::Index Index; - typedef typename Kernel::Scalar Scalar; - EIGEN_DEVICE_FUNC - static inline void run(Kernel &kernel) - { - for(Index j = 0; j < kernel.cols(); ++j) - { - Index maxi = (std::min)(j, kernel.rows()-1); - for(Index i = 0; i <= maxi; ++i) - kernel.assignCoeff(i, j); - if (ClearOpposite) - for(Index i = maxi+1; i < kernel.rows(); ++i) - kernel.assignCoeff(i, j, Scalar(0)); - } - } -}; - -template<typename Kernel, bool ClearOpposite> -struct triangular_assignment_loop<Kernel, Lower, Dynamic, ClearOpposite> -{ - typedef typename Kernel::Index Index; - typedef typename Kernel::Scalar Scalar; - EIGEN_DEVICE_FUNC - static inline void run(Kernel &kernel) - { - for(Index j = 0; j < kernel.cols(); ++j) - { - for(Index i = j; i < kernel.rows(); ++i) - kernel.assignCoeff(i, j); - Index maxi = (std::min)(j, kernel.rows()); - if (ClearOpposite) - for(Index i = 0; i < maxi; ++i) - kernel.assignCoeff(i, j, Scalar(0)); - } - } -}; - -template<typename Kernel, bool ClearOpposite> -struct triangular_assignment_loop<Kernel, StrictlyUpper, Dynamic, ClearOpposite> -{ - typedef typename Kernel::Index Index; - typedef typename Kernel::Scalar Scalar; - EIGEN_DEVICE_FUNC - static inline void run(Kernel &kernel) - { - for(Index j = 0; j < kernel.cols(); ++j) - { - Index maxi = (std::min)(j, kernel.rows()); - for(Index i = 0; i < maxi; ++i) - kernel.assignCoeff(i, j); - if (ClearOpposite) - for(Index i = maxi; i < kernel.rows(); ++i) - kernel.assignCoeff(i, j, Scalar(0)); - } - } -}; +// TODO: experiment with a recursive assignment procedure splitting the current +// triangular part into one rectangular and two triangular parts. -template<typename Kernel, bool ClearOpposite> -struct triangular_assignment_loop<Kernel, StrictlyLower, Dynamic, ClearOpposite> -{ - typedef typename Kernel::Index Index; - typedef typename Kernel::Scalar Scalar; - EIGEN_DEVICE_FUNC - static inline void run(Kernel &kernel) - { - for(Index j = 0; j < kernel.cols(); ++j) - { - for(Index i = j+1; i < kernel.rows(); ++i) - kernel.assignCoeff(i, j); - Index maxi = (std::min)(j, kernel.rows()-1); - if (ClearOpposite) - for(Index i = 0; i <= maxi; ++i) - kernel.assignCoeff(i, j, Scalar(0)); - } - } -}; -template<typename Kernel, bool ClearOpposite> -struct triangular_assignment_loop<Kernel, UnitUpper, Dynamic, ClearOpposite> -{ - typedef typename Kernel::Index Index; - typedef typename Kernel::Scalar Scalar; - EIGEN_DEVICE_FUNC - static inline void run(Kernel &kernel) - { - for(Index j = 0; j < kernel.cols(); ++j) - { - Index maxi = (std::min)(j, kernel.rows()); - Index i = 0; - for(; i < maxi; ++i) - kernel.assignCoeff(i, j); - - if(i<kernel.rows()) // then i==j - kernel.assignCoeff(i++, j, Scalar(1)); - - if (ClearOpposite) - { - for(; i < kernel.rows(); ++i) - kernel.assignCoeff(i, j, Scalar(0)); - } - } - } -}; -template<typename Kernel, bool ClearOpposite> -struct triangular_assignment_loop<Kernel, UnitLower, Dynamic, ClearOpposite> +template<typename Kernel, unsigned int Mode, bool ClearOpposite> +struct triangular_assignment_loop<Kernel, Mode, Dynamic, ClearOpposite> { typedef typename Kernel::Index Index; typedef typename Kernel::Scalar Scalar; @@ -1214,24 +1210,41 @@ struct triangular_assignment_loop<Kernel, UnitLower, Dynamic, ClearOpposite> { Index maxi = (std::min)(j, kernel.rows()); Index i = 0; - if (ClearOpposite) + if (((Mode&Lower) && ClearOpposite) || (Mode&Upper)) { for(; i < maxi; ++i) - kernel.assignCoeff(i, j, Scalar(0)); + if(Mode&Upper) kernel.assignCoeff(i, j); + else kernel.assignOppositeCoeff(i, j); } + else + i = maxi; if(i<kernel.rows()) // then i==j - kernel.assignCoeff(i++, j, Scalar(1)); + kernel.assignDiagonalCoeff(i++); - for(; i < kernel.rows(); ++i) - kernel.assignCoeff(i, j); + if (((Mode&Upper) && ClearOpposite) || (Mode&Lower)) + { + for(; i < kernel.rows(); ++i) + if(Mode&Lower) kernel.assignCoeff(i, j); + else kernel.assignOppositeCoeff(i, j); + } } } }; } // end namespace internal -#endif +/** Assigns a triangular or selfadjoint matrix to a dense matrix. + * If the matrix is triangular, the opposite part is set to zero. */ +template<typename Derived> +template<typename DenseDerived> +void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const +{ + other.derived().resize(this->rows(), this->cols()); + internal::call_triangular_assignment_loop<Derived::Mode,true /* ClearOpposite */>(other.derived(), derived().nestedExpression()); +} + +#endif // EIGEN_ENABLE_EVALUATORS } // end namespace Eigen |