diff options
author | 2009-11-09 09:08:03 -0500 | |
---|---|---|
committer | 2009-11-09 09:08:03 -0500 | |
commit | 92749eed11d000300cfa54654f1043cd52399ed8 (patch) | |
tree | ba227522582b2f9f4280ed1404e74c654e21ccb3 /Eigen | |
parent | 4b366b07be4e409239c61158a23d93e8ebf3811b (diff) | |
parent | 670651e2e0932c5edfe2a2da4b9f3c42af3b7dec (diff) |
* merge
* remove a ctor in QuaternionBase as it gives a strange error with GCC 4.4.2.
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/Eigenvalues | 1 | ||||
-rw-r--r-- | Eigen/src/Array/VectorwiseOp.h | 8 | ||||
-rw-r--r-- | Eigen/src/Core/Assign.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/Functors.h | 6 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 29 | ||||
-rw-r--r-- | Eigen/src/Core/NumTraits.h | 7 | ||||
-rw-r--r-- | Eigen/src/Core/Redux.h | 15 | ||||
-rw-r--r-- | Eigen/src/Core/Transpose.h | 3 | ||||
-rw-r--r-- | Eigen/src/Core/arch/SSE/PacketMath.h | 10 | ||||
-rw-r--r-- | Eigen/src/Core/util/Memory.h | 28 | ||||
-rw-r--r-- | Eigen/src/Core/util/Meta.h | 7 | ||||
-rw-r--r-- | Eigen/src/Core/util/StaticAssert.h | 1 | ||||
-rw-r--r-- | Eigen/src/Core/util/XprHelper.h | 2 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/ComplexSchur.h | 16 | ||||
-rw-r--r-- | Eigen/src/Geometry/Quaternion.h | 224 | ||||
-rw-r--r-- | Eigen/src/SVD/SVD.h | 11 | ||||
-rw-r--r-- | Eigen/src/Sparse/CholmodSupport.h | 5 |
17 files changed, 236 insertions, 139 deletions
diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues index 9a6443f39..8c6841549 100644 --- a/Eigen/Eigenvalues +++ b/Eigen/Eigenvalues @@ -8,6 +8,7 @@ #include "Cholesky" #include "Jacobi" #include "Householder" +#include "LU" // Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module #if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2) diff --git a/Eigen/src/Array/VectorwiseOp.h b/Eigen/src/Array/VectorwiseOp.h index 4cb0083fa..880567212 100644 --- a/Eigen/src/Array/VectorwiseOp.h +++ b/Eigen/src/Array/VectorwiseOp.h @@ -122,6 +122,7 @@ 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(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); EIGEN_MEMBER_FUNCTOR(maxCoeff, (Size-1)*NumTraits<Scalar>::AddCost); EIGEN_MEMBER_FUNCTOR(all, (Size-1)*NumTraits<Scalar>::AddCost); @@ -297,6 +298,13 @@ template<typename ExpressionType, int Direction> class VectorwiseOp const typename ReturnType<ei_member_sum>::Type sum() const { return _expression(); } + /** \returns a row (or column) vector expression of the mean + * of each column (or row) of the referenced expression. + * + * \sa MatrixBase::mean() */ + const typename ReturnType<ei_member_mean>::Type mean() const + { return _expression(); } + /** \returns a row (or column) vector expression representing * whether \b all coefficients of each respective column (or row) are \c true. * diff --git a/Eigen/src/Core/Assign.h b/Eigen/src/Core/Assign.h index 4bd1046a7..8dc015715 100644 --- a/Eigen/src/Core/Assign.h +++ b/Eigen/src/Core/Assign.h @@ -93,7 +93,7 @@ public: ? ( int(MayUnrollCompletely) && int(DstIsAligned) ? int(CompleteUnrolling) : int(NoUnrolling) ) : int(NoUnrolling) }; - + static void debug() { EIGEN_DEBUG_VAR(DstIsAligned) diff --git a/Eigen/src/Core/Functors.h b/Eigen/src/Core/Functors.h index cbaeb83e2..259f40244 100644 --- a/Eigen/src/Core/Functors.h +++ b/Eigen/src/Core/Functors.h @@ -350,7 +350,7 @@ struct ei_scalar_multiple_op { EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a * m_other; } EIGEN_STRONG_INLINE const PacketScalar packetOp(const PacketScalar& a) const { return ei_pmul(a, ei_pset1(m_other)); } - const Scalar m_other; + typename ei_makeconst<typename NumTraits<Scalar>::Nested>::type m_other; private: ei_scalar_multiple_op& operator=(const ei_scalar_multiple_op&); }; @@ -364,7 +364,7 @@ struct ei_scalar_multiple2_op { EIGEN_STRONG_INLINE ei_scalar_multiple2_op(const ei_scalar_multiple2_op& other) : m_other(other.m_other) { } EIGEN_STRONG_INLINE ei_scalar_multiple2_op(const Scalar2& other) : m_other(other) { } EIGEN_STRONG_INLINE result_type operator() (const Scalar1& a) const { return a * m_other; } - const Scalar2 m_other; + typename ei_makeconst<typename NumTraits<Scalar2>::Nested>::type m_other; }; template<typename Scalar1,typename Scalar2> struct ei_functor_traits<ei_scalar_multiple2_op<Scalar1,Scalar2> > @@ -393,7 +393,7 @@ struct ei_scalar_quotient1_impl<Scalar,false> { EIGEN_STRONG_INLINE ei_scalar_quotient1_impl(const ei_scalar_quotient1_impl& other) : m_other(other.m_other) { } EIGEN_STRONG_INLINE ei_scalar_quotient1_impl(const Scalar& other) : m_other(other) {} EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a / m_other; } - const Scalar m_other; + typename ei_makeconst<typename NumTraits<Scalar>::Nested>::type m_other; }; template<typename Scalar> struct ei_functor_traits<ei_scalar_quotient1_impl<Scalar,false> > diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 7b5310175..e5eed715b 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -145,12 +145,6 @@ template<typename Derived> class MatrixBase #endif }; - /** Default constructor. Just checks at compile-time for self-consistency of the flags. */ - MatrixBase() - { - ei_assert(ei_are_flags_consistent<Flags>::ret); - } - #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is the "real scalar" type; if the \a Scalar type is already real numbers * (e.g. int, float or double) then \a RealScalar is just the same as \a Scalar. If @@ -177,7 +171,7 @@ template<typename Derived> class MatrixBase inline int diagonalSize() const { return std::min(rows(),cols()); } /** \returns the number of nonzero coefficients which is in practice the number * of stored coefficients. */ - inline int nonZeros() const { return derived().nonZeros(); } + inline int nonZeros() const { return size(); } /** \returns true if either the number of rows or the number of columns is equal to 1. * In other words, this function returns * \code rows()==1 || cols()==1 \endcode @@ -645,8 +639,9 @@ template<typename Derived> class MatrixBase const CwiseBinaryOp<CustomBinaryOp, Derived, OtherDerived> binaryExpr(const MatrixBase<OtherDerived> &other, const CustomBinaryOp& func = CustomBinaryOp()) const; - + Scalar sum() const; + Scalar mean() const; Scalar trace() const; Scalar prod() const; @@ -811,6 +806,24 @@ template<typename Derived> class MatrixBase #ifdef EIGEN_MATRIXBASE_PLUGIN #include EIGEN_MATRIXBASE_PLUGIN #endif + + protected: + /** Default constructor. Do nothing. */ + MatrixBase() + { + /* Just checks for self-consistency of the flags. + * Only do it when debugging Eigen, as this borders on paranoiac and could slow compilation down + */ +#ifdef EIGEN_INTERNAL_DEBUGGING + EIGEN_STATIC_ASSERT(ei_are_flags_consistent<Flags>::ret, + INVALID_MATRIXBASE_TEMPLATE_PARAMETERS) +#endif + } + + private: + explicit MatrixBase(int); + MatrixBase(int,int); + template<typename OtherDerived> explicit MatrixBase(const MatrixBase<OtherDerived>&); }; #endif // EIGEN_MATRIXBASE_H diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h index 24afe54c5..304e2c1d6 100644 --- a/Eigen/src/Core/NumTraits.h +++ b/Eigen/src/Core/NumTraits.h @@ -52,6 +52,7 @@ template<> struct NumTraits<int> { typedef int Real; typedef double FloatingPoint; + typedef int Nested; enum { IsComplex = 0, HasFloatingPoint = 0, @@ -65,6 +66,7 @@ template<> struct NumTraits<float> { typedef float Real; typedef float FloatingPoint; + typedef float Nested; enum { IsComplex = 0, HasFloatingPoint = 1, @@ -78,6 +80,7 @@ template<> struct NumTraits<double> { typedef double Real; typedef double FloatingPoint; + typedef double Nested; enum { IsComplex = 0, HasFloatingPoint = 1, @@ -91,6 +94,7 @@ template<typename _Real> struct NumTraits<std::complex<_Real> > { typedef _Real Real; typedef std::complex<_Real> FloatingPoint; + typedef std::complex<_Real> Nested; enum { IsComplex = 1, HasFloatingPoint = NumTraits<Real>::HasFloatingPoint, @@ -104,6 +108,7 @@ template<> struct NumTraits<long long int> { typedef long long int Real; typedef long double FloatingPoint; + typedef long long int Nested; enum { IsComplex = 0, HasFloatingPoint = 0, @@ -117,6 +122,7 @@ template<> struct NumTraits<long double> { typedef long double Real; typedef long double FloatingPoint; + typedef long double Nested; enum { IsComplex = 0, HasFloatingPoint = 1, @@ -130,6 +136,7 @@ template<> struct NumTraits<bool> { typedef bool Real; typedef float FloatingPoint; + typedef bool Nested; enum { IsComplex = 0, HasFloatingPoint = 0, diff --git a/Eigen/src/Core/Redux.h b/Eigen/src/Core/Redux.h index 0df095750..9f796157a 100644 --- a/Eigen/src/Core/Redux.h +++ b/Eigen/src/Core/Redux.h @@ -342,7 +342,7 @@ MatrixBase<Derived>::maxCoeff() const /** \returns the sum of all coefficients of *this * - * \sa trace(), prod() + * \sa trace(), prod(), mean() */ template<typename Derived> EIGEN_STRONG_INLINE typename ei_traits<Derived>::Scalar @@ -351,12 +351,23 @@ MatrixBase<Derived>::sum() const return this->redux(Eigen::ei_scalar_sum_op<Scalar>()); } +/** \returns the mean of all coefficients of *this +* +* \sa trace(), prod(), sum() +*/ +template<typename Derived> +EIGEN_STRONG_INLINE typename ei_traits<Derived>::Scalar +MatrixBase<Derived>::mean() const +{ + return this->redux(Eigen::ei_scalar_sum_op<Scalar>()) / this->size(); +} + /** \returns the product of all coefficients of *this * * Example: \include MatrixBase_prod.cpp * Output: \verbinclude MatrixBase_prod.out * - * \sa sum() + * \sa sum(), mean(), trace() */ template<typename Derived> EIGEN_STRONG_INLINE typename ei_traits<Derived>::Scalar diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 8527edc2a..990aa3807 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -69,7 +69,6 @@ template<typename MatrixType> class Transpose inline int rows() const { return m_matrix.cols(); } inline int cols() const { return m_matrix.rows(); } - inline int nonZeros() const { return m_matrix.nonZeros(); } inline int stride() const { return m_matrix.stride(); } inline Scalar* data() { return m_matrix.data(); } inline const Scalar* data() const { return m_matrix.data(); } @@ -354,5 +353,5 @@ lazyAssign(const CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerivedA,CwiseUnaryOp<ei return lazyAssign(static_cast<const MatrixBase<CwiseBinaryOp<ei_scalar_sum_op<Scalar>,DerivedA,CwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, NestByValue<Eigen::Transpose<DerivedB> > > > >& >(other)); } #endif - + #endif // EIGEN_TRANSPOSE_H diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index eb1c2d311..60ccadc21 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -220,8 +220,14 @@ template<> EIGEN_STRONG_INLINE void ei_pstoreu<double>(double* to, const Packet2 template<> EIGEN_STRONG_INLINE void ei_pstoreu<float>(float* to, const Packet4f& from) { ei_pstoreu((double*)to, _mm_castps_pd(from)); } template<> EIGEN_STRONG_INLINE void ei_pstoreu<int>(int* to, const Packet4i& from) { ei_pstoreu((double*)to, _mm_castsi128_pd(from)); } -#ifdef _MSC_VER -// this fix internal compilation error +#if defined(_MSC_VER) && (_MSC_VER <= 1500) && defined(_WIN64) +// The temporary variable fixes an internal compilation error. +// Direct of the struct members fixed bug #62. +template<> EIGEN_STRONG_INLINE float ei_pfirst<Packet4f>(const Packet4f& a) { return a.m128_f32[0]; } +template<> EIGEN_STRONG_INLINE double ei_pfirst<Packet2d>(const Packet2d& a) { return a.m128d_f64[0]; } +template<> EIGEN_STRONG_INLINE int ei_pfirst<Packet4i>(const Packet4i& a) { int x = _mm_cvtsi128_si32(a); return x; } +#elif defined(_MSC_VER) && (_MSC_VER <= 1500) +// The temporary variable fixes an internal compilation error. template<> EIGEN_STRONG_INLINE float ei_pfirst<Packet4f>(const Packet4f& a) { float x = _mm_cvtss_f32(a); return x; } template<> EIGEN_STRONG_INLINE double ei_pfirst<Packet2d>(const Packet2d& a) { double x = _mm_cvtsd_f64(a); return x; } template<> EIGEN_STRONG_INLINE int ei_pfirst<Packet4i>(const Packet4i& a) { int x = _mm_cvtsi128_si32(a); return x; } diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index f8581eebc..bc5235582 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -83,7 +83,7 @@ inline void* ei_aligned_malloc(size_t size) ei_assert(false && "heap allocation is forbidden (EIGEN_NO_MALLOC is defined)"); #endif - void *result; + void *result; #if !EIGEN_ALIGN result = malloc(size); #elif EIGEN_MALLOC_ALREADY_ALIGNED @@ -97,7 +97,7 @@ inline void* ei_aligned_malloc(size_t size) #else result = ei_handmade_aligned_malloc(size); #endif - + #ifdef EIGEN_EXCEPTIONS if(result == 0) throw std::bad_alloc(); @@ -324,34 +324,34 @@ public: typedef aligned_allocator<U> other; }; - pointer address( reference value ) const + pointer address( reference value ) const { return &value; } - const_pointer address( const_reference value ) const + const_pointer address( const_reference value ) const { return &value; } - aligned_allocator() throw() + aligned_allocator() throw() { } - aligned_allocator( const aligned_allocator& ) throw() + aligned_allocator( const aligned_allocator& ) throw() { } template<class U> - aligned_allocator( const aligned_allocator<U>& ) throw() + aligned_allocator( const aligned_allocator<U>& ) throw() { } - ~aligned_allocator() throw() + ~aligned_allocator() throw() { } - size_type max_size() const throw() + size_type max_size() const throw() { return std::numeric_limits<size_type>::max(); } @@ -362,24 +362,24 @@ public: return static_cast<pointer>( ei_aligned_malloc( num * sizeof(T) ) ); } - void construct( pointer p, const T& value ) + void construct( pointer p, const T& value ) { ::new( p ) T( value ); } - void destroy( pointer p ) + void destroy( pointer p ) { p->~T(); } - void deallocate( pointer p, size_type /*num*/ ) + void deallocate( pointer p, size_type /*num*/ ) { ei_aligned_free( p ); } - + bool operator!=(const aligned_allocator<T>& other) const { return false; } - + bool operator==(const aligned_allocator<T>& other) const { return true; } }; diff --git a/Eigen/src/Core/util/Meta.h b/Eigen/src/Core/util/Meta.h index 3a960bea6..2fdfd93b5 100644 --- a/Eigen/src/Core/util/Meta.h +++ b/Eigen/src/Core/util/Meta.h @@ -64,6 +64,13 @@ template<typename T> struct ei_cleantype<T&> { typedef typename ei_cleant template<typename T> struct ei_cleantype<const T*> { typedef typename ei_cleantype<T>::type type; }; template<typename T> struct ei_cleantype<T*> { typedef typename ei_cleantype<T>::type type; }; +template<typename T> struct ei_makeconst { typedef const T type; }; +template<typename T> struct ei_makeconst<const T> { typedef const T type; }; +template<typename T> struct ei_makeconst<T&> { typedef const T& type; }; +template<typename T> struct ei_makeconst<const T&> { typedef const T& type; }; +template<typename T> struct ei_makeconst<T*> { typedef const T* type; }; +template<typename T> struct ei_makeconst<const T*> { typedef const T* type; }; + /** \internal Allows to enable/disable an overload * according to a compile time condition. */ diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h index 6210bc91c..56a57b706 100644 --- a/Eigen/src/Core/util/StaticAssert.h +++ b/Eigen/src/Core/util/StaticAssert.h @@ -76,6 +76,7 @@ THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES, THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES, INVALID_MATRIX_TEMPLATE_PARAMETERS, + INVALID_MATRIXBASE_TEMPLATE_PARAMETERS, BOTH_MATRICES_MUST_HAVE_THE_SAME_STORAGE_ORDER, THIS_METHOD_IS_ONLY_FOR_DIAGONAL_MATRIX, THE_MATRIX_OR_EXPRESSION_THAT_YOU_PASSED_DOES_NOT_HAVE_THE_EXPECTED_TYPE, diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index cea2faaa8..be4266f85 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -1,4 +1,4 @@ -// This file is part of Eigen, a lightweight C++ template library +// // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 0534715c4..a25af342d 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -167,10 +167,11 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) //locate the range in which to iterate while(iu > 0) { - d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1)); - sd = ei_norm1(m_matT.coeffRef(iu,iu-1)); + d = ei_norm1(m_matT.coeff(iu,iu)) + ei_norm1(m_matT.coeff(iu-1,iu-1)); + sd = ei_norm1(m_matT.coeff(iu,iu-1)); - if(sd >= eps * d) break; // FIXME : precision criterion ?? + if(!ei_isMuchSmallerThan(sd,d,eps)) + break; m_matT.coeffRef(iu,iu-1) = Complex(0); iter = 0; @@ -187,13 +188,14 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) } il = iu-1; - while( il > 0 ) + while(il > 0) { // check if the current 2x2 block on the diagonal is upper triangular - d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1)); - sd = ei_norm1(m_matT.coeffRef(il,il-1)); + d = ei_norm1(m_matT.coeff(il,il)) + ei_norm1(m_matT.coeff(il-1,il-1)); + sd = ei_norm1(m_matT.coeff(il,il-1)); - if(sd < eps * d) break; // FIXME : precision criterion ?? + if(ei_isMuchSmallerThan(sd,d,eps)) + break; --il; } diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 67b040165..b08a027c9 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -26,53 +26,27 @@ #ifndef EIGEN_QUATERNION_H #define EIGEN_QUATERNION_H -/** \geometry_module \ingroup Geometry_Module - * - * \class Quaternion - * - * \brief The quaternion class used to represent 3D orientations and rotations - * - * \param _Scalar the scalar type, i.e., the type of the coefficients - * - * This class represents a quaternion \f$ w+xi+yj+zk \f$ that is a convenient representation of - * orientations and rotations of objects in three dimensions. Compared to other representations - * like Euler angles or 3x3 matrices, quatertions offer the following advantages: - * \li \b compact storage (4 scalars) - * \li \b efficient to compose (28 flops), - * \li \b stable spherical interpolation - * - * The following two typedefs are provided for convenience: - * \li \c Quaternionf for \c float - * \li \c Quaterniond for \c double - * - * \sa class AngleAxis, class Transform - */ +/*************************************************************************** +* Definition of QuaternionBase<Derived> +* The implementation is at the end of the file +***************************************************************************/ template<typename Other, int OtherRows=Other::RowsAtCompileTime, int OtherCols=Other::ColsAtCompileTime> struct ei_quaternionbase_assign_impl; -template<typename Scalar> class Quaternion; // [XXX] => remove when Quaternion becomes Quaternion - -template<typename Derived> -struct ei_traits<QuaternionBase<Derived> > -{ - typedef typename ei_traits<Derived>::Scalar Scalar; - enum { - PacketAccess = ei_traits<Derived>::PacketAccess - }; -}; - template<class Derived> class QuaternionBase : public RotationBase<Derived, 3> { typedef RotationBase<Derived, 3> Base; public: using Base::operator*; + using Base::derived; - typedef typename ei_traits<QuaternionBase<Derived> >::Scalar Scalar; + typedef typename ei_traits<Derived>::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; + typedef typename ei_traits<Derived>::Coefficients Coefficients; // typedef typename Matrix<Scalar,4,1> Coefficients; /** the type of a 3D vector */ @@ -82,6 +56,8 @@ public: /** the equivalent angle-axis type */ typedef AngleAxis<Scalar> AngleAxisType; + + /** \returns the \c x coefficient */ inline Scalar x() const { return this->derived().coeffs().coeff(0); } /** \returns the \c y coefficient */ @@ -101,45 +77,52 @@ public: inline Scalar& w() { return this->derived().coeffs().coeffRef(3); } /** \returns a read-only vector expression of the imaginary part (x,y,z) */ - inline const VectorBlock<typename ei_traits<Derived>::Coefficients,3> vec() const { return this->derived().coeffs().template start<3>(); } + inline const VectorBlock<Coefficients,3> vec() const { return coeffs().template start<3>(); } /** \returns a vector expression of the imaginary part (x,y,z) */ - inline VectorBlock<typename ei_traits<Derived>::Coefficients,3> vec() { return this->derived().coeffs().template start<3>(); } + inline VectorBlock<Coefficients,3> vec() { return coeffs().template start<3>(); } /** \returns a read-only vector expression of the coefficients (x,y,z,w) */ - inline const typename ei_traits<Derived>::Coefficients& coeffs() const { return this->derived().coeffs(); } + inline const typename ei_traits<Derived>::Coefficients& coeffs() const { return derived().coeffs(); } /** \returns a vector expression of the coefficients (x,y,z,w) */ - inline typename ei_traits<Derived>::Coefficients& coeffs() { return this->derived().coeffs(); } + inline typename ei_traits<Derived>::Coefficients& coeffs() { return derived().coeffs(); } - template<class OtherDerived> QuaternionBase& operator=(const QuaternionBase<OtherDerived>& other); - QuaternionBase& operator=(const AngleAxisType& aa); - template<class OtherDerived> - QuaternionBase& operator=(const MatrixBase<OtherDerived>& m); + template<class OtherDerived> Derived& operator=(const QuaternionBase<OtherDerived>& other); + +// disabled this copy operator as it is giving very strange compilation errors when compiling +// test_stdvector with GCC 4.4.2. This looks like a GCC bug though, so feel free to re-enable it if it's +// useful; however notice that we already have the templated operator= above and e.g. in MatrixBase +// we didn't have to add, in addition to templated operator=, such a non-templated copy operator. +// Derived& operator=(const QuaternionBase& other) +// { return operator=<Derived>(other); } + + Derived& operator=(const AngleAxisType& aa); + template<class OtherDerived> Derived& operator=(const MatrixBase<OtherDerived>& m); /** \returns a quaternion representing an identity rotation * \sa MatrixBase::Identity() */ inline static Quaternion<Scalar> Identity() { return Quaternion<Scalar>(1, 0, 0, 0); } - /** \sa Quaternion2::Identity(), MatrixBase::setIdentity() + /** \sa QuaternionBase::Identity(), MatrixBase::setIdentity() */ inline QuaternionBase& setIdentity() { coeffs() << 0, 0, 0, 1; return *this; } /** \returns the squared norm of the quaternion's coefficients - * \sa Quaternion2::norm(), MatrixBase::squaredNorm() + * \sa QuaternionBase::norm(), MatrixBase::squaredNorm() */ inline Scalar squaredNorm() const { return coeffs().squaredNorm(); } /** \returns the norm of the quaternion's coefficients - * \sa Quaternion2::squaredNorm(), MatrixBase::norm() + * \sa QuaternionBase::squaredNorm(), MatrixBase::norm() */ inline Scalar norm() const { return coeffs().norm(); } /** Normalizes the quaternion \c *this * \sa normalized(), MatrixBase::normalize() */ inline void normalize() { coeffs().normalize(); } - /** \returns a normalized version of \c *this + /** \returns a normalized copy of \c *this * \sa normalize(), MatrixBase::normalized() */ inline Quaternion<Scalar> normalized() const { return Quaternion<Scalar>(coeffs().normalized()); } @@ -152,16 +135,16 @@ public: template<class OtherDerived> inline Scalar angularDistance(const QuaternionBase<OtherDerived>& other) const; - Matrix3 toRotationMatrix(void) const; + Matrix3 toRotationMatrix() const; template<typename Derived1, typename Derived2> - QuaternionBase& setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b); + Derived& setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b); template<class OtherDerived> inline Quaternion<Scalar> operator* (const QuaternionBase<OtherDerived>& q) const; - template<class OtherDerived> inline QuaternionBase& operator*= (const QuaternionBase<OtherDerived>& q); + template<class OtherDerived> inline Derived& operator*= (const QuaternionBase<OtherDerived>& q); - Quaternion<Scalar> inverse(void) const; - Quaternion<Scalar> conjugate(void) const; + Quaternion<Scalar> inverse() const; + Quaternion<Scalar> conjugate() const; template<class OtherDerived> Quaternion<Scalar> slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const; @@ -169,7 +152,8 @@ public: * determined by \a prec. * * \sa MatrixBase::isApprox() */ - bool isApprox(const QuaternionBase& other, RealScalar prec = precision<Scalar>()) const + template<class OtherDerived> + bool isApprox(const QuaternionBase<OtherDerived>& other, RealScalar prec = precision<Scalar>()) const { return coeffs().isApprox(other.coeffs(), prec); } Vector3 _transformVector(Vector3 v) const; @@ -187,13 +171,39 @@ public: } }; +/*************************************************************************** +* Definition/implementation of Quaternion<Scalar> +***************************************************************************/ + +/** \geometry_module \ingroup Geometry_Module + * + * \class Quaternion + * + * \brief The quaternion class used to represent 3D orientations and rotations + * + * \param _Scalar the scalar type, i.e., the type of the coefficients + * + * This class represents a quaternion \f$ w+xi+yj+zk \f$ that is a convenient representation of + * orientations and rotations of objects in three dimensions. Compared to other representations + * like Euler angles or 3x3 matrices, quatertions offer the following advantages: + * \li \b compact storage (4 scalars) + * \li \b efficient to compose (28 flops), + * \li \b stable spherical interpolation + * + * The following two typedefs are provided for convenience: + * \li \c Quaternionf for \c float + * \li \c Quaterniond for \c double + * + * \sa class AngleAxis, class Transform + */ + template<typename _Scalar> struct ei_traits<Quaternion<_Scalar> > { typedef _Scalar Scalar; typedef Matrix<_Scalar,4,1> Coefficients; enum{ - PacketAccess = Aligned + PacketAccess = Aligned }; }; @@ -239,7 +249,7 @@ public: explicit inline Quaternion(const MatrixBase<Derived>& other) { *this = other; } /** Copy constructor with scalar type conversion */ - template<class Derived> + template<typename Derived> inline explicit Quaternion(const QuaternionBase<Derived>& other) { m_coeffs = other.coeffs().template cast<Scalar>(); } @@ -250,16 +260,29 @@ protected: Coefficients m_coeffs; }; -/* ########### Map<Quaternion> */ +/** \ingroup Geometry_Module + * single precision quaternion type */ +typedef Quaternion<float> Quaternionf; +/** \ingroup Geometry_Module + * double precision quaternion type */ +typedef Quaternion<double> Quaterniond; + +/*************************************************************************** +* Specialization of Map<Quaternion<Scalar>> +***************************************************************************/ /** \class Map<Quaternion> * \nonstableyet * - * \brief Expression of a quaternion + * \brief Expression of a quaternion from a memory buffer * - * \param Scalar the type of the vector of diagonal coefficients + * \param _Scalar the type of the Quaternion coefficients + * \param PacketAccess see class Map * - * \sa class Quaternion, class QuaternionBase + * This is a specialization of class Map for Quaternion. This class allows to view + * a 4 scalar memory buffer as an Eigen's Quaternion object. + * + * \sa class Map, class Quaternion, class QuaternionBase */ template<typename _Scalar, int _PacketAccess> struct ei_traits<Map<Quaternion<_Scalar>, _PacketAccess> >: @@ -273,15 +296,23 @@ ei_traits<Quaternion<_Scalar> > }; template<typename _Scalar, int PacketAccess> -class Map<Quaternion<_Scalar>, PacketAccess > : public QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> >, ei_no_assignment_operator { +class Map<Quaternion<_Scalar>, PacketAccess > + : public QuaternionBase<Map<Quaternion<_Scalar>, PacketAccess> >, + ei_no_assignment_operator +{ public: - + typedef _Scalar Scalar; + typedef typename ei_traits<Map>::Coefficients Coefficients; - 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. */ + inline Map(const Scalar* coeffs) : m_coeffs(coeffs) {} - inline Map<Quaternion<Scalar>, PacketAccess >(const Scalar* coeffs) : m_coeffs(coeffs) {} - inline Coefficients& coeffs() { return m_coeffs;} inline const Coefficients& coeffs() const { return m_coeffs;} @@ -289,15 +320,20 @@ class Map<Quaternion<_Scalar>, PacketAccess > : public QuaternionBase<Map<Quater Coefficients m_coeffs; }; -typedef Map<Quaternion<double> > QuaternionMapd; -typedef Map<Quaternion<float> > QuaternionMapf; -typedef Map<Quaternion<double>, Aligned> QuaternionMapAlignedd; -typedef Map<Quaternion<float>, Aligned> QuaternionMapAlignedf; +typedef Map<Quaternion<double> > QuaternionMapd; +typedef Map<Quaternion<float> > QuaternionMapf; +typedef Map<Quaternion<double>, Aligned> QuaternionMapAlignedd; +typedef Map<Quaternion<float>, Aligned> QuaternionMapAlignedf; + +/*************************************************************************** +* Implementation of QuaternionBase methods +***************************************************************************/ // Generic Quaternion * Quaternion product -template<int Arch, class Derived, class OtherDerived, typename Scalar, int PacketAccess> struct ei_quat_product +// This product can be specialized for a given architecture via the Arch template argument. +template<int Arch, class Derived1, class Derived2, typename Scalar, int PacketAccess> struct ei_quat_product { - inline static Quaternion<Scalar> run(const QuaternionBase<Derived>& a, const QuaternionBase<OtherDerived>& b){ + inline static Quaternion<Scalar> run(const QuaternionBase<Derived1>& a, const QuaternionBase<Derived2>& b){ return Quaternion<Scalar> ( a.w() * b.w() - a.x() * b.x() - a.y() * b.y() - a.z() * b.z(), @@ -311,21 +347,22 @@ template<int Arch, class Derived, class OtherDerived, typename Scalar, int Packe /** \returns the concatenation of two rotations as a quaternion-quaternion product */ template <class Derived> template <class OtherDerived> -inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> QuaternionBase<Derived>::operator* (const QuaternionBase<OtherDerived>& other) const +inline Quaternion<typename ei_traits<Derived>::Scalar> +QuaternionBase<Derived>::operator* (const QuaternionBase<OtherDerived>& other) const { 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, - typename ei_traits<Derived>::Scalar, - ei_traits<Derived>::PacketAccess && ei_traits<OtherDerived>::PacketAccess>::run(*this, other); + return ei_quat_product<EiArch, Derived, OtherDerived, + typename ei_traits<Derived>::Scalar, + ei_traits<Derived>::PacketAccess && ei_traits<OtherDerived>::PacketAccess>::run(*this, other); } /** \sa operator*(Quaternion) */ template <class Derived> template <class OtherDerived> -inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator*= (const QuaternionBase<OtherDerived>& other) +inline Derived& QuaternionBase<Derived>::operator*= (const QuaternionBase<OtherDerived>& other) { - return (*this = *this * other); + return (derived() = derived() * other.derived()); } /** Rotation of a vector by a quaternion. @@ -350,21 +387,21 @@ QuaternionBase<Derived>::_transformVector(Vector3 v) const template<class Derived> template<class OtherDerived> -inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator=(const QuaternionBase<OtherDerived>& other) +inline Derived& QuaternionBase<Derived>::operator=(const QuaternionBase<OtherDerived>& other) { coeffs() = other.coeffs(); - return *this; + return derived(); } /** Set \c *this from an angle-axis \a aa and returns a reference to \c *this */ template<class Derived> -inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator=(const AngleAxisType& aa) +inline Derived& QuaternionBase<Derived>::operator=(const AngleAxisType& aa) { Scalar ha = Scalar(0.5)*aa.angle(); // Scalar(0.5) to suppress precision loss warnings this->w() = ei_cos(ha); this->vec() = ei_sin(ha) * aa.axis(); - return *this; + return derived(); } /** Set \c *this from the expression \a xpr: @@ -375,12 +412,12 @@ inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator=(const AngleAx template<class Derived> template<class MatrixDerived> -inline QuaternionBase<Derived>& QuaternionBase<Derived>::operator=(const MatrixBase<MatrixDerived>& xpr) +inline Derived& QuaternionBase<Derived>::operator=(const MatrixBase<MatrixDerived>& xpr) { EIGEN_STATIC_ASSERT((ei_is_same_type<typename Derived::Scalar, typename MatrixDerived::Scalar>::ret), YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) ei_quaternionbase_assign_impl<MatrixDerived>::run(*this, xpr.derived()); - return *this; + return derived(); } /** Convert the quaternion to a 3x3 rotation matrix. The quaternion is required to @@ -434,7 +471,7 @@ QuaternionBase<Derived>::toRotationMatrix(void) const */ template<class Derived> template<typename Derived1, typename Derived2> -inline QuaternionBase<Derived>& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b) +inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Derived1>& a, const MatrixBase<Derived2>& b) { Vector3 v0 = a.normalized(); Vector3 v1 = b.normalized(); @@ -458,7 +495,7 @@ inline QuaternionBase<Derived>& QuaternionBase<Derived>::setFromTwoVectors(const Scalar w2 = (Scalar(1)+c)*Scalar(0.5); this->w() = ei_sqrt(w2); this->vec() = axis * ei_sqrt(Scalar(1) - w2); - return *this; + return derived(); } Vector3 axis = v0.cross(v1); Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2)); @@ -466,17 +503,17 @@ inline QuaternionBase<Derived>& QuaternionBase<Derived>::setFromTwoVectors(const this->vec() = axis * invs; this->w() = s * Scalar(0.5); - return *this; + return derived(); } /** \returns the multiplicative inverse of \c *this * Note that in most cases, i.e., if you simply want the opposite rotation, * and/or the quaternion is normalized, then it is enough to use the conjugate. * - * \sa Quaternion2::conjugate() + * \sa QuaternionBase::conjugate() */ template <class Derived> -inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> QuaternionBase<Derived>::inverse() const +inline Quaternion<typename ei_traits<Derived>::Scalar> QuaternionBase<Derived>::inverse() const { // FIXME should this function be called multiplicativeInverse and conjugate() be called inverse() or opposite() ?? Scalar n2 = this->squaredNorm(); @@ -485,7 +522,7 @@ inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> Quaterni else { // return an invalid result to flag the error - return Quaternion<Scalar>(ei_traits<Derived>::Coefficients::Zero()); + return Quaternion<Scalar>(Coefficients::Zero()); } } @@ -496,7 +533,8 @@ inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> Quaterni * \sa Quaternion2::inverse() */ template <class Derived> -inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> QuaternionBase<Derived>::conjugate() const +inline Quaternion<typename ei_traits<Derived>::Scalar> +QuaternionBase<Derived>::conjugate() const { return Quaternion<Scalar>(this->w(),-this->x(),-this->y(),-this->z()); } @@ -506,11 +544,12 @@ inline Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> Quaterni */ template <class Derived> template <class OtherDerived> -inline typename ei_traits<QuaternionBase<Derived> >::Scalar QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& other) const +inline typename ei_traits<Derived>::Scalar +QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& other) const { double d = ei_abs(this->dot(other)); if (d>=1.0) - return 0; + return Scalar(0); return Scalar(2) * std::acos(d); } @@ -519,13 +558,14 @@ inline typename ei_traits<QuaternionBase<Derived> >::Scalar QuaternionBase<Deriv */ template <class Derived> template <class OtherDerived> -Quaternion<typename ei_traits<QuaternionBase<Derived> >::Scalar> QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const +Quaternion<typename ei_traits<Derived>::Scalar> +QuaternionBase<Derived>::slerp(Scalar t, const QuaternionBase<OtherDerived>& other) const { static const Scalar one = Scalar(1) - precision<Scalar>(); Scalar d = this->dot(other); Scalar absD = ei_abs(d); if (absD>=one) - return Quaternion<Scalar>(*this); + return Quaternion<Scalar>(derived()); // theta is the angle between the 2 quaternions Scalar theta = std::acos(absD); @@ -549,7 +589,7 @@ struct ei_quaternionbase_assign_impl<Other,3,3> // This algorithm comes from "Quaternion Calculus and Fast Animation", // Ken Shoemake, 1987 SIGGRAPH course notes Scalar t = mat.trace(); - if (t > 0) + if (t > Scalar(0)) { t = ei_sqrt(t + Scalar(1.0)); q.w() = Scalar(0.5)*t; diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 3b424be04..254885873 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -436,14 +436,13 @@ struct ei_solve_retval<SVD<_MatrixType>, Rhs> template<typename Dest> void evalTo(Dest& dst) const { - const int cols = this->cols(); ei_assert(rhs().rows() == dec().rows()); - for (int j=0; j<cols; ++j) + for (int j=0; j<cols(); ++j) { Matrix<Scalar,MatrixType::RowsAtCompileTime,1> aux = dec().matrixU().adjoint() * rhs().col(j); - for (int i = 0; i <dec().rows(); ++i) + for (int i = 0; i < dec().rows(); ++i) { Scalar si = dec().singularValues().coeff(i); if(si == RealScalar(0)) @@ -451,8 +450,10 @@ struct ei_solve_retval<SVD<_MatrixType>, Rhs> else aux.coeffRef(i) /= si; } - - dst.col(j) = dec().matrixV() * aux; + const int minsize = std::min(dec().rows(),dec().cols()); + dst.col(j).start(minsize) = aux.start(minsize); + if(dec().cols()>dec().rows()) dst.col(j).end(cols()-minsize).setZero(); + dst.col(j) = dec().matrixV() * dst.col(j); } } }; diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h index 30a33c3dc..f02374d7c 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/Eigen/src/Sparse/CholmodSupport.h @@ -126,6 +126,7 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType> typedef SparseLLT<MatrixType> Base; typedef typename Base::Scalar Scalar; typedef typename Base::RealScalar RealScalar; + typedef typename Base::CholMatrixType CholMatrixType; using Base::MatrixLIsDirty; using Base::SupernodalFactorIsDirty; using Base::m_flags; @@ -154,7 +155,7 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType> cholmod_finish(&m_cholmod); } - inline const typename Base::CholMatrixType& matrixL(void) const; + inline const CholMatrixType& matrixL() const; template<typename Derived> bool solveInPlace(MatrixBase<Derived> &b) const; @@ -198,7 +199,7 @@ void SparseLLT<MatrixType,Cholmod>::compute(const MatrixType& a) } template<typename MatrixType> -inline const typename SparseLLT<MatrixType>::CholMatrixType& +inline const typename SparseLLT<MatrixType,Cholmod>::CholMatrixType& SparseLLT<MatrixType,Cholmod>::matrixL() const { if (m_status & MatrixLIsDirty) |