diff options
author | Gael Guennebaud <g.gael@free.fr> | 2014-01-25 23:02:14 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2014-01-25 23:02:14 +0100 |
commit | 5fa7262e4c9fdbac2e5ddc4e2bb2241e04a10bef (patch) | |
tree | 2299d1899a1400668141cc0df18c59f4d292a3d6 /Eigen/src/Core/TriangularMatrix.h | |
parent | fef534f52e403e24637209f1ad843111c81122bd (diff) |
Refactor triangular assignment
Diffstat (limited to 'Eigen/src/Core/TriangularMatrix.h')
-rw-r--r-- | Eigen/src/Core/TriangularMatrix.h | 301 |
1 files changed, 157 insertions, 144 deletions
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 |