diff options
34 files changed, 235 insertions, 99 deletions
diff --git a/CTestConfig.cmake b/CTestConfig.cmake index 0557c491a..4c0027824 100644 --- a/CTestConfig.cmake +++ b/CTestConfig.cmake @@ -4,10 +4,14 @@ ## # The following are required to uses Dart and the Cdash dashboard ## ENABLE_TESTING() ## INCLUDE(CTest) -set(CTEST_PROJECT_NAME "Eigen3.2") +set(CTEST_PROJECT_NAME "Eigen") set(CTEST_NIGHTLY_START_TIME "00:00:00 UTC") set(CTEST_DROP_METHOD "http") set(CTEST_DROP_SITE "manao.inria.fr") -set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen3.2") +set(CTEST_DROP_LOCATION "/CDash/submit.php?project=Eigen") set(CTEST_DROP_SITE_CDASH TRUE) +set(CTEST_PROJECT_SUBPROJECTS +Official +Unsupported +) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index c5ae2c87e..d66110821 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -309,13 +309,6 @@ template<> struct ldlt_inplace<Lower> cutoff = abs(NumTraits<Scalar>::epsilon() * biggest_in_corner); } - // Finish early if the matrix is not full rank. - if(biggest_in_corner < cutoff) - { - for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i; - break; - } - transpositions.coeffRef(k) = index_of_biggest_in_corner; if(k != index_of_biggest_in_corner) { @@ -351,6 +344,7 @@ template<> struct ldlt_inplace<Lower> if(rs>0) A21.noalias() -= A20 * temp.head(k); } + if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) A21 /= mat.coeffRef(k,k); diff --git a/Eigen/src/Core/ArrayWrapper.h b/Eigen/src/Core/ArrayWrapper.h index 21830745b..4bb648024 100644 --- a/Eigen/src/Core/ArrayWrapper.h +++ b/Eigen/src/Core/ArrayWrapper.h @@ -49,7 +49,7 @@ class ArrayWrapper : public ArrayBase<ArrayWrapper<ExpressionType> > typedef typename internal::nested<ExpressionType>::type NestedExpressionType; EIGEN_DEVICE_FUNC - inline ArrayWrapper(ExpressionType& matrix) : m_expression(matrix) {} + EIGEN_STRONG_INLINE ArrayWrapper(ExpressionType& matrix) : m_expression(matrix) {} EIGEN_DEVICE_FUNC inline Index rows() const { return m_expression.rows(); } diff --git a/Eigen/src/Core/CommaInitializer.h b/Eigen/src/Core/CommaInitializer.h index 2bbf74b05..70cbfeff5 100644 --- a/Eigen/src/Core/CommaInitializer.h +++ b/Eigen/src/Core/CommaInitializer.h @@ -45,6 +45,18 @@ struct CommaInitializer m_xpr.block(0, 0, other.rows(), other.cols()) = other; } + /* Copy/Move constructor which transfers ownership. This is crucial in + * absence of return value optimization to avoid assertions during destruction. */ + // FIXME in C++11 mode this could be replaced by a proper RValue constructor + EIGEN_DEVICE_FUNC + inline CommaInitializer(const CommaInitializer& o) + : m_xpr(o.m_xpr), m_row(o.m_row), m_col(o.m_col), m_currentBlockRows(o.m_currentBlockRows) { + // Mark original object as finished. In absence of R-value references we need to const_cast: + const_cast<CommaInitializer&>(o).m_row = m_xpr.rows(); + const_cast<CommaInitializer&>(o).m_col = m_xpr.cols(); + const_cast<CommaInitializer&>(o).m_currentBlockRows = 0; + } + /* inserts a scalar value in the target matrix */ EIGEN_DEVICE_FUNC CommaInitializer& operator,(const Scalar& s) @@ -110,7 +122,7 @@ struct CommaInitializer EIGEN_DEVICE_FUNC inline XprType& finished() { return m_xpr; } - XprType& m_xpr; // target expression + XprType& m_xpr; // target expression Index m_row; // current row id Index m_col; // current col id Index m_currentBlockRows; // current block height diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 1f9f295c9..b95b3e6d8 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -378,8 +378,6 @@ template<typename Derived> class MatrixBase Scalar trace() const; -/////////// Array module /////////// - template<int p> EIGEN_DEVICE_FUNC RealScalar lpNorm() const; EIGEN_DEVICE_FUNC MatrixBase<Derived>& matrix() { return *this; } @@ -387,8 +385,10 @@ template<typename Derived> class MatrixBase /** \returns an \link Eigen::ArrayBase Array \endlink expression of this matrix * \sa ArrayBase::matrix() */ - EIGEN_DEVICE_FUNC ArrayWrapper<Derived> array() { return derived(); } - EIGEN_DEVICE_FUNC const ArrayWrapper<const Derived> array() const { return derived(); } + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ArrayWrapper<Derived> array() { return derived(); } + /** \returns a const \link Eigen::ArrayBase Array \endlink expression of this matrix + * \sa ArrayBase::matrix() */ + EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const ArrayWrapper<const Derived> array() const { return derived(); } /////////// LU module /////////// diff --git a/Eigen/src/Core/Transpose.h b/Eigen/src/Core/Transpose.h index 11b0e45a8..079f21fd9 100644 --- a/Eigen/src/Core/Transpose.h +++ b/Eigen/src/Core/Transpose.h @@ -300,7 +300,8 @@ struct inplace_transpose_selector<MatrixType,false> { // non square matrix * Notice however that this method is only useful if you want to replace a matrix by its own transpose. * If you just need the transpose of a matrix, use transpose(). * - * \note if the matrix is not square, then \c *this must be a resizable matrix. + * \note if the matrix is not square, then \c *this must be a resizable matrix. + * This excludes (non-square) fixed-size matrices, block-expressions and maps. * * \sa transpose(), adjoint(), adjointInPlace() */ template<typename Derived> @@ -331,6 +332,7 @@ inline void DenseBase<Derived>::transposeInPlace() * If you just need the adjoint of a matrix, use adjoint(). * * \note if the matrix is not square, then \c *this must be a resizable matrix. + * This excludes (non-square) fixed-size matrices, block-expressions and maps. * * \sa transpose(), adjoint(), transposeInPlace() */ template<typename Derived> diff --git a/Eigen/src/Core/products/GeneralBlockPanelKernel.h b/Eigen/src/Core/products/GeneralBlockPanelKernel.h index 780fa74d3..686ff84f1 100644 --- a/Eigen/src/Core/products/GeneralBlockPanelKernel.h +++ b/Eigen/src/Core/products/GeneralBlockPanelKernel.h @@ -1128,6 +1128,8 @@ EIGEN_DONT_INLINE void gemm_pack_lhs<Scalar, Index, Pack1, Pack2, StorageOrder, enum { PacketSize = packet_traits<Scalar>::size }; EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK LHS"); + EIGEN_UNUSED_VARIABLE(stride); + EIGEN_UNUSED_VARIABLE(offset); eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride)); eigen_assert( (StorageOrder==RowMajor) || ((Pack1%PacketSize)==0 && Pack1<=4*PacketSize) ); conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj; @@ -1215,6 +1217,8 @@ EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, nr, ColMajor, Conjugate, Pan ::operator()(Scalar* blockB, const Scalar* rhs, Index rhsStride, Index depth, Index cols, Index stride, Index offset) { EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS COLMAJOR"); + EIGEN_UNUSED_VARIABLE(stride); + EIGEN_UNUSED_VARIABLE(offset); eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride)); conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj; Index packet_cols = (cols/nr) * nr; @@ -1257,6 +1261,7 @@ EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, nr, ColMajor, Conjugate, Pan template<typename Scalar, typename Index, int nr, bool Conjugate, bool PanelMode> struct gemm_pack_rhs<Scalar, Index, nr, RowMajor, Conjugate, PanelMode> { + typedef typename packet_traits<Scalar>::type Packet; enum { PacketSize = packet_traits<Scalar>::size }; EIGEN_DONT_INLINE void operator()(Scalar* blockB, const Scalar* rhs, Index rhsStride, Index depth, Index cols, Index stride=0, Index offset=0); }; @@ -1266,6 +1271,8 @@ EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, nr, RowMajor, Conjugate, Pan ::operator()(Scalar* blockB, const Scalar* rhs, Index rhsStride, Index depth, Index cols, Index stride, Index offset) { EIGEN_ASM_COMMENT("EIGEN PRODUCT PACK RHS ROWMAJOR"); + EIGEN_UNUSED_VARIABLE(stride); + EIGEN_UNUSED_VARIABLE(offset); eigen_assert(((!PanelMode) && stride==0 && offset==0) || (PanelMode && stride>=depth && offset<=stride)); conj_if<NumTraits<Scalar>::IsComplex && Conjugate> cj; Index packet_cols = (cols/nr) * nr; @@ -1276,12 +1283,18 @@ EIGEN_DONT_INLINE void gemm_pack_rhs<Scalar, Index, nr, RowMajor, Conjugate, Pan if(PanelMode) count += nr * offset; for(Index k=0; k<depth; k++) { - const Scalar* b0 = &rhs[k*rhsStride + j2]; - blockB[count+0] = cj(b0[0]); - blockB[count+1] = cj(b0[1]); - if(nr==4) blockB[count+2] = cj(b0[2]); - if(nr==4) blockB[count+3] = cj(b0[3]); - count += nr; + if (nr == PacketSize) { + Packet A = ploadu<Packet>(&rhs[k*rhsStride + j2]); + pstoreu(blockB+count, cj.pconj(A)); + count += PacketSize; + } else { + const Scalar* b0 = &rhs[k*rhsStride + j2]; + blockB[count+0] = cj(b0[0]); + blockB[count+1] = cj(b0[1]); + if(nr==4) blockB[count+2] = cj(b0[2]); + if(nr==4) blockB[count+3] = cj(b0[3]); + count += nr; + } } // skip what we have after if(PanelMode) count += nr * (stride-offset-depth); diff --git a/Eigen/src/Core/products/GeneralMatrixVector.h b/Eigen/src/Core/products/GeneralMatrixVector.h index c8697c913..a73ce5ff0 100644 --- a/Eigen/src/Core/products/GeneralMatrixVector.h +++ b/Eigen/src/Core/products/GeneralMatrixVector.h @@ -80,11 +80,8 @@ EIGEN_DONT_INLINE static void run( Index rows, Index cols, const LhsScalar* lhs, Index lhsStride, const RhsScalar* rhs, Index rhsIncr, - ResScalar* res, Index - #ifdef EIGEN_INTERNAL_DEBUGGING - resIncr - #endif - , RhsScalar alpha); + ResScalar* res, Index resIncr, + RhsScalar alpha); }; template<typename Index, typename LhsScalar, bool ConjugateLhs, typename RhsScalar, bool ConjugateRhs, int Version> @@ -92,12 +89,10 @@ EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,ColMajor,Co Index rows, Index cols, const LhsScalar* lhs, Index lhsStride, const RhsScalar* rhs, Index rhsIncr, - ResScalar* res, Index - #ifdef EIGEN_INTERNAL_DEBUGGING - resIncr - #endif - , RhsScalar alpha) + ResScalar* res, Index resIncr, + RhsScalar alpha) { + EIGEN_UNUSED_VARIABLE(resIncr); eigen_internal_assert(resIncr==1); #ifdef _EIGEN_ACCUMULATE_PACKETS #error _EIGEN_ACCUMULATE_PACKETS has already been defined @@ -350,7 +345,7 @@ EIGEN_DONT_INLINE static void run( Index rows, Index cols, const LhsScalar* lhs, Index lhsStride, const RhsScalar* rhs, Index rhsIncr, - ResScalar* res, Index resIncr, + ResScalar* res, Index resIncr, ResScalar alpha); }; @@ -364,6 +359,7 @@ EIGEN_DONT_INLINE void general_matrix_vector_product<Index,LhsScalar,RowMajor,Co { EIGEN_UNUSED_VARIABLE(rhsIncr); eigen_internal_assert(rhsIncr==1); + #ifdef _EIGEN_ACCUMULATE_PACKETS #error _EIGEN_ACCUMULATE_PACKETS has already been defined #endif diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h index ddc07d535..20b5f1627 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -113,9 +113,9 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_product<Scalar,Index,StorageOrd for (size_t i=starti; i<alignedStart; ++i) { - res[i] += t0 * A0[i] + t1 * A1[i]; - t2 += numext::conj(A0[i]) * rhs[i]; - t3 += numext::conj(A1[i]) * rhs[i]; + res[i] += cj0.pmul(A0[i], t0) + cj0.pmul(A1[i],t1); + t2 += cj1.pmul(A0[i], rhs[i]); + t3 += cj1.pmul(A1[i], rhs[i]); } // Yes this an optimization for gcc 4.3 and 4.4 (=> huge speed up) // gcc 4.2 does this optimization automatically. diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 3a1dec129..6aa716698 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -252,7 +252,12 @@ #endif // Suppresses 'unused variable' warnings. -#define EIGEN_UNUSED_VARIABLE(var) (void)var; +namespace Eigen { + namespace internal { + template<typename T> void ignore_unused_variable(const T&) {} + } +} +#define EIGEN_UNUSED_VARIABLE(var) Eigen::internal::ignore_unused_variable(var); #if !defined(EIGEN_ASM_COMMENT) #if (defined __GNUC__) && ( defined(__i386__) || defined(__x86_64__) ) diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index bc472d720..286e963d2 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -274,12 +274,12 @@ inline void* aligned_realloc(void *ptr, size_t new_size, size_t old_size) // The defined(_mm_free) is just here to verify that this MSVC version // implements _mm_malloc/_mm_free based on the corresponding _aligned_ // functions. This may not always be the case and we just try to be safe. - #if defined(_MSC_VER) && defined(_mm_free) + #if defined(_MSC_VER) && (!defined(_WIN32_WCE)) && defined(_mm_free) result = _aligned_realloc(ptr,new_size,16); #else result = generic_aligned_realloc(ptr,new_size,old_size); #endif -#elif defined(_MSC_VER) +#elif defined(_MSC_VER) && (!defined(_WIN32_WCE)) result = _aligned_realloc(ptr,new_size,16); #else result = handmade_aligned_realloc(ptr,new_size,old_size); @@ -464,7 +464,7 @@ template<typename T, bool Align> inline void conditional_aligned_delete_auto(T * * There is also the variant first_aligned(const MatrixBase&) defined in DenseCoeffsBase.h. */ template<typename Scalar, typename Index> -static inline Index first_aligned(const Scalar* array, Index size) +inline Index first_aligned(const Scalar* array, Index size) { enum { PacketSize = packet_traits<Scalar>::size, PacketAlignedMask = PacketSize-1 @@ -492,7 +492,7 @@ static inline Index first_aligned(const Scalar* array, Index size) /** \internal Returns the smallest integer multiple of \a base and greater or equal to \a size */ template<typename Index> -inline static Index first_multiple(Index size, Index base) +inline Index first_multiple(Index size, Index base) { return ((size+base-1)/base)*base; } diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 803001272..a712692e5 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -34,8 +34,9 @@ struct quaternionbase_assign_impl; template<class Derived> class QuaternionBase : public RotationBase<Derived, 3> { + public: typedef RotationBase<Derived, 3> Base; -public: + using Base::operator*; using Base::derived; @@ -203,6 +204,8 @@ public: * \li \c Quaternionf for \c float * \li \c Quaterniond for \c double * + * \warning Operations interpreting the quaternion as rotation have undefined behavior if the quaternion is not normalized. + * * \sa class AngleAxis, class Transform */ @@ -223,10 +226,10 @@ struct traits<Quaternion<_Scalar,_Options> > template<typename _Scalar, int _Options> class Quaternion : public QuaternionBase<Quaternion<_Scalar,_Options> > { +public: typedef QuaternionBase<Quaternion<_Scalar,_Options> > Base; enum { IsAligned = internal::traits<Quaternion>::IsAligned }; -public: typedef _Scalar Scalar; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Quaternion) @@ -334,9 +337,9 @@ template<typename _Scalar, int _Options> class Map<const Quaternion<_Scalar>, _Options > : public QuaternionBase<Map<const Quaternion<_Scalar>, _Options> > { + public: typedef QuaternionBase<Map<const Quaternion<_Scalar>, _Options> > Base; - public: typedef _Scalar Scalar; typedef typename internal::traits<Map>::Coefficients Coefficients; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Map) @@ -344,7 +347,7 @@ class Map<const Quaternion<_Scalar>, _Options > /** 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: + * The pointer \a coeffs must reference the four coefficients of Quaternion in the following order: * \code *coeffs == {x, y, z, w} \endcode * * If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */ @@ -371,9 +374,9 @@ template<typename _Scalar, int _Options> class Map<Quaternion<_Scalar>, _Options > : public QuaternionBase<Map<Quaternion<_Scalar>, _Options> > { + public: typedef QuaternionBase<Map<Quaternion<_Scalar>, _Options> > Base; - public: typedef _Scalar Scalar; typedef typename internal::traits<Map>::Coefficients Coefficients; EIGEN_INHERIT_ASSIGNMENT_EQUAL_OPERATOR(Map) @@ -464,7 +467,7 @@ QuaternionBase<Derived>::_transformVector(Vector3 v) const // Note that this algorithm comes from the optimization by hand // of the conversion to a Matrix followed by a Matrix/Vector product. // It appears to be much faster than the common algorithm found - // in the litterature (30 versus 39 flops). It also requires two + // in the literature (30 versus 39 flops). It also requires two // Vector3 as temporaries. Vector3 uv = this->vec().cross(v); uv += uv; @@ -667,10 +670,10 @@ QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& oth { using std::acos; using std::abs; - double d = abs(this->dot(other)); - if (d>=1.0) + Scalar d = abs(this->dot(other)); + if (d>=Scalar(1)) return Scalar(0); - return static_cast<Scalar>(2 * acos(d)); + return Scalar(2) * acos(d); } diff --git a/Eigen/src/Geometry/Scaling.h b/Eigen/src/Geometry/Scaling.h index 1c25f36fe..023fba2ee 100644 --- a/Eigen/src/Geometry/Scaling.h +++ b/Eigen/src/Geometry/Scaling.h @@ -62,10 +62,10 @@ public: template<int Dim, int Mode, int Options> inline Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode)> operator* (const Transform<Scalar,Dim, Mode, Options>& t) const { - Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode)> res = t; - res.prescale(factor()); - return res; -} + Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode)> res = t; + res.prescale(factor()); + return res; + } /** Concatenates a uniform scaling and a linear transformation matrix */ // TODO returns an expression diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 0a5012cfe..f40644011 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -530,9 +530,9 @@ public: inline Transform& operator=(const UniformScaling<Scalar>& t); inline Transform& operator*=(const UniformScaling<Scalar>& s) { return scale(s.factor()); } - inline Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode)> operator*(const UniformScaling<Scalar>& s) const + inline TransformTimeDiagonalReturnType operator*(const UniformScaling<Scalar>& s) const { - Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode),Options> res = *this; + TransformTimeDiagonalReturnType res = *this; res.scale(s.factor()); return res; } diff --git a/Eigen/src/PaStiXSupport/PaStiXSupport.h b/Eigen/src/PaStiXSupport/PaStiXSupport.h index a955287d1..8a546dc2f 100644 --- a/Eigen/src/PaStiXSupport/PaStiXSupport.h +++ b/Eigen/src/PaStiXSupport/PaStiXSupport.h @@ -12,6 +12,14 @@ namespace Eigen { +#if defined(DCOMPLEX) + #define PASTIX_COMPLEX COMPLEX + #define PASTIX_DCOMPLEX DCOMPLEX +#else + #define PASTIX_COMPLEX std::complex<float> + #define PASTIX_DCOMPLEX std::complex<double> +#endif + /** \ingroup PaStiXSupport_Module * \brief Interface to the PaStix solver * @@ -74,14 +82,14 @@ namespace internal { if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } if (nbrhs == 0) {x = NULL; nbrhs=1;} - c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<COMPLEX*>(vals), perm, invp, reinterpret_cast<COMPLEX*>(x), nbrhs, iparm, dparm); + c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm); } void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<double> *vals, int *perm, int * invp, std::complex<double> *x, int nbrhs, int *iparm, double *dparm) { if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; } if (nbrhs == 0) {x = NULL; nbrhs=1;} - z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<DCOMPLEX*>(vals), perm, invp, reinterpret_cast<DCOMPLEX*>(x), nbrhs, iparm, dparm); + z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm); } // Convert the matrix to Fortran-style Numbering diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 1bf19d19e..96904c65f 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -378,7 +378,7 @@ template<typename _MatrixType> class ColPivHouseholderQR 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. - : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize(); + : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 8cdf14e4b..ea15da041 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -372,7 +372,7 @@ template<typename _MatrixType> class FullPivHouseholderQR 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. - : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize(); + : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); } /** \returns the number of nonzero pivots in the QR decomposition. @@ -449,7 +449,7 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons m_temp.resize(cols); - m_precision = NumTraits<Scalar>::epsilon() * size; + m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size); m_rows_transpositions.resize(size); m_cols_transpositions.resize(size); diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index 0e7e760fb..9ac18bcf6 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -223,7 +223,7 @@ class SparseMatrix if(isCompressed()) { - reserve(VectorXi::Constant(outerSize(), 2)); + reserve(Matrix<Index,Dynamic,1>::Constant(outerSize(), 2)); } return insertUncompressed(row,col); } @@ -939,12 +939,13 @@ void set_from_triplets(const InputIterator& begin, const InputIterator& end, Spa EIGEN_UNUSED_VARIABLE(Options); enum { IsRowMajor = SparseMatrixType::IsRowMajor }; typedef typename SparseMatrixType::Scalar Scalar; + typedef typename SparseMatrixType::Index Index; SparseMatrix<Scalar,IsRowMajor?ColMajor:RowMajor> trMat(mat.rows(),mat.cols()); if(begin!=end) { // pass 1: count the nnz per inner-vector - VectorXi wi(trMat.outerSize()); + Matrix<Index,Dynamic,1> wi(trMat.outerSize()); wi.setZero(); for(InputIterator it(begin); it!=end; ++it) { @@ -1018,7 +1019,7 @@ void SparseMatrix<Scalar,_Options,_Index>::sumupDuplicates() { eigen_assert(!isCompressed()); // TODO, in practice we should be able to use m_innerNonZeros for that task - VectorXi wi(innerSize()); + Matrix<Index,Dynamic,1> wi(innerSize()); wi.fill(-1); Index count = 0; // for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers @@ -1081,7 +1082,7 @@ EIGEN_DONT_INLINE SparseMatrix<Scalar,_Options,_Index>& SparseMatrix<Scalar,_Opt // prefix sum Index count = 0; - VectorXi positions(dest.outerSize()); + Matrix<Index,Dynamic,1> positions(dest.outerSize()); for (Index j=0; j<dest.outerSize(); ++j) { Index tmp = dest.m_outerIndex[j]; diff --git a/Eigen/src/SparseCore/SparsePermutation.h b/Eigen/src/SparseCore/SparsePermutation.h index b897b7595..b85be93f6 100644 --- a/Eigen/src/SparseCore/SparsePermutation.h +++ b/Eigen/src/SparseCore/SparsePermutation.h @@ -57,7 +57,7 @@ struct permut_sparsematrix_product_retval if(MoveOuter) { SparseMatrix<Scalar,SrcStorageOrder,Index> tmp(m_matrix.rows(), m_matrix.cols()); - VectorXi sizes(m_matrix.outerSize()); + Matrix<Index,Dynamic,1> sizes(m_matrix.outerSize()); for(Index j=0; j<m_matrix.outerSize(); ++j) { Index jp = m_permutation.indices().coeff(j); @@ -77,7 +77,7 @@ struct permut_sparsematrix_product_retval else { SparseMatrix<Scalar,int(SrcStorageOrder)==RowMajor?ColMajor:RowMajor,Index> tmp(m_matrix.rows(), m_matrix.cols()); - VectorXi sizes(tmp.outerSize()); + Matrix<Index,Dynamic,1> sizes(tmp.outerSize()); sizes.setZero(); PermutationMatrix<Dynamic,Dynamic,Index> perm; if((Side==OnTheLeft) ^ Transposed) diff --git a/bench/btl/CMakeLists.txt b/bench/btl/CMakeLists.txt index 119b470d9..0ac5c15fc 100644 --- a/bench/btl/CMakeLists.txt +++ b/bench/btl/CMakeLists.txt @@ -48,6 +48,12 @@ include_directories( # set(DEFAULT_LIBRARIES ${MKL_LIBRARIES}) # endif (MKL_FOUND) +find_library(EIGEN_BTL_RT_LIBRARY rt) +# if we cannot find it easily, then we don't need it! +if(NOT EIGEN_BTL_RT_LIBRARY) + set(EIGEN_BTL_RT_LIBRARY "") +endif() + MACRO(BTL_ADD_BENCH targetname) foreach(_current_var ${ARGN}) @@ -70,7 +76,7 @@ MACRO(BTL_ADD_BENCH targetname) IF(BUILD_${targetname}) ADD_EXECUTABLE(${targetname} ${_sources}) ADD_TEST(${targetname} "${targetname}") - target_link_libraries(${targetname} ${DEFAULT_LIBRARIES} rt) + target_link_libraries(${targetname} ${DEFAULT_LIBRARIES} ${EIGEN_BTL_RT_LIBRARY}) ENDIF(BUILD_${targetname}) ENDMACRO(BTL_ADD_BENCH) diff --git a/bench/btl/generic_bench/bench.hh b/bench/btl/generic_bench/bench.hh index 005c36395..7b7b951b5 100644 --- a/bench/btl/generic_bench/bench.hh +++ b/bench/btl/generic_bench/bench.hh @@ -102,8 +102,8 @@ BTL_DONT_INLINE void bench( int size_min, int size_max, int nb_point ) // merge the two data std::vector<int> newSizes; std::vector<double> newFlops; - int i=0; - int j=0; + unsigned int i=0; + unsigned int j=0; while (i<tab_sizes.size() && j<oldSizes.size()) { if (tab_sizes[i] == oldSizes[j]) diff --git a/bench/btl/generic_bench/btl.hh b/bench/btl/generic_bench/btl.hh index f1a88ff74..4670f5113 100644 --- a/bench/btl/generic_bench/btl.hh +++ b/bench/btl/generic_bench/btl.hh @@ -46,7 +46,7 @@ #if (defined __GNUC__) && (!defined __INTEL_COMPILER) && !defined(__arm__) && !defined(__powerpc__) #define BTL_DISABLE_SSE_EXCEPTIONS() { \ - int aux; \ + int aux = 0; \ asm( \ "stmxcsr %[aux] \n\t" \ "orl $32832, %[aux] \n\t" \ diff --git a/bench/btl/generic_bench/init/init_matrix.hh b/bench/btl/generic_bench/init/init_matrix.hh index 67cbd2073..6382d30c8 100644 --- a/bench/btl/generic_bench/init/init_matrix.hh +++ b/bench/btl/generic_bench/init/init_matrix.hh @@ -29,7 +29,7 @@ BTL_DONT_INLINE void init_row(Vector & X, int size, int row){ X.resize(size); - for (int j=0;j<X.size();j++){ + for (unsigned int j=0;j<X.size();j++){ X[j]=typename Vector::value_type(init_function(row,j)); } } @@ -42,7 +42,7 @@ BTL_DONT_INLINE void init_row(Vector & X, int size, int row){ template<double init_function(int,int),class Vector> BTL_DONT_INLINE void init_matrix(Vector & A, int size){ A.resize(size); - for (int row=0; row<A.size() ; row++){ + for (unsigned int row=0; row<A.size() ; row++){ init_row<init_function>(A[row],size,row); } } @@ -50,11 +50,11 @@ BTL_DONT_INLINE void init_matrix(Vector & A, int size){ template<double init_function(int,int),class Matrix> BTL_DONT_INLINE void init_matrix_symm(Matrix& A, int size){ A.resize(size); - for (int row=0; row<A.size() ; row++) + for (unsigned int row=0; row<A.size() ; row++) A[row].resize(size); - for (int row=0; row<A.size() ; row++){ + for (unsigned int row=0; row<A.size() ; row++){ A[row][row] = init_function(row,row); - for (int col=0; col<row ; col++){ + for (unsigned int col=0; col<row ; col++){ double x = init_function(row,col); A[row][col] = A[col][row] = x; } diff --git a/bench/btl/generic_bench/init/init_vector.hh b/bench/btl/generic_bench/init/init_vector.hh index efaf0c92e..518e87dbe 100644 --- a/bench/btl/generic_bench/init/init_vector.hh +++ b/bench/btl/generic_bench/init/init_vector.hh @@ -29,7 +29,7 @@ void init_vector(Vector & X, int size){ X.resize(size); - for (int i=0;i<X.size();i++){ + for (unsigned int i=0;i<X.size();i++){ X[i]=typename Vector::value_type(init_function(i)); } } diff --git a/bench/btl/generic_bench/timers/portable_perf_analyzer.hh b/bench/btl/generic_bench/timers/portable_perf_analyzer.hh index fc0f3168d..5e579fb49 100644 --- a/bench/btl/generic_bench/timers/portable_perf_analyzer.hh +++ b/bench/btl/generic_bench/timers/portable_perf_analyzer.hh @@ -78,7 +78,7 @@ public: // time measurement action.calculate(); _chronos.start(); - for (int ii=0;ii<_nb_calc;ii++) + for (unsigned int ii=0;ii<_nb_calc;ii++) { action.calculate(); } diff --git a/bench/btl/generic_bench/timers/portable_timer.hh b/bench/btl/generic_bench/timers/portable_timer.hh index e6ad309fe..c199811b6 100755 --- a/bench/btl/generic_bench/timers/portable_timer.hh +++ b/bench/btl/generic_bench/timers/portable_timer.hh @@ -34,7 +34,7 @@ // timer -------------------------------------------------------------------// // A timer object measures CPU time. -#ifdef _MSC_VER +#if defined(_MSC_VER) #define NOMINMAX #include <windows.h> @@ -87,6 +87,48 @@ }; // Portable_Timer +#elif defined(__APPLE__) +#include <CoreServices/CoreServices.h> +#include <mach/mach_time.h> + + +class Portable_Timer +{ + public: + + Portable_Timer() + { + } + + void start() + { + m_start_time = double(mach_absolute_time())*1e-9;; + + } + + void stop() + { + m_stop_time = double(mach_absolute_time())*1e-9;; + + } + + double elapsed() + { + return user_time(); + } + + double user_time() + { + return m_stop_time - m_start_time; + } + + +private: + + double m_stop_time, m_start_time; + +}; // Portable_Timer (Apple) + #else #include <sys/time.h> @@ -138,7 +180,7 @@ private: int m_clkid; double m_stop_time, m_start_time; -}; // Portable_Timer +}; // Portable_Timer (Linux) #endif diff --git a/bench/btl/libs/eigen3/eigen3_interface.hh b/bench/btl/libs/eigen3/eigen3_interface.hh index 31bcc1f93..0c8aa74da 100644 --- a/bench/btl/libs/eigen3/eigen3_interface.hh +++ b/bench/btl/libs/eigen3/eigen3_interface.hh @@ -52,8 +52,8 @@ public : static BTL_DONT_INLINE void matrix_from_stl(gene_matrix & A, stl_matrix & A_stl){ A.resize(A_stl[0].size(), A_stl.size()); - for (int j=0; j<A_stl.size() ; j++){ - for (int i=0; i<A_stl[j].size() ; i++){ + for (unsigned int j=0; j<A_stl.size() ; j++){ + for (unsigned int i=0; i<A_stl[j].size() ; i++){ A.coeffRef(i,j) = A_stl[j][i]; } } @@ -62,13 +62,13 @@ public : static BTL_DONT_INLINE void vector_from_stl(gene_vector & B, stl_vector & B_stl){ B.resize(B_stl.size(),1); - for (int i=0; i<B_stl.size() ; i++){ + for (unsigned int i=0; i<B_stl.size() ; i++){ B.coeffRef(i) = B_stl[i]; } } static BTL_DONT_INLINE void vector_to_stl(gene_vector & B, stl_vector & B_stl){ - for (int i=0; i<B_stl.size() ; i++){ + for (unsigned int i=0; i<B_stl.size() ; i++){ B_stl[i] = B.coeff(i); } } diff --git a/doc/TopicLinearAlgebraDecompositions.dox b/doc/TopicLinearAlgebraDecompositions.dox index 8649cc27b..77f2c92ab 100644 --- a/doc/TopicLinearAlgebraDecompositions.dox +++ b/doc/TopicLinearAlgebraDecompositions.dox @@ -249,7 +249,7 @@ For an introduction on linear solvers and decompositions, check this \link Tutor <dt><b>Implicit Multi Threading (MT)</b></dt> <dd>Means the algorithm can take advantage of multicore processors via OpenMP. "Implicit" means the algortihm itself is not parallelized, but that it relies on parallelized matrix-matrix product rountines.</dd> <dt><b>Explicit Multi Threading (MT)</b></dt> - <dd>Means the algorithm is explicitely parallelized to take advantage of multicore processors via OpenMP.</dd> + <dd>Means the algorithm is explicitly parallelized to take advantage of multicore processors via OpenMP.</dd> <dt><b>Meta-unroller</b></dt> <dd>Means the algorithm is automatically and explicitly unrolled for very small fixed size matrices.</dd> <dt><b></b></dt> diff --git a/doc/TopicMultithreading.dox b/doc/TopicMultithreading.dox index fb944af29..ba5e26290 100644 --- a/doc/TopicMultithreading.dox +++ b/doc/TopicMultithreading.dox @@ -39,7 +39,7 @@ int main(int argc, char** argv) } \endcode -\warning note that all functions generating random matrices are \b not re-entrant nor thread-safe. Those include DenseBase::Random(), and DenseBase::setRandom() despite a call to Eigen::initParallel(). This is because these functions are based on std::rand which is not re-entrant. For thread-safe random generator, we recommend the use of boost::random of c++11 random feature. +\warning note that all functions generating random matrices are \b not re-entrant nor thread-safe. Those include DenseBase::Random(), and DenseBase::setRandom() despite a call to Eigen::initParallel(). This is because these functions are based on std::rand which is not re-entrant. For thread-safe random generator, we recommend the use of boost::random or c++11 random feature. In the case your application is parallelized with OpenMP, you might want to disable Eigen's own parallization as detailed in the previous section. diff --git a/doc/snippets/EigenSolver_eigenvectors.cpp b/doc/snippets/EigenSolver_eigenvectors.cpp index 0fad4dadb..8355f76c9 100644 --- a/doc/snippets/EigenSolver_eigenvectors.cpp +++ b/doc/snippets/EigenSolver_eigenvectors.cpp @@ -1,4 +1,4 @@ MatrixXd ones = MatrixXd::Ones(3,3); EigenSolver<MatrixXd> es(ones); -cout << "The first eigenvector of the 3x3 matrix of ones is:" - << endl << es.eigenvectors().col(1) << endl; +cout << "The first eigenvector of the 3x3 matrix of ones is:" + << endl << es.eigenvectors().col(0) << endl; diff --git a/test/array.cpp b/test/array.cpp index f1deda7e3..5f49fc1ea 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -173,21 +173,14 @@ template<typename ArrayType> void array_real(const ArrayType& m) Scalar s1 = internal::random<Scalar>(); // these tests are mostly to check possible compilation issues. -// VERIFY_IS_APPROX(m1.sin(), std::sin(m1)); VERIFY_IS_APPROX(m1.sin(), sin(m1)); -// VERIFY_IS_APPROX(m1.cos(), std::cos(m1)); VERIFY_IS_APPROX(m1.cos(), cos(m1)); -// VERIFY_IS_APPROX(m1.asin(), std::asin(m1)); VERIFY_IS_APPROX(m1.asin(), asin(m1)); -// VERIFY_IS_APPROX(m1.acos(), std::acos(m1)); VERIFY_IS_APPROX(m1.acos(), acos(m1)); -// VERIFY_IS_APPROX(m1.tan(), std::tan(m1)); VERIFY_IS_APPROX(m1.tan(), tan(m1)); VERIFY_IS_APPROX(cos(m1+RealScalar(3)*m2), cos((m1+RealScalar(3)*m2).eval())); -// VERIFY_IS_APPROX(std::cos(m1+RealScalar(3)*m2), std::cos((m1+RealScalar(3)*m2).eval())); -// VERIFY_IS_APPROX(m1.abs().sqrt(), std::sqrt(std::abs(m1))); VERIFY_IS_APPROX(m1.abs().sqrt(), sqrt(abs(m1))); VERIFY_IS_APPROX(m1.abs(), sqrt(numext::abs2(m1))); @@ -196,9 +189,10 @@ template<typename ArrayType> void array_real(const ArrayType& m) if(!NumTraits<Scalar>::IsComplex) VERIFY_IS_APPROX(numext::real(m1), m1); - VERIFY_IS_APPROX(m1.abs().log() , log(abs(m1))); + // shift argument of logarithm so that it is not zero + Scalar smallNumber = NumTraits<Scalar>::dummy_precision(); + VERIFY_IS_APPROX((m1.abs() + smallNumber).log() , log(abs(m1) + smallNumber)); -// VERIFY_IS_APPROX(m1.exp(), std::exp(m1)); VERIFY_IS_APPROX(m1.exp() * m2.exp(), exp(m1+m2)); VERIFY_IS_APPROX(m1.exp(), exp(m1)); VERIFY_IS_APPROX(m1.exp() / m2.exp(),(m1-m2).exp()); @@ -242,7 +236,6 @@ template<typename ArrayType> void array_complex(const ArrayType& m) m2(i,j) = sqrt(m1(i,j)); VERIFY_IS_APPROX(m1.sqrt(), m2); -// VERIFY_IS_APPROX(m1.sqrt(), std::sqrt(m1)); VERIFY_IS_APPROX(m1.sqrt(), Eigen::sqrt(m1)); } diff --git a/test/block.cpp b/test/block.cpp index 189fa5aba..9ed5d7bc5 100644 --- a/test/block.cpp +++ b/test/block.cpp @@ -10,6 +10,26 @@ #define EIGEN_NO_STATIC_ASSERT // otherwise we fail at compile time on unused paths #include "main.h" +template<typename MatrixType, typename Index, typename Scalar> +typename Eigen::internal::enable_if<!NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type +block_real_only(const MatrixType &m1, Index r1, Index r2, Index c1, Index c2, const Scalar& s1) { + // check cwise-Functions: + VERIFY_IS_APPROX(m1.row(r1).cwiseMax(s1), m1.cwiseMax(s1).row(r1)); + VERIFY_IS_APPROX(m1.col(c1).cwiseMin(s1), m1.cwiseMin(s1).col(c1)); + + VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMin(s1), m1.cwiseMin(s1).block(r1,c1,r2-r1+1,c2-c1+1)); + VERIFY_IS_APPROX(m1.block(r1,c1,r2-r1+1,c2-c1+1).cwiseMax(s1), m1.cwiseMax(s1).block(r1,c1,r2-r1+1,c2-c1+1)); + + return Scalar(0); +} + +template<typename MatrixType, typename Index, typename Scalar> +typename Eigen::internal::enable_if<NumTraits<typename MatrixType::Scalar>::IsComplex,typename MatrixType::Scalar>::type +block_real_only(const MatrixType &, Index, Index, Index, Index, const Scalar&) { + return Scalar(0); +} + + template<typename MatrixType> void block(const MatrixType& m) { typedef typename MatrixType::Index Index; @@ -37,6 +57,8 @@ template<typename MatrixType> void block(const MatrixType& m) Index c1 = internal::random<Index>(0,cols-1); Index c2 = internal::random<Index>(c1,cols-1); + block_real_only(m1, r1, r2, c1, c1, s1); + //check row() and col() VERIFY_IS_EQUAL(m1.col(c1).transpose(), m1.transpose().row(c1)); //check operator(), both constant and non-constant, on row() and col() @@ -51,7 +73,8 @@ template<typename MatrixType> void block(const MatrixType& m) VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + s1 * m1_copy.col(c2)); m1.col(c1).col(0) += s1 * m1_copy.col(c2); VERIFY_IS_APPROX(m1.col(c1), m1_copy.col(c1) + Scalar(2) * s1 * m1_copy.col(c2)); - + + //check block() Matrix<Scalar,Dynamic,Dynamic> b1(1,1); b1(0,0) = m1(r1,c1); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index b980dc572..d4d90e467 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -179,6 +179,38 @@ template<typename MatrixType> void cholesky(const MatrixType& m) // restore if(sign == -1) symm = -symm; + + // check matrices coming from linear constraints with Lagrange multipliers + if(rows>=3) + { + SquareMatrixType A = symm; + int c = internal::random<int>(0,rows-2); + A.bottomRightCorner(c,c).setZero(); + // Make sure a solution exists: + vecX.setRandom(); + vecB = A * vecX; + vecX.setZero(); + ldltlo.compute(A); + VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); + vecX = ldltlo.solve(vecB); + VERIFY_IS_APPROX(A * vecX, vecB); + } + + // check non-full rank matrices + if(rows>=3) + { + int r = internal::random<int>(1,rows-1); + Matrix<Scalar,Dynamic,Dynamic> a = Matrix<Scalar,Dynamic,Dynamic>::Random(rows,r); + SquareMatrixType A = a * a.adjoint(); + // Make sure a solution exists: + vecX.setRandom(); + vecB = A * vecX; + vecX.setZero(); + ldltlo.compute(A); + VERIFY_IS_APPROX(A, ldltlo.reconstructedMatrix()); + vecX = ldltlo.solve(vecB); + VERIFY_IS_APPROX(A * vecX, vecB); + } } // update/downdate diff --git a/test/sizeof.cpp b/test/sizeof.cpp index d9ad35620..7044d2062 100644 --- a/test/sizeof.cpp +++ b/test/sizeof.cpp @@ -21,6 +21,8 @@ template<typename MatrixType> void verifySizeOf(const MatrixType&) void test_sizeof() { CALL_SUBTEST(verifySizeOf(Matrix<float, 1, 1>()) ); + CALL_SUBTEST(verifySizeOf(Vector2d()) ); + CALL_SUBTEST(verifySizeOf(Vector4f()) ); CALL_SUBTEST(verifySizeOf(Matrix4d()) ); CALL_SUBTEST(verifySizeOf(Matrix<double, 4, 2>()) ); CALL_SUBTEST(verifySizeOf(Matrix<bool, 7, 5>()) ); |