diff options
| -rw-r--r-- | Eigen/src/Core/TriangularMatrix.h | 317 | ||||
| -rw-r--r-- | 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<typename Derived> class TriangularBase : public EigenBase<Derived> RowsAtCompileTime = internal::traits<Derived>::RowsAtCompileTime, ColsAtCompileTime = internal::traits<Derived>::ColsAtCompileTime, MaxRowsAtCompileTime = internal::traits<Derived>::MaxRowsAtCompileTime, - MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime + MaxColsAtCompileTime = internal::traits<Derived>::MaxColsAtCompileTime, + + SizeAtCompileTime = (internal::size_at_compile_time<internal::traits<Derived>::RowsAtCompileTime, + internal::traits<Derived>::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<Derived>::Scalar Scalar; typedef typename internal::traits<Derived>::StorageKind StorageKind; @@ -408,15 +414,24 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView EIGEN_DEVICE_FUNC void swap(TriangularBase<OtherDerived> const & other) { +// #ifdef EIGEN_TEST_EVALUATORS +// swap_using_evaluator(this->derived(), other.derived()); +// #else TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); +// #endif } + // TODO: this overload is ambiguous and it should be deprecated (Gael) template<typename OtherDerived> EIGEN_DEVICE_FUNC void swap(MatrixBase<OtherDerived> const & other) { +// #ifdef EIGEN_TEST_EVALUATORS +// swap_using_evaluator(this->derived(), other.derived()); +// #else SwapWrapper<MatrixType> swaper(const_cast<MatrixType&>(m_matrix)); TriangularView<SwapWrapper<MatrixType>,Mode>(swaper).lazyAssign(other.derived()); +// #endif } EIGEN_DEVICE_FUNC @@ -895,6 +910,306 @@ bool MatrixBase<Derived>::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<TriangularView<.,.> > is valid. (currently TriangularBase::transpose() is overloaded to make it work) +template<typename MatrixType, unsigned int Mode> +struct evaluator_traits<TriangularView<MatrixType,Mode> > +{ + typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::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<typename MatrixType, unsigned int Mode, typename Kind> +struct evaluator<TriangularView<MatrixType,Mode>, Kind, typename MatrixType::Scalar> + : evaluator<MatrixType> +{ + typedef TriangularView<MatrixType,Mode> XprType; + typedef evaluator<MatrixType> Base; + typedef evaluator type; + evaluator(const XprType &xpr) : Base(xpr.nestedExpression()) {} +}; + +// Additional assignement kinds: +struct Triangular2Triangular {}; +struct Triangular2Dense {}; +struct Dense2Triangular {}; + + +template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> struct triangular_assignment_loop; + + +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) +{ +#ifdef EIGEN_DEBUG_ASSIGN + // TODO these traits should be computed from information provided by the evaluators + internal::copy_using_evaluator_traits<DstXprType, SrcXprType>::debug(); +#endif + eigen_assert(dst.rows() == src.rows() && dst.cols() == src.cols()); + + typedef typename evaluator<DstXprType>::type DstEvaluatorType; + typedef typename evaluator<SrcXprType>::type SrcEvaluatorType; + + DstEvaluatorType dstEvaluator(dst); + SrcEvaluatorType srcEvaluator(src); + + typedef generic_dense_assignment_kernel<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 + }; + + triangular_assignment_loop<Kernel, Mode, unroll ? int(DstXprType::SizeAtCompileTime) : Dynamic, ClearOpposite>::run(kernel); +} + +template<int Mode, bool ClearOpposite, typename DstXprType, typename SrcXprType> +void call_triangular_assignment_loop(const DstXprType& dst, const SrcXprType& src) +{ + call_triangular_assignment_loop<Mode,ClearOpposite>(dst, src, internal::assign_op<typename DstXprType::Scalar>()); +} + + +template<> struct AssignmentKind<TriangularShape,TriangularShape> { typedef Triangular2Triangular Kind; }; +template<> struct AssignmentKind<DenseShape,TriangularShape> { typedef Triangular2Dense Kind; }; +template<> struct AssignmentKind<TriangularShape,DenseShape> { typedef Dense2Triangular Kind; }; + + +template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> +struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Triangular, Scalar> +{ + static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) + { + eigen_assert(int(DstXprType::Mode) == int(SrcXprType::Mode)); + + call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func); + } +}; + +template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> +struct Assignment<DstXprType, SrcXprType, Functor, Triangular2Dense, Scalar> +{ + static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) + { + call_triangular_assignment_loop<SrcXprType::Mode, true>(dst, src, func); + } +}; + +template< typename DstXprType, typename SrcXprType, typename Functor, typename Scalar> +struct Assignment<DstXprType, SrcXprType, Functor, Dense2Triangular, Scalar> +{ + static void run(DstXprType &dst, const SrcXprType &src, const Functor &func) + { + call_triangular_assignment_loop<DstXprType::Mode, false>(dst, src, func); + } +}; + + +template<typename Kernel, unsigned int Mode, int UnrollCount, bool ClearOpposite> +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<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)); + } + } +}; + +// prevent buggy user code from causing an infinite recursion +template<typename Kernel, unsigned int Mode, bool ClearOpposite> +struct triangular_assignment_loop<Kernel, Mode, 0, ClearOpposite> +{ + 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<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); + } + } +}; + +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> +{ + 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<kernel.rows()) // then i==j + kernel.assignCoeff(i++, j, Scalar(1)); + + for(; i < kernel.rows(); ++i) + kernel.assignCoeff(i, j); + } + } +}; + +} // end namespace internal + +#endif + } // end namespace Eigen #endif // EIGEN_TRIANGULARMATRIX_H diff --git a/test/evaluators.cpp b/test/evaluators.cpp index 3fb0d9896..29192d7f9 100644 --- a/test/evaluators.cpp +++ b/test/evaluators.cpp @@ -3,6 +3,10 @@ #define EIGEN_ENABLE_EVALUATORS #endif +#ifdef EIGEN_TEST_EVALUATORS +#undef EIGEN_TEST_EVALUATORS +#endif + #include "main.h" namespace Eigen { @@ -101,7 +105,7 @@ void test_evaluators() copy_using_evaluator(w.transpose(), v_const); VERIFY_IS_APPROX(w,v_const.transpose().eval()); - +#if 0 // Testing Array evaluator { ArrayXXf a(2,3); @@ -401,4 +405,33 @@ void test_evaluators() arr_ref.row(1) /= (arr_ref.row(2) + 1); VERIFY_IS_APPROX(arr, arr_ref); } +#endif + { + // test triangular shapes + MatrixXd A = MatrixXd::Random(6,6), B(6,6), C(6,6); + A.setRandom();B.setRandom(); + VERIFY_IS_APPROX_EVALUATOR2(B, A.triangularView<Upper>(), MatrixXd(A.triangularView<Upper>())); + + A.setRandom();B.setRandom(); + VERIFY_IS_APPROX_EVALUATOR2(B, A.triangularView<UnitLower>(), MatrixXd(A.triangularView<UnitLower>())); + + A.setRandom();B.setRandom(); + VERIFY_IS_APPROX_EVALUATOR2(B, A.triangularView<UnitUpper>(), MatrixXd(A.triangularView<UnitUpper>())); + + A.setRandom();B.setRandom(); + C = B; C.triangularView<Upper>() = A; + copy_using_evaluator(B.triangularView<Upper>(), A); + VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView<Upper>(), A)"); + + A.setRandom();B.setRandom(); + C = B; C.triangularView<Lower>() = A.triangularView<Lower>(); + copy_using_evaluator(B.triangularView<Lower>(), A.triangularView<Lower>()); + VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView<Lower>(), A.triangularView<Lower>())"); + + + A.setRandom();B.setRandom(); + C = B; C.triangularView<Lower>() = A.triangularView<Upper>().transpose(); + copy_using_evaluator(B.triangularView<Lower>(), A.triangularView<Upper>().transpose()); + VERIFY(B.isApprox(C) && "copy_using_evaluator(B.triangularView<Lower>(), A.triangularView<Lower>().transpose())"); + } } |