diff options
| author |  Gael Guennebaud <g.gael@free.fr> | 2009-12-22 22:51:08 +0100 |
|---|---|---|
| committer | Gael Guennebaud <g.gael@free.fr> | 2009-12-22 22:51:08 +0100 |
| commit | eaaba30cacf078335ad41314bf1e4ce993f0b3b4 (patch) | |
| tree | 61156273b56dc19fc4e661a5aeb42d8e6112a307 /Eigen/src | |
| parent | e182e9c6166840b8db44e0be459a2e26d26fda0f (diff) | |
| parent | eeec7d842e05394703c42e7569230835f7842089 (diff) | |
merge with default branch
Diffstat (limited to 'Eigen/src')
88 files changed, 1354 insertions, 976 deletions
diff --git a/Eigen/src/Array/Functors.h b/Eigen/src/Array/Functors.h index fd259f7bc..975fd9c7f 100644 --- a/Eigen/src/Array/Functors.h +++ b/Eigen/src/Array/Functors.h @@ -56,7 +56,8 @@ struct ei_functor_traits<ei_scalar_add_op<Scalar> > * * \sa class CwiseUnaryOp, Cwise::sqrt() */ -template<typename Scalar> struct ei_scalar_sqrt_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_sqrt_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_sqrt_op) inline const Scalar operator() (const Scalar& a) const { return ei_sqrt(a); } typedef typename ei_packet_traits<Scalar>::type Packet; inline Packet packetOp(const Packet& a) const { return ei_psqrt(a); } @@ -77,7 +78,8 @@ struct ei_functor_traits<ei_scalar_sqrt_op<Scalar> > * * \sa class CwiseUnaryOp, Cwise::cos() */ -template<typename Scalar> struct ei_scalar_cos_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_cos_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_cos_op) inline Scalar operator() (const Scalar& a) const { return ei_cos(a); } typedef typename ei_packet_traits<Scalar>::type Packet; inline Packet packetOp(const Packet& a) const { return ei_pcos(a); } @@ -87,7 +89,7 @@ struct ei_functor_traits<ei_scalar_cos_op<Scalar> > { enum { Cost = 5 * NumTraits<Scalar>::MulCost, - PacketAccess = ei_packet_traits<Scalar>::HasCos && EIGEN_FAST_MATH + PacketAccess = ei_packet_traits<Scalar>::HasCos }; }; @@ -99,7 +101,8 @@ struct ei_functor_traits<ei_scalar_cos_op<Scalar> > * * \sa class CwiseUnaryOp, Cwise::sin() */ -template<typename Scalar> struct ei_scalar_sin_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_sin_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_sin_op) inline const Scalar operator() (const Scalar& a) const { return ei_sin(a); } typedef typename ei_packet_traits<Scalar>::type Packet; inline Packet packetOp(const Packet& a) const { return ei_psin(a); } @@ -109,7 +112,7 @@ struct ei_functor_traits<ei_scalar_sin_op<Scalar> > { enum { Cost = 5 * NumTraits<Scalar>::MulCost, - PacketAccess = ei_packet_traits<Scalar>::HasSin && EIGEN_FAST_MATH + PacketAccess = ei_packet_traits<Scalar>::HasSin }; }; @@ -143,6 +146,7 @@ struct ei_functor_traits<ei_scalar_pow_op<Scalar> > */ template<typename Scalar> struct ei_scalar_inverse_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_inverse_op) inline Scalar operator() (const Scalar& a) const { return Scalar(1)/a; } template<typename PacketScalar> inline const PacketScalar packetOp(const PacketScalar& a) const @@ -162,6 +166,7 @@ struct ei_functor_traits<ei_scalar_inverse_op<Scalar> > */ template<typename Scalar> struct ei_scalar_square_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_square_op) inline Scalar operator() (const Scalar& a) const { return a*a; } template<typename PacketScalar> inline const PacketScalar packetOp(const PacketScalar& a) const @@ -181,6 +186,7 @@ struct ei_functor_traits<ei_scalar_square_op<Scalar> > */ template<typename Scalar> struct ei_scalar_cube_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_cube_op) inline Scalar operator() (const Scalar& a) const { return a*a*a; } template<typename PacketScalar> inline const PacketScalar packetOp(const PacketScalar& a) const diff --git a/Eigen/src/Array/Random.h b/Eigen/src/Array/Random.h index 34b835776..97ca0fba3 100644 --- a/Eigen/src/Array/Random.h +++ b/Eigen/src/Array/Random.h @@ -25,8 +25,8 @@ #ifndef EIGEN_RANDOM_H #define EIGEN_RANDOM_H -template<typename Scalar> struct ei_scalar_random_op EIGEN_EMPTY_STRUCT { - inline ei_scalar_random_op(void) {} +template<typename Scalar> struct ei_scalar_random_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_random_op) inline const Scalar operator() (int, int) const { return ei_random<Scalar>(); } }; template<typename Scalar> diff --git a/Eigen/src/Array/Replicate.h b/Eigen/src/Array/Replicate.h index 5a9b10738..3f87e09fe 100644 --- a/Eigen/src/Array/Replicate.h +++ b/Eigen/src/Array/Replicate.h @@ -97,9 +97,6 @@ template<typename MatrixType,int RowFactor,int ColFactor> class Replicate const typename MatrixType::Nested m_matrix; const ei_int_if_dynamic<RowFactor> m_rowFactor; const ei_int_if_dynamic<ColFactor> m_colFactor; - - private: - Replicate& operator=(const Replicate&); }; /** \nonstableyet diff --git a/Eigen/src/Array/Select.h b/Eigen/src/Array/Select.h index 1983bd870..b1fab69c9 100644 --- a/Eigen/src/Array/Select.h +++ b/Eigen/src/Array/Select.h @@ -134,11 +134,11 @@ DenseBase<Derived>::select(const DenseBase<ThenDerived>& thenMatrix, */ template<typename Derived> template<typename ThenDerived> -inline const Select<Derived,ThenDerived, NestByValue<typename ThenDerived::ConstantReturnType> > +inline const Select<Derived,ThenDerived, typename ThenDerived::ConstantReturnType> DenseBase<Derived>::select(const DenseBase<ThenDerived>& thenMatrix, typename ThenDerived::Scalar elseScalar) const { - return Select<Derived,ThenDerived,NestByValue<typename ThenDerived::ConstantReturnType> >( + return Select<Derived,ThenDerived,typename ThenDerived::ConstantReturnType>( derived(), thenMatrix.derived(), ThenDerived::Constant(rows(),cols(),elseScalar)); } @@ -151,11 +151,11 @@ DenseBase<Derived>::select(const DenseBase<ThenDerived>& thenMatrix, */ template<typename Derived> template<typename ElseDerived> -inline const Select<Derived, NestByValue<typename ElseDerived::ConstantReturnType>, ElseDerived > +inline const Select<Derived, typename ElseDerived::ConstantReturnType, ElseDerived > DenseBase<Derived>::select(typename ElseDerived::Scalar thenScalar, const DenseBase<ElseDerived>& elseMatrix) const { - return Select<Derived,NestByValue<typename ElseDerived::ConstantReturnType>,ElseDerived>( + return Select<Derived,typename ElseDerived::ConstantReturnType,ElseDerived>( derived(), ElseDerived::Constant(rows(),cols(),thenScalar), elseMatrix.derived()); } diff --git a/Eigen/src/Array/VectorwiseOp.h b/Eigen/src/Array/VectorwiseOp.h index f6774a019..05dd69789 100644 --- a/Eigen/src/Array/VectorwiseOp.h +++ b/Eigen/src/Array/VectorwiseOp.h @@ -113,7 +113,8 @@ class PartialReduxExpr : ei_no_assignment_operator, #define EIGEN_MEMBER_FUNCTOR(MEMBER,COST) \ template <typename ResultType> \ - struct ei_member_##MEMBER EIGEN_EMPTY_STRUCT { \ + struct ei_member_##MEMBER { \ + EIGEN_EMPTY_STRUCT_CTOR(ei_member_##MEMBER) \ typedef ResultType result_type; \ template<typename Scalar, int Size> struct Cost \ { enum { value = COST }; }; \ @@ -124,6 +125,9 @@ class PartialReduxExpr : ei_no_assignment_operator, EIGEN_MEMBER_FUNCTOR(squaredNorm, Size * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(norm, (Size+5) * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost); +EIGEN_MEMBER_FUNCTOR(stableNorm, (Size+5) * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost); +EIGEN_MEMBER_FUNCTOR(blueNorm, (Size+5) * NumTraits<Scalar>::MulCost + (Size-1)*NumTraits<Scalar>::AddCost); +EIGEN_MEMBER_FUNCTOR(hypotNorm, (Size-1) * ei_functor_traits<ei_scalar_hypot_op<Scalar> >::Cost ); EIGEN_MEMBER_FUNCTOR(sum, (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(mean, (Size-1)*NumTraits<Scalar>::AddCost + NumTraits<Scalar>::MulCost); EIGEN_MEMBER_FUNCTOR(minCoeff, (Size-1)*NumTraits<Scalar>::AddCost); @@ -291,6 +295,33 @@ template<typename ExpressionType, int Direction> class VectorwiseOp const typename ReturnType<ei_member_norm>::Type norm() const { return _expression(); } + + /** \returns a row (or column) vector expression of the norm + * of each column (or row) of the referenced expression, using + * blue's algorithm. + * + * \sa MatrixBase::blueNorm() */ + const typename ReturnType<ei_member_blueNorm>::Type blueNorm() const + { return _expression(); } + + + /** \returns a row (or column) vector expression of the norm + * of each column (or row) of the referenced expression, avoiding + * underflow and overflow. + * + * \sa MatrixBase::stableNorm() */ + const typename ReturnType<ei_member_stableNorm>::Type stableNorm() const + { return _expression(); } + + + /** \returns a row (or column) vector expression of the norm + * of each column (or row) of the referenced expression, avoiding + * underflow and overflow using a concatenation of hypot() calls. + * + * \sa MatrixBase::hypotNorm() */ + const typename ReturnType<ei_member_hypotNorm>::Type hypotNorm() const + { return _expression(); } + /** \returns a row (or column) vector expression of the sum * of each column (or row) of the referenced expression. * @@ -409,22 +440,22 @@ template<typename ExpressionType, int Direction> class VectorwiseOp template<typename OtherDerived> CwiseBinaryOp<ei_scalar_sum_op<Scalar>, ExpressionType, - NestByValue<typename ExtendedType<OtherDerived>::Type> > + typename ExtendedType<OtherDerived>::Type> operator+(const MatrixBase<OtherDerived>& other) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); - return m_matrix + extendedTo(other).nestByValue(); + return m_matrix + extendedTo(other); } /** Returns the expression of the difference between each subvector of \c *this and the vector \a other */ template<typename OtherDerived> CwiseBinaryOp<ei_scalar_difference_op<Scalar>, ExpressionType, - NestByValue<typename ExtendedType<OtherDerived>::Type> > + typename ExtendedType<OtherDerived>::Type> operator-(const MatrixBase<OtherDerived>& other) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived); - return m_matrix - extendedTo(other).nestByValue(); + return m_matrix - extendedTo(other); } /////////// Geometry module /////////// @@ -451,19 +482,16 @@ template<typename ExpressionType, int Direction> class VectorwiseOp Direction==Horizontal ? 1 : int(ei_traits<ExpressionType>::ColsAtCompileTime)> HNormalized_Factors; typedef CwiseBinaryOp<ei_scalar_quotient_op<typename ei_traits<ExpressionType>::Scalar>, - NestByValue<HNormalized_Block>, - NestByValue<Replicate<NestByValue<HNormalized_Factors>, + HNormalized_Block, + Replicate<HNormalized_Factors, Direction==Vertical ? HNormalized_SizeMinusOne : 1, - Direction==Horizontal ? HNormalized_SizeMinusOne : 1> > > + Direction==Horizontal ? HNormalized_SizeMinusOne : 1> > HNormalizedReturnType; const HNormalizedReturnType hnormalized() const; protected: ExpressionTypeNested m_matrix; - - private: - VectorwiseOp& operator=(const VectorwiseOp&); }; /** \array_module diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index a1706e53e..ad737aaeb 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -224,18 +224,18 @@ template<> struct ei_llt_inplace<UpperTriangular> template<typename MatrixType> struct LLT_Traits<MatrixType,LowerTriangular> { typedef TriangularView<MatrixType, LowerTriangular> MatrixL; - typedef TriangularView<NestByValue<typename MatrixType::AdjointReturnType>, UpperTriangular> MatrixU; + typedef TriangularView<typename MatrixType::AdjointReturnType, UpperTriangular> MatrixU; inline static MatrixL getL(const MatrixType& m) { return m; } - inline static MatrixU getU(const MatrixType& m) { return m.adjoint().nestByValue(); } + inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } static bool inplace_decomposition(MatrixType& m) { return ei_llt_inplace<LowerTriangular>::blocked(m); } }; template<typename MatrixType> struct LLT_Traits<MatrixType,UpperTriangular> { - typedef TriangularView<NestByValue<typename MatrixType::AdjointReturnType>, LowerTriangular> MatrixL; + typedef TriangularView<typename MatrixType::AdjointReturnType, LowerTriangular> MatrixL; typedef TriangularView<MatrixType, UpperTriangular> MatrixU; - inline static MatrixL getL(const MatrixType& m) { return m.adjoint().nestByValue(); } + inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } inline static MatrixU getU(const MatrixType& m) { return m; } static bool inplace_decomposition(MatrixType& m) { return ei_llt_inplace<UpperTriangular>::blocked(m); } diff --git a/Eigen/src/Core/AnyMatrixBase.h b/Eigen/src/Core/AnyMatrixBase.h index 0522c1d5e..d1c16923b 100644 --- a/Eigen/src/Core/AnyMatrixBase.h +++ b/Eigen/src/Core/AnyMatrixBase.h @@ -57,7 +57,7 @@ template<typename Derived> struct AnyMatrixBase { derived().evalTo(dst); } /** \internal Don't use it, but do the equivalent: \code dst += *this; \endcode */ - template<typename Dest> inline void addToDense(Dest& dst) const + template<typename Dest> inline void addTo(Dest& dst) const { // This is the default implementation, // derived class can reimplement it in a more optimized way. @@ -67,7 +67,7 @@ template<typename Derived> struct AnyMatrixBase } /** \internal Don't use it, but do the equivalent: \code dst -= *this; \endcode */ - template<typename Dest> inline void subToDense(Dest& dst) const + template<typename Dest> inline void subTo(Dest& dst) const { // This is the default implementation, // derived class can reimplement it in a more optimized way. @@ -114,7 +114,7 @@ template<typename Derived> template<typename OtherDerived> Derived& DenseBase<Derived>::operator+=(const AnyMatrixBase<OtherDerived> &other) { - other.derived().addToDense(derived()); + other.derived().addTo(derived()); return derived(); } @@ -122,7 +122,7 @@ template<typename Derived> template<typename OtherDerived> Derived& DenseBase<Derived>::operator-=(const AnyMatrixBase<OtherDerived> &other) { - other.derived().subToDense(derived()); + other.derived().subTo(derived()); return derived(); } diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index 466205208..d6bf37c6e 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -28,7 +28,7 @@ #define EIGEN_ASSIGN_H /*************************************************************************** -* Part 1 : the logic deciding a strategy for vectorization and unrolling +* Part 1 : the logic deciding a strategy for traversal and unrolling * ***************************************************************************/ template <typename Derived, typename OtherDerived> @@ -53,44 +53,53 @@ private: }; enum { - MightVectorize = (int(Derived::Flags) & int(OtherDerived::Flags) & ActualPacketAccessBit) - && ((int(Derived::Flags)&RowMajorBit)==(int(OtherDerived::Flags)&RowMajorBit)), + StorageOrdersAgree = (int(Derived::Flags)&RowMajorBit)==(int(OtherDerived::Flags)&RowMajorBit), + MightVectorize = StorageOrdersAgree + && (int(Derived::Flags) & int(OtherDerived::Flags) & ActualPacketAccessBit), MayInnerVectorize = MightVectorize && int(InnerSize)!=Dynamic && int(InnerSize)%int(PacketSize)==0 && int(DstIsAligned) && int(SrcIsAligned), - MayLinearVectorize = MightVectorize && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit) - && (DstIsAligned || InnerMaxSize == Dynamic),/* If the destination isn't aligned, - we have to do runtime checks and we don't unroll, so it's only good for large enough sizes. See remark below - about InnerMaxSize. */ - MaySliceVectorize = MightVectorize && int(InnerMaxSize)>=3*PacketSize /* slice vectorization can be slow, so we only - want it if the slices are big, which is indicated by InnerMaxSize rather than InnerSize, think of the case - of a dynamic block in a fixed-size matrix */ + MayLinearize = StorageOrdersAgree && (int(Derived::Flags) & int(OtherDerived::Flags) & LinearAccessBit), + MayLinearVectorize = MightVectorize && MayLinearize + && (DstIsAligned || InnerMaxSize == Dynamic), + /* If the destination isn't aligned, we have to do runtime checks and we don't unroll, + so it's only good for large enough sizes. See remark below about InnerMaxSize. */ + MaySliceVectorize = MightVectorize && int(InnerMaxSize)>=3*PacketSize + /* slice vectorization can be slow, so we only want it if the slices are big, which is + indicated by InnerMaxSize rather than InnerSize, think of the case of a dynamic block + in a fixed-size matrix */ }; public: enum { - Vectorization = int(MayInnerVectorize) ? int(InnerVectorization) - : int(MayLinearVectorize) ? int(LinearVectorization) - : int(MaySliceVectorize) ? int(SliceVectorization) - : int(NoVectorization) + Traversal = int(MayInnerVectorize) ? int(InnerVectorizedTraversal) + : int(MayLinearVectorize) ? int(LinearVectorizedTraversal) + : int(MaySliceVectorize) ? int(SliceVectorizedTraversal) + : int(MayLinearize) ? int(LinearTraversal) + : int(DefaultTraversal), + Vectorized = int(Traversal) == InnerVectorizedTraversal + || int(Traversal) == LinearVectorizedTraversal + || int(Traversal) == SliceVectorizedTraversal }; private: enum { - UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Vectorization) == int(NoVectorization) ? 1 : int(PacketSize)), + UnrollingLimit = EIGEN_UNROLLING_LIMIT * (Vectorized ? int(PacketSize) : 1), MayUnrollCompletely = int(Derived::SizeAtCompileTime) * int(OtherDerived::CoeffReadCost) <= int(UnrollingLimit), - MayUnrollInner = int(InnerSize * OtherDerived::CoeffReadCost) <= int(UnrollingLimit) + MayUnrollInner = int(InnerSize) * int(OtherDerived::CoeffReadCost) <= int(UnrollingLimit) }; public: enum { - Unrolling = (int(Vectorization) == int(InnerVectorization) || int(Vectorization) == int(NoVectorization)) - ? ( - int(MayUnrollCompletely) ? int(CompleteUnrolling) - : int(MayUnrollInner) ? int(InnerUnrolling) - : int(NoUnrolling) - ) - : int(Vectorization) == int(LinearVectorization) - ? ( int(MayUnrollCompletely) && int(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) ) + Unrolling = (int(Traversal) == int(InnerVectorizedTraversal) || int(Traversal) == int(DefaultTraversal)) + ? ( + int(MayUnrollCompletely) ? int(CompleteUnrolling) + : int(MayUnrollInner) ? int(InnerUnrolling) + : int(NoUnrolling) + ) + : int(Traversal) == int(LinearVectorizedTraversal) + ? ( int(MayUnrollCompletely) && int(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) ) + : int(Traversal) == int(LinearTraversal) + ? ( int(MayUnrollCompletely) ? int(CompleteUnrolling) : int(NoUnrolling) ) : int(NoUnrolling) }; @@ -102,11 +111,12 @@ public: EIGEN_DEBUG_VAR(InnerSize) EIGEN_DEBUG_VAR(InnerMaxSize) EIGEN_DEBUG_VAR(PacketSize) + EIGEN_DEBUG_VAR(StorageOrdersAgree) EIGEN_DEBUG_VAR(MightVectorize) EIGEN_DEBUG_VAR(MayInnerVectorize) EIGEN_DEBUG_VAR(MayLinearVectorize) EIGEN_DEBUG_VAR(MaySliceVectorize) - EIGEN_DEBUG_VAR(Vectorization) + EIGEN_DEBUG_VAR(Traversal) EIGEN_DEBUG_VAR(UnrollingLimit) EIGEN_DEBUG_VAR(MayUnrollCompletely) EIGEN_DEBUG_VAR(MayUnrollInner) @@ -118,12 +128,12 @@ public: * Part 2 : meta-unrollers ***************************************************************************/ -/*********************** -*** No vectorization *** -***********************/ +/************************ +*** Default traversal *** +************************/ template<typename Derived1, typename Derived2, int Index, int Stop> -struct ei_assign_novec_CompleteUnrolling +struct ei_assign_DefaultTraversal_CompleteUnrolling { enum { row = int(Derived1::Flags)&RowMajorBit @@ -137,18 +147,18 @@ struct ei_assign_novec_CompleteUnrolling EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src) { dst.copyCoeff(row, col, src); - ei_assign_novec_CompleteUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src); + ei_assign_DefaultTraversal_CompleteUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src); } }; template<typename Derived1, typename Derived2, int Stop> -struct ei_assign_novec_CompleteUnrolling<Derived1, Derived2, Stop, Stop> +struct ei_assign_DefaultTraversal_CompleteUnrolling<Derived1, Derived2, Stop, Stop> { EIGEN_STRONG_INLINE static void run(Derived1 &, const Derived2 &) {} }; template<typename Derived1, typename Derived2, int Index, int Stop> -struct ei_assign_novec_InnerUnrolling +struct ei_assign_DefaultTraversal_InnerUnrolling { EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src, int row_or_col) { @@ -156,16 +166,36 @@ struct ei_assign_novec_InnerUnrolling const int row = rowMajor ? row_or_col : Index; const int col = rowMajor ? Index : row_or_col; dst.copyCoeff(row, col, src); - ei_assign_novec_InnerUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src, row_or_col); + ei_assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src, row_or_col); } }; template<typename Derived1, typename Derived2, int Stop> -struct ei_assign_novec_InnerUnrolling<Derived1, Derived2, Stop, Stop> +struct ei_assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, Stop, Stop> { EIGEN_STRONG_INLINE static void run(Derived1 &, const Derived2 &, int) {} }; +/*********************** +*** Linear traversal *** +***********************/ + +template<typename Derived1, typename Derived2, int Index, int Stop> +struct ei_assign_LinearTraversal_CompleteUnrolling +{ + EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src) + { + dst.copyCoeff(Index, src); + ei_assign_LinearTraversal_CompleteUnrolling<Derived1, Derived2, Index+1, Stop>::run(dst, src); + } +}; + +template<typename Derived1, typename Derived2, int Stop> +struct ei_assign_LinearTraversal_CompleteUnrolling<Derived1, Derived2, Stop, Stop> +{ + EIGEN_STRONG_INLINE static void run(Derived1 &, const Derived2 &) {} +}; + /************************** *** Inner vectorization *** **************************/ @@ -221,16 +251,16 @@ struct ei_assign_innervec_InnerUnrolling<Derived1, Derived2, Stop, Stop> ***************************************************************************/ template<typename Derived1, typename Derived2, - int Vectorization = ei_assign_traits<Derived1, Derived2>::Vectorization, + int Traversal = ei_assign_traits<Derived1, Derived2>::Traversal, int Unrolling = ei_assign_traits<Derived1, Derived2>::Unrolling> struct ei_assign_impl; -/*********************** -*** No vectorization *** -***********************/ +/************************ +*** Default traversal *** +************************/ template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, NoVectorization, NoUnrolling> +struct ei_assign_impl<Derived1, Derived2, DefaultTraversal, NoUnrolling> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -248,17 +278,17 @@ struct ei_assign_impl<Derived1, Derived2, NoVectorization, NoUnrolling> }; template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, NoVectorization, CompleteUnrolling> +struct ei_assign_impl<Derived1, Derived2, DefaultTraversal, CompleteUnrolling> { EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src) { - ei_assign_novec_CompleteUnrolling<Derived1, Derived2, 0, Derived1::SizeAtCompileTime> + ei_assign_DefaultTraversal_CompleteUnrolling<Derived1, Derived2, 0, Derived1::SizeAtCompileTime> ::run(dst, src); } }; template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, NoVectorization, InnerUnrolling> +struct ei_assign_impl<Derived1, Derived2, DefaultTraversal, InnerUnrolling> { EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src) { @@ -266,17 +296,42 @@ struct ei_assign_impl<Derived1, Derived2, NoVectorization, InnerUnrolling> const int innerSize = rowMajor ? Derived1::ColsAtCompileTime : Derived1::RowsAtCompileTime; const int outerSize = dst.outerSize(); for(int j = 0; j < outerSize; ++j) - ei_assign_novec_InnerUnrolling<Derived1, Derived2, 0, innerSize> + ei_assign_DefaultTraversal_InnerUnrolling<Derived1, Derived2, 0, innerSize> ::run(dst, src, j); } }; +/*********************** +*** Linear traversal *** +***********************/ + +template<typename Derived1, typename Derived2> +struct ei_assign_impl<Derived1, Derived2, LinearTraversal, NoUnrolling> +{ + inline static void run(Derived1 &dst, const Derived2 &src) + { + const int size = dst.size(); + for(int i = 0; i < size; ++i) + dst.copyCoeff(i, src); + } +}; + +template<typename Derived1, typename Derived2> +struct ei_assign_impl<Derived1, Derived2, LinearTraversal, CompleteUnrolling> +{ + EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src) + { + ei_assign_LinearTraversal_CompleteUnrolling<Derived1, Derived2, 0, Derived1::SizeAtCompileTime> + ::run(dst, src); + } +}; + /************************** *** Inner vectorization *** **************************/ template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, InnerVectorization, NoUnrolling> +struct ei_assign_impl<Derived1, Derived2, InnerVectorizedTraversal, NoUnrolling> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -295,7 +350,7 @@ struct ei_assign_impl<Derived1, Derived2, InnerVectorization, NoUnrolling> }; template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, InnerVectorization, CompleteUnrolling> +struct ei_assign_impl<Derived1, Derived2, InnerVectorizedTraversal, CompleteUnrolling> { EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src) { @@ -305,7 +360,7 @@ struct ei_assign_impl<Derived1, Derived2, InnerVectorization, CompleteUnrolling> }; template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, InnerVectorization, InnerUnrolling> +struct ei_assign_impl<Derived1, Derived2, InnerVectorizedTraversal, InnerUnrolling> { EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src) { @@ -323,7 +378,7 @@ struct ei_assign_impl<Derived1, Derived2, InnerVectorization, InnerUnrolling> ***************************/ template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, LinearVectorization, NoUnrolling> +struct ei_assign_impl<Derived1, Derived2, LinearVectorizedTraversal, NoUnrolling> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -347,7 +402,7 @@ struct ei_assign_impl<Derived1, Derived2, LinearVectorization, NoUnrolling> }; template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling> +struct ei_assign_impl<Derived1, Derived2, LinearVectorizedTraversal, CompleteUnrolling> { EIGEN_STRONG_INLINE static void run(Derived1 &dst, const Derived2 &src) { @@ -356,7 +411,7 @@ struct ei_assign_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling const int alignedSize = (size/packetSize)*packetSize; ei_assign_innervec_CompleteUnrolling<Derived1, Derived2, 0, alignedSize>::run(dst, src); - ei_assign_novec_CompleteUnrolling<Derived1, Derived2, alignedSize, size>::run(dst, src); + ei_assign_DefaultTraversal_CompleteUnrolling<Derived1, Derived2, alignedSize, size>::run(dst, src); } }; @@ -365,7 +420,7 @@ struct ei_assign_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling ***************************/ template<typename Derived1, typename Derived2> -struct ei_assign_impl<Derived1, Derived2, SliceVectorization, NoUnrolling> +struct ei_assign_impl<Derived1, Derived2, SliceVectorizedTraversal, NoUnrolling> { inline static void run(Derived1 &dst, const Derived2 &src) { @@ -430,6 +485,9 @@ EIGEN_STRONG_INLINE Derived& DenseBase<Derived> #ifdef EIGEN_DEBUG_ASSIGN ei_assign_traits<Derived, OtherDerived>::debug(); #endif +#ifndef EIGEN_NO_DEBUG + checkTransposeAliasing(other.derived()); +#endif return derived(); } diff --git a/Eigen/src/Core/BandMatrix.h b/Eigen/src/Core/BandMatrix.h index 841d1786a..7943e6280 100644 --- a/Eigen/src/Core/BandMatrix.h +++ b/Eigen/src/Core/BandMatrix.h @@ -70,7 +70,7 @@ class BandMatrix : public AnyMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs MaxColsAtCompileTime = ei_traits<BandMatrix>::MaxColsAtCompileTime }; typedef typename ei_traits<BandMatrix>::Scalar Scalar; - typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> PlainMatrixType; + typedef Matrix<Scalar,RowsAtCompileTime,ColsAtCompileTime> DenseMatrixType; protected: enum { @@ -87,7 +87,7 @@ class BandMatrix : public AnyMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs : m_data(1+supers+subs,cols), m_rows(rows), m_supers(supers), m_subs(subs) { - m_data.setConstant(666); + //m_data.setConstant(666); } /** \returns the number of columns */ @@ -141,7 +141,7 @@ class BandMatrix : public AnyMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs }; typedef Block<DataType,1, DiagonalSize> BuildType; typedef typename ei_meta_if<Conjugate, - CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>,NestByValue<BuildType> >, + CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>,BuildType >, BuildType>::ret Type; }; @@ -171,9 +171,9 @@ class BandMatrix : public AnyMatrixBase<BandMatrix<_Scalar,Rows,Cols,Supers,Subs return Block<DataType,1,Dynamic>(m_data, supers()-i, std::max(0,i), 1, diagonalLength(i)); } - PlainMatrixType toDense() const + DenseMatrixType toDenseMatrix() const { - PlainMatrixType res(rows(),cols()); + DenseMatrixType res(rows(),cols()); res.setZero(); res.diagonal() = diagonal(); for (int i=1; i<=supers();++i) diff --git a/Eigen/src/Core/Coeffs.h b/Eigen/src/Core/Coeffs.h index c8bc9db85..b8af2531e 100644 --- a/Eigen/src/Core/Coeffs.h +++ b/Eigen/src/Core/Coeffs.h @@ -379,6 +379,38 @@ EIGEN_STRONG_INLINE void DenseBase<Derived>::copyPacket(int index, const DenseBa other.derived().template packet<LoadMode>(index)); } + +template<typename Derived, typename Integer, bool JustReturnZero> +struct ei_alignmentOffset_impl +{ + inline static Integer run(const DenseBase<Derived>&, Integer) + { return 0; } +}; + +template<typename Derived, typename Integer> +struct ei_alignmentOffset_impl<Derived, Integer, false> +{ + inline static Integer run(const DenseBase<Derived>& m, Integer maxOffset) + { + return ei_alignmentOffset(&m.const_cast_derived().coeffRef(0,0), maxOffset); + } +}; + +/** \internal \returns the number of elements which have to be skipped, starting + * from the address of coeffRef(0,0), to find the first 16-byte aligned element. + * + * \note If the expression doesn't have the DirectAccessBit, this function returns 0. + * + * There is also the variant ei_alignmentOffset(const Scalar*, Integer) defined in Memory.h. + */ +template<typename Derived, typename Integer> +inline static Integer ei_alignmentOffset(const DenseBase<Derived>& m, Integer maxOffset) +{ + return ei_alignmentOffset_impl<Derived, Integer, + (Derived::Flags & AlignedBit) || !(Derived::Flags & DirectAccessBit)> + ::run(m, maxOffset); +} + #endif #endif // EIGEN_COEFFS_H diff --git a/Eigen/src/Core/CommaInitializer.h b/Eigen/src/Core/CommaInitializer.h index 428e4d82d..a312b40f5 100644 --- a/Eigen/src/Core/CommaInitializer.h +++ b/Eigen/src/Core/CommaInitializer.h @@ -116,9 +116,6 @@ struct CommaInitializer int m_row; // current row id int m_col; // current col id int m_currentBlockRows; // current block height - -private: - CommaInitializer& operator=(const CommaInitializer&); }; /** \anchor MatrixBaseCommaInitRef diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index 954ea9275..4d7ab6598 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -65,11 +65,16 @@ struct ei_traits<CwiseBinaryOp<BinaryOp, Lhs, Rhs> > : ei_traits<Lhs> RhsCoeffReadCost = _RhsNested::CoeffReadCost, LhsFlags = _LhsNested::Flags, RhsFlags = _RhsNested::Flags, + StorageOrdersAgree = (int(Lhs::Flags)&RowMajorBit)==(int(Rhs::Flags)&RowMajorBit), Flags = (int(LhsFlags) | int(RhsFlags)) & ( HereditaryBits - | (int(LhsFlags) & int(RhsFlags) & (LinearAccessBit | AlignedBit)) - | (ei_functor_traits<BinaryOp>::PacketAccess && ((int(LhsFlags) & RowMajorBit)==(int(RhsFlags) & RowMajorBit)) - ? (int(LhsFlags) & int(RhsFlags) & PacketAccessBit) : 0)), + | (int(LhsFlags) & int(RhsFlags) & + ( AlignedBit + | (StorageOrdersAgree ? LinearAccessBit : 0) + | (ei_functor_traits<BinaryOp>::PacketAccess && StorageOrdersAgree ? PacketAccessBit : 0) + ) + ) + ), CoeffReadCost = LhsCoeffReadCost + RhsCoeffReadCost + ei_functor_traits<BinaryOp>::Cost }; }; diff --git a/Eigen/src/Core/DenseBase.h b/Eigen/src/Core/DenseBase.h index ea6383998..a7e74b81d 100644 --- a/Eigen/src/Core/DenseBase.h +++ b/Eigen/src/Core/DenseBase.h @@ -274,21 +274,12 @@ template<typename Derived> class DenseBase Eigen::Transpose<Derived> transpose(); const Eigen::Transpose<Derived> transpose() const; void transposeInPlace(); - #ifndef EIGEN_NO_DEBUG - template<typename OtherDerived> - Derived& lazyAssign(const Transpose<OtherDerived>& other); - template<typename DerivedA, typename DerivedB> - Derived& lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,Transpose<DerivedA>,DerivedB>& other); - template<typename DerivedA, typename DerivedB> - Derived& lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerivedA,Transpose<DerivedB> >& other); - +#ifndef EIGEN_NO_DEBUG + protected: template<typename OtherDerived> - Derived& lazyAssign(const CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<OtherDerived> > >& other); - template<typename DerivedA, typename DerivedB> - Derived& lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<DerivedA> > >,DerivedB>& other); - template<typename DerivedA, typename DerivedB> - Derived& lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerivedA,CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<DerivedB> > > >& other); - #endif + void checkTransposeAliasing(const OtherDerived& other) const; + public: +#endif RowXpr row(int i); const RowXpr row(int i) const; @@ -382,18 +373,19 @@ template<typename Derived> class DenseBase template<typename OtherDerived> bool isApprox(const DenseBase<OtherDerived>& other, - RealScalar prec = precision<Scalar>()) const; + RealScalar prec = dummy_precision<Scalar>()) const; bool isMuchSmallerThan(const RealScalar& other, - RealScalar prec = precision<Scalar>()) const; + RealScalar prec = dummy_precision<Scalar>()) const; template<typename OtherDerived> bool isMuchSmallerThan(const DenseBase<OtherDerived>& other, - RealScalar prec = precision<Scalar>()) const; + RealScalar prec = dummy_precision<Scalar>()) const; - bool isApproxToConstant(const Scalar& value, RealScalar prec = precision<Scalar>()) const; - bool isConstant(const Scalar& value, RealScalar prec = precision<Scalar>()) const; - bool isZero(RealScalar prec = precision<Scalar>()) const; - bool isOnes(RealScalar prec = precision<Scalar>()) const; + bool isApproxToConstant(const Scalar& value, RealScalar prec = dummy_precision<Scalar>()) const; + bool isConstant(const Scalar& value, RealScalar prec = dummy_precision<Scalar>()) const; + bool isZero(RealScalar prec = dummy_precision<Scalar>()) const; + bool isOnes(RealScalar prec = dummy_precision<Scalar>()) const; + // FIXME EIGEN_STRONG_INLINE Derived& operator*=(const Scalar& other) { SelfCwiseBinaryOp<ei_scalar_product_op<Scalar>, Derived> tmp(derived()); @@ -409,6 +401,7 @@ template<typename Derived> class DenseBase return derived(); } + // FIXME // template<typename OtherDerived> // inline bool operator==(const DenseBase<OtherDerived>& other) const // { return cwiseEqual(other).all(); } @@ -487,11 +480,11 @@ template<typename Derived> class DenseBase const DenseBase<ElseDerived>& elseMatrix) const; template<typename ThenDerived> - inline const Select<Derived,ThenDerived, NestByValue<typename ThenDerived::ConstantReturnType> > + inline const Select<Derived,ThenDerived, typename ThenDerived::ConstantReturnType> select(const DenseBase<ThenDerived>& thenMatrix, typename ThenDerived::Scalar elseScalar) const; template<typename ElseDerived> - inline const Select<Derived, NestByValue<typename ElseDerived::ConstantReturnType>, ElseDerived > + inline const Select<Derived, typename ElseDerived::ConstantReturnType, ElseDerived > select(typename ElseDerived::Scalar thenScalar, const DenseBase<ElseDerived>& elseMatrix) const; template<int p> RealScalar lpNorm() const; diff --git a/Eigen/src/Core/DenseStorageBase.h b/Eigen/src/Core/DenseStorageBase.h index ed116d708..7dc8d38c4 100644 --- a/Eigen/src/Core/DenseStorageBase.h +++ b/Eigen/src/Core/DenseStorageBase.h @@ -175,8 +175,14 @@ class DenseStorageBase : public _Base<Derived> && (RowsAtCompileTime == Dynamic || RowsAtCompileTime == rows) && (MaxColsAtCompileTime == Dynamic || MaxColsAtCompileTime >= cols) && (ColsAtCompileTime == Dynamic || ColsAtCompileTime == cols)); - m_storage.resize(rows * cols, rows, cols); - EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED + #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO + int size = rows*cols; + bool size_changed = size != this->size(); + m_storage.resize(size, rows, cols); + if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED + #else + m_storage.resize(rows*cols, rows, cols); + #endif } /** Resizes \c *this to a vector of length \a size @@ -194,11 +200,16 @@ class DenseStorageBase : public _Base<Derived> { EIGEN_STATIC_ASSERT_VECTOR_ONLY(DenseStorageBase) ei_assert(SizeAtCompileTime == Dynamic || SizeAtCompileTime == size); + #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO + bool size_changed = size != this->size(); + #endif if(RowsAtCompileTime == 1) m_storage.resize(size, 1, size); else m_storage.resize(size, size, 1); - EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED + #ifdef EIGEN_INITIALIZE_MATRICES_BY_ZERO + if(size_changed) EIGEN_INITIALIZE_BY_ZERO_IF_THAT_OPTION_IS_ENABLED + #endif } /** Resizes the matrix, changing only the number of columns. For the parameter of type NoChange_t, just pass the special value \c NoChange diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h index 6f93737ff..bd23b2e09 100644 --- a/Eigen/src/Core/DiagonalMatrix.h +++ b/Eigen/src/Core/DiagonalMatrix.h @@ -26,6 +26,7 @@ #ifndef EIGEN_DIAGONALMATRIX_H #define EIGEN_DIAGONALMATRIX_H +#ifndef EIGEN_PARSED_BY_DOXYGEN template<typename Derived> class DiagonalBase : public AnyMatrixBase<Derived> { @@ -44,19 +45,17 @@ class DiagonalBase : public AnyMatrixBase<Derived> typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType; - #ifndef EIGEN_PARSED_BY_DOXYGEN inline const Derived& derived() const { return *static_cast<const Derived*>(this); } inline Derived& derived() { return *static_cast<Derived*>(this); } - #endif // not EIGEN_PARSED_BY_DOXYGEN DenseMatrixType toDenseMatrix() const { return derived(); } template<typename DenseDerived> void evalTo(MatrixBase<DenseDerived> &other) const; template<typename DenseDerived> - void addToDense(MatrixBase<DenseDerived> &other) const + void addTo(MatrixBase<DenseDerived> &other) const { other.diagonal() += diagonal(); } template<typename DenseDerived> - void subToDense(MatrixBase<DenseDerived> &other) const + void subTo(MatrixBase<DenseDerived> &other) const { other.diagonal() -= diagonal(); } inline const DiagonalVectorType& diagonal() const { return derived().diagonal(); } @@ -83,16 +82,18 @@ void DiagonalBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const other.setZero(); other.diagonal() = diagonal(); } +#endif /** \class DiagonalMatrix - * \nonstableyet * * \brief Represents a diagonal matrix with its storage * * \param _Scalar the type of coefficients - * \param _Size the dimension of the matrix, or Dynamic + * \param SizeAtCompileTime the dimension of the matrix, or Dynamic + * \param MaxSizeAtCompileTime the dimension of the matrix, or Dynamic. This parameter is optional and defaults + * to SizeAtCompileTime. Most of the time, you do not need to specify it. * - * \sa class Matrix + * \sa class DiagonalWrapper */ template<typename _Scalar, int SizeAtCompileTime, int MaxSizeAtCompileTime> struct ei_traits<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> > @@ -106,10 +107,11 @@ class DiagonalMatrix : public DiagonalBase<DiagonalMatrix<_Scalar,SizeAtCompileTime,MaxSizeAtCompileTime> > { public: - + #ifndef EIGEN_PARSED_BY_DOXYGEN typedef typename ei_traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType; typedef const DiagonalMatrix& Nested; typedef _Scalar Scalar; + #endif protected: @@ -117,7 +119,9 @@ class DiagonalMatrix public: + /** const version of diagonal(). */ inline const DiagonalVectorType& diagonal() const { return m_diagonal; } + /** \returns a reference to the stored vector of diagonal coefficients. */ inline DiagonalVectorType& diagonal() { return m_diagonal; } /** Default constructor without initialization */ @@ -126,23 +130,27 @@ class DiagonalMatrix /** Constructs a diagonal matrix with given dimension */ inline DiagonalMatrix(int dim) : m_diagonal(dim) {} - /** 2D only */ + /** 2D constructor. */ inline DiagonalMatrix(const Scalar& x, const Scalar& y) : m_diagonal(x,y) {} - /** 3D only */ + /** 3D constructor. */ inline DiagonalMatrix(const Scalar& x, const Scalar& y, const Scalar& z) : m_diagonal(x,y,z) {} + /** Copy constructor. */ template<typename OtherDerived> inline DiagonalMatrix(const DiagonalBase<OtherDerived>& other) : m_diagonal(other.diagonal()) {} + #ifndef EIGEN_PARSED_BY_DOXYGEN /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */ inline DiagonalMatrix(const DiagonalMatrix& other) : m_diagonal(other.diagonal()) {} + #endif /** generic constructor from expression of the diagonal coefficients */ template<typename OtherDerived> explicit inline DiagonalMatrix(const MatrixBase<OtherDerived>& other) : m_diagonal(other) {} + /** Copy operator. */ template<typename OtherDerived> DiagonalMatrix& operator=(const DiagonalBase<OtherDerived>& other) { @@ -150,6 +158,7 @@ class DiagonalMatrix return *this; } + #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is a special case of the templated operator=. Its purpose is to * prevent a default operator= from hiding the templated operator=. */ @@ -158,23 +167,28 @@ class DiagonalMatrix m_diagonal = other.m_diagonal(); return *this; } + #endif + /** Resizes to given size. */ inline void resize(int size) { m_diagonal.resize(size); } + /** Sets all coefficients to zero. */ inline void setZero() { m_diagonal.setZero(); } + /** Resizes and sets all coefficients to zero. */ inline void setZero(int size) { m_diagonal.setZero(size); } + /** Sets this matrix to be the identity matrix of the current size. */ inline void setIdentity() { m_diagonal.setOnes(); } + /** Sets this matrix to be the identity matrix of the given size. */ inline void setIdentity(int size) { m_diagonal.setOnes(size); } }; /** \class DiagonalWrapper - * \nonstableyet * * \brief Expression of a diagonal matrix * * \param _DiagonalVectorType the type of the vector of diagonal coefficients * - * This class is an expression of a diagonal matrix with given vector of diagonal - * coefficients. It is the return type of MatrixBase::asDiagonal() + * This class is an expression of a diagonal matrix, but not storing its own vector of diagonal coefficients, + * instead wrapping an existing vector expression. It is the return type of MatrixBase::asDiagonal() * and most of the time this is the only way that it is used. * * \sa class DiagonalMatrix, class DiagonalBase, MatrixBase::asDiagonal() @@ -198,18 +212,22 @@ class DiagonalWrapper : public DiagonalBase<DiagonalWrapper<_DiagonalVectorType> >, ei_no_assignment_operator { public: + #ifndef EIGEN_PARSED_BY_DOXYGEN typedef _DiagonalVectorType DiagonalVectorType; typedef DiagonalWrapper Nested; + #endif + /** Constructor from expression of diagonal coefficients to wrap. */ inline DiagonalWrapper(const DiagonalVectorType& diagonal) : m_diagonal(diagonal) {} + + /** \returns a const reference to the wrapped expression of diagonal coefficients. */ const DiagonalVectorType& diagonal() const { return m_diagonal; } protected: const typename DiagonalVectorType::Nested m_diagonal; }; -/** \nonstableyet - * \returns a pseudo-expression of a diagonal matrix with *this as vector of diagonal coefficients +/** \returns a pseudo-expression of a diagonal matrix with *this as vector of diagonal coefficients * * \only_for_vectors * @@ -225,8 +243,7 @@ MatrixBase<Derived>::asDiagonal() const return derived(); } -/** \nonstableyet - * \returns true if *this is approximately equal to a diagonal matrix, +/** \returns true if *this is approximately equal to a diagonal matrix, * within the precision given by \a prec. * * Example: \include MatrixBase_isDiagonal.cpp diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 4a164a99e..f0c520b1f 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -34,10 +34,10 @@ struct ei_dot_traits { public: enum { - Vectorization = (int(Derived1::Flags)&int(Derived2::Flags)&ActualPacketAccessBit) + Traversal = (int(Derived1::Flags)&int(Derived2::Flags)&ActualPacketAccessBit) && (int(Derived1::Flags)&int(Derived2::Flags)&LinearAccessBit) - ? LinearVectorization - : NoVectorization + ? LinearVectorizedTraversal + : DefaultTraversal }; private: @@ -46,7 +46,7 @@ private: PacketSize = ei_packet_traits<Scalar>::size, Cost = Derived1::SizeAtCompileTime * (Derived1::CoeffReadCost + Derived2::CoeffReadCost + NumTraits<Scalar>::MulCost) + (Derived1::SizeAtCompileTime-1) * NumTraits<Scalar>::AddCost, - UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Vectorization) == int(NoVectorization) ? 1 : int(PacketSize)) + UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Traversal) == int(DefaultTraversal) ? 1 : int(PacketSize)) }; public: @@ -142,13 +142,13 @@ struct ei_dot_vec_unroller<Derived1, Derived2, Index, Stop, true> ***************************************************************************/ template<typename Derived1, typename Derived2, - int Vectorization = ei_dot_traits<Derived1, Derived2>::Vectorization, + int Traversal = ei_dot_traits<Derived1, Derived2>::Traversal, int Unrolling = ei_dot_traits<Derived1, Derived2>::Unrolling > struct ei_dot_impl; template<typename Derived1, typename Derived2> -struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling> +struct ei_dot_impl<Derived1, Derived2, DefaultTraversal, NoUnrolling> { typedef typename Derived1::Scalar Scalar; static Scalar run(const Derived1& v1, const Derived2& v2) @@ -163,12 +163,12 @@ struct ei_dot_impl<Derived1, Derived2, NoVectorization, NoUnrolling> }; template<typename Derived1, typename Derived2> -struct ei_dot_impl<Derived1, Derived2, NoVectorization, CompleteUnrolling> +struct ei_dot_impl<Derived1, Derived2, DefaultTraversal, CompleteUnrolling> : public ei_dot_novec_unroller<Derived1, Derived2, 0, Derived1::SizeAtCompileTime> {}; template<typename Derived1, typename Derived2> -struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling> +struct ei_dot_impl<Derived1, Derived2, LinearVectorizedTraversal, NoUnrolling> { typedef typename Derived1::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; @@ -221,20 +221,20 @@ struct ei_dot_impl<Derived1, Derived2, LinearVectorization, NoUnrolling> }; template<typename Derived1, typename Derived2> -struct ei_dot_impl<Derived1, Derived2, LinearVectorization, CompleteUnrolling> +struct ei_dot_impl<Derived1, Derived2, LinearVectorizedTraversal, CompleteUnrolling> { typedef typename Derived1::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; enum { PacketSize = ei_packet_traits<Scalar>::size, Size = Derived1::SizeAtCompileTime, - VectorizationSize = (Size / PacketSize) * PacketSize + VectorizedSize = (Size / PacketSize) * PacketSize }; static Scalar run(const Derived1& v1, const Derived2& v2) { - Scalar res = ei_predux(ei_dot_vec_unroller<Derived1, Derived2, 0, VectorizationSize>::run(v1, v2)); - if (VectorizationSize != Size) - res += ei_dot_novec_unroller<Derived1, Derived2, VectorizationSize, Size-VectorizationSize>::run(v1, v2); + Scalar res = ei_predux(ei_dot_vec_unroller<Derived1, Derived2, 0, VectorizedSize>::run(v1, v2)); + if (VectorizedSize != Size) + res += ei_dot_novec_unroller<Derived1, Derived2, VectorizedSize, Size-VectorizedSize>::run(v1, v2); return res; } }; diff --git a/Eigen/src/Core/ExpressionMaker.h b/Eigen/src/Core/ExpressionMaker.h index 1d265b63c..7e2b81d4a 100644 --- a/Eigen/src/Core/ExpressionMaker.h +++ b/Eigen/src/Core/ExpressionMaker.h @@ -37,12 +37,6 @@ template<typename XprType> struct ei_shape_of // matrix. Unless we change the overall design, here is a workaround. // There is an example in unsuported/Eigen/src/AutoDiff/AutoDiffScalar. -template<typename XprType, int Shape = ei_shape_of<XprType>::ret> -struct MakeNestByValue -{ - typedef NestByValue<XprType> Type; -}; - template<typename Func, typename XprType, int Shape = ei_shape_of<XprType>::ret> struct MakeCwiseUnaryOp { diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index 259f40244..aa3eba5cc 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -32,7 +32,8 @@ * * \sa class CwiseBinaryOp, MatrixBase::operator+, class VectorwiseOp, MatrixBase::sum() */ -template<typename Scalar> struct ei_scalar_sum_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_sum_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_sum_op) EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return a + b; } template<typename PacketScalar> EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const @@ -54,7 +55,8 @@ struct ei_functor_traits<ei_scalar_sum_op<Scalar> > { * * \sa class CwiseBinaryOp, Cwise::operator*(), class VectorwiseOp, MatrixBase::redux() */ -template<typename Scalar> struct ei_scalar_product_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_product_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_product_op) EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return a * b; } template<typename PacketScalar> EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const @@ -76,7 +78,8 @@ struct ei_functor_traits<ei_scalar_product_op<Scalar> > { * * \sa class CwiseBinaryOp, MatrixBase::cwiseMin, class VectorwiseOp, MatrixBase::minCoeff() */ -template<typename Scalar> struct ei_scalar_min_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_min_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_min_op) EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return std::min(a, b); } template<typename PacketScalar> EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const @@ -98,7 +101,8 @@ struct ei_functor_traits<ei_scalar_min_op<Scalar> > { * * \sa class CwiseBinaryOp, MatrixBase::cwiseMax, class VectorwiseOp, MatrixBase::maxCoeff() */ -template<typename Scalar> struct ei_scalar_max_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_max_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_max_op) EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return std::max(a, b); } template<typename PacketScalar> EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const @@ -120,7 +124,8 @@ struct ei_functor_traits<ei_scalar_max_op<Scalar> > { * * \sa MatrixBase::stableNorm(), class Redux */ -template<typename Scalar> struct ei_scalar_hypot_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_hypot_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_hypot_op) // typedef typename NumTraits<Scalar>::Real result_type; EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& _x, const Scalar& _y) const { @@ -142,7 +147,8 @@ struct ei_functor_traits<ei_scalar_hypot_op<Scalar> > { * * \sa class CwiseBinaryOp, MatrixBase::operator- */ -template<typename Scalar> struct ei_scalar_difference_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_difference_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_difference_op) EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return a - b; } template<typename PacketScalar> EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const @@ -161,7 +167,8 @@ struct ei_functor_traits<ei_scalar_difference_op<Scalar> > { * * \sa class CwiseBinaryOp, Cwise::operator/() */ -template<typename Scalar> struct ei_scalar_quotient_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_quotient_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_quotient_op) EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a, const Scalar& b) const { return a / b; } template<typename PacketScalar> EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a, const PacketScalar& b) const @@ -185,7 +192,8 @@ struct ei_functor_traits<ei_scalar_quotient_op<Scalar> > { * * \sa class CwiseUnaryOp, MatrixBase::operator- */ -template<typename Scalar> struct ei_scalar_opposite_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_opposite_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_opposite_op) EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const { return -a; } template<typename PacketScalar> EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const @@ -203,7 +211,8 @@ struct ei_functor_traits<ei_scalar_opposite_op<Scalar> > * * \sa class CwiseUnaryOp, Cwise::abs */ -template<typename Scalar> struct ei_scalar_abs_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_abs_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_abs_op) typedef typename NumTraits<Scalar>::Real result_type; EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return ei_abs(a); } template<typename PacketScalar> @@ -224,7 +233,8 @@ struct ei_functor_traits<ei_scalar_abs_op<Scalar> > * * \sa class CwiseUnaryOp, Cwise::abs2 */ -template<typename Scalar> struct ei_scalar_abs2_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_abs2_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_abs2_op) typedef typename NumTraits<Scalar>::Real result_type; EIGEN_STRONG_INLINE const result_type operator() (const Scalar& a) const { return ei_abs2(a); } template<typename PacketScalar> @@ -240,7 +250,8 @@ struct ei_functor_traits<ei_scalar_abs2_op<Scalar> > * * \sa class CwiseUnaryOp, MatrixBase::conjugate() */ -template<typename Scalar> struct ei_scalar_conjugate_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_conjugate_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_conjugate_op) EIGEN_STRONG_INLINE const Scalar operator() (const Scalar& a) const { return ei_conj(a); } template<typename PacketScalar> EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const { return a; } @@ -260,7 +271,8 @@ struct ei_functor_traits<ei_scalar_conjugate_op<Scalar> > * \sa class CwiseUnaryOp, MatrixBase::cast() */ template<typename Scalar, typename NewType> -struct ei_scalar_cast_op EIGEN_EMPTY_STRUCT { +struct ei_scalar_cast_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_cast_op) typedef NewType result_type; EIGEN_STRONG_INLINE const NewType operator() (const Scalar& a) const { return static_cast<NewType>(a); } }; @@ -274,7 +286,8 @@ struct ei_functor_traits<ei_scalar_cast_op<Scalar,NewType> > * \sa class CwiseUnaryOp, MatrixBase::real() */ template<typename Scalar> -struct ei_scalar_real_op EIGEN_EMPTY_STRUCT { +struct ei_scalar_real_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_real_op) typedef typename NumTraits<Scalar>::Real result_type; EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return ei_real(a); } EIGEN_STRONG_INLINE result_type& operator() (Scalar& a) const { return ei_real_ref(a); } @@ -289,7 +302,8 @@ struct ei_functor_traits<ei_scalar_real_op<Scalar> > * \sa class CwiseUnaryOp, MatrixBase::imag() */ template<typename Scalar> -struct ei_scalar_imag_op EIGEN_EMPTY_STRUCT { +struct ei_scalar_imag_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_imag_op) typedef typename NumTraits<Scalar>::Real result_type; EIGEN_STRONG_INLINE result_type operator() (const Scalar& a) const { return ei_imag(a); } EIGEN_STRONG_INLINE result_type& operator() (Scalar& a) const { return ei_imag_ref(a); } @@ -304,7 +318,8 @@ struct ei_functor_traits<ei_scalar_imag_op<Scalar> > * * \sa class CwiseUnaryOp, Cwise::exp() */ -template<typename Scalar> struct ei_scalar_exp_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_exp_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_exp_op) inline const Scalar operator() (const Scalar& a) const { return ei_exp(a); } typedef typename ei_packet_traits<Scalar>::type Packet; inline Packet packetOp(const Packet& a) const { return ei_pexp(a); } @@ -319,7 +334,8 @@ struct ei_functor_traits<ei_scalar_exp_op<Scalar> > * * \sa class CwiseUnaryOp, Cwise::log() */ -template<typename Scalar> struct ei_scalar_log_op EIGEN_EMPTY_STRUCT { +template<typename Scalar> struct ei_scalar_log_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_log_op) inline const Scalar operator() (const Scalar& a) const { return ei_log(a); } typedef typename ei_packet_traits<Scalar>::type Packet; inline Packet packetOp(const Packet& a) const { return ei_plog(a); } @@ -351,8 +367,6 @@ struct ei_scalar_multiple_op { EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const { return ei_pmul(a, ei_pset1(m_other)); } typename ei_makeconst<typename NumTraits<Scalar>::Nested>::type m_other; -private: - ei_scalar_multiple_op& operator=(const ei_scalar_multiple_op&); }; template<typename Scalar> struct ei_functor_traits<ei_scalar_multiple_op<Scalar> > @@ -380,8 +394,6 @@ struct ei_scalar_quotient1_impl { EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const { return ei_pmul(a, ei_pset1(m_other)); } const Scalar m_other; -private: - ei_scalar_quotient1_impl& operator=(const ei_scalar_quotient1_impl&); }; template<typename Scalar> struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,true> > @@ -427,15 +439,13 @@ struct ei_scalar_constant_op { EIGEN_STRONG_INLINE const Scalar operator() (int, int = 0) const { return m_other; } EIGEN_STRONG_INLINE const PacketScalar packetOp() const { return ei_pset1(m_other); } const Scalar m_other; -private: - ei_scalar_constant_op& operator=(const ei_scalar_constant_op&); }; template<typename Scalar> struct ei_functor_traits<ei_scalar_constant_op<Scalar> > { enum { Cost = 1, PacketAccess = ei_packet_traits<Scalar>::size>1, IsRepeatable = true }; }; -template<typename Scalar> struct ei_scalar_identity_op EIGEN_EMPTY_STRUCT { - EIGEN_STRONG_INLINE ei_scalar_identity_op(void) {} +template<typename Scalar> struct ei_scalar_identity_op { + EIGEN_EMPTY_STRUCT_CTOR(ei_scalar_identity_op) EIGEN_STRONG_INLINE const Scalar operator() (int row, int col) const { return row==col ? Scalar(1) : Scalar(0); } }; template<typename Scalar> diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 05469b340..7ffddcbf8 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -30,7 +30,7 @@ template<typename T> inline typename NumTraits<T>::Real epsilon() return std::numeric_limits<typename NumTraits<T>::Real>::epsilon(); } -template<typename T> inline typename NumTraits<T>::Real precision(); +template<typename T> inline typename NumTraits<T>::Real dummy_precision(); template<typename T> inline T ei_random(T a, T b); template<typename T> inline T ei_random(); @@ -55,7 +55,7 @@ template<typename T> inline typename NumTraits<T>::Real ei_hypot(T x, T y) *** int *** **************/ -template<> inline int precision<int>() { return 0; } +template<> inline int dummy_precision<int>() { return 0; } inline int ei_real(int x) { return x; } inline int& ei_real_ref(int& x) { return x; } inline int ei_imag(int) { return 0; } @@ -92,15 +92,15 @@ template<> inline int ei_random() { return ei_random<int>(-ei_random_amplitude<int>(), ei_random_amplitude<int>()); } -inline bool ei_isMuchSmallerThan(int a, int, int = precision<int>()) +inline bool ei_isMuchSmallerThan(int a, int, int = dummy_precision<int>()) { return a == 0; } -inline bool ei_isApprox(int a, int b, int = precision<int>()) +inline bool ei_isApprox(int a, int b, int = dummy_precision<int>()) { return a == b; } -inline bool ei_isApproxOrLessThan(int a, int b, int = precision<int>()) +inline bool ei_isApproxOrLessThan(int a, int b, int = dummy_precision<int>()) { return a <= b; } @@ -109,7 +109,7 @@ inline bool ei_isApproxOrLessThan(int a, int b, int = precision<int>()) *** float *** **************/ -template<> inline float precision<float>() { return 1e-5f; } +template<> inline float dummy_precision<float>() { return 1e-5f; } inline float ei_real(float x) { return x; } inline float& ei_real_ref(float& x) { return x; } inline float ei_imag(float) { return 0.f; } @@ -140,15 +140,15 @@ template<> inline float ei_random() { return ei_random<float>(-ei_random_amplitude<float>(), ei_random_amplitude<float>()); } -inline bool ei_isMuchSmallerThan(float a, float b, float prec = precision<float>()) +inline bool ei_isMuchSmallerThan(float a, float b, float prec = dummy_precision<float>()) { return ei_abs(a) <= ei_abs(b) * prec; } -inline bool ei_isApprox(float a, float b, float prec = precision<float>()) +inline bool ei_isApprox(float a, float b, float prec = dummy_precision<float>()) { return ei_abs(a - b) <= std::min(ei_abs(a), ei_abs(b)) * prec; } -inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float>()) +inline bool ei_isApproxOrLessThan(float a, float b, float prec = dummy_precision<float>()) { return a <= b || ei_isApprox(a, b, prec); } @@ -157,7 +157,7 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float *** double *** **************/ -template<> inline double precision<double>() { return 1e-12; } +template<> inline double dummy_precision<double>() { return 1e-12; } inline double ei_real(double x) { return x; } inline double& ei_real_ref(double& x) { return x; } @@ -189,15 +189,15 @@ template<> inline double ei_random() { return ei_random<double>(-ei_random_amplitude<double>(), ei_random_amplitude<double>()); } -inline bool ei_isMuchSmallerThan(double a, double b, double prec = precision<double>()) +inline bool ei_isMuchSmallerThan(double a, double b, double prec = dummy_precision<double>()) { return ei_abs(a) <= ei_abs(b) * prec; } -inline bool ei_isApprox(double a, double b, double prec = precision<double>()) +inline bool ei_isApprox(double a, double b, double prec = dummy_precision<double>()) { return ei_abs(a - b) <= std::min(ei_abs(a), ei_abs(b)) * prec; } -inline bool ei_isApproxOrLessThan(double a, double b, double prec = precision<double>()) +inline bool ei_isApproxOrLessThan(double a, double b, double prec = dummy_precision<double>()) { return a <= b || ei_isApprox(a, b, prec); } @@ -206,7 +206,7 @@ inline bool ei_isApproxOrLessThan(double a, double b, double prec = precision<do *** complex<float> *** *********************/ -template<> inline float precision<std::complex<float> >() { return precision<float>(); } +template<> inline float dummy_precision<std::complex<float> >() { return dummy_precision<float>(); } inline float ei_real(const std::complex<float>& x) { return std::real(x); } inline float ei_imag(const std::complex<float>& x) { return std::imag(x); } inline float& ei_real_ref(std::complex<float>& x) { return reinterpret_cast<float*>(&x)[0]; } @@ -224,15 +224,15 @@ template<> inline std::complex<float> ei_random() { return std::complex<float>(ei_random<float>(), ei_random<float>()); } -inline bool ei_isMuchSmallerThan(const std::complex<float>& a, const std::complex<float>& b, float prec = precision<float>()) +inline bool ei_isMuchSmallerThan(const std::complex<float>& a, const std::complex<float>& b, float prec = dummy_precision<float>()) { return ei_abs2(a) <= ei_abs2(b) * prec * prec; } -inline bool ei_isMuchSmallerThan(const std::complex<float>& a, float b, float prec = precision<float>()) +inline bool ei_isMuchSmallerThan(const std::complex<float>& a, float b, float prec = dummy_precision<float>()) { return ei_abs2(a) <= ei_abs2(b) * prec * prec; } -inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>& b, float prec = precision<float>()) +inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>& b, float prec = dummy_precision<float>()) { return ei_isApprox(ei_real(a), ei_real(b), prec) && ei_isApprox(ei_imag(a), ei_imag(b), prec); @@ -243,7 +243,7 @@ inline bool ei_isApprox(const std::complex<float>& a, const std::complex<float>& *** complex<double> *** **********************/ -template<> inline double precision<std::complex<double> >() { return precision<double>(); } +template<> inline double dummy_precision<std::complex<double> >() { return dummy_precision<double>(); } inline double ei_real(const std::complex<double>& x) { return std::real(x); } inline double ei_imag(const std::complex<double>& x) { return std::imag(x); } inline double& ei_real_ref(std::complex<double>& x) { return reinterpret_cast<double*>(&x)[0]; } @@ -261,15 +261,15 @@ template<> inline std::complex<double> ei_random() { return std::complex<double>(ei_random<double>(), ei_random<double>()); } -inline bool ei_isMuchSmallerThan(const std::complex<double>& a, const std::complex<double>& b, double prec = precision<double>()) +inline bool ei_isMuchSmallerThan(const std::complex<double>& a, const std::complex<double>& b, double prec = dummy_precision<double>()) { return ei_abs2(a) <= ei_abs2(b) * prec * prec; } -inline bool ei_isMuchSmallerThan(const std::complex<double>& a, double b, double prec = precision<double>()) +inline bool ei_isMuchSmallerThan(const std::complex<double>& a, double b, double prec = dummy_precision<double>()) { return ei_abs2(a) <= ei_abs2(b) * prec * prec; } -inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double>& b, double prec = precision<double>()) +inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double>& b, double prec = dummy_precision<double>()) { return ei_isApprox(ei_real(a), ei_real(b), prec) && ei_isApprox(ei_imag(a), ei_imag(b), prec); @@ -281,7 +281,7 @@ inline bool ei_isApprox(const std::complex<double>& a, const std::complex<double *** long double *** ******************/ -template<> inline long double precision<long double>() { return precision<double>(); } +template<> inline long double dummy_precision<long double>() { return dummy_precision<double>(); } inline long double ei_real(long double x) { return x; } inline long double& ei_real_ref(long double& x) { return x; } inline long double ei_imag(long double) { return 0.; } @@ -304,15 +304,15 @@ template<> inline long double ei_random() { return ei_random<double>(-ei_random_amplitude<double>(), ei_random_amplitude<double>()); } -inline bool ei_isMuchSmallerThan(long double a, long double b, long double prec = precision<long double>()) +inline bool ei_isMuchSmallerThan(long double a, long double b, long double prec = dummy_precision<long double>()) { return ei_abs(a) <= ei_abs(b) * prec; } -inline bool ei_isApprox(long double a, long double b, long double prec = precision<long double>()) +inline bool ei_isApprox(long double a, long double b, long double prec = dummy_precision<long double>()) { return ei_abs(a - b) <= std::min(ei_abs(a), ei_abs(b)) * prec; } -inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec = precision<long double>()) +inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec = dummy_precision<long double>()) { return a <= b || ei_isApprox(a, b, prec); } @@ -321,7 +321,7 @@ inline bool ei_isApproxOrLessThan(long double a, long double b, long double prec *** bool *** **************/ -template<> inline bool precision<bool>() { return 0; } +template<> inline bool dummy_precision<bool>() { return 0; } inline bool ei_real(bool x) { return x; } inline bool& ei_real_ref(bool& x) { return x; } inline bool ei_imag(bool) { return 0; } @@ -334,15 +334,15 @@ template<> inline bool ei_random() { return (ei_random<int>(0,1) == 1); } -inline bool ei_isMuchSmallerThan(bool a, bool, bool = precision<bool>()) +inline bool ei_isMuchSmallerThan(bool a, bool, bool = dummy_precision<bool>()) { return !a; } -inline bool ei_isApprox(bool a, bool b, bool = precision<bool>()) +inline bool ei_isApprox(bool a, bool b, bool = dummy_precision<bool>()) { return a == b; } -inline bool ei_isApproxOrLessThan(bool a, bool b, bool = precision<bool>()) +inline bool ei_isApproxOrLessThan(bool a, bool b, bool = dummy_precision<bool>()) { return int(a) <= int(b); } diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 75c0a73b5..24c618549 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -152,7 +152,6 @@ class Matrix using Base::coeff; using Base::coeffRef; - /** Copies the value of the expression \a other into \c *this with automatic resizing. * * *this might be resized to match the dimensions of \a other. If *this was a null matrix (not already initialized), @@ -286,9 +285,6 @@ class Matrix other.evalTo(*this); } - /** Destructor */ - inline ~Matrix() {} - /** \sa MatrixBase::operator=(const AnyMatrixBase<OtherDerived>&) */ template<typename OtherDerived> EIGEN_STRONG_INLINE Matrix(const AnyMatrixBase<OtherDerived> &other) diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 53dd6f991..f7b32a650 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -144,7 +144,7 @@ template<typename Derived> class MatrixBase typedef CwiseNullaryOp<ei_scalar_constant_op<Scalar>,Derived> ConstantReturnType; /** \internal the return type of MatrixBase::adjoint() */ typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex, - CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<Derived> > >, + CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Eigen::Transpose<Derived> >, Transpose<Derived> >::ret AdjointReturnType; /** \internal the return type of MatrixBase::eigenvalues() */ @@ -259,16 +259,16 @@ template<typename Derived> class MatrixBase Derived& setIdentity(); - bool isIdentity(RealScalar prec = precision<Scalar>()) const; - bool isDiagonal(RealScalar prec = precision<Scalar>()) const; + bool isIdentity(RealScalar prec = dummy_precision<Scalar>()) const; + bool isDiagonal(RealScalar prec = dummy_precision<Scalar>()) const; - bool isUpperTriangular(RealScalar prec = precision<Scalar>()) const; - bool isLowerTriangular(RealScalar prec = precision<Scalar>()) const; + bool isUpperTriangular(RealScalar prec = dummy_precision<Scalar>()) const; + bool isLowerTriangular(RealScalar prec = dummy_precision<Scalar>()) const; template<typename OtherDerived> bool isOrthogonal(const MatrixBase<OtherDerived>& other, - RealScalar prec = precision<Scalar>()) const; - bool isUnitary(RealScalar prec = precision<Scalar>()) const; + RealScalar prec = dummy_precision<Scalar>()) const; + bool isUnitary(RealScalar prec = dummy_precision<Scalar>()) const; /** \returns true if each coefficients of \c *this and \a other are all exactly equal. * \warning When using floating point scalar values you probably should rather use a @@ -332,13 +332,13 @@ template<typename Derived> class MatrixBase ResultType& inverse, typename ResultType::Scalar& determinant, bool& invertible, - const RealScalar& absDeterminantThreshold = precision<Scalar>() + const RealScalar& absDeterminantThreshold = dummy_precision<Scalar>() ) const; template<typename ResultType> void computeInverseWithCheck( ResultType& inverse, bool& invertible, - const RealScalar& absDeterminantThreshold = precision<Scalar>() + const RealScalar& absDeterminantThreshold = dummy_precision<Scalar>() ) const; Scalar determinant() const; @@ -376,7 +376,7 @@ template<typename Derived> class MatrixBase ei_traits<Derived>::ColsAtCompileTime==1 ? SizeMinusOne : 1, ei_traits<Derived>::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne; typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, - NestByValue<StartMinusOne> > HNormalizedReturnType; + StartMinusOne > HNormalizedReturnType; const HNormalizedReturnType hnormalized() const; typedef Homogeneous<Derived,MatrixBase<Derived>::ColsAtCompileTime==1?Vertical:Horizontal> HomogeneousReturnType; diff --git a/Eigen/src/Core/MatrixStorage.h b/Eigen/src/Core/MatrixStorage.h index 73b17e63e..8bfa728b6 100644 --- a/Eigen/src/Core/MatrixStorage.h +++ b/Eigen/src/Core/MatrixStorage.h @@ -117,7 +117,6 @@ template<typename T, int Size, int _Options> class ei_matrix_storage<T, Size, Dy inline ei_matrix_storage(ei_constructor_without_unaligned_array_assert) : m_data(ei_constructor_without_unaligned_array_assert()), m_rows(0), m_cols(0) {} inline ei_matrix_storage(int, int rows, int cols) : m_rows(rows), m_cols(cols) {} - inline ~ei_matrix_storage() {} inline void swap(ei_matrix_storage& other) { std::swap(m_data,other.m_data); std::swap(m_rows,other.m_rows); std::swap(m_cols,other.m_cols); } inline int rows(void) const {return m_rows;} @@ -141,7 +140,6 @@ template<typename T, int Size, int _Cols, int _Options> class ei_matrix_storage< inline ei_matrix_storage(ei_constructor_without_unaligned_array_assert) : m_data(ei_constructor_without_unaligned_array_assert()), m_rows(0) {} inline ei_matrix_storage(int, int rows, int) : m_rows(rows) {} - inline ~ei_matrix_storage() {} inline void swap(ei_matrix_storage& other) { std::swap(m_data,other.m_data); std::swap(m_rows,other.m_rows); } inline int rows(void) const {return m_rows;} inline int cols(void) const {return _Cols;} @@ -163,7 +161,6 @@ template<typename T, int Size, int _Rows, int _Options> class ei_matrix_storage< inline ei_matrix_storage(ei_constructor_without_unaligned_array_assert) : m_data(ei_constructor_without_unaligned_array_assert()), m_cols(0) {} inline ei_matrix_storage(int, int, int cols) : m_cols(cols) {} - inline ~ei_matrix_storage() {} inline void swap(ei_matrix_storage& other) { std::swap(m_data,other.m_data); std::swap(m_cols,other.m_cols); } inline int rows(void) const {return _Rows;} inline int cols(void) const {return m_cols;} diff --git a/Eigen/src/Core/Minor.h b/Eigen/src/Core/Minor.h index 3ed0fec6b..629dbe609 100644 --- a/Eigen/src/Core/Minor.h +++ b/Eigen/src/Core/Minor.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -54,7 +54,8 @@ struct ei_traits<Minor<MatrixType> > MaxColsAtCompileTime = (MatrixType::MaxColsAtCompileTime != Dynamic) ? int(MatrixType::MaxColsAtCompileTime) - 1 : Dynamic, Flags = _MatrixTypeNested::Flags & HereditaryBits, - CoeffReadCost = _MatrixTypeNested::CoeffReadCost + CoeffReadCost = _MatrixTypeNested::CoeffReadCost // minor is used typically on tiny matrices, + // where loops are unrolled and the 'if' evaluates at compile time }; }; diff --git a/Eigen/src/Core/NestByValue.h b/Eigen/src/Core/NestByValue.h index 94a8f8078..85a672779 100644 --- a/Eigen/src/Core/NestByValue.h +++ b/Eigen/src/Core/NestByValue.h @@ -102,9 +102,6 @@ template<typename ExpressionType> class NestByValue protected: const ExpressionType m_expression; - - private: - NestByValue& operator=(const NestByValue&); }; /** \returns an expression of the temporary version of *this. diff --git a/Eigen/src/Core/NoAlias.h b/Eigen/src/Core/NoAlias.h index 7ed848bce..bfea5c91d 100644 --- a/Eigen/src/Core/NoAlias.h +++ b/Eigen/src/Core/NoAlias.h @@ -73,9 +73,6 @@ class NoAlias protected: ExpressionType& m_expression; - - private: - NoAlias& operator=(const NoAlias&); }; /** \returns a pseudo expression of \c *this with an operator= assuming diff --git a/Eigen/src/Core/PermutationMatrix.h b/Eigen/src/Core/PermutationMatrix.h index aaccb4e7b..284baf678 100644 --- a/Eigen/src/Core/PermutationMatrix.h +++ b/Eigen/src/Core/PermutationMatrix.h @@ -25,18 +25,24 @@ #ifndef EIGEN_PERMUTATIONMATRIX_H #define EIGEN_PERMUTATIONMATRIX_H -/** \nonstableyet - * \class PermutationMatrix +/** \class PermutationMatrix * * \brief Permutation matrix * * \param SizeAtCompileTime the number of rows/cols, or Dynamic - * \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. + * \param MaxSizeAtCompileTime the maximum number of rows/cols, or Dynamic. This optional parameter defaults to SizeAtCompileTime. Most of the time, you should not have to specify it. * * This class represents a permutation matrix, internally stored as a vector of integers. - * The convention followed here is the same as on <a href="http://en.wikipedia.org/wiki/Permutation_matrix">Wikipedia</a>, - * namely: the matrix of permutation \a p is the matrix such that on each row \a i, the only nonzero coefficient is - * in column p(i). + * The convention followed here is that if \f$ \sigma \f$ is a permutation, the corresponding permutation matrix + * \f$ P_\sigma \f$ is such that if \f$ (e_1,\ldots,e_p) \f$ is the canonical basis, we have: + * \f[ P_\sigma(e_i) = e_{\sigma(i)}. \f] + * This convention ensures that for any two permutations \f$ \sigma, \tau \f$, we have: + * \f[ P_{\sigma\circ\tau} = P_\sigma P_\tau. \f] + * + * Permutation matrices are square and invertible. + * + * Notice that in addition to the member functions and operators listed here, there also are non-member + * operator* to multiply a PermutationMatrix with any kind of matrix expression (MatrixBase) on either side. * * \sa class DiagonalMatrix */ @@ -53,6 +59,7 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi { public: + #ifndef EIGEN_PARSED_BY_DOXYGEN typedef ei_traits<PermutationMatrix> Traits; typedef Matrix<int,SizeAtCompileTime,SizeAtCompileTime,0,MaxSizeAtCompileTime,MaxSizeAtCompileTime> DenseMatrixType; @@ -65,25 +72,37 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi MaxColsAtCompileTime = Traits::MaxColsAtCompileTime }; typedef typename Traits::Scalar Scalar; + #endif - typedef Matrix<int, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> IndicesType; + typedef Matrix<int, SizeAtCompileTime, 1, 0, MaxSizeAtCompileTime, 1> IndicesType; inline PermutationMatrix() { } + /** Copy constructor. */ template<int OtherSize, int OtherMaxSize> inline PermutationMatrix(const PermutationMatrix<OtherSize, OtherMaxSize>& other) : m_indices(other.indices()) {} - /** copy constructor. prevent a default copy constructor from hiding the other templated constructor */ + #ifndef EIGEN_PARSED_BY_DOXYGEN + /** Standard copy constructor. Defined only to prevent a default copy constructor + * from hiding the other templated constructor */ inline PermutationMatrix(const PermutationMatrix& other) : m_indices(other.indices()) {} + #endif - /** generic constructor from expression of the indices */ + /** Generic constructor from expression of the indices. The indices + * array has the meaning that the permutations sends each integer i to indices[i]. + * + * \warning It is your responsibility to check that the indices array that you passes actually + * describes a permutation, i.e., each value between 0 and n-1 occurs exactly once, where n is the + * array's size. + */ template<typename Other> - explicit inline PermutationMatrix(const MatrixBase<Other>& other) : m_indices(other) + explicit inline PermutationMatrix(const MatrixBase<Other>& indices) : m_indices(indices) {} + /** Copies the other permutation into *this */ template<int OtherSize, int OtherMaxSize> PermutationMatrix& operator=(const PermutationMatrix<OtherSize, OtherMaxSize>& other) { @@ -91,42 +110,143 @@ class PermutationMatrix : public AnyMatrixBase<PermutationMatrix<SizeAtCompileTi return *this; } + #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is a special case of the templated operator=. Its purpose is to * prevent a default operator= from hiding the templated operator=. */ PermutationMatrix& operator=(const PermutationMatrix& other) { - m_indices = other.m_indices(); + m_indices = other.m_indices; return *this; } + #endif - inline PermutationMatrix(int rows, int cols) : m_indices(rows) - { - ei_assert(rows == cols); - } + /** Constructs an uninitialized permutation matrix of given size. + */ + inline PermutationMatrix(int size) : m_indices(size) + {} - /** \returns the number of columns */ + /** \returns the number of rows */ inline int rows() const { return m_indices.size(); } - /** \returns the number of rows */ + /** \returns the number of columns */ inline int cols() const { return m_indices.size(); } + #ifndef EIGEN_PARSED_BY_DOXYGEN template<typename DenseDerived> void evalTo(MatrixBase<DenseDerived>& other) const { other.setZero(); for (int i=0; i<rows();++i) - other.coeffRef(i,m_indices.coeff(i)) = typename DenseDerived::Scalar(1); + other.coeffRef(m_indices.coeff(i),i) = typename DenseDerived::Scalar(1); } + #endif + /** \returns a Matrix object initialized from this permutation matrix. Notice that it + * is inefficient to return this Matrix object by value. For efficiency, favor using + * the Matrix constructor taking AnyMatrixBase objects. + */ DenseMatrixType toDenseMatrix() const { return *this; } + /** const version of indices(). */ const IndicesType& indices() const { return m_indices; } + /** \returns a reference to the stored array representing the permutation. */ IndicesType& indices() { return m_indices; } - + + /** Resizes to given size. + */ + inline void resize(int size) + { + m_indices.resize(size); + } + + /** Sets *this to be the identity permutation matrix */ + void setIdentity() + { + for(int i = 0; i < m_indices.size(); ++i) + m_indices.coeffRef(i) = i; + } + + /** Sets *this to be the identity permutation matrix of given size. + */ + void setIdentity(int size) + { + resize(size); + setIdentity(); + } + + /** Multiplies *this by the transposition \f$(ij)\f$ on the left. + * + * \returns a reference to *this. + * + * \warning This is much slower than applyTranspositionOnTheRight(int,int): + * this has linear complexity and requires a lot of branching. + * + * \sa applyTranspositionOnTheRight(int,int) + */ + PermutationMatrix& applyTranspositionOnTheLeft(int i, int j) + { + ei_assert(i>=0 && j>=0 && i<m_indices.size() && j<m_indices.size()); + for(int k = 0; k < m_indices.size(); ++k) + { + if(m_indices.coeff(k) == i) m_indices.coeffRef(k) = j; + else if(m_indices.coeff(k) == j) m_indices.coeffRef(k) = i; + } + return *this; + } + + /** Multiplies *this by the transposition \f$(ij)\f$ on the right. + * + * \returns a reference to *this. + * + * This is a fast operation, it only consists in swapping two indices. + * + * \sa applyTranspositionOnTheLeft(int,int) + */ + PermutationMatrix& applyTranspositionOnTheRight(int i, int j) + { + ei_assert(i>=0 && j>=0 && i<m_indices.size() && j<m_indices.size()); + std::swap(m_indices.coeffRef(i), m_indices.coeffRef(j)); + return *this; + } + + /**** inversion and multiplication helpers to hopefully get RVO ****/ + +#ifndef EIGEN_PARSED_BY_DOXYGEN + protected: + enum Inverse_t {Inverse}; + PermutationMatrix(Inverse_t, const PermutationMatrix& other) + : m_indices(other.m_indices.size()) + { + for (int i=0; i<rows();++i) m_indices.coeffRef(other.m_indices.coeff(i)) = i; + } + enum Product_t {Product}; + PermutationMatrix(Product_t, const PermutationMatrix& lhs, const PermutationMatrix& rhs) + : m_indices(lhs.m_indices.size()) + { + ei_assert(lhs.cols() == rhs.rows()); + for (int i=0; i<rows();++i) m_indices.coeffRef(i) = lhs.m_indices.coeff(rhs.m_indices.coeff(i)); + } +#endif + + public: + /** \returns the inverse permutation matrix. + * + * \note \note_try_to_help_rvo + */ + inline PermutationMatrix inverse() const + { return PermutationMatrix(Inverse, *this); } + /** \returns the product permutation matrix. + * + * \note \note_try_to_help_rvo + */ + template<int OtherSize, int OtherMaxSize> + inline PermutationMatrix operator*(const PermutationMatrix<OtherSize, OtherMaxSize>& other) const + { return PermutationMatrix(Product, *this, other); } + protected: IndicesType m_indices; @@ -185,7 +305,7 @@ struct ei_permut_matrix_product_retval Dest, Side==OnTheLeft ? 1 : Dest::RowsAtCompileTime, Side==OnTheRight ? 1 : Dest::ColsAtCompileTime - >(dst, Side==OnTheRight ? m_permutation.indices().coeff(i) : i) + >(dst, Side==OnTheLeft ? m_permutation.indices().coeff(i) : i) = @@ -193,7 +313,7 @@ struct ei_permut_matrix_product_retval MatrixTypeNestedCleaned, Side==OnTheLeft ? 1 : MatrixType::RowsAtCompileTime, Side==OnTheRight ? 1 : MatrixType::ColsAtCompileTime - >(m_matrix, Side==OnTheLeft ? m_permutation.indices().coeff(i) : i); + >(m_matrix, Side==OnTheRight ? m_permutation.indices().coeff(i) : i); } } diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 7db35eaad..15f44de89 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -58,10 +58,12 @@ enum { OuterProduct, InnerProduct, UnrolledProduct, GemvProduct, GemmProduct }; template<typename Lhs, typename Rhs> struct ei_product_type { + typedef typename ei_cleantype<Lhs>::type _Lhs; + typedef typename ei_cleantype<Rhs>::type _Rhs; enum { - Rows = Lhs::RowsAtCompileTime, - Cols = Rhs::ColsAtCompileTime, - Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime) + Rows = _Lhs::RowsAtCompileTime, + Cols = _Rhs::ColsAtCompileTime, + Depth = EIGEN_ENUM_MIN(_Lhs::ColsAtCompileTime,_Rhs::RowsAtCompileTime) }; // the splitting into different lines of code here, introducing the _select enums and the typedef below, @@ -211,9 +213,6 @@ class GeneralProduct<Lhs, Rhs, OuterProduct> { ei_outer_product_selector<(int(Dest::Flags)&RowMajorBit) ? RowMajor : ColMajor>::run(*this, dest, alpha); } - - private: - GeneralProduct& operator=(const GeneralProduct&); }; template<> struct ei_outer_product_selector<ColMajor> { @@ -279,9 +278,6 @@ class GeneralProduct<Lhs, Rhs, GemvProduct> ei_gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor, bool(ei_blas_traits<MatrixType>::ActualAccess)>::run(*this, dst, alpha); } - -private: - GeneralProduct& operator=(const GeneralProduct&); }; // The vector is on the left => transposition diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index d835f16b7..3248e7293 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -114,14 +114,6 @@ class ProductBase : public MatrixBase<Derived> template<typename Dest> inline void scaleAndAddTo(Dest& dst,Scalar alpha) const { derived().scaleAndAddTo(dst,alpha); } - PlainMatrixType eval() const - { - PlainMatrixType res(rows(), cols()); - res.setZero(); - derived().evalTo(res); - return res; - } - EIGEN_DEPRECATED const Flagged<ProductBase, 0, EvalBeforeAssigningBit> lazy() const { return *this; } @@ -140,8 +132,6 @@ class ProductBase : public MatrixBase<Derived> void coeffRef(int,int); void coeff(int) const; void coeffRef(int); - - ProductBase& operator=(const ProductBase&); }; template<typename NestedProduct> diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h index b6b28c00b..92522f86c 100644 --- a/Eigen/src/Core/Redux.h +++ b/Eigen/src/Core/Redux.h @@ -54,16 +54,16 @@ private: public: enum { - Vectorization = int(MayLinearVectorize) ? int(LinearVectorization) - : int(MaySliceVectorize) ? int(SliceVectorization) - : int(NoVectorization) + Traversal = int(MayLinearVectorize) ? int(LinearVectorizedTraversal) + : int(MaySliceVectorize) ? int(SliceVectorizedTraversal) + : int(DefaultTraversal) }; private: enum { Cost = Derived::SizeAtCompileTime * Derived::CoeffReadCost + (Derived::SizeAtCompileTime-1) * NumTraits<typename Derived::Scalar>::AddCost, - UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Vectorization) == int(NoVectorization) ? 1 : int(PacketSize)) + UnrollingLimit = EIGEN_UNROLLING_LIMIT * (int(Traversal) == int(DefaultTraversal) ? 1 : int(PacketSize)) }; public: @@ -171,13 +171,13 @@ struct ei_redux_vec_unroller<Func, Derived, Start, 1> ***************************************************************************/ template<typename Func, typename Derived, - int Vectorization = ei_redux_traits<Func, Derived>::Vectorization, + int Traversal = ei_redux_traits<Func, Derived>::Traversal, int Unrolling = ei_redux_traits<Func, Derived>::Unrolling > struct ei_redux_impl; template<typename Func, typename Derived> -struct ei_redux_impl<Func, Derived, NoVectorization, NoUnrolling> +struct ei_redux_impl<Func, Derived, DefaultTraversal, NoUnrolling> { typedef typename Derived::Scalar Scalar; static Scalar run(const Derived& mat, const Func& func) @@ -195,12 +195,12 @@ struct ei_redux_impl<Func, Derived, NoVectorization, NoUnrolling> }; template<typename Func, typename Derived> -struct ei_redux_impl<Func,Derived, NoVectorization, CompleteUnrolling> +struct ei_redux_impl<Func,Derived, DefaultTraversal, CompleteUnrolling> : public ei_redux_novec_unroller<Func,Derived, 0, Derived::SizeAtCompileTime> {}; template<typename Func, typename Derived> -struct ei_redux_impl<Func, Derived, LinearVectorization, NoUnrolling> +struct ei_redux_impl<Func, Derived, LinearVectorizedTraversal, NoUnrolling> { typedef typename Derived::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; @@ -209,10 +209,7 @@ struct ei_redux_impl<Func, Derived, LinearVectorization, NoUnrolling> { const int size = mat.size(); const int packetSize = ei_packet_traits<Scalar>::size; - const int alignedStart = (Derived::Flags & AlignedBit) - || !(Derived::Flags & DirectAccessBit) - ? 0 - : ei_alignmentOffset(&mat.const_cast_derived().coeffRef(0), size); + const int alignedStart = ei_alignmentOffset(mat,size); enum { alignment = (Derived::Flags & DirectAccessBit) || (Derived::Flags & AlignedBit) ? Aligned : Unaligned @@ -246,7 +243,7 @@ struct ei_redux_impl<Func, Derived, LinearVectorization, NoUnrolling> }; template<typename Func, typename Derived> -struct ei_redux_impl<Func, Derived, SliceVectorization, NoUnrolling> +struct ei_redux_impl<Func, Derived, SliceVectorizedTraversal, NoUnrolling> { typedef typename Derived::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; @@ -277,7 +274,7 @@ struct ei_redux_impl<Func, Derived, SliceVectorization, NoUnrolling> else // too small to vectorize anything. // since this is dynamic-size hence inefficient anyway for such small sizes, don't try to optimize. { - res = ei_redux_impl<Func, Derived, NoVectorization, NoUnrolling>::run(mat, func); + res = ei_redux_impl<Func, Derived, DefaultTraversal, NoUnrolling>::run(mat, func); } return res; @@ -285,20 +282,20 @@ struct ei_redux_impl<Func, Derived, SliceVectorization, NoUnrolling> }; template<typename Func, typename Derived> -struct ei_redux_impl<Func, Derived, LinearVectorization, CompleteUnrolling> +struct ei_redux_impl<Func, Derived, LinearVectorizedTraversal, CompleteUnrolling> { typedef typename Derived::Scalar Scalar; typedef typename ei_packet_traits<Scalar>::type PacketScalar; enum { PacketSize = ei_packet_traits<Scalar>::size, Size = Derived::SizeAtCompileTime, - VectorizationSize = (Size / PacketSize) * PacketSize + VectorizedSize = (Size / PacketSize) * PacketSize }; EIGEN_STRONG_INLINE static Scalar run(const Derived& mat, const Func& func) { Scalar res = func.predux(ei_redux_vec_unroller<Func, Derived, 0, Size / PacketSize>::run(mat,func)); - if (VectorizationSize != Size) - res = func(res,ei_redux_novec_unroller<Func, Derived, VectorizationSize, Size-VectorizationSize>::run(mat,func)); + if (VectorizedSize != Size) + res = func(res,ei_redux_novec_unroller<Func, Derived, VectorizedSize, Size-VectorizedSize>::run(mat,func)); return res; } }; diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index 3e435853c..9ab8cb2c1 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -150,7 +150,6 @@ template<typename MatrixType, unsigned int UpLo> class SelfAdjointView const LDLT<PlainMatrixType> ldlt() const; protected: - const typename MatrixType::Nested m_matrix; }; diff --git a/Eigen/src/Core/Swap.h b/Eigen/src/Core/Swap.h index 2af78485b..010d8bb8b 100644 --- a/Eigen/src/Core/Swap.h +++ b/Eigen/src/Core/Swap.h @@ -106,9 +106,6 @@ template<typename ExpressionType> class SwapWrapper protected: ExpressionType& m_expression; - - private: - SwapWrapper& operator=(const SwapWrapper&); }; /** swaps *this with the expression \a other. diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 513ce522a..63d77f134 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -215,7 +215,7 @@ template<typename Derived> inline const typename MatrixBase<Derived>::AdjointReturnType MatrixBase<Derived>::adjoint() const { - return this->transpose().nestByValue(); + return this->transpose(); } /*************************************************************************** @@ -297,13 +297,7 @@ inline void MatrixBase<Derived>::adjointInPlace() #ifndef EIGEN_NO_DEBUG -// The following is to detect aliasing problems in the following common cases: -// a = a.transpose() -// a = a.transpose() + X -// a = X + a.transpose() -// a = a.adjoint() -// a = a.adjoint() + X -// a = X + a.adjoint() +// The following is to detect aliasing problems in most common cases. template<typename T, int Access=ei_blas_traits<T>::ActualAccess> struct ei_extract_data_selector { @@ -323,63 +317,31 @@ template<typename T> typename T::Scalar* ei_extract_data(const T& m) return ei_extract_data_selector<T>::run(m); } -template<typename Derived> -template<typename OtherDerived> -Derived& DenseBase<Derived>::lazyAssign(const Transpose<OtherDerived>& other) -{ - ei_assert(ei_extract_data(other) != ei_extract_data(derived()) - && "aliasing detected during tranposition, please use transposeInPlace()"); - return lazyAssign(static_cast<const DenseBase<Transpose<OtherDerived> >& >(other)); -} - -template<typename Derived> -template<typename DerivedA, typename DerivedB> -Derived& DenseBase<Derived>:: -lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,Transpose<DerivedA>,DerivedB>& other) -{ - ei_assert(ei_extract_data(derived()) != ei_extract_data(other.lhs()) - && "aliasing detected during tranposition, please evaluate your expression"); - return lazyAssign(static_cast<const DenseBase<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,Transpose<DerivedA>,DerivedB> >& >(other)); -} - -template<typename Derived> -template<typename DerivedA, typename DerivedB> -Derived& DenseBase<Derived>:: -lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerivedA,Transpose<DerivedB> >& other) -{ - ei_assert(ei_extract_data(derived()) != ei_extract_data(other.rhs()) - && "aliasing detected during tranposition, please evaluate your expression"); - return lazyAssign(static_cast<const DenseBase<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerivedA,Transpose<DerivedB> > >& >(other)); -} - -template<typename Derived> -template<typename OtherDerived> Derived& -DenseBase<Derived>:: -lazyAssign(const CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<OtherDerived> > >& other) +template<typename Scalar, bool DestIsTranposed, typename OtherDerived> +struct ei_check_transpose_aliasing_selector { - ei_assert(ei_extract_data(other) != ei_extract_data(derived()) - && "aliasing detected during tranposition, please use adjointInPlace()"); - return lazyAssign(static_cast<const DenseBase<CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<OtherDerived> > > >& >(other)); -} + static bool run(const Scalar* dest, const OtherDerived& src) + { + return (ei_blas_traits<OtherDerived>::IsTransposed != DestIsTranposed) && (dest==(Scalar*)ei_extract_data(src)); + } +}; -template<typename Derived> -template<typename DerivedA, typename DerivedB> -Derived& DenseBase<Derived>:: -lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<DerivedA> > >,DerivedB>& other) +template<typename Scalar, bool DestIsTranposed, typename BinOp, typename DerivedA, typename DerivedB> +struct ei_check_transpose_aliasing_selector<Scalar,DestIsTranposed,CwiseBinaryOp<BinOp,DerivedA,DerivedB> > { - ei_assert(ei_extract_data(derived()) != ei_extract_data(other.lhs()) - && "aliasing detected during tranposition, please evaluate your expression"); - return lazyAssign(static_cast<const DenseBase<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<DerivedA> > >,DerivedB> >& >(other)); -} + static bool run(const Scalar* dest, const CwiseBinaryOp<BinOp,DerivedA,DerivedB>& src) + { + return ((ei_blas_traits<DerivedA>::IsTransposed != DestIsTranposed) && dest==(Scalar*)ei_extract_data(src.lhs())) + || ((ei_blas_traits<DerivedB>::IsTransposed != DestIsTranposed) && dest==(Scalar*)ei_extract_data(src.rhs())); + } +}; template<typename Derived> -template<typename DerivedA, typename DerivedB> -Derived& DenseBase<Derived>:: -lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerivedA,CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<DerivedB> > > >& other) +template<typename OtherDerived> +void DenseBase<Derived>::checkTransposeAliasing(const OtherDerived& other) const { - ei_assert(ei_extract_data(derived()) != ei_extract_data(other.rhs()) - && "aliasing detected during tranposition, please evaluate your expression"); - return lazyAssign(static_cast<const DenseBase<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerivedA,CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<DerivedB> > > > >& >(other)); + ei_assert((!ei_check_transpose_aliasing_selector<Scalar,ei_blas_traits<Derived>::IsTransposed,OtherDerived>::run(ei_extract_data(derived()), other)) + && "aliasing detected during tranposition, use transposeInPlace() or evaluate the rhs into a temporary using .eval()"); } #endif diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index e60d57e70..aaf781d1f 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -26,22 +26,11 @@ #ifndef EIGEN_TRIANGULARMATRIX_H #define EIGEN_TRIANGULARMATRIX_H -/** \nonstableyet - * \class TriangularBase - * - * \brief Expression of a triangular matrix extracted from a given matrix - * - * \param MatrixType the type of the object in which we are taking the triangular part - * \param Mode the kind of triangular matrix expression to construct. Can be UpperTriangular, - * LowerTriangular, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field; - * it must have either UpperBit or LowerBit, and additionnaly it may have either - * TraingularBit or SelfadjointBit. +/** \internal * - * This class represents an expression of the upper or lower triangular part of - * a square matrix, possibly with a further assumption on the diagonal. It is the return type - * of MatrixBase::part() and most of the time this is the only way it is used. + * \class TriangularBase * - * \sa MatrixBase::part() + * \brief Base class for triangular part in a matrix */ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived> { @@ -99,11 +88,11 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived> void check_coordinates(int row, int col) { - ei_assert(col>0 && col<cols() && row>0 && row<rows()); + ei_assert(col>=0 && col<cols() && row>=0 && row<rows()); ei_assert( (Mode==UpperTriangular && col>=row) || (Mode==LowerTriangular && col<=row) - || (Mode==StrictlyUpperTriangular && col>row) - || (Mode==StrictlyLowerTriangular && col<row)); + || ((Mode==StrictlyUpperTriangular || Mode==UnitUpperTriangular) && col>row) + || ((Mode==StrictlyLowerTriangular || Mode==UnitLowerTriangular) && col<row)); } void check_coordinates_internal(int row, int col) @@ -115,19 +104,21 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived> }; - /** \class TriangularView - * \nonstableyet * - * \brief Expression of a triangular part of a dense matrix + * \brief Base class for triangular part in a matrix * - * \param MatrixType the type of the dense matrix storing the coefficients + * \param MatrixType the type of the object in which we are taking the triangular part + * \param Mode the kind of triangular matrix expression to construct. Can be UpperTriangular, + * LowerTriangular, UpperSelfadjoint, or LowerSelfadjoint. This is in fact a bit field; + * it must have either UpperBit or LowerBit, and additionnaly it may have either + * TraingularBit or SelfadjointBit. * - * This class is an expression of a triangular part of a matrix with given dense - * storage of the coefficients. It is the return type of MatrixBase::triangularPart() - * and most of the time this is the only way that it is used. + * This class represents a triangular part of a matrix, not necessarily square. Strictly speaking, for rectangular + * matrices one should speak ok "trapezoid" parts. This class is the return type + * of MatrixBase::triangularView() and most of the time this is the only way it is used. * - * \sa class TriangularBase, MatrixBase::triangularPart(), class DiagonalWrapper + * \sa MatrixBase::triangularView() */ template<typename MatrixType, unsigned int _Mode> struct ei_traits<TriangularView<MatrixType, _Mode> > : ei_traits<MatrixType> @@ -155,7 +146,7 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView typedef TriangularBase<TriangularView> Base; typedef typename ei_traits<TriangularView>::Scalar Scalar; typedef _MatrixType MatrixType; - typedef typename MatrixType::PlainMatrixType PlainMatrixType; + typedef typename MatrixType::PlainMatrixType DenseMatrixType; typedef typename MatrixType::Nested MatrixTypeNested; typedef typename ei_cleantype<MatrixTypeNested>::type _MatrixTypeNested; @@ -231,23 +222,23 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView /** \sa MatrixBase::adjoint() */ - inline TriangularView<NestByValue<typename MatrixType::AdjointReturnType>,TransposeMode> adjoint() - { return m_matrix.adjoint().nestByValue(); } + inline TriangularView<typename MatrixType::AdjointReturnType,TransposeMode> adjoint() + { return m_matrix.adjoint(); } /** \sa MatrixBase::adjoint() const */ - inline const TriangularView<NestByValue<typename MatrixType::AdjointReturnType>,TransposeMode> adjoint() const - { return m_matrix.adjoint().nestByValue(); } + inline const TriangularView<typename MatrixType::AdjointReturnType,TransposeMode> adjoint() const + { return m_matrix.adjoint(); } /** \sa MatrixBase::transpose() */ - inline TriangularView<NestByValue<Transpose<MatrixType> >,TransposeMode> transpose() - { return m_matrix.transpose().nestByValue(); } + inline TriangularView<Transpose<MatrixType>,TransposeMode> transpose() + { return m_matrix.transpose(); } /** \sa MatrixBase::transpose() const */ - inline const TriangularView<NestByValue<Transpose<MatrixType> >,TransposeMode> transpose() const - { return m_matrix.transpose().nestByValue(); } + inline const TriangularView<Transpose<MatrixType>,TransposeMode> transpose() const + { return m_matrix.transpose(); } - PlainMatrixType toDense() const + DenseMatrixType toDenseMatrix() const { - PlainMatrixType res(rows(), cols()); - res = *this; + DenseMatrixType res(rows(), cols()); + evalToLazy(res); return res; } @@ -351,20 +342,7 @@ struct ei_triangular_assignment_selector } } }; -template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite> -struct ei_triangular_assignment_selector<Derived1, Derived2, Mode, 1, ClearOpposite> -{ - inline static void run(Derived1 &dst, const Derived2 &src) - { - if(Mode&UnitDiagBit) - { - if(ClearOpposite) - dst.coeffRef(0, 0) = 1; - } - else if(!(Mode & ZeroDiagBit)) - dst.copyCoeff(0, 0, src); - } -}; + // prevent buggy user code from causing an infinite recursion template<typename Derived1, typename Derived2, unsigned int Mode, bool ClearOpposite> struct ei_triangular_assignment_selector<Derived1, Derived2, Mode, 0, ClearOpposite> @@ -379,14 +357,16 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UpperTriangular, Dy { for(int j = 0; j < dst.cols(); ++j) { - for(int i = 0; i <= j; ++i) + int maxi = std::min(j, dst.rows()-1); + for(int i = 0; i <= maxi; ++i) dst.copyCoeff(i, j, src); if (ClearOpposite) - for(int i = j+1; i < dst.rows(); ++i) + for(int i = maxi+1; i < dst.rows(); ++i) dst.coeffRef(i, j) = 0; } } }; + template<typename Derived1, typename Derived2, bool ClearOpposite> struct ei_triangular_assignment_selector<Derived1, Derived2, LowerTriangular, Dynamic, ClearOpposite> { @@ -396,8 +376,9 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, LowerTriangular, Dy { for(int i = j; i < dst.rows(); ++i) dst.copyCoeff(i, j, src); + int maxi = std::min(j, dst.rows()); if (ClearOpposite) - for(int i = 0; i < j; ++i) + for(int i = 0; i < maxi; ++i) dst.coeffRef(i, j) = 0; } } @@ -410,14 +391,16 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyUpperTriang { for(int j = 0; j < dst.cols(); ++j) { - for(int i = 0; i < j; ++i) + int maxi = std::min(j, dst.rows()); + for(int i = 0; i < maxi; ++i) dst.copyCoeff(i, j, src); if (ClearOpposite) - for(int i = j; i < dst.rows(); ++i) + for(int i = maxi; i < dst.rows(); ++i) dst.coeffRef(i, j) = 0; } } }; + template<typename Derived1, typename Derived2, bool ClearOpposite> struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyLowerTriangular, Dynamic, ClearOpposite> { @@ -427,8 +410,9 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, StrictlyLowerTriang { for(int i = j+1; i < dst.rows(); ++i) dst.copyCoeff(i, j, src); + int maxi = std::min(j, dst.rows()-1); if (ClearOpposite) - for(int i = 0; i <= j; ++i) + for(int i = 0; i <= maxi; ++i) dst.coeffRef(i, j) = 0; } } @@ -441,15 +425,16 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UnitUpperTriangular { for(int j = 0; j < dst.cols(); ++j) { - for(int i = 0; i < j; ++i) + int maxi = std::min(j, dst.rows()); + for(int i = 0; i < maxi; ++i) dst.copyCoeff(i, j, src); if (ClearOpposite) { - for(int i = j+1; i < dst.rows(); ++i) + for(int i = maxi+1; i < dst.rows(); ++i) dst.coeffRef(i, j) = 0; - dst.coeffRef(j, j) = 1; } } + dst.diagonal().setOnes(); } }; template<typename Derived1, typename Derived2, bool ClearOpposite> @@ -459,15 +444,16 @@ struct ei_triangular_assignment_selector<Derived1, Derived2, UnitLowerTriangular { for(int j = 0; j < dst.cols(); ++j) { - for(int i = j+1; i < dst.rows(); ++i) + int maxi = std::min(j, dst.rows()); + for(int i = maxi+1; i < dst.rows(); ++i) dst.copyCoeff(i, j, src); if (ClearOpposite) { - for(int i = 0; i < j; ++i) + for(int i = 0; i < maxi; ++i) dst.coeffRef(i, j) = 0; - dst.coeffRef(j, j) = 1; } } + dst.diagonal().setOnes(); } }; @@ -514,7 +500,7 @@ TriangularView<MatrixType, Mode>::operator=(const TriangularBase<OtherDerived>& ei_assert(Mode == OtherDerived::Mode); if(ei_traits<OtherDerived>::Flags & EvalBeforeAssigningBit) { - typename OtherDerived::PlainMatrixType other_evaluated(other.rows(), other.cols()); + typename OtherDerived::DenseMatrixType other_evaluated(other.rows(), other.cols()); other_evaluated.template triangularView<Mode>().lazyAssign(other.derived()); lazyAssign(other_evaluated); } @@ -633,17 +619,20 @@ const TriangularView<Derived, Mode> MatrixBase<Derived>::triangularView() const template<typename Derived> bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const { - if(cols() != rows()) return false; RealScalar maxAbsOnUpperTriangularPart = static_cast<RealScalar>(-1); for(int j = 0; j < cols(); ++j) - for(int i = 0; i <= j; ++i) + { + int maxi = std::min(j, rows()-1); + for(int i = 0; i <= maxi; ++i) { RealScalar absValue = ei_abs(coeff(i,j)); if(absValue > maxAbsOnUpperTriangularPart) maxAbsOnUpperTriangularPart = absValue; } - for(int j = 0; j < cols()-1; ++j) + } + RealScalar threshold = maxAbsOnUpperTriangularPart * prec; + for(int j = 0; j < cols(); ++j) for(int i = j+1; i < rows(); ++i) - if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnUpperTriangularPart, prec)) return false; + if(ei_abs(coeff(i, j)) > threshold) return false; return true; } @@ -655,7 +644,6 @@ bool MatrixBase<Derived>::isUpperTriangular(RealScalar prec) const template<typename Derived> bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const { - if(cols() != rows()) return false; RealScalar maxAbsOnLowerTriangularPart = static_cast<RealScalar>(-1); for(int j = 0; j < cols(); ++j) for(int i = j; i < rows(); ++i) @@ -663,9 +651,13 @@ bool MatrixBase<Derived>::isLowerTriangular(RealScalar prec) const RealScalar absValue = ei_abs(coeff(i,j)); if(absValue > maxAbsOnLowerTriangularPart) maxAbsOnLowerTriangularPart = absValue; } + RealScalar threshold = maxAbsOnLowerTriangularPart * prec; for(int j = 1; j < cols(); ++j) - for(int i = 0; i < j; ++i) - if(!ei_isMuchSmallerThan(coeff(i, j), maxAbsOnLowerTriangularPart, prec)) return false; + { + int maxi = std::min(j, rows()-1); + for(int i = 0; i < maxi; ++i) + if(ei_abs(coeff(i, j)) > threshold) return false; + } return true; } diff --git a/Eigen/src/Core/arch/SSE/MathFunctions.h b/Eigen/src/Core/arch/SSE/MathFunctions.h index 8399aac30..3c0020248 100644 --- a/Eigen/src/Core/arch/SSE/MathFunctions.h +++ b/Eigen/src/Core/arch/SSE/MathFunctions.h @@ -365,6 +365,8 @@ static EIGEN_DONT_INLINE EIGEN_UNUSED Packet4f ei_pcos(Packet4f x) return _mm_xor_ps(y, sign_bit); } +// This is Quake3's fast inverse square root. +// For detail see here: http://www.beyond3d.com/content/articles/8/ static EIGEN_UNUSED Packet4f ei_psqrt(Packet4f _x) { Packet4f half = ei_pmul(_x, ei_pset1(.5f)); diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index 29c89c310..a633c7b7c 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -58,8 +58,8 @@ template<> struct ei_packet_traits<float> : ei_default_packet_traits { typedef Packet4f type; enum {size=4}; enum { - HasSin = 1, - HasCos = 1, + HasSin = EIGEN_FAST_MATH, + HasCos = EIGEN_FAST_MATH, HasLog = 1, HasExp = 1, HasSqrt = 1 @@ -118,6 +118,9 @@ template<> EIGEN_STRONG_INLINE Packet4f ei_pmul<Packet4f>(const Packet4f& a, con template<> EIGEN_STRONG_INLINE Packet2d ei_pmul<Packet2d>(const Packet2d& a, const Packet2d& b) { return _mm_mul_pd(a,b); } template<> EIGEN_STRONG_INLINE Packet4i ei_pmul<Packet4i>(const Packet4i& a, const Packet4i& b) { +#ifdef __SSE4_1__ + return _mm_mullo_epi32(a,b); +#else // this version is slightly faster than 4 scalar products return ei_vec4i_swizzle1( ei_vec4i_swizzle2( @@ -126,6 +129,7 @@ template<> EIGEN_STRONG_INLINE Packet4i ei_pmul<Packet4i>(const Packet4i& a, con ei_vec4i_swizzle1(b,1,0,3,2)), 0,2,0,2), 0,2,1,3); +#endif } template<> EIGEN_STRONG_INLINE Packet4f ei_pdiv<Packet4f>(const Packet4f& a, const Packet4f& b) { return _mm_div_ps(a,b); } diff --git a/Eigen/src/Core/products/GeneralUnrolled.h b/Eigen/src/Core/products/GeneralUnrolled.h index 7e2eebf08..f04c27a95 100644 --- a/Eigen/src/Core/products/GeneralUnrolled.h +++ b/Eigen/src/Core/products/GeneralUnrolled.h @@ -36,7 +36,7 @@ * Note that here the inner-loops should always be unrolled. */ -template<int VectorizationMode, int Index, typename Lhs, typename Rhs, typename RetScalar> +template<int Traversal, int Index, typename Lhs, typename Rhs, typename RetScalar> struct ei_product_coeff_impl; template<int StorageOrder, int Index, typename Lhs, typename Rhs, typename PacketScalar, int LoadMode> @@ -118,7 +118,7 @@ template<typename LhsNested, typename RhsNested> class GeneralProduct<LhsNested, CanVectorizeInner = ei_traits<GeneralProduct>::CanVectorizeInner }; - typedef ei_product_coeff_impl<CanVectorizeInner ? InnerVectorization : NoVectorization, + typedef ei_product_coeff_impl<CanVectorizeInner ? InnerVectorizedTraversal : DefaultTraversal, Unroll ? InnerSize-1 : Dynamic, _LhsNested, _RhsNested, Scalar> ScalarCoeffImpl; @@ -185,17 +185,17 @@ template<typename LhsNested, typename RhsNested> class GeneralProduct<LhsNested, **************************************/ template<int Index, typename Lhs, typename Rhs, typename RetScalar> -struct ei_product_coeff_impl<NoVectorization, Index, Lhs, Rhs, RetScalar> +struct ei_product_coeff_impl<DefaultTraversal, Index, Lhs, Rhs, RetScalar> { EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, RetScalar &res) { - ei_product_coeff_impl<NoVectorization, Index-1, Lhs, Rhs, RetScalar>::run(row, col, lhs, rhs, res); + ei_product_coeff_impl<DefaultTraversal, Index-1, Lhs, Rhs, RetScalar>::run(row, col, lhs, rhs, res); res += lhs.coeff(row, Index) * rhs.coeff(Index, col); } }; template<typename Lhs, typename Rhs, typename RetScalar> -struct ei_product_coeff_impl<NoVectorization, 0, Lhs, Rhs, RetScalar> +struct ei_product_coeff_impl<DefaultTraversal, 0, Lhs, Rhs, RetScalar> { EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, RetScalar &res) { @@ -204,7 +204,7 @@ struct ei_product_coeff_impl<NoVectorization, 0, Lhs, Rhs, RetScalar> }; template<typename Lhs, typename Rhs, typename RetScalar> -struct ei_product_coeff_impl<NoVectorization, Dynamic, Lhs, Rhs, RetScalar> +struct ei_product_coeff_impl<DefaultTraversal, Dynamic, Lhs, Rhs, RetScalar> { EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, RetScalar& res) { @@ -217,7 +217,7 @@ struct ei_product_coeff_impl<NoVectorization, Dynamic, Lhs, Rhs, RetScalar> // prevent buggy user code from causing an infinite recursion template<typename Lhs, typename Rhs, typename RetScalar> -struct ei_product_coeff_impl<NoVectorization, -1, Lhs, Rhs, RetScalar> +struct ei_product_coeff_impl<DefaultTraversal, -1, Lhs, Rhs, RetScalar> { EIGEN_STRONG_INLINE static void run(int, int, const Lhs&, const Rhs&, RetScalar&) {} }; @@ -247,7 +247,7 @@ struct ei_product_coeff_vectorized_unroller<0, Lhs, Rhs, PacketScalar> }; template<int Index, typename Lhs, typename Rhs, typename RetScalar> -struct ei_product_coeff_impl<InnerVectorization, Index, Lhs, Rhs, RetScalar> +struct ei_product_coeff_impl<InnerVectorizedTraversal, Index, Lhs, Rhs, RetScalar> { typedef typename Lhs::PacketScalar PacketScalar; enum { PacketSize = ei_packet_traits<typename Lhs::Scalar>::size }; @@ -255,7 +255,7 @@ struct ei_product_coeff_impl<InnerVectorization, Index, Lhs, Rhs, RetScalar> { PacketScalar pres; ei_product_coeff_vectorized_unroller<Index+1-PacketSize, Lhs, Rhs, PacketScalar>::run(row, col, lhs, rhs, pres); - ei_product_coeff_impl<NoVectorization,Index,Lhs,Rhs,RetScalar>::run(row, col, lhs, rhs, res); + ei_product_coeff_impl<DefaultTraversal,Index,Lhs,Rhs,RetScalar>::run(row, col, lhs, rhs, res); res = ei_predux(pres); } }; @@ -268,7 +268,7 @@ struct ei_product_coeff_vectorized_dyn_selector res = ei_dot_impl< Block<Lhs, 1, ei_traits<Lhs>::ColsAtCompileTime>, Block<Rhs, ei_traits<Rhs>::RowsAtCompileTime, 1>, - LinearVectorization, NoUnrolling>::run(lhs.row(row), rhs.col(col)); + LinearVectorizedTraversal, NoUnrolling>::run(lhs.row(row), rhs.col(col)); } }; @@ -282,7 +282,7 @@ struct ei_product_coeff_vectorized_dyn_selector<Lhs,Rhs,1,RhsCols> res = ei_dot_impl< Lhs, Block<Rhs, ei_traits<Rhs>::RowsAtCompileTime, 1>, - LinearVectorization, NoUnrolling>::run(lhs, rhs.col(col)); + LinearVectorizedTraversal, NoUnrolling>::run(lhs, rhs.col(col)); } }; @@ -294,7 +294,7 @@ struct ei_product_coeff_vectorized_dyn_selector<Lhs,Rhs,LhsRows,1> res = ei_dot_impl< Block<Lhs, 1, ei_traits<Lhs>::ColsAtCompileTime>, Rhs, - LinearVectorization, NoUnrolling>::run(lhs.row(row), rhs); + LinearVectorizedTraversal, NoUnrolling>::run(lhs.row(row), rhs); } }; @@ -306,12 +306,12 @@ struct ei_product_coeff_vectorized_dyn_selector<Lhs,Rhs,1,1> res = ei_dot_impl< Lhs, Rhs, - LinearVectorization, NoUnrolling>::run(lhs, rhs); + LinearVectorizedTraversal, NoUnrolling>::run(lhs, rhs); } }; template<typename Lhs, typename Rhs, typename RetScalar> -struct ei_product_coeff_impl<InnerVectorization, Dynamic, Lhs, Rhs, RetScalar> +struct ei_product_coeff_impl<InnerVectorizedTraversal, Dynamic, Lhs, Rhs, RetScalar> { EIGEN_STRONG_INLINE static void run(int row, int col, const Lhs& lhs, const Rhs& rhs, typename Lhs::Scalar &res) { diff --git a/Eigen/src/Core/util/BlasUtil.h b/Eigen/src/Core/util/BlasUtil.h index 94154108c..fa21ceebb 100644 --- a/Eigen/src/Core/util/BlasUtil.h +++ b/Eigen/src/Core/util/BlasUtil.h @@ -160,6 +160,7 @@ template<typename XprType> struct ei_blas_traits typedef XprType _ExtractType; enum { IsComplex = NumTraits<Scalar>::IsComplex, + IsTransposed = false, NeedToConjugate = false, ActualAccess = int(ei_traits<XprType>::Flags)&DirectAccessBit ? HasDirectAccess : NoDirectAccess }; @@ -214,20 +215,6 @@ struct ei_blas_traits<CwiseUnaryOp<ei_scalar_opposite_op<Scalar>, NestedXpr> > { return - Base::extractScalarFactor(x._expression()); } }; -// pop NestByValue -template<typename NestedXpr> -struct ei_blas_traits<NestByValue<NestedXpr> > - : ei_blas_traits<NestedXpr> -{ - typedef typename NestedXpr::Scalar Scalar; - typedef ei_blas_traits<NestedXpr> Base; - typedef NestByValue<NestedXpr> XprType; - typedef typename Base::ExtractType ExtractType; - static inline ExtractType extract(const XprType& x) { return Base::extract(static_cast<const NestedXpr&>(x)); } - static inline Scalar extractScalarFactor(const XprType& x) - { return Base::extractScalarFactor(static_cast<const NestedXpr&>(x)); } -}; - // pop/push transpose template<typename NestedXpr> struct ei_blas_traits<Transpose<NestedXpr> > @@ -241,6 +228,9 @@ struct ei_blas_traits<Transpose<NestedXpr> > ExtractType, typename ExtractType::PlainMatrixType >::ret DirectLinearAccessType; + enum { + IsTransposed = Base::IsTransposed ? 0 : 1 + }; static inline const ExtractType extract(const XprType& x) { return Base::extract(x._expression()); } static inline Scalar extractScalarFactor(const XprType& x) { return Base::extractScalarFactor(x._expression()); } }; diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index beb19cdaa..8483772df 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -195,16 +195,19 @@ enum DirectionType { Vertical, Horizontal, BothDirections }; enum ProductEvaluationMode { NormalProduct, CacheFriendlyProduct }; enum { + /** \internal Default traversal, no vectorization, no index-based access */ + DefaultTraversal, + /** \internal No vectorization, use index-based access to have only one for loop instead of 2 nested loops */ + LinearTraversal, /** \internal Equivalent to a slice vectorization for fixed-size matrices having good alignment * and good size */ - InnerVectorization, + InnerVectorizedTraversal, /** \internal Vectorization path using a single loop plus scalar loops for the * unaligned boundaries */ - LinearVectorization, + LinearVectorizedTraversal, /** \internal Generic vectorization path using one vectorized loop per row/column with some * scalar loops to handle the unaligned boundaries */ - SliceVectorization, - NoVectorization + SliceVectorizedTraversal }; enum { @@ -218,8 +221,7 @@ enum { RowMajor = 0x1, // it is only a coincidence that this is equal to RowMajorBit -- don't rely on that /** \internal Align the matrix itself if it is vectorizable fixed-size */ AutoAlign = 0, - /** \internal Don't require alignment for the matrix itself (the array of coefficients, if dynamically allocated, may still be - requested to be aligned) */ + /** \internal Don't require alignment for the matrix itself (the array of coefficients, if dynamically allocated, may still be requested to be aligned) */ // FIXME --- clarify the situation DontAlign = 0x2 }; @@ -267,19 +269,23 @@ enum TransformTraits { Projective = 0x20 }; -const int EiArch_Generic = 0x0; -const int EiArch_SSE = 0x1; -const int EiArch_AltiVec = 0x2; - -enum DenseStorageMatrix {}; -enum DenseStorageArray {}; - +namespace Architecture +{ + enum Type { + Generic = 0x0, + SSE = 0x1, + AltiVec = 0x2, #if defined EIGEN_VECTORIZE_SSE - const int EiArch = EiArch_SSE; + Target = SSE #elif defined EIGEN_VECTORIZE_ALTIVEC - const int EiArch = EiArch_AltiVec; + Target = AltiVec #else - const int EiArch = EiArch_Generic; + Target = Generic #endif + }; +} + +enum DenseStorageMatrix {}; +enum DenseStorageArray {}; #endif // EIGEN_CONSTANTS_H diff --git a/Eigen/src/Core/util/DisableMSVCWarnings.h b/Eigen/src/Core/util/DisableMSVCWarnings.h index c08d0389f..dcc71143d 100644 --- a/Eigen/src/Core/util/DisableMSVCWarnings.h +++ b/Eigen/src/Core/util/DisableMSVCWarnings.h @@ -1,6 +1,9 @@ #ifdef _MSC_VER - // 4273 - QtAlignedMalloc, inconsistent dll linkage + // 4273 - QtAlignedMalloc, inconsistent DLL linkage + // 4100 - unreferenced formal parameter (occurred e.g. in aligned_allocator::destroy(pointer p)) + // 4101 - unreferenced local variable + // 4512 - assignment operator could not be generated #pragma warning( push ) - #pragma warning( disable : 4181 4244 4127 4211 4273 4522 4717 ) + #pragma warning( disable : 4100 4101 4181 4244 4127 4211 4273 4512 4522 4717 ) #endif diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index b3a1955fb..4ea5013f7 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -29,8 +29,8 @@ #undef minor #define EIGEN_WORLD_VERSION 2 -#define EIGEN_MAJOR_VERSION 90 -#define EIGEN_MINOR_VERSION 1 +#define EIGEN_MAJOR_VERSION 91 +#define EIGEN_MINOR_VERSION 0 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index bc5235582..524bec2fc 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -209,16 +209,20 @@ template<typename T, bool Align> inline void ei_conditional_aligned_delete(T *pt ei_conditional_aligned_free<Align>(ptr); } -/** \internal \returns the number of elements which have to be skipped such that data are 16 bytes aligned */ -template<typename Scalar> -inline static int ei_alignmentOffset(const Scalar* ptr, int maxOffset) +/** \internal \returns the number of elements which have to be skipped to + * find the first 16-byte aligned element + * + * There is also the variant ei_alignmentOffset(const MatrixBase&, Integer) defined in Coeffs.h. + */ +template<typename Scalar, typename Integer> +inline static Integer ei_alignmentOffset(const Scalar* ptr, Integer maxOffset) { typedef typename ei_packet_traits<Scalar>::type Packet; - const int PacketSize = ei_packet_traits<Scalar>::size; - const int PacketAlignedMask = PacketSize-1; + const Integer PacketSize = ei_packet_traits<Scalar>::size; + const Integer PacketAlignedMask = PacketSize-1; const bool Vectorized = PacketSize>1; return Vectorized - ? std::min<int>( (PacketSize - (int((size_t(ptr)/sizeof(Scalar))) & PacketAlignedMask)) + ? std::min<Integer>( (PacketSize - (Integer((size_t(ptr)/sizeof(Scalar))) & PacketAlignedMask)) & PacketAlignedMask, maxOffset) : 0; } diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 6d6540acc..6d4a5c7bc 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -28,10 +28,11 @@ // just a workaround because GCC seems to not really like empty structs #ifdef __GNUG__ - struct ei_empty_struct{char _ei_dummy_;}; - #define EIGEN_EMPTY_STRUCT : Eigen::ei_empty_struct + #define EIGEN_EMPTY_STRUCT_CTOR(X) \ + EIGEN_STRONG_INLINE X() {} \ + EIGEN_STRONG_INLINE X(const X&) {} #else - #define EIGEN_EMPTY_STRUCT + #define EIGEN_EMPTY_STRUCT_CTOR(X) #endif //classes inheriting ei_no_assignment_operator don't generate a default operator=. @@ -45,10 +46,10 @@ class ei_no_assignment_operator * can be accessed using value() and setValue(). * Otherwise, this class is an empty structure and value() just returns the template parameter Value. */ -template<int Value> class ei_int_if_dynamic EIGEN_EMPTY_STRUCT +template<int Value> class ei_int_if_dynamic { public: - ei_int_if_dynamic() {} + EIGEN_EMPTY_STRUCT_CTOR(ei_int_if_dynamic) explicit ei_int_if_dynamic(int) {} static int value() { return Value; } void setValue(int) {} @@ -214,8 +215,35 @@ template<typename T> struct ei_plain_matrix_type_row_major > type; }; +// we should be able to get rid of this one too template<typename T> struct ei_must_nest_by_value { enum { ret = false }; }; -template<typename T> struct ei_must_nest_by_value<NestByValue<T> > { enum { ret = true }; }; + +/** +* The reference selector for template expressions. The idea is that we don't +* need to use references for expressions since they are light weight proxy +* objects which should generate no copying overhead. +**/ +template <typename T> +struct ei_ref_selector +{ + typedef T type; +}; + +/** +* Matrices on the other hand side should only be copied, when it is sure +* we gain by copying (see arithmetic cost check and eval before nesting flag). +* Note: This is an optimization measure that comprises potential (though little) +* to create erroneous code. Any user, utilizing ei_nested outside of +* Eigen needs to take care that no references to temporaries are +* stored or that this potential danger is at least communicated +* to the user. +**/ +template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols> +struct ei_ref_selector< Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> > +{ + typedef Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> MatrixType; + typedef MatrixType const& type; +}; /** \internal Determines how a given expression should be nested into another one. * For example, when you do a * (b+c), Eigen will determine how the expression b+c should be @@ -241,15 +269,12 @@ template<typename T, int n=1, typename PlainMatrixType = typename ei_eval<T>::ty CostEval = (n+1) * int(NumTraits<typename ei_traits<T>::Scalar>::ReadCost), CostNoEval = (n-1) * int(ei_traits<T>::CoeffReadCost) }; + typedef typename ei_meta_if< - ei_must_nest_by_value<T>::ret, - T, - typename ei_meta_if< - (int(ei_traits<T>::Flags) & EvalBeforeNestingBit) - || ( int(CostEval) <= int(CostNoEval) ), + ( int(ei_traits<T>::Flags) & EvalBeforeNestingBit ) || + ( int(CostEval) <= int(CostNoEval) ), PlainMatrixType, - const T& - >::ret + typename ei_ref_selector<T>::type >::ret type; }; @@ -302,7 +327,7 @@ template<typename ExpressionType> struct HNormalizedReturnType { ei_traits<ExpressionType>::ColsAtCompileTime==1 ? SizeMinusOne : 1, ei_traits<ExpressionType>::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne; typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<ExpressionType>::Scalar>, - NestByValue<StartMinusOne> > Type; + StartMinusOne > Type; }; template<typename XprType, typename CastType> struct ei_cast_return_type diff --git a/Eigen/src/Eigen2Support/Cwise.h b/Eigen/src/Eigen2Support/Cwise.h index c8d470228..ab7056a1d 100644 --- a/Eigen/src/Eigen2Support/Cwise.h +++ b/Eigen/src/Eigen2Support/Cwise.h @@ -40,7 +40,7 @@ * convenient macro to defined the return type of a cwise comparison to a scalar */ #define EIGEN_CWISE_COMP_TO_SCALAR_RETURN_TYPE(OP) \ CwiseBinaryOp<OP<typename ei_traits<ExpressionType>::Scalar>, ExpressionType, \ - NestByValue<typename ExpressionType::ConstantReturnType> > + typename ExpressionType::ConstantReturnType > /** \class Cwise * @@ -166,9 +166,6 @@ template<typename ExpressionType> class Cwise protected: ExpressionTypeNested m_matrix; - - private: - Cwise& operator=(const Cwise&); }; diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index cce76b288..9f7df49bc 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -58,10 +58,10 @@ template<typename _MatrixType> class HessenbergDecomposition typedef Matrix<RealScalar, Size, 1> DiagonalType; typedef Matrix<RealScalar, SizeMinusOne, 1> SubDiagonalType; - typedef typename NestByValue<Diagonal<MatrixType,0> >::RealReturnType DiagonalReturnType; + typedef typename Diagonal<MatrixType,0>::RealReturnType DiagonalReturnType; - typedef typename NestByValue<Diagonal< - NestByValue<Block<MatrixType,SizeMinusOne,SizeMinusOne> >,0 > >::RealReturnType SubDiagonalReturnType; + typedef typename Diagonal< + Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >::RealReturnType SubDiagonalReturnType; /** This constructor initializes a HessenbergDecomposition object for * further use with HessenbergDecomposition::compute() diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index 04b9f72d7..d8dcfb047 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -61,15 +61,15 @@ template<typename _MatrixType> class Tridiagonalization typedef Matrix<RealScalar, SizeMinusOne, 1> SubDiagonalType; typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex, - typename NestByValue<Diagonal<MatrixType,0> >::RealReturnType, + typename Diagonal<MatrixType,0>::RealReturnType, Diagonal<MatrixType,0> >::ret DiagonalReturnType; typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex, - typename NestByValue<Diagonal< - NestByValue<Block<MatrixType,SizeMinusOne,SizeMinusOne> >,0 > >::RealReturnType, + typename Diagonal< + Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 >::RealReturnType, Diagonal< - NestByValue<Block<MatrixType,SizeMinusOne,SizeMinusOne> >,0 > + Block<MatrixType,SizeMinusOne,SizeMinusOne>,0 > >::ret SubDiagonalReturnType; /** This constructor initializes a Tridiagonalization object for @@ -144,7 +144,7 @@ template<typename MatrixType> const typename Tridiagonalization<MatrixType>::DiagonalReturnType Tridiagonalization<MatrixType>::diagonal(void) const { - return m_matrix.diagonal().nestByValue(); + return m_matrix.diagonal(); } /** \returns an expression of the sub-diagonal vector */ @@ -153,8 +153,7 @@ const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType Tridiagonalization<MatrixType>::subDiagonal(void) const { int n = m_matrix.rows(); - return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1) - .nestByValue().diagonal().nestByValue(); + return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal(); } /** constructs and returns the tridiagonal matrix T. diff --git a/Eigen/src/Geometry/AlignedBox.h b/Eigen/src/Geometry/AlignedBox.h index 707cfeb4f..d47adb7ab 100644 --- a/Eigen/src/Geometry/AlignedBox.h +++ b/Eigen/src/Geometry/AlignedBox.h @@ -171,7 +171,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const AlignedBox& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + bool isApprox(const AlignedBox& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const { return m_min.isApprox(other.m_min, prec) && m_max.isApprox(other.m_max, prec); } protected: diff --git a/Eigen/src/Geometry/AngleAxis.h b/Eigen/src/Geometry/AngleAxis.h index d948ea40b..da698bbfe 100644 --- a/Eigen/src/Geometry/AngleAxis.h +++ b/Eigen/src/Geometry/AngleAxis.h @@ -146,7 +146,7 @@ public: * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const AngleAxis& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + bool isApprox(const AngleAxis& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const { return m_axis.isApprox(other.m_axis, prec) && ei_isApprox(m_angle,other.m_angle, prec); } }; @@ -165,7 +165,7 @@ template<typename QuatDerived> AngleAxis<Scalar>& AngleAxis<Scalar>::operator=(const QuaternionBase<QuatDerived>& q) { Scalar n2 = q.vec().squaredNorm(); - if (n2 < precision<Scalar>()*precision<Scalar>()) + if (n2 < dummy_precision<Scalar>()*dummy_precision<Scalar>()) { m_angle = 0; m_axis << 1, 0, 0; diff --git a/Eigen/src/Geometry/EulerAngles.h b/Eigen/src/Geometry/EulerAngles.h index dbcc7ae89..b6dbf8ae9 100644 --- a/Eigen/src/Geometry/EulerAngles.h +++ b/Eigen/src/Geometry/EulerAngles.h @@ -50,7 +50,7 @@ MatrixBase<Derived>::eulerAngles(int a0, int a1, int a2) const Matrix<Scalar,3,1> res; typedef Matrix<typename Derived::Scalar,2,1> Vector2; - const Scalar epsilon = precision<Scalar>(); + const Scalar epsilon = dummy_precision<Scalar>(); const int odd = ((a0+1)%3 == a1) ? 0 : 1; const int i = a0; diff --git a/Eigen/src/Geometry/Homogeneous.h b/Eigen/src/Geometry/Homogeneous.h index 6cf6a381d..a601e29cf 100644 --- a/Eigen/src/Geometry/Homogeneous.h +++ b/Eigen/src/Geometry/Homogeneous.h @@ -176,7 +176,7 @@ MatrixBase<Derived>::hnormalized() const EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived); return StartMinusOne(derived(),0,0, ColsAtCompileTime==1?size()-1:1, - ColsAtCompileTime==1?1:size()-1).nestByValue() / coeff(size()-1); + ColsAtCompileTime==1?1:size()-1) / coeff(size()-1); } /** \geometry_module @@ -193,8 +193,7 @@ VectorwiseOp<ExpressionType,Direction>::hnormalized() const { return HNormalized_Block(_expression(),0,0, Direction==Vertical ? _expression().rows()-1 : _expression().rows(), - Direction==Horizontal ? _expression().cols()-1 : _expression().cols()).nestByValue() - .cwiseQuotient( + Direction==Horizontal ? _expression().cols()-1 : _expression().cols()).cwiseQuotient( Replicate<NestByValue<HNormalized_Factors>, Direction==Vertical ? HNormalized_SizeMinusOne : 1, Direction==Horizontal ? HNormalized_SizeMinusOne : 1> @@ -202,9 +201,9 @@ VectorwiseOp<ExpressionType,Direction>::hnormalized() const Direction==Vertical ? _expression().rows()-1:0, Direction==Horizontal ? _expression().cols()-1:0, Direction==Vertical ? 1 : _expression().rows(), - Direction==Horizontal ? 1 : _expression().cols()).nestByValue(), + Direction==Horizontal ? 1 : _expression().cols()), Direction==Vertical ? _expression().rows()-1 : 1, - Direction==Horizontal ? _expression().cols()-1 : 1).nestByValue()); + Direction==Horizontal ? _expression().cols()-1 : 1)); } template<typename MatrixType,typename Lhs> @@ -281,7 +280,6 @@ struct ei_homogeneous_right_product_impl<Homogeneous<MatrixType,Horizontal>,Rhs> const typename MatrixType::Nested m_lhs; const typename Rhs::Nested m_rhs; - }; #endif // EIGEN_HOMOGENEOUS_H diff --git a/Eigen/src/Geometry/Hyperplane.h b/Eigen/src/Geometry/Hyperplane.h index ab65ca2ae..aab3d5b35 100644 --- a/Eigen/src/Geometry/Hyperplane.h +++ b/Eigen/src/Geometry/Hyperplane.h @@ -52,9 +52,9 @@ public: typedef _Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; typedef Matrix<Scalar,AmbientDimAtCompileTime,1> VectorType; - typedef Matrix<Scalar,AmbientDimAtCompileTime==Dynamic + typedef Matrix<Scalar,int(AmbientDimAtCompileTime)==Dynamic ? Dynamic - : AmbientDimAtCompileTime+1,1> Coefficients; + : int(AmbientDimAtCompileTime)+1,1> Coefficients; typedef Block<Coefficients,AmbientDimAtCompileTime,1> NormalReturnType; /** Default constructor without initialization */ @@ -257,7 +257,7 @@ public: * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const Hyperplane& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + bool isApprox(const Hyperplane& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const { return m_coeffs.isApprox(other.m_coeffs, prec); } protected: diff --git a/Eigen/src/Geometry/OrthoMethods.h b/Eigen/src/Geometry/OrthoMethods.h index 631b82bf5..79baeadd5 100644 --- a/Eigen/src/Geometry/OrthoMethods.h +++ b/Eigen/src/Geometry/OrthoMethods.h @@ -90,8 +90,9 @@ MatrixBase<Derived>::cross3(const MatrixBase<OtherDerived>& other) const const DerivedNested lhs(derived()); const OtherDerivedNested rhs(other.derived()); - return ei_cross3_impl<EiArch,typename ei_cleantype<DerivedNested>::type, - typename ei_cleantype<OtherDerivedNested>::type>::run(lhs,rhs); + return ei_cross3_impl<Architecture::Target, + typename ei_cleantype<DerivedNested>::type, + typename ei_cleantype<OtherDerivedNested>::type>::run(lhs,rhs); } /** \returns a matrix expression of the cross product of each column or row diff --git a/Eigen/src/Geometry/ParametrizedLine.h b/Eigen/src/Geometry/ParametrizedLine.h index 6a4fcb1c5..21a5595b9 100644 --- a/Eigen/src/Geometry/ParametrizedLine.h +++ b/Eigen/src/Geometry/ParametrizedLine.h @@ -123,7 +123,7 @@ public: * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const ParametrizedLine& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + bool isApprox(const ParametrizedLine& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const { return m_origin.isApprox(other.m_origin, prec) && m_direction.isApprox(other.m_direction, prec); } protected: diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 660e10d1e..861eff19c 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -152,9 +152,9 @@ public: /** \returns the conjugated quaternion */ Quaternion<Scalar> conjugate() const; - /** \returns an interpolation for a constant motion between \a other and \c *this + /** \returns an interpolation for a constant motion between \a other and \c *this * \a t in [0;1] - * see http://en.wikipedia.org/wiki/Slerp + * see http://en.wikipedia.org/wiki/Slerp */ template<class OtherDerived> Quaternion<Scalar> slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const; @@ -163,7 +163,7 @@ public: * * \sa MatrixBase::isApprox() */ template<class OtherDerived> - bool isApprox(const QuaternionBase<OtherDerived>& other, RealScalar prec = precision<Scalar>()) const + bool isApprox(const QuaternionBase<OtherDerived>& other, RealScalar prec = dummy_precision<Scalar>()) const { return coeffs().isApprox(other.coeffs(), prec); } /** return the result vector of \a v through the rotation*/ @@ -221,7 +221,7 @@ struct ei_traits<Quaternion<_Scalar> > template<typename _Scalar> class Quaternion : public QuaternionBase<Quaternion<_Scalar> >{ typedef QuaternionBase<Quaternion<_Scalar> > Base; -public: +public: typedef _Scalar Scalar; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Quaternion<Scalar>) @@ -304,27 +304,20 @@ template<typename _Scalar, int PacketAccess> class Map<Quaternion<_Scalar>, PacketAccess > : public QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> > { - public: - typedef _Scalar Scalar; - typedef Map<Quaternion<Scalar>, PacketAccess > MapQuat; - - private: - Map<Quaternion<Scalar>, PacketAccess >(); - Map<Quaternion<Scalar>, PacketAccess >(const Map<Quaternion<Scalar>, PacketAccess>&); - typedef QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> > Base; + public: - EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(MapQuat) + typedef _Scalar Scalar; + typedef typename ei_traits<Map>::Coefficients Coefficients; + EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Map) using Base::operator*=; - typedef typename ei_traits<Map<Quaternion<Scalar>, PacketAccess> >::Coefficients Coefficients; - /** Constructs a Mapped Quaternion object from the pointer \a coeffs * * The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order: * \code *coeffs == {x, y, z, w} \endcode * - * If the template paramter PacketAccess is set to Aligned, then the pointer coeffs must be aligned. */ + * If the template parameter PacketAccess is set to Aligned, then the pointer coeffs must be aligned. */ EIGEN_STRONG_INLINE Map(const Scalar* coeffs) : m_coeffs(coeffs) {} inline Coefficients& coeffs() { return m_coeffs;} @@ -374,7 +367,7 @@ QuaternionBase<Derived>::operator* (const QuaternionBase<OtherDerived>& other) c { EIGEN_STATIC_ASSERT((ei_is_same_type<typename Derived::Scalar, typename OtherDerived::Scalar>::ret), YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) - return ei_quat_product<EiArch, Derived, OtherDerived, + return ei_quat_product<Architecture::Target, Derived, OtherDerived, typename ei_traits<Derived>::Scalar, ei_traits<Derived>::PacketAccess && ei_traits<OtherDerived>::PacketAccess>::run(*this, other); } @@ -514,7 +507,7 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri // under the constraint: // ||x|| = 1 // which yields a singular value problem - if (c < Scalar(-1)+precision<Scalar>()) + if (c < Scalar(-1)+dummy_precision<Scalar>()) { c = std::max<Scalar>(c,-1); Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose(); @@ -590,20 +583,29 @@ template <class OtherDerived> Quaternion<typename ei_traits<Derived>::Scalar> QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const { - static const Scalar one = Scalar(1) - precision<Scalar>(); + static const Scalar one = Scalar(1) - epsilon<Scalar>(); Scalar d = this->dot(other); Scalar absD = ei_abs(d); - if (absD>=one) - return Quaternion<Scalar>(derived()); - // theta is the angle between the 2 quaternions - Scalar theta = std::acos(absD); - Scalar sinTheta = ei_sin(theta); + Scalar scale0; + Scalar scale1; - Scalar scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta; - Scalar scale1 = ei_sin( ( t * theta) ) / sinTheta; - if (d<0) - scale1 = -scale1; + if (absD>=one) + { + scale0 = Scalar(1) - t; + scale1 = t; + } + else + { + // theta is the angle between the 2 quaternions + Scalar theta = std::acos(absD); + Scalar sinTheta = ei_sin(theta); + + scale0 = ei_sin( ( Scalar(1) - t ) * theta) / sinTheta; + scale1 = ei_sin( ( t * theta) ) / sinTheta; + if (d<0) + scale1 = -scale1; + } return Quaternion<Scalar>(scale0 * coeffs() + scale1 * other.coeffs()); } diff --git a/Eigen/src/Geometry/Rotation2D.h b/Eigen/src/Geometry/Rotation2D.h index b80fcace0..7f24a3eae 100644 --- a/Eigen/src/Geometry/Rotation2D.h +++ b/Eigen/src/Geometry/Rotation2D.h @@ -121,7 +121,7 @@ public: * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const Rotation2D& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + bool isApprox(const Rotation2D& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const { return ei_isApprox(m_angle,other.m_angle, prec); } }; diff --git a/Eigen/src/Geometry/Scaling.h b/Eigen/src/Geometry/Scaling.h index ce191b5da..4695914fd 100644 --- a/Eigen/src/Geometry/Scaling.h +++ b/Eigen/src/Geometry/Scaling.h @@ -107,7 +107,7 @@ public: * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const UniformScaling& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + bool isApprox(const UniformScaling& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const { return ei_isApprox(m_factor, other.factor(), prec); } }; diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index cf961b98b..89df73505 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -187,7 +187,7 @@ public: /** type of read/write reference to the affine part of the transformation */ typedef typename ei_meta_if<int(Mode)==int(AffineCompact), MatrixType&, - NestByValue<Block<MatrixType,Dim,HDim> > >::ret AffinePartNested; + Block<MatrixType,Dim,HDim> >::ret AffinePartNested; /** type of a vector */ typedef Matrix<Scalar,Dim,1> VectorType; /** type of a read/write reference to the translation part of the rotation */ @@ -424,7 +424,7 @@ public: * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const Transform& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + bool isApprox(const Transform& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const { return m_matrix.isApprox(other.m_matrix, prec); } /** Sets the last row to [0 ... 0 1] diff --git a/Eigen/src/Geometry/Translation.h b/Eigen/src/Geometry/Translation.h index 1fff03810..b7477df9f 100644 --- a/Eigen/src/Geometry/Translation.h +++ b/Eigen/src/Geometry/Translation.h @@ -154,7 +154,7 @@ public: * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const Translation& other, typename NumTraits<Scalar>::Real prec = precision<Scalar>()) const + bool isApprox(const Translation& other, typename NumTraits<Scalar>::Real prec = dummy_precision<Scalar>()) const { return m_coeffs.isApprox(other.m_coeffs, prec); } }; diff --git a/Eigen/src/Geometry/arch/Geometry_SSE.h b/Eigen/src/Geometry/arch/Geometry_SSE.h index a6ed10d82..297932f92 100644 --- a/Eigen/src/Geometry/arch/Geometry_SSE.h +++ b/Eigen/src/Geometry/arch/Geometry_SSE.h @@ -26,7 +26,8 @@ #ifndef EIGEN_GEOMETRY_SSE_H #define EIGEN_GEOMETRY_SSE_H -template<class Derived, class OtherDerived> struct ei_quat_product<EiArch_SSE, Derived, OtherDerived, float, Aligned> +template<class Derived, class OtherDerived> +struct ei_quat_product<Architecture::SSE, Derived, OtherDerived, float, Aligned> { inline static Quaternion<float> run(const QuaternionBase<Derived>& _a, const QuaternionBase<OtherDerived>& _b) { @@ -48,7 +49,8 @@ template<class Derived, class OtherDerived> struct ei_quat_product<EiArch_SSE, D }; template<typename VectorLhs,typename VectorRhs> -struct ei_cross3_impl<EiArch_SSE,VectorLhs,VectorRhs,float,true> { +struct ei_cross3_impl<Architecture::SSE,VectorLhs,VectorRhs,float,true> +{ inline static typename ei_plain_matrix_type<VectorLhs>::type run(const VectorLhs& lhs, const VectorRhs& rhs) { diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h index 85aa90362..25e962001 100644 --- a/Eigen/src/Householder/HouseholderSequence.h +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -73,27 +73,31 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence CoeffsType>::ret> ConjugateReturnType; HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false) - : m_vectors(v), m_coeffs(h), m_trans(trans) + : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()) + {} + + HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors) + : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors) {} int rows() const { return m_vectors.rows(); } int cols() const { return m_vectors.rows(); } HouseholderSequence transpose() const - { return HouseholderSequence(m_vectors, m_coeffs, !m_trans); } + { return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors); } ConjugateReturnType conjugate() const - { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); } + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors); } ConjugateReturnType adjoint() const - { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); } + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors); } ConjugateReturnType inverse() const { return adjoint(); } /** \internal */ template<typename DestType> void evalTo(DestType& dst) const { - int vecs = std::min(m_vectors.cols(),m_vectors.rows()); + int vecs = m_actualVectors; int length = m_vectors.rows(); dst.setIdentity(); Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows()); @@ -111,22 +115,22 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence /** \internal */ template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const { - int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors - int length = m_vectors.rows(); // size of the largest householder vector - Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.rows()); + int vecs = m_actualVectors; // number of householder vectors + int length = m_vectors.rows(); // size of the largest householder vector + Matrix<Scalar,1,Dest::RowsAtCompileTime> temp(dst.rows()); for(int k = 0; k < vecs; ++k) { int actual_k = m_trans ? vecs-k-1 : k; - dst.corner(BottomRight, dst.rows(), length-k) - .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0)); + dst.corner(BottomRight, dst.rows(), length-actual_k) + .applyHouseholderOnTheRight(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); } } /** \internal */ template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const { - int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors - int length = m_vectors.rows(); // size of the largest householder vector + int vecs = m_actualVectors; // number of householder vectors + int length = m_vectors.rows(); // size of the largest householder vector Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols()); for(int k = 0; k < vecs; ++k) { @@ -156,9 +160,7 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence typename VectorsType::Nested m_vectors; typename CoeffsType::Nested m_coeffs; bool m_trans; - -private: - HouseholderSequence& operator=(const HouseholderSequence&); + int m_actualVectors; }; template<typename VectorsType, typename CoeffsType> @@ -167,4 +169,10 @@ HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsTyp return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans); } +template<typename VectorsType, typename CoeffsType> +HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors) +{ + return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors); +} + #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H diff --git a/Eigen/src/LU/CMakeLists.txt b/Eigen/src/LU/CMakeLists.txt index 2788595b0..e0d8d78c1 100644 --- a/Eigen/src/LU/CMakeLists.txt +++ b/Eigen/src/LU/CMakeLists.txt @@ -4,3 +4,5 @@ INSTALL(FILES ${Eigen_LU_SRCS} DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LU COMPONENT Devel ) + +ADD_SUBDIRECTORY(arch) diff --git a/Eigen/src/LU/Determinant.h b/Eigen/src/LU/Determinant.h index 8870d9f20..27ad6abe9 100644 --- a/Eigen/src/LU/Determinant.h +++ b/Eigen/src/LU/Determinant.h @@ -118,7 +118,9 @@ template<typename Derived> inline typename ei_traits<Derived>::Scalar MatrixBase<Derived>::determinant() const { assert(rows() == cols()); - return ei_determinant_impl<Derived>::run(derived()); + typedef typename ei_nested<Derived,RowsAtCompileTime>::type Nested; + Nested nested(derived()); + return ei_determinant_impl<typename ei_cleantype<Nested>::type>::run(nested); } #endif // EIGEN_DETERMINANT_H diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 89eec38c7..149af64bc 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -14,7 +14,7 @@ // the License, or (at your option) any later version. // // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the // GNU General Public License for more details. // @@ -63,8 +63,8 @@ template<typename _MatrixType> class FullPivLU typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; - typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType; + typedef PermutationMatrix<MatrixType::ColsAtCompileTime> PermutationQType; + typedef PermutationMatrix<MatrixType::RowsAtCompileTime> PermutationPType; /** * \brief Default Constructor. @@ -119,26 +119,22 @@ template<typename _MatrixType> class FullPivLU * diagonal coefficient of U. */ RealScalar maxPivot() const { return m_maxpivot; } - - /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, - * representing the P permutation i.e. the permutation of the rows. For its precise meaning, - * see the examples given in the documentation of class FullPivLU. + + /** \returns the permutation matrix P * * \sa permutationQ() */ - inline const IntColVectorType& permutationP() const + inline const PermutationPType& permutationP() const { ei_assert(m_isInitialized && "LU is not initialized."); return m_p; } - /** \returns a vector of integers, whose size is the number of columns of the matrix being - * decomposed, representing the Q permutation i.e. the permutation of the columns. - * For its precise meaning, see the examples given in the documentation of class FullPivLU. + /** \returns the permutation matrix Q * * \sa permutationP() */ - inline const IntRowVectorType& permutationQ() const + inline const PermutationQType& permutationQ() const { ei_assert(m_isInitialized && "LU is not initialized."); return m_q; @@ -238,8 +234,9 @@ template<typename _MatrixType> class FullPivLU * who need to determine when pivots are to be considered nonzero. This is not used for the * LU decomposition itself. * - * When it needs to get the threshold value, Eigen calls threshold(). By default, this calls - * defaultThreshold(). Once you have called the present method setThreshold(const RealScalar&), + * When it needs to get the threshold value, Eigen calls threshold(). By default, this + * uses a formula to automatically determine a reasonable threshold. + * Once you have called the present method setThreshold(const RealScalar&), * your value is used instead. * * \param threshold The new value to use as the threshold. @@ -307,7 +304,7 @@ template<typename _MatrixType> class FullPivLU inline int dimensionOfKernel() const { ei_assert(m_isInitialized && "LU is not initialized."); - return m_lu.cols() - rank(); + return cols() - rank(); } /** \returns true if the matrix of which *this is the LU decomposition represents an injective @@ -320,7 +317,7 @@ template<typename _MatrixType> class FullPivLU inline bool isInjective() const { ei_assert(m_isInitialized && "LU is not initialized."); - return rank() == m_lu.cols(); + return rank() == cols(); } /** \returns true if the matrix of which *this is the LU decomposition represents a surjective @@ -333,7 +330,7 @@ template<typename _MatrixType> class FullPivLU inline bool isSurjective() const { ei_assert(m_isInitialized && "LU is not initialized."); - return rank() == m_lu.rows(); + return rank() == rows(); } /** \returns true if the matrix of which *this is the LU decomposition is invertible. @@ -355,12 +352,12 @@ template<typename _MatrixType> class FullPivLU * * \sa MatrixBase::inverse() */ - inline const ei_solve_retval<FullPivLU,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const + inline const ei_solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const { ei_assert(m_isInitialized && "LU is not initialized."); ei_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!"); - return ei_solve_retval<FullPivLU,NestByValue<typename MatrixType::IdentityReturnType> > - (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); + return ei_solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> + (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols())); } inline int rows() const { return m_lu.rows(); } @@ -368,8 +365,8 @@ template<typename _MatrixType> class FullPivLU protected: MatrixType m_lu; - IntColVectorType m_p; - IntRowVectorType m_q; + PermutationPType m_p; + PermutationQType m_q; int m_det_pq, m_nonzero_pivots; RealScalar m_maxpivot, m_prescribedThreshold; bool m_isInitialized, m_usePrescribedThreshold; @@ -393,8 +390,6 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix) { m_isInitialized = true; m_lu = matrix; - m_p.resize(matrix.rows()); - m_q.resize(matrix.cols()); const int size = matrix.diagonalSize(); const int rows = matrix.rows(); @@ -408,6 +403,7 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix) m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) m_maxpivot = RealScalar(0); + for(int k = 0; k < size; ++k) { // First, we need to find the pivot. @@ -424,10 +420,10 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix) // if the pivot (hence the corner) is exactly zero, terminate to avoid generating nan/inf values if(biggest_in_corner == RealScalar(0)) { - // before exiting, make sure to initialize the still uninitialized row_transpositions + // before exiting, make sure to initialize the still uninitialized transpositions // in a sane state without destroying what we already have. m_nonzero_pivots = k; - for(int i = k; i < size; i++) + for(int i = k; i < size; ++i) { rows_transpositions.coeffRef(i) = i; cols_transpositions.coeffRef(i) = i; @@ -463,13 +459,13 @@ FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix) // the main loop is over, we still have to accumulate the transpositions to find the // permutations P and Q - for(int k = 0; k < matrix.rows(); ++k) m_p.coeffRef(k) = k; + m_p.setIdentity(rows); for(int k = size-1; k >= 0; --k) - std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k))); + m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k)); - for(int k = 0; k < matrix.cols(); ++k) m_q.coeffRef(k) = k; + m_q.setIdentity(cols); for(int k = 0; k < size; ++k) - std::swap(m_q.coeffRef(k), m_q.coeffRef(cols_transpositions.coeff(k))); + m_q.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k)); m_det_pq = (number_of_transpositions%2) ? -1 : 1; return *this; @@ -562,9 +558,9 @@ struct ei_kernel_retval<FullPivLU<_MatrixType> > m.col(i).swap(m.col(pivots.coeff(i))); // see the negative sign in the next line, that's what we were talking about above. - for(int i = 0; i < rank(); ++i) dst.row(dec().permutationQ().coeff(i)) = -m.row(i).end(dimker); - for(int i = rank(); i < cols; ++i) dst.row(dec().permutationQ().coeff(i)).setZero(); - for(int k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().coeff(rank()+k), k) = Scalar(1); + for(int i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).end(dimker); + for(int i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero(); + for(int k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1); } }; @@ -601,7 +597,7 @@ struct ei_image_retval<FullPivLU<_MatrixType> > ei_internal_assert(p == rank()); for(int i = 0; i < rank(); ++i) - dst.col(i) = originalMatrix().col(dec().permutationQ().coeff(pivots.coeff(i))); + dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i))); } }; @@ -623,7 +619,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> * Step 4: result = Q * c; */ - const int rows = dec().matrixLU().rows(), cols = dec().matrixLU().cols(), + const int rows = dec().rows(), cols = dec().cols(), nonzero_pivots = dec().nonzeroPivots(); ei_assert(rhs().rows() == rows); const int smalldim = std::min(rows, cols); @@ -637,8 +633,7 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> typename Rhs::PlainMatrixType c(rhs().rows(), rhs().cols()); // Step 1 - for(int i = 0; i < rows; ++i) - c.row(dec().permutationP().coeff(i)) = rhs().row(i); + c = dec().permutationP() * rhs(); // Step 2 dec().matrixLU() @@ -660,9 +655,9 @@ struct ei_solve_retval<FullPivLU<_MatrixType>, Rhs> // Step 4 for(int i = 0; i < nonzero_pivots; ++i) - dst.row(dec().permutationQ().coeff(i)) = c.row(i); + dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i); for(int i = nonzero_pivots; i < dec().matrixLU().cols(); ++i) - dst.row(dec().permutationQ().coeff(i)).setZero(); + dst.row(dec().permutationQ().indices().coeff(i)).setZero(); } }; diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index 13de19ace..ea2330da3 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -182,91 +182,37 @@ struct ei_compute_inverse_and_det_with_check<MatrixType, ResultType, 3> *** Size 4 implementation *** ****************************/ -template<typename MatrixType, typename ResultType> -void ei_compute_inverse_size4_helper(const MatrixType& matrix, ResultType& result) +template<int Arch, typename Scalar, typename MatrixType, typename ResultType> +struct ei_compute_inverse_size4 { - /* Let's split M into four 2x2 blocks: - * (P Q) - * (R S) - * If P is invertible, with inverse denoted by P_inverse, and if - * (S - R*P_inverse*Q) is also invertible, then the inverse of M is - * (P' Q') - * (R' S') - * where - * S' = (S - R*P_inverse*Q)^(-1) - * P' = P1 + (P1*Q) * S' *(R*P_inverse) - * Q' = -(P_inverse*Q) * S' - * R' = -S' * (R*P_inverse) - */ - typedef Block<ResultType,2,2> XprBlock22; - typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22; - Block22 P_inverse; - ei_compute_inverse<XprBlock22, Block22>::run(matrix.template block<2,2>(0,0), P_inverse); - const Block22 Q = matrix.template block<2,2>(0,2); - const Block22 P_inverse_times_Q = P_inverse * Q; - const XprBlock22 R = matrix.template block<2,2>(2,0); - const Block22 R_times_P_inverse = R * P_inverse; - const Block22 R_times_P_inverse_times_Q = R_times_P_inverse * Q; - const XprBlock22 S = matrix.template block<2,2>(2,2); - const Block22 X = S - R_times_P_inverse_times_Q; - Block22 Y; - ei_compute_inverse<Block22, Block22>::run(X, Y); - result.template block<2,2>(2,2) = Y; - result.template block<2,2>(2,0) = - Y * R_times_P_inverse; - const Block22 Z = P_inverse_times_Q * Y; - result.template block<2,2>(0,2) = - Z; - result.template block<2,2>(0,0) = P_inverse + Z * R_times_P_inverse; -} + static void run(const MatrixType& matrix, ResultType& result) + { + result.coeffRef(0,0) = matrix.minor(0,0).determinant(); + result.coeffRef(1,0) = -matrix.minor(0,1).determinant(); + result.coeffRef(2,0) = matrix.minor(0,2).determinant(); + result.coeffRef(3,0) = -matrix.minor(0,3).determinant(); + result.coeffRef(0,2) = matrix.minor(2,0).determinant(); + result.coeffRef(1,2) = -matrix.minor(2,1).determinant(); + result.coeffRef(2,2) = matrix.minor(2,2).determinant(); + result.coeffRef(3,2) = -matrix.minor(2,3).determinant(); + result.coeffRef(0,1) = -matrix.minor(1,0).determinant(); + result.coeffRef(1,1) = matrix.minor(1,1).determinant(); + result.coeffRef(2,1) = -matrix.minor(1,2).determinant(); + result.coeffRef(3,1) = matrix.minor(1,3).determinant(); + result.coeffRef(0,3) = -matrix.minor(3,0).determinant(); + result.coeffRef(1,3) = matrix.minor(3,1).determinant(); + result.coeffRef(2,3) = -matrix.minor(3,2).determinant(); + result.coeffRef(3,3) = matrix.minor(3,3).determinant(); + result /= (matrix.col(0).cwise()*result.row(0).transpose()).sum(); + } +}; template<typename MatrixType, typename ResultType> struct ei_compute_inverse<MatrixType, ResultType, 4> + : ei_compute_inverse_size4<Architecture::Target, typename MatrixType::Scalar, + MatrixType, ResultType> { - static inline void run(const MatrixType& _matrix, ResultType& result) - { - typedef typename ResultType::Scalar Scalar; - typedef typename MatrixType::RealScalar RealScalar; - - // we will do row permutations on the matrix. This copy should have negligible cost. - // if not, consider working in-place on the matrix (const-cast it, but then undo the permutations - // to nevertheless honor constness) - typename MatrixType::PlainMatrixType matrix(_matrix); - - // let's extract from the 2 first colums a 2x2 block whose determinant is as big as possible. - int good_row0, good_row1, good_i; - Matrix<RealScalar,6,1> absdet; - - // any 2x2 block with determinant above this threshold will be considered good enough - RealScalar d = (matrix.col(0).squaredNorm()+matrix.col(1).squaredNorm()) * RealScalar(1e-2); - #define ei_inv_size4_helper_macro(i,row0,row1) \ - absdet[i] = ei_abs(matrix.coeff(row0,0)*matrix.coeff(row1,1) \ - - matrix.coeff(row0,1)*matrix.coeff(row1,0)); \ - if(absdet[i] > d) { good_row0=row0; good_row1=row1; goto good; } - ei_inv_size4_helper_macro(0,0,1) - ei_inv_size4_helper_macro(1,0,2) - ei_inv_size4_helper_macro(2,0,3) - ei_inv_size4_helper_macro(3,1,2) - ei_inv_size4_helper_macro(4,1,3) - ei_inv_size4_helper_macro(5,2,3) - - // no 2x2 block has determinant bigger than the threshold. So just take the one that - // has the biggest determinant - absdet.maxCoeff(&good_i); - good_row0 = good_i <= 2 ? 0 : good_i <= 4 ? 1 : 2; - good_row1 = good_i <= 2 ? good_i+1 : good_i <= 4 ? good_i-1 : 3; - - // now good_row0 and good_row1 are correctly set - good: - - // do row permutations to move this 2x2 block to the top - matrix.row(0).swap(matrix.row(good_row0)); - matrix.row(1).swap(matrix.row(good_row1)); - // now applying our helper function is numerically stable - ei_compute_inverse_size4_helper(matrix, result); - // Since we did row permutations on the original matrix, we need to do column permutations - // in the reverse order on the inverse - result.col(1).swap(result.col(good_row1)); - result.col(0).swap(result.col(good_row0)); - } + // FIXME empty? }; template<typename MatrixType, typename ResultType> @@ -300,7 +246,8 @@ template<typename MatrixType> struct ei_inverse_impl : public ReturnByValue<ei_inverse_impl<MatrixType> > { // for 2x2, it's worth giving a chance to avoid evaluating. - // for larger sizes, evaluating has negligible cost and limits code size. + // for larger sizes, evaluating has negligible cost, limits code size, + // and allows for vectorized paths. typedef typename ei_meta_if< MatrixType::RowsAtCompileTime == 2, typename ei_nested<MatrixType,2>::type, @@ -326,7 +273,7 @@ struct ei_inverse_impl : public ReturnByValue<ei_inverse_impl<MatrixType> > * * \returns the matrix inverse of this matrix. * - * For small fixed sizes up to 4x4, this method uses ad-hoc methods (cofactors up to 3x3, Euler's trick for 4x4). + * For small fixed sizes up to 4x4, this method uses cofactors. * In the general case, this method uses class PartialPivLU. * * \note This matrix must be invertible, otherwise the result is undefined. If you need an diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h index 37cffd7a1..deb29b5c6 100644 --- a/Eigen/src/LU/PartialPivLU.h +++ b/Eigen/src/LU/PartialPivLU.h @@ -64,10 +64,8 @@ template<typename _MatrixType> class PartialPivLU typedef _MatrixType MatrixType; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType; - typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType; - typedef Matrix<Scalar, 1, MatrixType::ColsAtCompileTime> RowVectorType; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVectorType; + typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> PermutationVectorType; + typedef PermutationMatrix<MatrixType::RowsAtCompileTime> PermutationType; enum { MaxSmallDimAtCompileTime = EIGEN_ENUM_MIN( MatrixType::MaxColsAtCompileTime, @@ -105,11 +103,9 @@ template<typename _MatrixType> class PartialPivLU return m_lu; } - /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, - * representing the P permutation i.e. the permutation of the rows. For its precise meaning, - * see the examples given in the documentation of class FullPivLU. + /** \returns the permutation matrix P. */ - inline const IntColVectorType& permutationP() const + inline const PermutationType& permutationP() const { ei_assert(m_isInitialized && "PartialPivLU is not initialized."); return m_p; @@ -147,11 +143,11 @@ template<typename _MatrixType> class PartialPivLU * * \sa MatrixBase::inverse(), LU::inverse() */ - inline const ei_solve_retval<PartialPivLU,NestByValue<typename MatrixType::IdentityReturnType> > inverse() const + inline const ei_solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const { ei_assert(m_isInitialized && "PartialPivLU is not initialized."); - return ei_solve_retval<PartialPivLU,NestByValue<typename MatrixType::IdentityReturnType> > - (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()).nestByValue()); + return ei_solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> + (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols())); } /** \returns the determinant of the matrix of which @@ -174,7 +170,7 @@ template<typename _MatrixType> class PartialPivLU protected: MatrixType m_lu; - IntColVectorType m_p; + PermutationType m_p; int m_det_p; bool m_isInitialized; }; @@ -379,20 +375,19 @@ template<typename MatrixType> PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix) { m_lu = matrix; - m_p.resize(matrix.rows()); ei_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); const int size = matrix.rows(); - IntColVectorType rows_transpositions(size); + PermutationVectorType rows_transpositions(size); int nb_transpositions; ei_partial_lu_inplace(m_lu, rows_transpositions, nb_transpositions); m_det_p = (nb_transpositions%2) ? -1 : 1; - for(int k = 0; k < size; ++k) m_p.coeffRef(k) = k; + m_p.setIdentity(size); for(int k = size-1; k >= 0; --k) - std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k))); + m_p.applyTranspositionOnTheRight(k, rows_transpositions.coeff(k)); m_isInitialized = true; return *this; @@ -428,7 +423,7 @@ struct ei_solve_retval<PartialPivLU<_MatrixType>, Rhs> dst.resize(size, rhs().cols()); // Step 1 - for(int i = 0; i < size; ++i) dst.row(dec().permutationP().coeff(i)) = rhs().row(i); + dst = dec().permutationP() * rhs(); // Step 2 dec().matrixLU().template triangularView<UnitLowerTriangular>().solveInPlace(dst); diff --git a/Eigen/src/LU/arch/CMakeLists.txt b/Eigen/src/LU/arch/CMakeLists.txt new file mode 100644 index 000000000..f6b7ed9ec --- /dev/null +++ b/Eigen/src/LU/arch/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_LU_arch_SRCS "*.h") + +INSTALL(FILES + ${Eigen_LU_arch_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LU/arch COMPONENT Devel + ) diff --git a/Eigen/src/LU/arch/Inverse_SSE.h b/Eigen/src/LU/arch/Inverse_SSE.h new file mode 100644 index 000000000..cded9195c --- /dev/null +++ b/Eigen/src/LU/arch/Inverse_SSE.h @@ -0,0 +1,151 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 1999 Intel Corporation +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +// The SSE code for the 4x4 float matrix inverse in this file comes from the file +// ftp://download.intel.com/design/PentiumIII/sml/24504301.pdf +// See page ii of that document for legal stuff. Not being lawyers, we just assume +// here that if Intel makes this document publically available, with source code +// and detailed explanations, it's because they want their CPUs to be fed with +// good code, and therefore they presumably don't mind us using it in Eigen. + +#ifndef EIGEN_INVERSE_SSE_H +#define EIGEN_INVERSE_SSE_H + +template<typename MatrixType, typename ResultType> +struct ei_compute_inverse_size4<Architecture::SSE, float, MatrixType, ResultType> +{ + static void run(const MatrixType& matrix, ResultType& result) + { + // Variables (Streaming SIMD Extensions registers) which will contain cofactors and, later, the + // lines of the inverted matrix. + __m128 minor0, minor1, minor2, minor3; + + // Variables which will contain the lines of the reference matrix and, later (after the transposition), + // the columns of the original matrix. + __m128 row0, row1, row2, row3; + + // Temporary variables and the variable that will contain the matrix determinant. + __m128 det, tmp1; + + // Matrix transposition + const float *src = matrix.data(); + tmp1 = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)src)), (__m64*)(src+ 4)); + row1 = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+8))), (__m64*)(src+12)); + row0 = _mm_shuffle_ps(tmp1, row1, 0x88); + row1 = _mm_shuffle_ps(row1, tmp1, 0xDD); + tmp1 = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+ 2))), (__m64*)(src+ 6)); + row3 = _mm_loadh_pi(_mm_castpd_ps(_mm_load_sd((double*)(src+10))), (__m64*)(src+14)); + row2 = _mm_shuffle_ps(tmp1, row3, 0x88); + row3 = _mm_shuffle_ps(row3, tmp1, 0xDD); + + + // Cofactors calculation. Because in the process of cofactor computation some pairs in three- + // element products are repeated, it is not reasonable to load these pairs anew every time. The + // values in the registers with these pairs are formed using shuffle instruction. Cofactors are + // calculated row by row (4 elements are placed in 1 SP FP SIMD floating point register). + + tmp1 = _mm_mul_ps(row2, row3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + minor0 = _mm_mul_ps(row1, tmp1); + minor1 = _mm_mul_ps(row0, tmp1); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + minor0 = _mm_sub_ps(_mm_mul_ps(row1, tmp1), minor0); + minor1 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor1); + minor1 = _mm_shuffle_ps(minor1, minor1, 0x4E); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row1, row2); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + minor0 = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor0); + minor3 = _mm_mul_ps(row0, tmp1); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + minor0 = _mm_sub_ps(minor0, _mm_mul_ps(row3, tmp1)); + minor3 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor3); + minor3 = _mm_shuffle_ps(minor3, minor3, 0x4E); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(_mm_shuffle_ps(row1, row1, 0x4E), row3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + row2 = _mm_shuffle_ps(row2, row2, 0x4E); + minor0 = _mm_add_ps(_mm_mul_ps(row2, tmp1), minor0); + minor2 = _mm_mul_ps(row0, tmp1); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + minor0 = _mm_sub_ps(minor0, _mm_mul_ps(row2, tmp1)); + minor2 = _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor2); + minor2 = _mm_shuffle_ps(minor2, minor2, 0x4E); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row1); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + minor2 = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor2); + minor3 = _mm_sub_ps(_mm_mul_ps(row2, tmp1), minor3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + minor2 = _mm_sub_ps(_mm_mul_ps(row3, tmp1), minor2); + minor3 = _mm_sub_ps(minor3, _mm_mul_ps(row2, tmp1)); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row3); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + minor1 = _mm_sub_ps(minor1, _mm_mul_ps(row2, tmp1)); + minor2 = _mm_add_ps(_mm_mul_ps(row1, tmp1), minor2); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + minor1 = _mm_add_ps(_mm_mul_ps(row2, tmp1), minor1); + minor2 = _mm_sub_ps(minor2, _mm_mul_ps(row1, tmp1)); + // ----------------------------------------------- + tmp1 = _mm_mul_ps(row0, row2); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0xB1); + minor1 = _mm_add_ps(_mm_mul_ps(row3, tmp1), minor1); + minor3 = _mm_sub_ps(minor3, _mm_mul_ps(row1, tmp1)); + tmp1 = _mm_shuffle_ps(tmp1, tmp1, 0x4E); + minor1 = _mm_sub_ps(minor1, _mm_mul_ps(row3, tmp1)); + minor3 = _mm_add_ps(_mm_mul_ps(row1, tmp1), minor3); + + // Evaluation of determinant and its reciprocal value. In the original Intel document, + // 1/det was evaluated using a fast rcpps command with subsequent approximation using + // the Newton-Raphson algorithm. Here, we go for a IEEE-compliant division instead, + // so as to not compromise precision at all. + det = _mm_mul_ps(row0, minor0); + det = _mm_add_ps(_mm_shuffle_ps(det, det, 0x4E), det); + det = _mm_add_ss(_mm_shuffle_ps(det, det, 0xB1), det); +// tmp1= _mm_rcp_ss(det); +// det= _mm_sub_ss(_mm_add_ss(tmp1, tmp1), _mm_mul_ss(det, _mm_mul_ss(tmp1, tmp1))); + det = _mm_div_ss(_mm_set_ss(1.0f), det); // <--- yay, one original line not copied from Intel + det = _mm_shuffle_ps(det, det, 0x00); + // warning, Intel's variable naming is very confusing: now 'det' is 1/det ! + + // Multiplication of cofactors by 1/det. Storing the inverse matrix to the address in pointer src. + minor0 = _mm_mul_ps(det, minor0); + float *dst = result.data(); + _mm_storel_pi((__m64*)(dst), minor0); + _mm_storeh_pi((__m64*)(dst+2), minor0); + minor1 = _mm_mul_ps(det, minor1); + _mm_storel_pi((__m64*)(dst+4), minor1); + _mm_storeh_pi((__m64*)(dst+6), minor1); + minor2 = _mm_mul_ps(det, minor2); + _mm_storel_pi((__m64*)(dst+ 8), minor2); + _mm_storeh_pi((__m64*)(dst+10), minor2); + minor3 = _mm_mul_ps(det, minor3); + _mm_storel_pi((__m64*)(dst+12), minor3); + _mm_storeh_pi((__m64*)(dst+14), minor3); + } +}; + +#endif // EIGEN_INVERSE_SSE_H diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 02864caa5..b4c1a5fcc 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -57,10 +57,9 @@ template<typename _MatrixType> class ColPivHouseholderQR typedef typename MatrixType::RealScalar RealScalar; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType; + typedef PermutationMatrix<ColsAtCompileTime> PermutationType; typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType; - typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType; typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; - typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType; typedef typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; @@ -75,7 +74,8 @@ template<typename _MatrixType> class ColPivHouseholderQR ColPivHouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs(std::min(matrix.rows(),matrix.cols())), - m_isInitialized(false) + m_isInitialized(false), + m_usePrescribedThreshold(false) { compute(matrix); } @@ -105,7 +105,7 @@ template<typename _MatrixType> class ColPivHouseholderQR return ei_solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived()); } - HouseholderSequenceType matrixQ(void) const; + HouseholderSequenceType householderQ(void) const; /** \returns a reference to the matrix where the Householder QR decomposition is stored */ @@ -117,7 +117,7 @@ template<typename _MatrixType> class ColPivHouseholderQR ColPivHouseholderQR& compute(const MatrixType& matrix); - const IntRowVectorType& colsPermutation() const + const PermutationType& colsPermutation() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); return m_cols_permutation; @@ -154,54 +154,63 @@ template<typename _MatrixType> class ColPivHouseholderQR /** \returns the rank of the matrix of which *this is the QR decomposition. * - * \note This is computed at the time of the construction of the QR decomposition. This - * method does not perform any further computation. + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). */ inline int rank() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return m_rank; + RealScalar premultiplied_threshold = ei_abs(m_maxpivot) * threshold(); + int result = 0; + for(int i = 0; i < m_nonzero_pivots; ++i) + result += (ei_abs(m_qr.coeff(i,i)) > premultiplied_threshold); + return result; } /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. * - * \note Since the rank is computed at the time of the construction of the QR decomposition, this - * method almost does not perform any further computation. + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). */ inline int dimensionOfKernel() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return m_qr.cols() - m_rank; + return cols() - rank(); } /** \returns true if the matrix of which *this is the QR decomposition represents an injective * linear map, i.e. has trivial kernel; false otherwise. * - * \note Since the rank is computed at the time of the construction of the QR decomposition, this - * method almost does not perform any further computation. + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). */ inline bool isInjective() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return m_rank == m_qr.cols(); + return rank() == cols(); } /** \returns true if the matrix of which *this is the QR decomposition represents a surjective * linear map; false otherwise. * - * \note Since the rank is computed at the time of the construction of the QR decomposition, this - * method almost does not perform any further computation. + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). */ inline bool isSurjective() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return m_rank == m_qr.rows(); + return rank() == rows(); } /** \returns true if the matrix of which *this is the QR decomposition is invertible. * - * \note Since the rank is computed at the time of the construction of the QR decomposition, this - * method almost does not perform any further computation. + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). */ inline bool isInvertible() const { @@ -215,25 +224,92 @@ template<typename _MatrixType> class ColPivHouseholderQR * Use isInvertible() to first determine whether this matrix is invertible. */ inline const - ei_solve_retval<ColPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> > + ei_solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType> inverse() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return ei_solve_retval<ColPivHouseholderQR,NestByValue<typename MatrixType::IdentityReturnType> > - (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()).nestByValue()); + return ei_solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType> + (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols())); } inline int rows() const { return m_qr.rows(); } inline int cols() const { return m_qr.cols(); } const HCoeffsType& hCoeffs() const { return m_hCoeffs; } + /** Allows to prescribe a threshold to be used by certain methods, such as rank(), + * who need to determine when pivots are to be considered nonzero. This is not used for the + * QR decomposition itself. + * + * When it needs to get the threshold value, Eigen calls threshold(). By default, this + * uses a formula to automatically determine a reasonable threshold. + * Once you have called the present method setThreshold(const RealScalar&), + * your value is used instead. + * + * \param threshold The new value to use as the threshold. + * + * A pivot will be considered nonzero if its absolute value is strictly greater than + * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ + * where maxpivot is the biggest pivot. + * + * If you want to come back to the default behavior, call setThreshold(Default_t) + */ + ColPivHouseholderQR& setThreshold(const RealScalar& threshold) + { + m_usePrescribedThreshold = true; + m_prescribedThreshold = threshold; + } + + /** Allows to come back to the default behavior, letting Eigen use its default formula for + * determining the threshold. + * + * You should pass the special object Eigen::Default as parameter here. + * \code qr.setThreshold(Eigen::Default); \endcode + * + * See the documentation of setThreshold(const RealScalar&). + */ + ColPivHouseholderQR& setThreshold(Default_t) + { + m_usePrescribedThreshold = false; + } + + /** Returns the threshold that will be used by certain methods such as rank(). + * + * See the documentation of setThreshold(const RealScalar&). + */ + RealScalar threshold() const + { + ei_assert(m_isInitialized || m_usePrescribedThreshold); + return m_usePrescribedThreshold ? m_prescribedThreshold + // this formula comes from experimenting (see "LU precision tuning" thread on the list) + // and turns out to be identical to Higham's formula used already in LDLt. + : epsilon<Scalar>() * m_qr.diagonalSize(); + } + + /** \returns the number of nonzero pivots in the QR decomposition. + * Here nonzero is meant in the exact sense, not in a fuzzy sense. + * So that notion isn't really intrinsically interesting, but it is + * still useful when implementing algorithms. + * + * \sa rank() + */ + inline int nonzeroPivots() const + { + ei_assert(m_isInitialized && "LU is not initialized."); + return m_nonzero_pivots; + } + + /** \returns the absolute value of the biggest pivot, i.e. the biggest + * diagonal coefficient of U. + */ + RealScalar maxPivot() const { return m_maxpivot; } + protected: MatrixType m_qr; HCoeffsType m_hCoeffs; - IntRowVectorType m_cols_permutation; - bool m_isInitialized; - RealScalar m_precision; - int m_rank; + PermutationType m_cols_permutation; + bool m_isInitialized, m_usePrescribedThreshold; + RealScalar m_prescribedThreshold, m_maxpivot; + int m_nonzero_pivots; int m_det_pq; }; @@ -260,63 +336,85 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const { int rows = matrix.rows(); int cols = matrix.cols(); - int size = std::min(rows,cols); - m_rank = size; + int size = matrix.diagonalSize(); m_qr = matrix; m_hCoeffs.resize(size); RowVectorType temp(cols); - m_precision = epsilon<Scalar>() * size; - IntRowVectorType cols_transpositions(matrix.cols()); - m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; RealRowVectorType colSqNorms(cols); for(int k = 0; k < cols; ++k) colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); - RealScalar biggestColSqNorm = colSqNorms.maxCoeff(); - for (int k = 0; k < size; ++k) - { - int biggest_col_in_corner; - RealScalar biggestColSqNormInCorner = colSqNorms.end(cols-k).maxCoeff(&biggest_col_in_corner); - biggest_col_in_corner += k; + RealScalar threshold_helper = colSqNorms.maxCoeff() * ei_abs2(epsilon<Scalar>()) / rows; - // if the corner is negligible, then we have less than full rank, and we can finish early - if(ei_isMuchSmallerThan(biggestColSqNormInCorner, biggestColSqNorm, m_precision)) + m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) + m_maxpivot = RealScalar(0); + + for(int k = 0; k < size; ++k) + { + // first, we look up in our table colSqNorms which column has the biggest squared norm + int biggest_col_index; + RealScalar biggest_col_sq_norm = colSqNorms.end(cols-k).maxCoeff(&biggest_col_index); + biggest_col_index += k; + + // since our table colSqNorms accumulates imprecision at every step, we must now recompute + // the actual squared norm of the selected column. + // Note that not doing so does result in solve() sometimes returning inf/nan values + // when running the unit test with 1000 repetitions. + biggest_col_sq_norm = m_qr.col(biggest_col_index).end(rows-k).squaredNorm(); + + // we store that back into our table: it can't hurt to correct our table. + colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm; + + // if the current biggest column is smaller than epsilon times the initial biggest column, + // terminate to avoid generating nan/inf values. + // Note that here, if we test instead for "biggest == 0", we get a failure every 1000 (or so) + // repetitions of the unit test, with the result of solve() filled with large values of the order + // of 1/(size*epsilon). + if(biggest_col_sq_norm < threshold_helper * (rows-k)) { - m_rank = k; - for(int i = k; i < size; i++) - { - cols_transpositions.coeffRef(i) = i; - m_hCoeffs.coeffRef(i) = Scalar(0); - } + m_nonzero_pivots = k; + m_hCoeffs.end(size-k).setZero(); + m_qr.corner(BottomRight,rows-k,cols-k) + .template triangularView<StrictlyLowerTriangular>() + .setZero(); break; } - cols_transpositions.coeffRef(k) = biggest_col_in_corner; - if(k != biggest_col_in_corner) { - m_qr.col(k).swap(m_qr.col(biggest_col_in_corner)); - std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_in_corner)); + // apply the transposition to the columns + cols_transpositions.coeffRef(k) = biggest_col_index; + if(k != biggest_col_index) { + m_qr.col(k).swap(m_qr.col(biggest_col_index)); + std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_index)); ++number_of_transpositions; } + // generate the householder vector, store it below the diagonal RealScalar beta; m_qr.col(k).end(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); + + // apply the householder transformation to the diagonal coefficient m_qr.coeffRef(k,k) = beta; + // remember the maximum absolute value of diagonal coefficients + if(ei_abs(beta) > m_maxpivot) m_maxpivot = ei_abs(beta); + + // apply the householder transformation m_qr.corner(BottomRight, rows-k, cols-k-1) .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); + // update our table of squared norms of the columns colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwiseAbs2(); } - for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k; - for(int k = 0; k < size; ++k) - std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k))); + m_cols_permutation.setIdentity(cols); + for(int k = 0; k < m_nonzero_pivots; ++k) + m_cols_permutation.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k)); m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; @@ -332,13 +430,12 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - const int rows = dec().rows(), cols = dec().cols(); + const int rows = dec().rows(), cols = dec().cols(), + nonzero_pivots = dec().nonzeroPivots(); dst.resize(cols, rhs().cols()); ei_assert(rhs().rows() == rows); - // FIXME introduce nonzeroPivots() and use it here. and more generally, - // make the same improvements in this dec as in FullPivLU. - if(dec().rank()==0) + if(nonzero_pivots == 0) { dst.setZero(); return; @@ -348,37 +445,37 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs> // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T c.applyOnTheLeft(householderSequence( - dec().matrixQR().corner(TopLeft,rows,dec().rank()), - dec().hCoeffs().start(dec().rank())).transpose() - ); - - if(!dec().isSurjective()) - { - // is c is in the image of R ? - RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec().rank(), c.cols()).cwiseAbs().maxCoeff(); - RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec().rank(), c.cols()).cwiseAbs().maxCoeff(); - // FIXME brain dead - const RealScalar m_precision = epsilon<Scalar>() * std::min(rows,cols); - if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4)) - return; - } + dec().matrixQR(), + dec().hCoeffs(), + true, + dec().nonzeroPivots() + )); dec().matrixQR() - .corner(TopLeft, dec().rank(), dec().rank()) + .corner(TopLeft, nonzero_pivots, nonzero_pivots) + .template triangularView<UpperTriangular>() + .solveInPlace(c.corner(TopLeft, nonzero_pivots, c.cols())); + + + typename Rhs::PlainMatrixType d(c); + d.corner(TopLeft, nonzero_pivots, c.cols()) + = dec().matrixQR() + .corner(TopLeft, nonzero_pivots, nonzero_pivots) .template triangularView<UpperTriangular>() - .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols())); + * c.corner(TopLeft, nonzero_pivots, c.cols()); - for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().coeff(i)) = c.row(i); - for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero(); + for(int i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); + for(int i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); } }; /** \returns the matrix Q as a sequence of householder transformations */ template<typename MatrixType> -typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>::matrixQ() const +typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType> + ::householderQ() const { ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots); } #endif // EIGEN_HIDE_HEAVY_CODE diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index a1eaabf23..51609ca1a 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -58,6 +58,7 @@ template<typename _MatrixType> class FullPivHouseholderQR typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType; typedef Matrix<int, 1, ColsAtCompileTime> IntRowVectorType; + typedef PermutationMatrix<ColsAtCompileTime> PermutationType; typedef Matrix<int, RowsAtCompileTime, 1> IntColVectorType; typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; @@ -112,7 +113,7 @@ template<typename _MatrixType> class FullPivHouseholderQR FullPivHouseholderQR& compute(const MatrixType& matrix); - const IntRowVectorType& colsPermutation() const + const PermutationType& colsPermutation() const { ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); return m_cols_permutation; @@ -215,12 +216,12 @@ template<typename _MatrixType> class FullPivHouseholderQR * \note If this matrix is not invertible, the returned matrix has undefined coefficients. * Use isInvertible() to first determine whether this matrix is invertible. */ inline const - ei_solve_retval<FullPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> > + ei_solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType> inverse() const { ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); - return ei_solve_retval<FullPivHouseholderQR,NestByValue<typename MatrixType::IdentityReturnType> > - (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()).nestByValue()); + return ei_solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType> + (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols())); } inline int rows() const { return m_qr.rows(); } @@ -231,7 +232,7 @@ template<typename _MatrixType> class FullPivHouseholderQR MatrixType m_qr; HCoeffsType m_hCoeffs; IntColVectorType m_rows_transpositions; - IntRowVectorType m_cols_permutation; + PermutationType m_cols_permutation; bool m_isInitialized; RealScalar m_precision; int m_rank; @@ -273,7 +274,6 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons m_rows_transpositions.resize(matrix.rows()); IntRowVectorType cols_transpositions(matrix.cols()); - m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; RealScalar biggest(0); @@ -322,9 +322,9 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); } - for(int k = 0; k < matrix.cols(); ++k) m_cols_permutation.coeffRef(k) = k; + m_cols_permutation.setIdentity(cols); for(int k = 0; k < size; ++k) - std::swap(m_cols_permutation.coeffRef(k), m_cols_permutation.coeffRef(cols_transpositions.coeff(k))); + m_cols_permutation.applyTranspositionOnTheRight(k, cols_transpositions.coeff(k)); m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; @@ -379,8 +379,8 @@ struct ei_solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs> .template triangularView<UpperTriangular>() .solveInPlace(c.corner(TopLeft, dec().rank(), c.cols())); - for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().coeff(i)) = c.row(i); - for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().coeff(i)).setZero(); + for(int i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i); + for(int i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero(); } }; diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index 9d8a33610..895ae046a 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -104,11 +104,10 @@ template<typename _MatrixType> class HouseholderQR ei_assert(m_isInitialized && "HouseholderQR is not initialized."); return ei_solve_retval<HouseholderQR, Rhs>(*this, b.derived()); } - - MatrixQType matrixQ() const; - HouseholderSequenceType matrixQAsHouseholderSequence() const + HouseholderSequenceType householderQ() const { + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); } @@ -240,14 +239,6 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs> } }; -/** \returns the matrix Q */ -template<typename MatrixType> -typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const -{ - ei_assert(m_isInitialized && "HouseholderQR is not initialized."); - return matrixQAsHouseholderSequence(); -} - #endif // EIGEN_HIDE_HEAVY_CODE /** \return the Householder QR decomposition of \c *this. diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 927ef6591..5792c5767 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -236,9 +236,8 @@ struct ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options, true> FullPivHouseholderQR<MatrixType> qr(matrix); work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>(); if(ComputeU) svd.m_matrixU = qr.matrixQ(); - if(ComputeV) - for(int i = 0; i < cols; i++) - svd.m_matrixV.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); + if(ComputeV) svd.m_matrixV = qr.colsPermutation(); + return true; } else return false; @@ -281,9 +280,7 @@ struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true> FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint()); work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>().adjoint(); if(ComputeV) svd.m_matrixV = qr.matrixQ(); - if(ComputeU) - for(int i = 0; i < rows; i++) - svd.m_matrixU.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); + if(ComputeU) svd.m_matrixU = qr.colsPermutation(); return true; } else return false; @@ -297,8 +294,6 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma int rows = matrix.rows(); int cols = matrix.cols(); int diagSize = std::min(rows, cols); - if(ComputeU) m_matrixU = MatrixUType::Zero(rows,rows); - if(ComputeV) m_matrixV = MatrixVType::Zero(cols,cols); m_singularValues.resize(diagSize); const RealScalar precision = 2 * epsilon<Scalar>(); @@ -306,8 +301,8 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const Ma && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, work_matrix, *this)) { work_matrix = matrix.block(0,0,diagSize,diagSize); - if(ComputeU) m_matrixU.diagonal().setOnes(); - if(ComputeV) m_matrixV.diagonal().setOnes(); + if(ComputeU) m_matrixU.setIdentity(rows,rows); + if(ComputeV) m_matrixV.setIdentity(cols,cols); } bool finished = false; diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 6ed9db20a..7a8e4f312 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -190,11 +190,10 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) SingularValuesType& W = m_sigma; bool flag; - int i,its,j,k,nm; - int l=0; + int i=0,its=0,j=0,k=0,l=0,nm=0; Scalar anorm, c, f, g, h, s, scale, x, y, z; bool convergence = true; - Scalar eps = precision<Scalar>(); + Scalar eps = dummy_precision<Scalar>(); Matrix<Scalar,Dynamic,1> rv1(n); g = scale = anorm = 0; diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h index 474626848..2988999d6 100644 --- a/Eigen/src/Sparse/AmbiVector.h +++ b/Eigen/src/Sparse/AmbiVector.h @@ -126,10 +126,6 @@ template<typename _Scalar> class AmbiVector int m_llStart; int m_llCurrent; int m_llSize; - - private: - AmbiVector(const AmbiVector&); - }; /** \returns the number of non zeros in the current sub vector */ @@ -300,7 +296,7 @@ class AmbiVector<_Scalar>::Iterator * In practice, all coefficients having a magnitude smaller than \a epsilon * are skipped. */ - Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*precision<RealScalar>()) + Iterator(const AmbiVector& vec, RealScalar epsilon = RealScalar(0.1)*dummy_precision<RealScalar>()) : m_vector(vec) { m_epsilon = epsilon; diff --git a/Eigen/src/Sparse/CompressedStorage.h b/Eigen/src/Sparse/CompressedStorage.h index abf42914f..b25b05e91 100644 --- a/Eigen/src/Sparse/CompressedStorage.h +++ b/Eigen/src/Sparse/CompressedStorage.h @@ -93,7 +93,7 @@ class CompressedStorage void append(const Scalar& v, int i) { - int id = m_size; + int id = static_cast<int>(m_size); resize(m_size+1, 1); m_values[id] = v; m_indices[id] = i; @@ -135,7 +135,7 @@ class CompressedStorage else end = mid; } - return start; + return static_cast<int>(start); } /** \returns the stored value at index \a key @@ -185,7 +185,7 @@ class CompressedStorage return m_values[id]; } - void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) + void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar>()) { size_t k = 0; size_t n = size(); diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h index 9644dbcd5..4373e1cda 100644 --- a/Eigen/src/Sparse/DynamicSparseMatrix.h +++ b/Eigen/src/Sparse/DynamicSparseMatrix.h @@ -58,6 +58,13 @@ struct ei_traits<DynamicSparseMatrix<_Scalar, _Flags> > }; }; +template<typename _Scalar, int _Options> +struct ei_ref_selector< DynamicSparseMatrix<_Scalar, _Options> > +{ + typedef DynamicSparseMatrix<_Scalar, _Options> MatrixType; + typedef MatrixType const& type; +}; + template<typename _Scalar, int _Flags> class DynamicSparseMatrix : public SparseMatrixBase<DynamicSparseMatrix<_Scalar, _Flags> > @@ -82,7 +89,7 @@ class DynamicSparseMatrix inline int rows() const { return IsRowMajor ? outerSize() : m_innerSize; } inline int cols() const { return IsRowMajor ? m_innerSize : outerSize(); } inline int innerSize() const { return m_innerSize; } - inline int outerSize() const { return m_data.size(); } + inline int outerSize() const { return static_cast<int>(m_data.size()); } inline int innerNonZeros(int j) const { return m_data[j].size(); } std::vector<CompressedStorage<Scalar> >& _data() { return m_data; } @@ -122,7 +129,7 @@ class DynamicSparseMatrix { int res = 0; for (int j=0; j<outerSize(); ++j) - res += m_data[j].size(); + res += static_cast<int>(m_data[j].size()); return res; } @@ -189,7 +196,7 @@ class DynamicSparseMatrix const int inner = IsRowMajor ? col : row; int startId = 0; - int id = m_data[outer].size() - 1; + int id = static_cast<int>(m_data[outer].size()) - 1; m_data[outer].resize(id+2,1); while ( (id >= startId) && (m_data[outer].index(id) > inner) ) @@ -209,7 +216,7 @@ class DynamicSparseMatrix inline void finalize() {} - void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) + void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar>()) { for (int j=0; j<outerSize(); ++j) m_data[j].prune(reference,epsilon); @@ -314,7 +321,6 @@ class DynamicSparseMatrix<Scalar,_Flags>::InnerIterator : public SparseVector<Sc inline int row() const { return IsRowMajor ? m_outer : Base::index(); } inline int col() const { return IsRowMajor ? Base::index() : m_outer; } - protected: const int m_outer; }; diff --git a/Eigen/src/Sparse/RandomSetter.h b/Eigen/src/Sparse/RandomSetter.h index 50824eba1..b34ca19a8 100644 --- a/Eigen/src/Sparse/RandomSetter.h +++ b/Eigen/src/Sparse/RandomSetter.h @@ -322,7 +322,7 @@ class RandomSetter { int nz = 0; for (int k=0; k<m_outerPackets; ++k) - nz += m_hashmaps[k].size(); + nz += static_cast<int>(m_hashmaps[k].size()); return nz; } diff --git a/Eigen/src/Sparse/SparseDiagonalProduct.h b/Eigen/src/Sparse/SparseDiagonalProduct.h index 642db5789..ceb4d6576 100644 --- a/Eigen/src/Sparse/SparseDiagonalProduct.h +++ b/Eigen/src/Sparse/SparseDiagonalProduct.h @@ -86,8 +86,7 @@ class SparseDiagonalProduct typedef ei_sparse_diagonal_product_inner_iterator_selector <_LhsNested,_RhsNested,SparseDiagonalProduct,LhsMode,RhsMode> InnerIterator; - template<typename _Lhs, typename _Rhs> - EIGEN_STRONG_INLINE SparseDiagonalProduct(const _Lhs& lhs, const _Rhs& rhs) + EIGEN_STRONG_INLINE SparseDiagonalProduct(const Lhs& lhs, const Rhs& rhs) : m_lhs(lhs), m_rhs(rhs) { ei_assert(lhs.cols() == rhs.rows() && "invalid sparse matrix * diagonal matrix product"); @@ -156,16 +155,16 @@ class ei_sparse_diagonal_product_inner_iterator_selector : public CwiseBinaryOp< ei_scalar_product_op<typename Rhs::Scalar>, SparseInnerVectorSet<Lhs,1>, - NestByValue<Transpose<typename Rhs::DiagonalVectorType> > >::InnerIterator + Transpose<typename Rhs::DiagonalVectorType> >::InnerIterator { typedef typename CwiseBinaryOp< ei_scalar_product_op<typename Rhs::Scalar>, SparseInnerVectorSet<Lhs,1>, - NestByValue<Transpose<typename Rhs::DiagonalVectorType> > >::InnerIterator Base; + Transpose<typename Rhs::DiagonalVectorType> >::InnerIterator Base; public: inline ei_sparse_diagonal_product_inner_iterator_selector( const SparseDiagonalProductType& expr, int outer) - : Base(expr.lhs().innerVector(outer) .cwiseProduct(expr.rhs().diagonal().transpose().nestByValue()), 0) + : Base(expr.lhs().innerVector(outer) .cwiseProduct(expr.rhs().diagonal().transpose()), 0) {} }; diff --git a/Eigen/src/Sparse/SparseExpressionMaker.h b/Eigen/src/Sparse/SparseExpressionMaker.h new file mode 100644 index 000000000..8e31d55ef --- /dev/null +++ b/Eigen/src/Sparse/SparseExpressionMaker.h @@ -0,0 +1,42 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_SPARSE_EXPRESSIONMAKER_H +#define EIGEN_SPARSE_EXPRESSIONMAKER_H + +template<typename Func, typename XprType> +struct MakeCwiseUnaryOp<Func,XprType,IsSparse> +{ + typedef SparseCwiseUnaryOp<Func,XprType> Type; +}; + +template<typename Func, typename A, typename B> +struct MakeCwiseBinaryOp<Func,A,B,IsSparse> +{ + typedef SparseCwiseBinaryOp<Func,A,B> Type; +}; + +// TODO complete the list + +#endif // EIGEN_SPARSE_EXPRESSIONMAKER_H diff --git a/Eigen/src/Sparse/SparseLDLT.h b/Eigen/src/Sparse/SparseLDLT.h index a30bd5ccd..625357e2b 100644 --- a/Eigen/src/Sparse/SparseLDLT.h +++ b/Eigen/src/Sparse/SparseLDLT.h @@ -94,7 +94,7 @@ class SparseLDLT : m_flags(flags), m_status(0) { ei_assert((MatrixType::Flags&RowMajorBit)==0); - m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>(); + m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>(); } /** Creates a LDLT object and compute the respective factorization of \a matrix using @@ -103,7 +103,7 @@ class SparseLDLT : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0) { ei_assert((MatrixType::Flags&RowMajorBit)==0); - m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>(); + m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>(); compute(matrix); } diff --git a/Eigen/src/Sparse/SparseLLT.h b/Eigen/src/Sparse/SparseLLT.h index 6307b2493..68ae43022 100644 --- a/Eigen/src/Sparse/SparseLLT.h +++ b/Eigen/src/Sparse/SparseLLT.h @@ -54,7 +54,7 @@ class SparseLLT SparseLLT(int flags = 0) : m_flags(flags), m_status(0) { - m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>(); + m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>(); } /** Creates a LLT object and compute the respective factorization of \a matrix using @@ -62,7 +62,7 @@ class SparseLLT SparseLLT(const MatrixType& matrix, int flags = 0) : m_matrix(matrix.rows(), matrix.cols()), m_flags(flags), m_status(0) { - m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>(); + m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>(); compute(matrix); } diff --git a/Eigen/src/Sparse/SparseLU.h b/Eigen/src/Sparse/SparseLU.h index 3f8d0f8db..2ec6d0e74 100644 --- a/Eigen/src/Sparse/SparseLU.h +++ b/Eigen/src/Sparse/SparseLU.h @@ -59,7 +59,7 @@ class SparseLU SparseLU(int flags = 0) : m_flags(flags), m_status(0) { - m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>(); + m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>(); } /** Creates a LU object and compute the respective factorization of \a matrix using @@ -67,7 +67,7 @@ class SparseLU SparseLU(const MatrixType& matrix, int flags = 0) : /*m_matrix(matrix.rows(), matrix.cols()),*/ m_flags(flags), m_status(0) { - m_precision = RealScalar(0.1) * Eigen::precision<RealScalar>(); + m_precision = RealScalar(0.1) * Eigen::dummy_precision<RealScalar>(); compute(matrix); } diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index 5e89d3f53..7feb2aff8 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -58,6 +58,13 @@ struct ei_traits<SparseMatrix<_Scalar, _Options> > }; template<typename _Scalar, int _Options> +struct ei_ref_selector<SparseMatrix<_Scalar, _Options> > +{ + typedef SparseMatrix<_Scalar, _Options> MatrixType; + typedef MatrixType const& type; +}; + +template<typename _Scalar, int _Options> class SparseMatrix : public SparseMatrixBase<SparseMatrix<_Scalar, _Options> > { @@ -132,7 +139,7 @@ class SparseMatrix } /** \returns the number of non zero coefficients */ - inline int nonZeros() const { return m_data.size(); } + inline int nonZeros() const { return static_cast<int>(m_data.size()); } /** \deprecated use setZero() and reserve() * Initializes the filling process of \c *this. @@ -230,7 +237,7 @@ class SparseMatrix // we start a new inner vector while (previousOuter>=0 && m_outerIndex[previousOuter]==0) { - m_outerIndex[previousOuter] = m_data.size(); + m_outerIndex[previousOuter] = static_cast<int>(m_data.size()); --previousOuter; } m_outerIndex[outer+1] = m_outerIndex[outer]; @@ -329,7 +336,7 @@ class SparseMatrix */ inline void finalize() { - int size = m_data.size(); + int size = static_cast<int>(m_data.size()); int i = m_outerSize; // find the last filled column while (i>=0 && m_outerIndex[i]==0) @@ -342,7 +349,7 @@ class SparseMatrix } } - void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) + void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar>()) { int k = 0; for (int j=0; j<m_outerSize; ++j) diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index 4ea1023af..80cb7620c 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -105,7 +105,7 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived // typedef SparseCwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType; /** \internal the return type of MatrixBase::adjoint() */ typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex, - CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, SparseNestByValue<Eigen::Transpose<Derived> > >, + CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, Eigen::Transpose<Derived> >, Transpose<Derived> >::ret AdjointReturnType; @@ -399,7 +399,7 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived Transpose<Derived> transpose() { return derived(); } const Transpose<Derived> transpose() const { return derived(); } // void transposeInPlace(); - const AdjointReturnType adjoint() const { return transpose().nestByValue(); } + const AdjointReturnType adjoint() const { return transpose(); } // sub-vector SparseInnerVectorSet<Derived,1> row(int i); @@ -510,32 +510,32 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived template<typename OtherDerived> bool isApprox(const SparseMatrixBase<OtherDerived>& other, - RealScalar prec = precision<Scalar>()) const + RealScalar prec = dummy_precision<Scalar>()) const { return toDense().isApprox(other.toDense(),prec); } template<typename OtherDerived> bool isApprox(const MatrixBase<OtherDerived>& other, - RealScalar prec = precision<Scalar>()) const + RealScalar prec = dummy_precision<Scalar>()) const { return toDense().isApprox(other,prec); } // bool isMuchSmallerThan(const RealScalar& other, -// RealScalar prec = precision<Scalar>()) const; +// RealScalar prec = dummy_precision<Scalar>()) const; // template<typename OtherDerived> // bool isMuchSmallerThan(const MatrixBase<OtherDerived>& other, -// RealScalar prec = precision<Scalar>()) const; +// RealScalar prec = dummy_precision<Scalar>()) const; -// bool isApproxToConstant(const Scalar& value, RealScalar prec = precision<Scalar>()) const; -// bool isZero(RealScalar prec = precision<Scalar>()) const; -// bool isOnes(RealScalar prec = precision<Scalar>()) const; -// bool isIdentity(RealScalar prec = precision<Scalar>()) const; -// bool isDiagonal(RealScalar prec = precision<Scalar>()) const; +// bool isApproxToConstant(const Scalar& value, RealScalar prec = dummy_precision<Scalar>()) const; +// bool isZero(RealScalar prec = dummy_precision<Scalar>()) const; +// bool isOnes(RealScalar prec = dummy_precision<Scalar>()) const; +// bool isIdentity(RealScalar prec = dummy_precision<Scalar>()) const; +// bool isDiagonal(RealScalar prec = dummy_precision<Scalar>()) const; -// bool isUpperTriangular(RealScalar prec = precision<Scalar>()) const; -// bool isLowerTriangular(RealScalar prec = precision<Scalar>()) const; +// bool isUpperTriangular(RealScalar prec = dummy_precision<Scalar>()) const; +// bool isLowerTriangular(RealScalar prec = dummy_precision<Scalar>()) const; // template<typename OtherDerived> // bool isOrthogonal(const MatrixBase<OtherDerived>& other, -// RealScalar prec = precision<Scalar>()) const; -// bool isUnitary(RealScalar prec = precision<Scalar>()) const; +// RealScalar prec = dummy_precision<Scalar>()) const; +// bool isUnitary(RealScalar prec = dummy_precision<Scalar>()) const; // template<typename OtherDerived> // inline bool operator==(const MatrixBase<OtherDerived>& other) const @@ -571,9 +571,7 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived */ // inline int stride(void) const { return derived().stride(); } - inline const SparseNestByValue<Derived> nestByValue() const; - - +// FIXME // ConjugateReturnType conjugate() const; // const RealReturnType real() const; // const ImagReturnType imag() const; @@ -626,11 +624,11 @@ template<typename Derived> class SparseMatrixBase : public AnyMatrixBase<Derived const MatrixBase<ElseDerived>& elseMatrix) const; template<typename ThenDerived> - inline const Select<Derived,ThenDerived, NestByValue<typename ThenDerived::ConstantReturnType> > + inline const Select<Derived,ThenDerived, typename ThenDerived::ConstantReturnType> select(const MatrixBase<ThenDerived>& thenMatrix, typename ThenDerived::Scalar elseScalar) const; template<typename ElseDerived> - inline const Select<Derived, NestByValue<typename ElseDerived::ConstantReturnType>, ElseDerived > + inline const Select<Derived, typename ElseDerived::ConstantReturnType, ElseDerived > select(typename ElseDerived::Scalar thenScalar, const MatrixBase<ElseDerived>& elseMatrix) const; template<int p> RealScalar lpNorm() const; diff --git a/Eigen/src/Sparse/SparseNestByValue.h b/Eigen/src/Sparse/SparseNestByValue.h deleted file mode 100644 index b48277232..000000000 --- a/Eigen/src/Sparse/SparseNestByValue.h +++ /dev/null @@ -1,84 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr> -// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_SPARSENESTBYVALUE_H -#define EIGEN_SPARSENESTBYVALUE_H - -/** \class SparseNestByValue - * - * \brief Expression which must be nested by value - * - * \param ExpressionType the type of the object of which we are requiring nesting-by-value - * - * This class is the return type of MatrixBase::nestByValue() - * and most of the time this is the only way it is used. - * - * \sa SparseMatrixBase::nestByValue(), class NestByValue - */ -template<typename ExpressionType> -struct ei_traits<SparseNestByValue<ExpressionType> > : public ei_traits<ExpressionType> -{}; - -template<typename ExpressionType> class SparseNestByValue - : public SparseMatrixBase<SparseNestByValue<ExpressionType> > -{ - public: - - typedef typename ExpressionType::InnerIterator InnerIterator; - - EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseNestByValue) - - inline SparseNestByValue(const ExpressionType& matrix) : m_expression(matrix) {} - - EIGEN_STRONG_INLINE int rows() const { return m_expression.rows(); } - EIGEN_STRONG_INLINE int cols() const { return m_expression.cols(); } - - operator const ExpressionType&() const { return m_expression; } - - protected: - const ExpressionType m_expression; -}; - -/** \returns an expression of the temporary version of *this. - */ -template<typename Derived> -inline const SparseNestByValue<Derived> -SparseMatrixBase<Derived>::nestByValue() const -{ - return SparseNestByValue<Derived>(derived()); -} - -// template<typename MatrixType> -// class SparseNestByValue<MatrixType>::InnerIterator : public MatrixType::InnerIterator -// { -// typedef typename MatrixType::InnerIterator Base; -// public: -// -// EIGEN_STRONG_INLINE InnerIterator(const SparseNestByValue& expr, int outer) -// : Base(expr.m_expression, outer) -// {} -// }; - -#endif // EIGEN_SPARSENESTBYVALUE_H diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index 6da1bee25..965565b73 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -123,11 +123,10 @@ template<typename _Scalar, int _Flags = 0> class DynamicSparseMatrix; template<typename _Scalar, int _Flags = 0> class SparseVector; template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix; -template<typename MatrixType> class SparseNestByValue; -template<typename MatrixType, int Size> class SparseInnerVectorSet; -template<typename MatrixType, int Mode> class SparseTriangularView; -template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView; -template<typename Lhs, typename Rhs> class SparseDiagonalProduct; +template<typename MatrixType, int Size> class SparseInnerVectorSet; +template<typename MatrixType, int Mode> class SparseTriangularView; +template<typename MatrixType, unsigned int UpLo> class SparseSelfAdjointView; +template<typename Lhs, typename Rhs> class SparseDiagonalProduct; template<typename Lhs, typename Rhs> class SparseProduct; @@ -156,6 +155,4 @@ template<typename T> class ei_eval<T,Sparse> typedef SparseMatrix<_Scalar, _Flags> type; }; -template<typename T> struct ei_must_nest_by_value<SparseNestByValue<T> > { enum { ret = true }; }; - #endif // EIGEN_SPARSEUTIL_H diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h index 68e984fae..5ac3f6be7 100644 --- a/Eigen/src/Sparse/SparseVector.h +++ b/Eigen/src/Sparse/SparseVector.h @@ -53,6 +53,13 @@ struct ei_traits<SparseVector<_Scalar, _Options> > }; template<typename _Scalar, int _Options> +struct ei_ref_selector< SparseVector<_Scalar, _Options> > +{ + typedef SparseVector<_Scalar, _Options> MatrixType; + typedef MatrixType const& type; +}; + +template<typename _Scalar, int _Options> class SparseVector : public SparseMatrixBase<SparseVector<_Scalar, _Options> > { @@ -119,7 +126,7 @@ class SparseVector inline void setZero() { m_data.clear(); } /** \returns the number of non zero coefficients */ - inline int nonZeros() const { return m_data.size(); } + inline int nonZeros() const { return static_cast<int>(m_data.size()); } inline void startVec(int outer) { @@ -202,7 +209,7 @@ class SparseVector EIGEN_DEPRECATED void endFill() {} inline void finalize() {} - void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>()) + void prune(Scalar reference, RealScalar epsilon = dummy_precision<RealScalar>()) { m_data.prune(reference,epsilon); } @@ -368,13 +375,13 @@ class SparseVector<Scalar,_Options>::InnerIterator { public: InnerIterator(const SparseVector& vec, int outer=0) - : m_data(vec.m_data), m_id(0), m_end(m_data.size()) + : m_data(vec.m_data), m_id(0), m_end(static_cast<int>(m_data.size())) { ei_assert(outer==0); } InnerIterator(const CompressedStorage<Scalar>& data) - : m_data(data), m_id(0), m_end(m_data.size()) + : m_data(data), m_id(0), m_end(static_cast<int>(m_data.size())) {} template<unsigned int Added, unsigned int Removed> diff --git a/Eigen/src/misc/Image.h b/Eigen/src/misc/Image.h index 9ed5d5f70..2f39d6b9d 100644 --- a/Eigen/src/misc/Image.h +++ b/Eigen/src/misc/Image.h @@ -67,9 +67,9 @@ template<typename _DecompositionType> struct ei_image_retval_base } protected: - const DecompositionType& m_dec; - int m_rank, m_cols; - const MatrixType& m_originalMatrix; + const DecompositionType& m_dec; + int m_rank, m_cols; + const MatrixType& m_originalMatrix; }; #define EIGEN_MAKE_IMAGE_HELPERS(DecompositionType) \ diff --git a/Eigen/src/misc/Kernel.h b/Eigen/src/misc/Kernel.h index 717eef450..908c408e9 100644 --- a/Eigen/src/misc/Kernel.h +++ b/Eigen/src/misc/Kernel.h @@ -67,8 +67,8 @@ template<typename _DecompositionType> struct ei_kernel_retval_base } protected: - const DecompositionType& m_dec; - int m_rank, m_cols; + const DecompositionType& m_dec; + int m_rank, m_cols; }; #define EIGEN_MAKE_KERNEL_HELPERS(DecompositionType) \ diff --git a/Eigen/src/misc/Solve.h b/Eigen/src/misc/Solve.h index d93869121..4ab0775fc 100644 --- a/Eigen/src/misc/Solve.h +++ b/Eigen/src/misc/Solve.h @@ -61,8 +61,8 @@ template<typename _DecompositionType, typename Rhs> struct ei_solve_retval_base } protected: - const DecompositionType& m_dec; - const typename Rhs::Nested m_rhs; + const DecompositionType& m_dec; + const typename Rhs::Nested m_rhs; }; #define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType,Rhs) \ |