diff options
author | Gael Guennebaud <g.gael@free.fr> | 2013-08-01 16:26:57 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2013-08-01 16:26:57 +0200 |
commit | ddf775363147fc7ee778b42c21b642f085193f55 (patch) | |
tree | dc08320f7a4dd5812599a84f0c8fa2c31a70e210 /Eigen | |
parent | 55b57fcba6e56bea5c084cc756b50a447985e5c2 (diff) |
Add nvcc support for small eigenvalues decompositions and workaround lack of support for std::swap and std::numeric_limits
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/Core | 2 | ||||
-rw-r--r-- | Eigen/src/Core/Assign.h | 8 | ||||
-rw-r--r-- | Eigen/src/Core/CwiseNullaryOp.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/DenseBase.h | 6 | ||||
-rw-r--r-- | Eigen/src/Core/Matrix.h | 1 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 40 | ||||
-rw-r--r-- | Eigen/src/Core/NumTraits.h | 12 | ||||
-rw-r--r-- | Eigen/src/Core/PlainObjectBase.h | 11 | ||||
-rw-r--r-- | Eigen/src/Core/SelfAdjointView.h | 24 | ||||
-rw-r--r-- | Eigen/src/Core/Swap.h | 14 | ||||
-rw-r--r-- | Eigen/src/Core/TriangularMatrix.h | 76 | ||||
-rw-r--r-- | Eigen/src/Core/util/Meta.h | 24 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 40 | ||||
-rw-r--r-- | Eigen/src/Geometry/OrthoMethods.h | 3 |
14 files changed, 234 insertions, 29 deletions
diff --git a/Eigen/Core b/Eigen/Core index d15795aed..1aaac5ae5 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -30,7 +30,7 @@ #endif #if defined(__CUDA_ARCH__) - #define EIGEN_USING_STD_MATH(FUNC) + #define EIGEN_USING_STD_MATH(FUNC) using ::FUNC; #else #define EIGEN_USING_STD_MATH(FUNC) using std::FUNC; #endif diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index b5a5a9fdb..906adcf82 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -532,6 +532,7 @@ struct assign_selector<Derived,OtherDerived,false,false> { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Derived& run(Derived& dst, const OtherDerived& other) { return dst.lazyAssign(other.derived()); } template<typename ActualDerived, typename ActualOtherDerived> + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Derived& evalTo(ActualDerived& dst, const ActualOtherDerived& other) { other.evalTo(dst); return dst; } }; template<typename Derived, typename OtherDerived> @@ -544,6 +545,7 @@ struct assign_selector<Derived,OtherDerived,false,true> { EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Derived& run(Derived& dst, const OtherDerived& other) { return dst.lazyAssign(other.transpose()); } template<typename ActualDerived, typename ActualOtherDerived> + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Derived& evalTo(ActualDerived& dst, const ActualOtherDerived& other) { Transpose<ActualDerived> dstTrans(dst); other.evalTo(dstTrans); return dst; } }; template<typename Derived, typename OtherDerived> @@ -556,18 +558,21 @@ struct assign_selector<Derived,OtherDerived,true,true> { template<typename Derived> template<typename OtherDerived> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::operator=(const DenseBase<OtherDerived>& other) { return internal::assign_selector<Derived,OtherDerived>::run(derived(), other.derived()); } template<typename Derived> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& DenseBase<Derived>::operator=(const DenseBase& other) { return internal::assign_selector<Derived,Derived>::run(derived(), other.derived()); } template<typename Derived> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const MatrixBase& other) { return internal::assign_selector<Derived,Derived>::run(derived(), other.derived()); @@ -575,6 +580,7 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const MatrixBase& ot template<typename Derived> template <typename OtherDerived> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const DenseBase<OtherDerived>& other) { return internal::assign_selector<Derived,OtherDerived>::run(derived(), other.derived()); @@ -582,6 +588,7 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const DenseBase<Othe template<typename Derived> template <typename OtherDerived> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const EigenBase<OtherDerived>& other) { return internal::assign_selector<Derived,OtherDerived,false>::evalTo(derived(), other.derived()); @@ -589,6 +596,7 @@ EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const EigenBase<Othe template<typename Derived> template<typename OtherDerived> +EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Derived& MatrixBase<Derived>::operator=(const ReturnByValue<OtherDerived>& other) { return internal::assign_selector<Derived,OtherDerived,false>::evalTo(derived(), other.derived()); diff --git a/Eigen/src/Core/CwiseNullaryOp.h b/Eigen/src/Core/CwiseNullaryOp.h index 14c78cb85..1d4ee50a8 100644 --- a/Eigen/src/Core/CwiseNullaryOp.h +++ b/Eigen/src/Core/CwiseNullaryOp.h @@ -746,6 +746,7 @@ namespace internal { template<typename Derived, bool Big = (Derived::SizeAtCompileTime>=16)> struct setIdentity_impl { + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Derived& run(Derived& m) { return m = Derived::Identity(m.rows(), m.cols()); @@ -756,6 +757,7 @@ template<typename Derived> struct setIdentity_impl<Derived, true> { typedef typename Derived::Index Index; + EIGEN_DEVICE_FUNC static EIGEN_STRONG_INLINE Derived& run(Derived& m) { m.setZero(); diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h index 097717075..889e2eafb 100644 --- a/Eigen/src/Core/DenseBase.h +++ b/Eigen/src/Core/DenseBase.h @@ -391,7 +391,8 @@ template<typename Derived> class DenseBase /** swaps *this with the expression \a other. * */ - template<typename OtherDerived> EIGEN_DEVICE_FUNC + template<typename OtherDerived> + EIGEN_DEVICE_FUNC void swap(const DenseBase<OtherDerived>& other, int = OtherDerived::ThisConstantIsPrivateInPlainObjectBase) { @@ -401,7 +402,8 @@ template<typename Derived> class DenseBase /** swaps *this with the matrix or array \a other. * */ - template<typename OtherDerived> EIGEN_DEVICE_FUNC + template<typename OtherDerived> + EIGEN_DEVICE_FUNC void swap(PlainObjectBase<OtherDerived>& other) { SwapWrapper<Derived>(derived()).lazyAssign(other.derived()); diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 8c3780ca2..b80d97c8c 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -331,6 +331,7 @@ class Matrix * of same type it is enough to swap the data pointers. */ template<typename OtherDerived> + EIGEN_DEVICE_FUNC void swap(MatrixBase<OtherDerived> const & other) { this->_swap(other.derived()); } diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 8b8a868b7..ff4526d75 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -278,31 +278,41 @@ template<typename Derived> class MatrixBase template<unsigned int Mode> struct TriangularViewReturnType { typedef TriangularView<Derived, Mode> Type; }; template<unsigned int Mode> struct ConstTriangularViewReturnType { typedef const TriangularView<const Derived, Mode> Type; }; - template<unsigned int Mode> typename TriangularViewReturnType<Mode>::Type triangularView(); - template<unsigned int Mode> typename ConstTriangularViewReturnType<Mode>::Type triangularView() const; + template<unsigned int Mode> + EIGEN_DEVICE_FUNC + typename TriangularViewReturnType<Mode>::Type triangularView(); + template<unsigned int Mode> + EIGEN_DEVICE_FUNC + typename ConstTriangularViewReturnType<Mode>::Type triangularView() const; template<unsigned int UpLo> struct SelfAdjointViewReturnType { typedef SelfAdjointView<Derived, UpLo> Type; }; template<unsigned int UpLo> struct ConstSelfAdjointViewReturnType { typedef const SelfAdjointView<const Derived, UpLo> Type; }; - template<unsigned int UpLo> typename SelfAdjointViewReturnType<UpLo>::Type selfadjointView(); - template<unsigned int UpLo> typename ConstSelfAdjointViewReturnType<UpLo>::Type selfadjointView() const; + template<unsigned int UpLo> + EIGEN_DEVICE_FUNC + typename SelfAdjointViewReturnType<UpLo>::Type selfadjointView(); + template<unsigned int UpLo> + EIGEN_DEVICE_FUNC + typename ConstSelfAdjointViewReturnType<UpLo>::Type selfadjointView() const; const SparseView<Derived> sparseView(const Scalar& m_reference = Scalar(0), const typename NumTraits<Scalar>::Real& m_epsilon = NumTraits<Scalar>::dummy_precision()) const; - static const IdentityReturnType Identity(); - static const IdentityReturnType Identity(Index rows, Index cols); - static const BasisReturnType Unit(Index size, Index i); - static const BasisReturnType Unit(Index i); - static const BasisReturnType UnitX(); - static const BasisReturnType UnitY(); - static const BasisReturnType UnitZ(); - static const BasisReturnType UnitW(); + EIGEN_DEVICE_FUNC static const IdentityReturnType Identity(); + EIGEN_DEVICE_FUNC static const IdentityReturnType Identity(Index rows, Index cols); + EIGEN_DEVICE_FUNC static const BasisReturnType Unit(Index size, Index i); + EIGEN_DEVICE_FUNC static const BasisReturnType Unit(Index i); + EIGEN_DEVICE_FUNC static const BasisReturnType UnitX(); + EIGEN_DEVICE_FUNC static const BasisReturnType UnitY(); + EIGEN_DEVICE_FUNC static const BasisReturnType UnitZ(); + EIGEN_DEVICE_FUNC static const BasisReturnType UnitW(); EIGEN_DEVICE_FUNC const DiagonalWrapper<const Derived> asDiagonal() const; const PermutationWrapper<const Derived> asPermutation() const; + EIGEN_DEVICE_FUNC Derived& setIdentity(); + EIGEN_DEVICE_FUNC Derived& setIdentity(Index rows, Index cols); bool isIdentity(const RealScalar& prec = NumTraits<Scalar>::dummy_precision()) const; @@ -430,11 +440,17 @@ template<typename Derived> class MatrixBase }; #endif // EIGEN_PARSED_BY_DOXYGEN template<typename OtherDerived> + EIGEN_DEVICE_FUNC typename cross_product_return_type<OtherDerived>::type cross(const MatrixBase<OtherDerived>& other) const; + template<typename OtherDerived> + EIGEN_DEVICE_FUNC PlainObject cross3(const MatrixBase<OtherDerived>& other) const; + + EIGEN_DEVICE_FUNC PlainObject unitOrthogonal(void) const; + Matrix<Scalar,3,1> eulerAngles(Index a0, Index a1, Index a2) const; #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h index eb5306bf2..2b6633c9c 100644 --- a/Eigen/src/Core/NumTraits.h +++ b/Eigen/src/Core/NumTraits.h @@ -68,16 +68,22 @@ template<typename T> struct GenericNumTraits >::type NonInteger; typedef T Nested; - static inline Real epsilon() { return std::numeric_limits<T>::epsilon(); } + EIGEN_DEVICE_FUNC + static inline Real epsilon() + { + #if defined(__CUDA_ARCH__) + return internal::device::numeric_limits<T>::epsilon(); + #else + return std::numeric_limits<T>::epsilon(); + #endif + } EIGEN_DEVICE_FUNC static inline Real dummy_precision() { // make sure to override this for floating-point types return Real(0); } - EIGEN_DEVICE_FUNC static inline T highest() { return (std::numeric_limits<T>::max)(); } - EIGEN_DEVICE_FUNC static inline T lowest() { return IsInteger ? (std::numeric_limits<T>::min)() : (-(std::numeric_limits<T>::max)()); } #ifdef EIGEN2_SUPPORT diff --git a/Eigen/src/Core/PlainObjectBase.h b/Eigen/src/Core/PlainObjectBase.h index 39d422f78..618e00e02 100644 --- a/Eigen/src/Core/PlainObjectBase.h +++ b/Eigen/src/Core/PlainObjectBase.h @@ -128,7 +128,9 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type enum { NeedsToAlign = SizeAtCompileTime != Dynamic && (internal::traits<Derived>::Flags & AlignedBit) != 0 }; EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) + EIGEN_DEVICE_FUNC Base& base() { return *static_cast<Base*>(this); } + EIGEN_DEVICE_FUNC const Base& base() const { return *static_cast<const Base*>(this); } EIGEN_DEVICE_FUNC @@ -351,6 +353,7 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type * Matrices are resized relative to the top-left element. In case values need to be * appended to the matrix they will be uninitialized. */ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void conservativeResize(Index nbRows, Index nbCols) { internal::conservative_resize_like_impl<Derived>::run(*this, nbRows, nbCols); @@ -363,6 +366,7 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type * * In case the matrix is growing, new rows will be uninitialized. */ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void conservativeResize(Index nbRows, NoChange_t) { // Note: see the comment in conservativeResize(Index,Index) @@ -376,6 +380,7 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type * * In case the matrix is growing, new columns will be uninitialized. */ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, Index nbCols) { // Note: see the comment in conservativeResize(Index,Index) @@ -390,6 +395,7 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type * * When values are appended, they will be uninitialized. */ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void conservativeResize(Index size) { internal::conservative_resize_like_impl<Derived>::run(*this, size); @@ -405,6 +411,7 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type * appended to the matrix they will copied from \c other. */ template<typename OtherDerived> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void conservativeResizeLike(const DenseBase<OtherDerived>& other) { internal::conservative_resize_like_impl<Derived,OtherDerived>::run(*this, other); @@ -647,6 +654,7 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type } template<typename T0, typename T1> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void _init2(Index nbRows, Index nbCols, typename internal::enable_if<Base::SizeAtCompileTime!=2,T0>::type* = 0) { EIGEN_STATIC_ASSERT(bool(NumTraits<T0>::IsInteger) && @@ -670,6 +678,7 @@ class PlainObjectBase : public internal::dense_xpr_base<Derived>::type * data pointers. */ template<typename OtherDerived> + EIGEN_DEVICE_FUNC void _swap(DenseBase<OtherDerived> const & other) { enum { SwapPointers = internal::is_same<Derived, OtherDerived>::value && Base::SizeAtCompileTime==Dynamic }; @@ -790,6 +799,7 @@ struct conservative_resize_like_impl<Derived,OtherDerived,true> template<typename MatrixTypeA, typename MatrixTypeB, bool SwapPointers> struct matrix_swap_impl { + EIGEN_DEVICE_FUNC static inline void run(MatrixTypeA& a, MatrixTypeB& b) { a.base().swap(b); @@ -799,6 +809,7 @@ struct matrix_swap_impl template<typename MatrixTypeA, typename MatrixTypeB> struct matrix_swap_impl<MatrixTypeA, MatrixTypeB, true> { + EIGEN_DEVICE_FUNC static inline void run(MatrixTypeA& a, MatrixTypeB& b) { static_cast<typename MatrixTypeA::Base&>(a).m_storage.swap(static_cast<typename MatrixTypeB::Base&>(b).m_storage); diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 6fa7cd15e..8231e3f5c 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -69,17 +69,23 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView }; typedef typename MatrixType::PlainObject PlainObject; + EIGEN_DEVICE_FUNC inline SelfAdjointView(MatrixType& matrix) : m_matrix(matrix) {} + EIGEN_DEVICE_FUNC inline Index rows() const { return m_matrix.rows(); } + EIGEN_DEVICE_FUNC inline Index cols() const { return m_matrix.cols(); } + EIGEN_DEVICE_FUNC inline Index outerStride() const { return m_matrix.outerStride(); } + EIGEN_DEVICE_FUNC inline Index innerStride() const { return m_matrix.innerStride(); } /** \sa MatrixBase::coeff() * \warning the coordinates must fit into the referenced triangular part */ + EIGEN_DEVICE_FUNC inline Scalar coeff(Index row, Index col) const { Base::check_coordinates_internal(row, col); @@ -89,6 +95,7 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView /** \sa MatrixBase::coeffRef() * \warning the coordinates must fit into the referenced triangular part */ + EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index row, Index col) { Base::check_coordinates_internal(row, col); @@ -96,13 +103,17 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView } /** \internal */ + EIGEN_DEVICE_FUNC const MatrixTypeNestedCleaned& _expression() const { return m_matrix; } + EIGEN_DEVICE_FUNC const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } + EIGEN_DEVICE_FUNC MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); } /** Efficient self-adjoint matrix times vector/matrix product */ template<typename OtherDerived> + EIGEN_DEVICE_FUNC SelfadjointProductMatrix<MatrixType,Mode,false,OtherDerived,0,OtherDerived::IsVectorAtCompileTime> operator*(const MatrixBase<OtherDerived>& rhs) const { @@ -113,6 +124,7 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView /** Efficient vector/matrix times self-adjoint matrix product */ template<typename OtherDerived> friend + EIGEN_DEVICE_FUNC SelfadjointProductMatrix<OtherDerived,0,OtherDerived::IsVectorAtCompileTime,MatrixType,Mode,false> operator*(const MatrixBase<OtherDerived>& lhs, const SelfAdjointView& rhs) { @@ -132,6 +144,7 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView * \sa rankUpdate(const MatrixBase<DerivedU>&, Scalar) */ template<typename DerivedU, typename DerivedV> + EIGEN_DEVICE_FUNC SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const MatrixBase<DerivedV>& v, const Scalar& alpha = Scalar(1)); /** Perform a symmetric rank K update of the selfadjoint matrix \c *this: @@ -145,6 +158,7 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView * \sa rankUpdate(const MatrixBase<DerivedU>&, const MatrixBase<DerivedV>&, Scalar) */ template<typename DerivedU> + EIGEN_DEVICE_FUNC SelfAdjointView& rankUpdate(const MatrixBase<DerivedU>& u, const Scalar& alpha = Scalar(1)); /////////// Cholesky module /////////// @@ -159,11 +173,14 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView /** Return type of eigenvalues() */ typedef Matrix<RealScalar, internal::traits<MatrixType>::ColsAtCompileTime, 1> EigenvaluesReturnType; + EIGEN_DEVICE_FUNC EigenvaluesReturnType eigenvalues() const; + EIGEN_DEVICE_FUNC RealScalar operatorNorm() const; #ifdef EIGEN2_SUPPORT template<typename OtherDerived> + EIGEN_DEVICE_FUNC SelfAdjointView& operator=(const MatrixBase<OtherDerived>& other) { enum { @@ -174,6 +191,7 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView return *this; } template<typename OtherMatrixType, unsigned int OtherMode> + EIGEN_DEVICE_FUNC SelfAdjointView& operator=(const TriangularView<OtherMatrixType, OtherMode>& other) { enum { @@ -209,6 +227,7 @@ struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), U row = (UnrollCount-1) % Derived1::RowsAtCompileTime }; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), UnrollCount-1, ClearOpposite>::run(dst, src); @@ -223,6 +242,7 @@ struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Upper), U template<typename Derived1, typename Derived2, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, 0, ClearOpposite> { + EIGEN_DEVICE_FUNC static inline void run(Derived1 &, const Derived2 &) {} }; @@ -234,6 +254,7 @@ struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), U row = (UnrollCount-1) % Derived1::RowsAtCompileTime }; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), UnrollCount-1, ClearOpposite>::run(dst, src); @@ -248,6 +269,7 @@ struct triangular_assignment_selector<Derived1, Derived2, (SelfAdjoint|Lower), U template<typename Derived1, typename Derived2, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, 0, ClearOpposite> { + EIGEN_DEVICE_FUNC static inline void run(Derived1 &, const Derived2 &) {} }; @@ -255,6 +277,7 @@ template<typename Derived1, typename Derived2, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dynamic, ClearOpposite> { typedef typename Derived1::Index Index; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { for(Index j = 0; j < dst.cols(); ++j) @@ -272,6 +295,7 @@ struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Upper, Dyn template<typename Derived1, typename Derived2, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, SelfAdjoint|Lower, Dynamic, ClearOpposite> { + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { typedef typename Derived1::Index Index; diff --git a/Eigen/src/Core/Swap.h b/Eigen/src/Core/Swap.h index bf58bd599..d602fba65 100644 --- a/Eigen/src/Core/Swap.h +++ b/Eigen/src/Core/Swap.h @@ -33,11 +33,16 @@ template<typename ExpressionType> class SwapWrapper EIGEN_DENSE_PUBLIC_INTERFACE(SwapWrapper) typedef typename internal::packet_traits<Scalar>::type Packet; + EIGEN_DEVICE_FUNC inline SwapWrapper(ExpressionType& xpr) : m_expression(xpr) {} + EIGEN_DEVICE_FUNC inline Index rows() const { return m_expression.rows(); } + EIGEN_DEVICE_FUNC inline Index cols() const { return m_expression.cols(); } + EIGEN_DEVICE_FUNC inline Index outerStride() const { return m_expression.outerStride(); } + EIGEN_DEVICE_FUNC inline Index innerStride() const { return m_expression.innerStride(); } typedef typename internal::conditional< @@ -46,30 +51,37 @@ template<typename ExpressionType> class SwapWrapper const Scalar >::type ScalarWithConstIfNotLvalue; + EIGEN_DEVICE_FUNC inline ScalarWithConstIfNotLvalue* data() { return m_expression.data(); } + EIGEN_DEVICE_FUNC inline const Scalar* data() const { return m_expression.data(); } + EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index rowId, Index colId) { return m_expression.const_cast_derived().coeffRef(rowId, colId); } + EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index index) { return m_expression.const_cast_derived().coeffRef(index); } + EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index rowId, Index colId) const { return m_expression.coeffRef(rowId, colId); } + EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index index) const { return m_expression.coeffRef(index); } template<typename OtherDerived> + EIGEN_DEVICE_FUNC void copyCoeff(Index rowId, Index colId, const DenseBase<OtherDerived>& other) { OtherDerived& _other = other.const_cast_derived(); @@ -81,6 +93,7 @@ template<typename ExpressionType> class SwapWrapper } template<typename OtherDerived> + EIGEN_DEVICE_FUNC void copyCoeff(Index index, const DenseBase<OtherDerived>& other) { OtherDerived& _other = other.const_cast_derived(); @@ -115,6 +128,7 @@ template<typename ExpressionType> class SwapWrapper _other.template writePacket<LoadMode>(index, tmp); } + EIGEN_DEVICE_FUNC ExpressionType& expression() const { return m_expression; } protected: diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index fba07365f..1d6e34650 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -44,29 +44,39 @@ template<typename Derived> class TriangularBase : public EigenBase<Derived> typedef typename internal::traits<Derived>::DenseMatrixType DenseMatrixType; typedef DenseMatrixType DenseType; + EIGEN_DEVICE_FUNC inline TriangularBase() { eigen_assert(!((Mode&UnitDiag) && (Mode&ZeroDiag))); } + EIGEN_DEVICE_FUNC inline Index rows() const { return derived().rows(); } + EIGEN_DEVICE_FUNC inline Index cols() const { return derived().cols(); } + EIGEN_DEVICE_FUNC inline Index outerStride() const { return derived().outerStride(); } + EIGEN_DEVICE_FUNC inline Index innerStride() const { return derived().innerStride(); } + EIGEN_DEVICE_FUNC inline Scalar coeff(Index row, Index col) const { return derived().coeff(row,col); } + EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index row, Index col) { return derived().coeffRef(row,col); } /** \see MatrixBase::copyCoeff(row,col) */ template<typename Other> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void copyCoeff(Index row, Index col, Other& other) { derived().coeffRef(row, col) = other.coeff(row, col); } + EIGEN_DEVICE_FUNC inline Scalar operator()(Index row, Index col) const { check_coordinates(row, col); return coeff(row,col); } + EIGEN_DEVICE_FUNC inline Scalar& operator()(Index row, Index col) { check_coordinates(row, col); @@ -74,15 +84,20 @@ template<typename Derived> class TriangularBase : public EigenBase<Derived> } #ifndef EIGEN_PARSED_BY_DOXYGEN + EIGEN_DEVICE_FUNC inline const Derived& derived() const { return *static_cast<const Derived*>(this); } + EIGEN_DEVICE_FUNC inline Derived& derived() { return *static_cast<Derived*>(this); } #endif // not EIGEN_PARSED_BY_DOXYGEN template<typename DenseDerived> + EIGEN_DEVICE_FUNC void evalTo(MatrixBase<DenseDerived> &other) const; template<typename DenseDerived> + EIGEN_DEVICE_FUNC void evalToLazy(MatrixBase<DenseDerived> &other) const; + EIGEN_DEVICE_FUNC DenseMatrixType toDenseMatrix() const { DenseMatrixType res(rows(), cols()); @@ -189,36 +204,52 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView | (Mode & (ZeroDiag)) }; + EIGEN_DEVICE_FUNC inline TriangularView(const MatrixType& matrix) : m_matrix(matrix) {} + EIGEN_DEVICE_FUNC inline Index rows() const { return m_matrix.rows(); } + EIGEN_DEVICE_FUNC inline Index cols() const { return m_matrix.cols(); } + EIGEN_DEVICE_FUNC inline Index outerStride() const { return m_matrix.outerStride(); } + EIGEN_DEVICE_FUNC inline Index innerStride() const { return m_matrix.innerStride(); } - /** \sa MatrixBase::operator+=() */ - template<typename Other> TriangularView& operator+=(const DenseBase<Other>& other) { return *this = m_matrix + other.derived(); } + /** \sa MatrixBase::operator+=() */ + template<typename Other> + EIGEN_DEVICE_FUNC + TriangularView& operator+=(const DenseBase<Other>& other) { return *this = m_matrix + other.derived(); } /** \sa MatrixBase::operator-=() */ - template<typename Other> TriangularView& operator-=(const DenseBase<Other>& other) { return *this = m_matrix - other.derived(); } + template<typename Other> + EIGEN_DEVICE_FUNC + TriangularView& operator-=(const DenseBase<Other>& other) { return *this = m_matrix - other.derived(); } /** \sa MatrixBase::operator*=() */ + EIGEN_DEVICE_FUNC TriangularView& operator*=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix * other; } /** \sa MatrixBase::operator/=() */ + EIGEN_DEVICE_FUNC TriangularView& operator/=(const typename internal::traits<MatrixType>::Scalar& other) { return *this = m_matrix / other; } /** \sa MatrixBase::fill() */ + EIGEN_DEVICE_FUNC void fill(const Scalar& value) { setConstant(value); } /** \sa MatrixBase::setConstant() */ + EIGEN_DEVICE_FUNC TriangularView& setConstant(const Scalar& value) { return *this = MatrixType::Constant(rows(), cols(), value); } /** \sa MatrixBase::setZero() */ + EIGEN_DEVICE_FUNC TriangularView& setZero() { return setConstant(Scalar(0)); } /** \sa MatrixBase::setOnes() */ + EIGEN_DEVICE_FUNC TriangularView& setOnes() { return setConstant(Scalar(1)); } /** \sa MatrixBase::coeff() * \warning the coordinates must fit into the referenced triangular part */ + EIGEN_DEVICE_FUNC inline Scalar coeff(Index row, Index col) const { Base::check_coordinates_internal(row, col); @@ -228,49 +259,62 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView /** \sa MatrixBase::coeffRef() * \warning the coordinates must fit into the referenced triangular part */ + EIGEN_DEVICE_FUNC inline Scalar& coeffRef(Index row, Index col) { Base::check_coordinates_internal(row, col); return m_matrix.const_cast_derived().coeffRef(row, col); } + EIGEN_DEVICE_FUNC const MatrixTypeNestedCleaned& nestedExpression() const { return m_matrix; } + EIGEN_DEVICE_FUNC MatrixTypeNestedCleaned& nestedExpression() { return *const_cast<MatrixTypeNestedCleaned*>(&m_matrix); } /** Assigns a triangular matrix to a triangular part of a dense matrix */ template<typename OtherDerived> + EIGEN_DEVICE_FUNC TriangularView& operator=(const TriangularBase<OtherDerived>& other); template<typename OtherDerived> + EIGEN_DEVICE_FUNC TriangularView& operator=(const MatrixBase<OtherDerived>& other); + EIGEN_DEVICE_FUNC TriangularView& operator=(const TriangularView& other) { return *this = other.nestedExpression(); } template<typename OtherDerived> + EIGEN_DEVICE_FUNC void lazyAssign(const TriangularBase<OtherDerived>& other); template<typename OtherDerived> + EIGEN_DEVICE_FUNC void lazyAssign(const MatrixBase<OtherDerived>& other); /** \sa MatrixBase::conjugate() */ + EIGEN_DEVICE_FUNC inline TriangularView<MatrixConjugateReturnType,Mode> conjugate() { return m_matrix.conjugate(); } /** \sa MatrixBase::conjugate() const */ + EIGEN_DEVICE_FUNC inline const TriangularView<MatrixConjugateReturnType,Mode> conjugate() const { return m_matrix.conjugate(); } /** \sa MatrixBase::adjoint() const */ + EIGEN_DEVICE_FUNC inline const TriangularView<const typename MatrixType::AdjointReturnType,TransposeMode> adjoint() const { return m_matrix.adjoint(); } /** \sa MatrixBase::transpose() */ + EIGEN_DEVICE_FUNC inline TriangularView<Transpose<MatrixType>,TransposeMode> transpose() { EIGEN_STATIC_ASSERT_LVALUE(MatrixType) return m_matrix.const_cast_derived().transpose(); } /** \sa MatrixBase::transpose() const */ + EIGEN_DEVICE_FUNC inline const TriangularView<Transpose<MatrixType>,TransposeMode> transpose() const { return m_matrix.transpose(); @@ -278,6 +322,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView /** Efficient triangular matrix times vector/matrix product */ template<typename OtherDerived> + EIGEN_DEVICE_FUNC TriangularProduct<Mode,true,MatrixType,false,OtherDerived, OtherDerived::IsVectorAtCompileTime> operator*(const MatrixBase<OtherDerived>& rhs) const { @@ -288,6 +333,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView /** Efficient vector/matrix times triangular matrix product */ template<typename OtherDerived> friend + EIGEN_DEVICE_FUNC TriangularProduct<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false> operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs) { @@ -326,26 +372,32 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView #endif // EIGEN2_SUPPORT template<int Side, typename Other> + EIGEN_DEVICE_FUNC inline const internal::triangular_solve_retval<Side,TriangularView, Other> solve(const MatrixBase<Other>& other) const; template<int Side, typename OtherDerived> + EIGEN_DEVICE_FUNC void solveInPlace(const MatrixBase<OtherDerived>& other) const; template<typename Other> + EIGEN_DEVICE_FUNC inline const internal::triangular_solve_retval<OnTheLeft,TriangularView, Other> solve(const MatrixBase<Other>& other) const { return solve<OnTheLeft>(other); } template<typename OtherDerived> + EIGEN_DEVICE_FUNC void solveInPlace(const MatrixBase<OtherDerived>& other) const { return solveInPlace<OnTheLeft>(other); } + EIGEN_DEVICE_FUNC const SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() const { EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR); return SelfAdjointView<MatrixTypeNestedNonRef,Mode>(m_matrix); } + EIGEN_DEVICE_FUNC SelfAdjointView<MatrixTypeNestedNonRef,Mode> selfadjointView() { EIGEN_STATIC_ASSERT((Mode&UnitDiag)==0,PROGRAMMING_ERROR); @@ -353,18 +405,21 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView } template<typename OtherDerived> + EIGEN_DEVICE_FUNC void swap(TriangularBase<OtherDerived> const & other) { TriangularView<SwapWrapper<MatrixType>,Mode>(const_cast<MatrixType&>(m_matrix)).lazyAssign(other.derived()); } template<typename OtherDerived> + EIGEN_DEVICE_FUNC void swap(MatrixBase<OtherDerived> const & other) { SwapWrapper<MatrixType> swaper(const_cast<MatrixType&>(m_matrix)); TriangularView<SwapWrapper<MatrixType>,Mode>(swaper).lazyAssign(other.derived()); } + EIGEN_DEVICE_FUNC Scalar determinant() const { if (Mode & UnitDiag) @@ -377,6 +432,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView // TODO simplify the following: template<typename ProductDerived, typename Lhs, typename Rhs> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other) { setZero(); @@ -384,12 +440,14 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView } template<typename ProductDerived, typename Lhs, typename Rhs> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other) { return assignProduct(other,1); } template<typename ProductDerived, typename Lhs, typename Rhs> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other) { return assignProduct(other,-1); @@ -397,6 +455,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView template<typename ProductDerived> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other) { setZero(); @@ -404,12 +463,14 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView } template<typename ProductDerived> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other) { return assignProduct(other,other.alpha()); } template<typename ProductDerived> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other) { return assignProduct(other,-other.alpha()); @@ -418,6 +479,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView protected: template<typename ProductDerived, typename Lhs, typename Rhs> + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha); MatrixTypeNested m_matrix; @@ -439,6 +501,7 @@ struct triangular_assignment_selector typedef typename Derived1::Scalar Scalar; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { triangular_assignment_selector<Derived1, Derived2, Mode, UnrollCount-1, ClearOpposite>::run(dst, src); @@ -467,6 +530,7 @@ struct triangular_assignment_selector template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite> { + EIGEN_DEVICE_FUNC static inline void run(Derived1 &, const Derived2 &) {} }; @@ -475,6 +539,7 @@ struct triangular_assignment_selector<Derived1, Derived2, Upper, Dynamic, ClearO { typedef typename Derived1::Index Index; typedef typename Derived1::Scalar Scalar; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { for(Index j = 0; j < dst.cols(); ++j) @@ -493,6 +558,7 @@ template<typename Derived1, typename Derived2, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, Lower, Dynamic, ClearOpposite> { typedef typename Derived1::Index Index; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { for(Index j = 0; j < dst.cols(); ++j) @@ -512,6 +578,7 @@ struct triangular_assignment_selector<Derived1, Derived2, StrictlyUpper, Dynamic { typedef typename Derived1::Index Index; typedef typename Derived1::Scalar Scalar; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { for(Index j = 0; j < dst.cols(); ++j) @@ -530,6 +597,7 @@ template<typename Derived1, typename Derived2, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, StrictlyLower, Dynamic, ClearOpposite> { typedef typename Derived1::Index Index; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { for(Index j = 0; j < dst.cols(); ++j) @@ -548,6 +616,7 @@ template<typename Derived1, typename Derived2, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, UnitUpper, Dynamic, ClearOpposite> { typedef typename Derived1::Index Index; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { for(Index j = 0; j < dst.cols(); ++j) @@ -568,6 +637,7 @@ template<typename Derived1, typename Derived2, bool ClearOpposite> struct triangular_assignment_selector<Derived1, Derived2, UnitLower, Dynamic, ClearOpposite> { typedef typename Derived1::Index Index; + EIGEN_DEVICE_FUNC static inline void run(Derived1 &dst, const Derived2 &src) { for(Index j = 0; j < dst.cols(); ++j) diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 71d587108..aea168b46 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -88,7 +88,31 @@ template<bool Condition, typename T> struct enable_if; template<typename T> struct enable_if<true,T> { typedef T type; }; +#if defined(__CUDA_ARCH__) +template<typename T> EIGEN_DEVICE_FUNC void swap(T &a, T &b) { T tmp = b; b = a; a = tmp; } +namespace device { +template<typename T> struct numeric_limits +{ + EIGEN_DEVICE_FUNC + static T epsilon() { return 0; } +}; +template<> struct numeric_limits<float> +{ + EIGEN_DEVICE_FUNC + static float epsilon() { return __FLT_EPSILON__; } +}; +template<> struct numeric_limits<double> +{ + EIGEN_DEVICE_FUNC + static double epsilon() { return __DBL_EPSILON__; } +}; + +} + +#else +template<typename T> EIGEN_STRONG_INLINE void swap(T &a, T &b) { std::swap(a,b); } +#endif /** \internal * A base class do disable default copy ctor and copy assignement operator. diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 3993046a8..23a334b5c 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -109,6 +109,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver() : m_eivec(), m_eivalues(), @@ -128,6 +129,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa compute() for an example */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver(Index size) : m_eivec(size, size), m_eivalues(size), @@ -150,6 +152,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa compute(const MatrixType&, int) */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors) : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), @@ -189,6 +192,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa SelfAdjointEigenSolver(const MatrixType&, int) */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); /** \brief Computes eigendecomposition of given matrix using a direct algorithm @@ -205,6 +209,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa compute(const MatrixType&, int options) */ + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors); /** \brief Returns the eigenvectors of given matrix. @@ -225,6 +230,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa eigenvalues() */ + EIGEN_DEVICE_FUNC const MatrixType& eigenvectors() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -247,6 +253,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \sa eigenvectors(), MatrixBase::eigenvalues() */ + EIGEN_DEVICE_FUNC const RealVectorType& eigenvalues() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -271,6 +278,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * \sa operatorInverseSqrt(), * \ref MatrixFunctions_Module "MatrixFunctions Module" */ + EIGEN_DEVICE_FUNC MatrixType operatorSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -296,6 +304,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * \sa operatorSqrt(), MatrixBase::inverse(), * \ref MatrixFunctions_Module "MatrixFunctions Module" */ + EIGEN_DEVICE_FUNC MatrixType operatorInverseSqrt() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -307,6 +316,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver * * \returns \c Success if computation was succesful, \c NoConvergence otherwise. */ + EIGEN_DEVICE_FUNC ComputationInfo info() const { eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); @@ -321,6 +331,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver static const int m_maxIterations = 30; #ifdef EIGEN2_SUPPORT + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors) : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), @@ -330,6 +341,7 @@ template<typename _MatrixType> class SelfAdjointEigenSolver compute(matrix, computeEigenvectors); } + EIGEN_DEVICE_FUNC SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) : m_eivec(matA.cols(), matA.cols()), m_eivalues(matA.cols()), @@ -339,11 +351,13 @@ template<typename _MatrixType> class SelfAdjointEigenSolver static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); } + EIGEN_DEVICE_FUNC void compute(const MatrixType& matrix, bool computeEigenvectors) { compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); } + EIGEN_DEVICE_FUNC void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) { compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); @@ -377,10 +391,12 @@ template<typename _MatrixType> class SelfAdjointEigenSolver */ namespace internal { template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> +EIGEN_DEVICE_FUNC static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n); } template<typename MatrixType> +EIGEN_DEVICE_FUNC SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> ::compute(const MatrixType& matrix, int options) { @@ -481,6 +497,7 @@ namespace internal { template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues { + EIGEN_DEVICE_FUNC static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options) { eig.compute(A,options); } }; @@ -491,12 +508,13 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 typedef typename SolverType::RealVectorType VectorType; typedef typename SolverType::Scalar Scalar; + EIGEN_DEVICE_FUNC static inline void computeRoots(const MatrixType& m, VectorType& roots) { - using std::sqrt; - using std::atan2; - using std::cos; - using std::sin; + EIGEN_USING_STD_MATH(sqrt) + EIGEN_USING_STD_MATH(atan2) + EIGEN_USING_STD_MATH(cos) + EIGEN_USING_STD_MATH(sin) const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0); const Scalar s_sqrt3 = sqrt(Scalar(3.0)); @@ -531,15 +549,16 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 // Sort in increasing order. if (roots(0) >= roots(1)) - std::swap(roots(0),roots(1)); + internal::swap(roots(0),roots(1)); if (roots(1) >= roots(2)) { - std::swap(roots(1),roots(2)); + internal::swap(roots(1),roots(2)); if (roots(0) >= roots(1)) - std::swap(roots(0),roots(1)); + internal::swap(roots(0),roots(1)); } } + EIGEN_DEVICE_FUNC static inline void run(SolverType& solver, const MatrixType& mat, int options) { using std::sqrt; @@ -660,12 +679,14 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 }; // 2x2 direct eigenvalues decomposition, code from Hauke Heibel -template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false> +template<typename SolverType> +struct direct_selfadjoint_eigenvalues<SolverType,2,false> { typedef typename SolverType::MatrixType MatrixType; typedef typename SolverType::RealVectorType VectorType; typedef typename SolverType::Scalar Scalar; + EIGEN_DEVICE_FUNC static inline void computeRoots(const MatrixType& m, VectorType& roots) { using std::sqrt; @@ -675,6 +696,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2 roots(1) = t1 + t0; } + EIGEN_DEVICE_FUNC static inline void run(SolverType& solver, const MatrixType& mat, int options) { using std::sqrt; @@ -728,6 +750,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2 } template<typename MatrixType> +EIGEN_DEVICE_FUNC SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> ::computeDirect(const MatrixType& matrix, int options) { @@ -737,6 +760,7 @@ SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> namespace internal { template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> +EIGEN_DEVICE_FUNC static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) { using std::abs; diff --git a/Eigen/src/Geometry/OrthoMethods.h b/Eigen/src/Geometry/OrthoMethods.h index 556bc8160..26be3ee5b 100644 --- a/Eigen/src/Geometry/OrthoMethods.h +++ b/Eigen/src/Geometry/OrthoMethods.h @@ -132,6 +132,7 @@ struct unitOrthogonal_selector typedef typename NumTraits<Scalar>::Real RealScalar; typedef typename Derived::Index Index; typedef Matrix<Scalar,2,1> Vector2; + EIGEN_DEVICE_FUNC static inline VectorType run(const Derived& src) { VectorType perp = VectorType::Zero(src.size()); @@ -154,6 +155,7 @@ struct unitOrthogonal_selector<Derived,3> typedef typename plain_matrix_type<Derived>::type VectorType; typedef typename traits<Derived>::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; + EIGEN_DEVICE_FUNC static inline VectorType run(const Derived& src) { VectorType perp; @@ -192,6 +194,7 @@ template<typename Derived> struct unitOrthogonal_selector<Derived,2> { typedef typename plain_matrix_type<Derived>::type VectorType; + EIGEN_DEVICE_FUNC static inline VectorType run(const Derived& src) { return VectorType(-numext::conj(src.y()), numext::conj(src.x())).normalized(); } }; |