From 626821b0e34a624e8fa8980339b771e155722ace Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Mon, 2 Dec 2013 14:06:17 +0100 Subject: [PATCH] Add evaluator/assignment to TriangularView expressions --- Eigen/src/Core/TriangularMatrix.h | 317 +++++++++++++++++++++++++++++- test/evaluators.cpp | 35 +++- 2 files changed, 350 insertions(+), 2 deletions(-) diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index 1d6e34650..96ed9cef9 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -36,7 +36,13 @@ template class TriangularBase : public EigenBase RowsAtCompileTime = internal::traits::RowsAtCompileTime, ColsAtCompileTime = internal::traits::ColsAtCompileTime, MaxRowsAtCompileTime = internal::traits::MaxRowsAtCompileTime, - MaxColsAtCompileTime = internal::traits::MaxColsAtCompileTime + MaxColsAtCompileTime = internal::traits::MaxColsAtCompileTime, + + SizeAtCompileTime = (internal::size_at_compile_time::RowsAtCompileTime, + internal::traits::ColsAtCompileTime>::ret) + /**< This is equal to the number of coefficients, i.e. the number of + * rows times the number of columns, or to \a Dynamic if this is not + * known at compile-time. \sa RowsAtCompileTime, ColsAtCompileTime */ }; typedef typename internal::traits::Scalar Scalar; typedef typename internal::traits::StorageKind StorageKind; @@ -408,15 +414,24 @@ template class TriangularView EIGEN_DEVICE_FUNC void swap(TriangularBase const & other) { +// #ifdef EIGEN_TEST_EVALUATORS +// swap_using_evaluator(this->derived(), other.derived()); +// #else TriangularView,Mode>(const_cast(m_matrix)).lazyAssign(other.derived()); +// #endif } + // TODO: this overload is ambiguous and it should be deprecated (Gael) template EIGEN_DEVICE_FUNC void swap(MatrixBase const & other) { +// #ifdef EIGEN_TEST_EVALUATORS +// swap_using_evaluator(this->derived(), other.derived()); +// #else SwapWrapper swaper(const_cast(m_matrix)); TriangularView,Mode>(swaper).lazyAssign(other.derived()); +// #endif } EIGEN_DEVICE_FUNC @@ -895,6 +910,306 @@ bool MatrixBase::isLowerTriangular(const RealScalar& prec) const return true; } + +#ifdef EIGEN_ENABLE_EVALUATORS + +/*************************************************************************** +**************************************************************************** +* Evaluators and Assignment of triangular expressions +*************************************************************************** +***************************************************************************/ + +namespace internal { + + +// TODO currently a triangular expression has the form TriangularView<.,.> +// in the future triangular-ness should be defined by the expression traits +// such that Transpose > is valid. (currently TriangularBase::transpose() is overloaded to make it work) +template +struct evaluator_traits > +{ + typedef typename storage_kind_to_evaluator_kind::Kind Kind; + typedef TriangularShape Shape; + + // 1 if assignment A = B assumes aliasing when B is of type T and thus B needs to be evaluated into a + // temporary; 0 if not. + static const int AssumeAliasing = 0; +}; + +template +struct evaluator, Kind, typename MatrixType::Scalar> + : evaluator +{ + typedef TriangularView XprType; + typedef evaluator Base; + typedef evaluator type; + evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {} +}; + +// Additional assignement kinds: +struct Triangular2Triangular {}; +struct Triangular2Dense {}; +struct Dense2Triangular {}; + + +template struct triangular_assignment_loop; + + +template +void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src, const Functor &func) +{ +#ifdef EIGEN_DEBUG_ASSIGN + // TODO these traits should be computed from information provided by the evaluators + internal::copy_using_evaluator_traits::debug(); +#endif + eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); + + typedef typename evaluator::type DstEvaluatorType; + typedef typename evaluator::type SrcEvaluatorType; + + DstEvaluatorType dstEvaluator(dst); + SrcEvaluatorType srcEvaluator(src); + + typedef generic_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 + }; + + triangular_assignment_loop::run(kernel); +} + +template +void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src) +{ + call_triangular_assignment_loop(dst, src, internal::assign_op()); +} + + +template<> struct AssignmentKind { typedef Triangular2Triangular Kind; }; +template<> struct AssignmentKind { typedef Triangular2Dense Kind; }; +template<> struct AssignmentKind { typedef Dense2Triangular Kind; }; + + +template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> +struct Assignment +{ + static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) + { + eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode)); + + call_triangular_assignment_loop(dst, src, func); + } +}; + +template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> +struct Assignment +{ + static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) + { + call_triangular_assignment_loop(dst, src, func); + } +}; + +template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> +struct Assignment +{ + static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) + { + call_triangular_assignment_loop(dst, src, func); + } +}; + + +template +struct triangular_assignment_loop +{ + enum { + col = (UnrollCount-1) / Kernel::RowsAtCompileTime, + row = (UnrollCount-1) % Kernel::RowsAtCompileTime + }; + + typedef typename Kernel::Scalar Scalar; + + EIGEN_DEVICE_FUNC + 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)); + } + } +}; + +// prevent buggy user code from causing an infinite recursion +template +struct triangular_assignment_loop +{ + EIGEN_DEVICE_FUNC + 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 +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 +{ + 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 +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) + { + for(; i < maxi; ++i) + kernel.assignCoeff(i, j, Scalar(0)); + } + + if(i(), MatrixXd(A.triangularView())); + + A.setRandom();B.setRandom(); + VERIFY_IS_APPROX_EVALUATOR2(B, A.triangularView(), MatrixXd(A.triangularView())); + + A.setRandom();B.setRandom(); + VERIFY_IS_APPROX_EVALUATOR2(B, A.triangularView(), MatrixXd(A.triangularView())); + + A.setRandom();B.setRandom(); + C = B; C.triangularView() = A; + copy_using_evaluator(B.triangularView(), A); + VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView(), A)"); + + A.setRandom();B.setRandom(); + C = B; C.triangularView() = A.triangularView(); + copy_using_evaluator(B.triangularView(), A.triangularView()); + VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView(), A.triangularView())"); + + + A.setRandom();B.setRandom(); + C = B; C.triangularView() = A.triangularView().transpose(); + copy_using_evaluator(B.triangularView(), A.triangularView().transpose()); + VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView(), A.triangularView().transpose())"); + } }