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 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 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(outer, Index); - typedef typename Kernel::DstEvaluatorType::XprType DstXprType; - enum { NextIndex = Index + packet_traits::size }; + enum { NextIndex = Index + packet_traits::size }; copy_using_evaluator_innervec_InnerUnrolling::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 class SelfAdjointView namespace internal { +#ifndef EIGEN_TEST_EVALUATORS + template struct triangular_assignment_selector { @@ -336,6 +338,8 @@ struct triangular_assignment_selector // in the future selfadjoint-ness should be defined by the expression traits @@ -348,6 +352,7 @@ struct evaluator_traits > 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 class generic_dense_assignment_kernel, Specialized> : public generic_dense_assignment_kernel, BuiltIn> { +protected: typedef generic_dense_assignment_kernel, 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 class TriangularView template EIGEN_DEVICE_FUNC - void lazyAssign(const MatrixBase& other); - + void lazyAssign(const MatrixBase& other); + /** \sa MatrixBase::conjugate() */ EIGEN_DEVICE_FUNC inline TriangularView conjugate() @@ -469,6 +469,8 @@ template class TriangularView return m_matrix.diagonal().prod(); } +#ifndef EIGEN_TEST_EVALUATORS + // TODO simplify the following: template EIGEN_DEVICE_FUNC @@ -514,7 +516,9 @@ template class TriangularView { return assignProduct(other,-other.alpha()); } - + +#endif // EIGEN_TEST_EVALUATORS + protected: template @@ -530,6 +534,8 @@ template class TriangularView namespace internal { +#ifndef EIGEN_TEST_EVALUATORS + template struct triangular_assignment_selector { @@ -694,8 +700,52 @@ struct triangular_assignment_selector +template +inline TriangularView& +TriangularView::operator=(const MatrixBase& other) +{ + internal::call_assignment(*this, other.template triangularView(), internal::assign_op()); + return *this; +} + +// FIXME should we keep that possibility +template +template +void TriangularView::lazyAssign(const MatrixBase& other) +{ + internal::call_assignment(this->noalias(), other.template triangularView()); +} + + + +template +template +inline TriangularView& +TriangularView::operator=(const TriangularBase& other) +{ + eigen_assert(Mode == int(OtherDerived::Mode)); + internal::call_assignment(*this, other.derived()); + return *this; +} + +template +template +void TriangularView::lazyAssign(const TriangularBase& other) +{ + eigen_assert(Mode == int(OtherDerived::Mode)); + internal::call_assignment(this->noalias(), other.derived()); +} + +#else + // FIXME should we keep that possibility template template @@ -770,6 +820,8 @@ void TriangularView::lazyAssign(const TriangularBase::run(m_matrix.const_cast_derived(), other.derived().nestedExpression()); } +#endif // EIGEN_TEST_EVALUATORS + /*************************************************************************** * Implementation of TriangularBase methods ***************************************************************************/ @@ -790,6 +842,8 @@ void TriangularBase::evalTo(MatrixBase &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 @@ -811,6 +865,8 @@ void TriangularBase::evalToLazy(MatrixBase &other) const >::run(other.derived(), derived().nestedExpression()); } +#endif // EIGEN_TEST_EVALUATORS + /*************************************************************************** * Implementation of TriangularView methods ***************************************************************************/ @@ -962,15 +1018,15 @@ struct evaluator_traits > template struct evaluator, Kind, typename MatrixType::Scalar> - : evaluator + : evaluator::type> { typedef TriangularView XprType; - typedef evaluator Base; + typedef evaluator::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 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 +class triangular_dense_assignment_kernel : public generic_dense_assignment_kernel +{ +protected: + typedef generic_dense_assignment_kernel 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 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 Kernel; + typedef triangular_dense_assignment_kernel Kernel; Kernel kernel(dstEvaluator, srcEvaluator, func, dst.const_cast_derived()); enum { unroll = DstXprType::SizeAtCompileTime != Dynamic - && internal::traits::CoeffReadCost != Dynamic - && DstXprType::SizeAtCompileTime * internal::traits::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT +// && internal::traits::CoeffReadCost != Dynamic +// && DstXprType::SizeAtCompileTime * internal::traits::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT }; triangular_assignment_loop::run(kernel); @@ -1050,9 +1159,13 @@ struct Assignment template 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::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 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 + +// TODO: experiment with a recursive assignment procedure splitting the current // triangular part into one rectangular and two triangular parts. -template -struct triangular_assignment_loop -{ - 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 -struct triangular_assignment_loop -{ - 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 -struct triangular_assignment_loop -{ - 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)); - } - } -}; - -template -struct triangular_assignment_loop -{ - 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 -struct triangular_assignment_loop +template +struct triangular_assignment_loop { typedef typename Kernel::Index Index; typedef typename Kernel::Scalar Scalar; @@ -1188,50 +1210,41 @@ struct triangular_assignment_loop { Index maxi = (std::min)(j, kernel.rows()); Index i = 0; - for(; i < maxi; ++i) - kernel.assignCoeff(i, j); - - if(i -struct triangular_assignment_loop -{ - 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; - 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 +template +void TriangularBase::evalToLazy(MatrixBase &other) const +{ + other.derived().resize(this->rows(), this->cols()); + internal::call_triangular_assignment_loop(other.derived(), derived().nestedExpression()); +} + +#endif // EIGEN_ENABLE_EVALUATORS } // end namespace Eigen