From 17f61b27741fdbf19d2c73648a75bc37a9637d14 Mon Sep 17 00:00:00 2001 From: cignoni Date: Wed, 27 Feb 2013 21:07:14 +0000 Subject: [PATCH] Updated Eigen to the latest version. --- eigenlib/Eigen/src/Core/ArrayWrapper.h | 14 ++++++ eigenlib/Eigen/src/Core/CommaInitializer.h | 2 + eigenlib/Eigen/src/Core/DiagonalMatrix.h | 12 +++++ eigenlib/Eigen/src/Core/Functors.h | 44 +++++-------------- eigenlib/Eigen/src/Core/MatrixBase.h | 2 +- eigenlib/Eigen/src/Core/StableNorm.h | 1 - eigenlib/Eigen/src/Core/TriangularMatrix.h | 3 +- .../src/Core/products/GeneralMatrixVector.h | 22 ++++++---- .../products/TriangularMatrixMatrix_MKL.h | 20 ++++----- .../products/TriangularMatrixVector_MKL.h | 8 ++-- eigenlib/Eigen/src/Core/util/Macros.h | 2 +- eigenlib/Eigen/src/Core/util/Memory.h | 4 +- eigenlib/Eigen/src/Core/util/XprHelper.h | 4 +- eigenlib/Eigen/src/Eigenvalues/ComplexSchur.h | 8 ++-- .../Eigen/src/Eigenvalues/ComplexSchur_MKL.h | 2 +- eigenlib/Eigen/src/Eigenvalues/RealSchur.h | 10 +++-- .../src/Eigenvalues/SelfAdjointEigenSolver.h | 11 ++++- .../Eigenvalues/SelfAdjointEigenSolver_MKL.h | 2 +- eigenlib/Eigen/src/Geometry/AlignedBox.h | 2 +- .../Eigen/src/PardisoSupport/PardisoSupport.h | 5 ++- .../Eigen/src/QR/ColPivHouseholderQR_MKL.h | 2 +- eigenlib/Eigen/src/SVD/JacobiSVD_MKL.h | 2 +- .../src/SparseCore/SparseDiagonalProduct.h | 12 ++++- eigenlib/Eigen/src/SparseCore/SparseMatrix.h | 42 +++++++++++++----- eigenlib/Eigen/src/SparseCore/SparseUtil.h | 5 ++- .../Eigen/src/SuperLUSupport/SuperLUSupport.h | 6 +-- .../Eigen/src/plugins/ArrayCwiseBinaryOps.h | 6 ++- .../Eigen/src/plugins/ArrayCwiseUnaryOps.h | 1 + eigenlib/howto.txt | 2 +- .../src/MatrixFunctions/MatrixLogarithm.h | 34 +++++++------- 30 files changed, 172 insertions(+), 118 deletions(-) diff --git a/eigenlib/Eigen/src/Core/ArrayWrapper.h b/eigenlib/Eigen/src/Core/ArrayWrapper.h index 87af7fda..65ffd64c 100644 --- a/eigenlib/Eigen/src/Core/ArrayWrapper.h +++ b/eigenlib/Eigen/src/Core/ArrayWrapper.h @@ -121,6 +121,13 @@ class ArrayWrapper : public ArrayBase > return m_expression; } + /** Forwards the resizing request to the nested expression + * \sa DenseBase::resize(Index) */ + void resize(Index newSize) { m_expression.const_cast_derived().resize(newSize); } + /** Forwards the resizing request to the nested expression + * \sa DenseBase::resize(Index,Index)*/ + void resize(Index nbRows, Index nbCols) { m_expression.const_cast_derived().resize(nbRows,nbCols); } + protected: NestedExpressionType m_expression; }; @@ -231,6 +238,13 @@ class MatrixWrapper : public MatrixBase > return m_expression; } + /** Forwards the resizing request to the nested expression + * \sa DenseBase::resize(Index) */ + void resize(Index newSize) { m_expression.const_cast_derived().resize(newSize); } + /** Forwards the resizing request to the nested expression + * \sa DenseBase::resize(Index,Index)*/ + void resize(Index nbRows, Index nbCols) { m_expression.const_cast_derived().resize(nbRows,nbCols); } + protected: NestedExpressionType m_expression; }; diff --git a/eigenlib/Eigen/src/Core/CommaInitializer.h b/eigenlib/Eigen/src/Core/CommaInitializer.h index 4adce641..f20c1774 100644 --- a/eigenlib/Eigen/src/Core/CommaInitializer.h +++ b/eigenlib/Eigen/src/Core/CommaInitializer.h @@ -65,6 +65,8 @@ struct CommaInitializer template CommaInitializer& operator,(const DenseBase& other) { + if(other.cols()==0 || other.rows()==0) + return *this; if (m_col==m_xpr.cols()) { m_row+=m_currentBlockRows; diff --git a/eigenlib/Eigen/src/Core/DiagonalMatrix.h b/eigenlib/Eigen/src/Core/DiagonalMatrix.h index 88190da6..6e8b50fa 100644 --- a/eigenlib/Eigen/src/Core/DiagonalMatrix.h +++ b/eigenlib/Eigen/src/Core/DiagonalMatrix.h @@ -20,6 +20,7 @@ class DiagonalBase : public EigenBase public: typedef typename internal::traits::DiagonalVectorType DiagonalVectorType; typedef typename DiagonalVectorType::Scalar Scalar; + typedef typename DiagonalVectorType::RealScalar RealScalar; typedef typename internal::traits::StorageKind StorageKind; typedef typename internal::traits::Index Index; @@ -65,6 +66,17 @@ class DiagonalBase : public EigenBase return diagonal().cwiseInverse(); } + inline const DiagonalWrapper, const DiagonalVectorType> > + operator*(const Scalar& scalar) const + { + return diagonal() * scalar; + } + friend inline const DiagonalWrapper, const DiagonalVectorType> > + operator*(const Scalar& scalar, const DiagonalBase& other) + { + return other.diagonal() * scalar; + } + #ifdef EIGEN2_SUPPORT template bool isApprox(const DiagonalBase& other, typename NumTraits::Real precision = NumTraits::dummy_precision()) const diff --git a/eigenlib/Eigen/src/Core/Functors.h b/eigenlib/Eigen/src/Core/Functors.h index 278c46c6..2f46abfd 100644 --- a/eigenlib/Eigen/src/Core/Functors.h +++ b/eigenlib/Eigen/src/Core/Functors.h @@ -447,7 +447,7 @@ struct functor_traits > * indeed it seems better to declare m_other as a Packet and do the pset1() once * in the constructor. However, in practice: * - GCC does not like m_other as a Packet and generate a load every time it needs it - * - on the other hand GCC is able to moves the pset1() away the loop :) + * - on the other hand GCC is able to moves the pset1() outside the loop :) * - simpler code ;) * (ICC and gcc 4.4 seems to perform well in both cases, the issue is visible with y = a*x + b*y) */ @@ -478,33 +478,6 @@ template struct functor_traits > { enum { Cost = NumTraits::MulCost, PacketAccess = false }; }; -template -struct scalar_quotient1_impl { - typedef typename packet_traits::type Packet; - // FIXME default copy constructors seems bugged with std::complex<> - EIGEN_STRONG_INLINE scalar_quotient1_impl(const scalar_quotient1_impl& other) : m_other(other.m_other) { } - EIGEN_STRONG_INLINE scalar_quotient1_impl(const Scalar& other) : m_other(static_cast(1) / other) {} - EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a * m_other; } - EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const - { return internal::pmul(a, pset1(m_other)); } - const Scalar m_other; -}; -template -struct functor_traits > -{ enum { Cost = NumTraits::MulCost, PacketAccess = packet_traits::HasMul }; }; - -template -struct scalar_quotient1_impl { - // FIXME default copy constructors seems bugged with std::complex<> - EIGEN_STRONG_INLINE scalar_quotient1_impl(const scalar_quotient1_impl& other) : m_other(other.m_other) { } - EIGEN_STRONG_INLINE scalar_quotient1_impl(const Scalar& other) : m_other(other) {} - EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a / m_other; } - typename add_const_on_value_type::Nested>::type m_other; -}; -template -struct functor_traits > -{ enum { Cost = 2 * NumTraits::MulCost, PacketAccess = false }; }; - /** \internal * \brief Template functor to divide a scalar by a fixed other one * @@ -514,14 +487,19 @@ struct functor_traits > * \sa class CwiseUnaryOp, MatrixBase::operator/ */ template -struct scalar_quotient1_op : scalar_quotient1_impl::IsInteger > { - EIGEN_STRONG_INLINE scalar_quotient1_op(const Scalar& other) - : scalar_quotient1_impl::IsInteger >(other) {} +struct scalar_quotient1_op { + typedef typename packet_traits::type Packet; + // FIXME default copy constructors seems bugged with std::complex<> + EIGEN_STRONG_INLINE scalar_quotient1_op(const scalar_quotient1_op& other) : m_other(other.m_other) { } + EIGEN_STRONG_INLINE scalar_quotient1_op(const Scalar& other) : m_other(other) {} + EIGEN_STRONG_INLINE Scalar operator() (const Scalar& a) const { return a / m_other; } + EIGEN_STRONG_INLINE const Packet packetOp(const Packet& a) const + { return internal::pdiv(a, pset1(m_other)); } + typename add_const_on_value_type::Nested>::type m_other; }; template struct functor_traits > -: functor_traits::IsInteger> > -{}; +{ enum { Cost = 2 * NumTraits::MulCost, PacketAccess = packet_traits::HasDiv }; }; // nullary functors diff --git a/eigenlib/Eigen/src/Core/MatrixBase.h b/eigenlib/Eigen/src/Core/MatrixBase.h index c1e0ed13..36ea2cee 100644 --- a/eigenlib/Eigen/src/Core/MatrixBase.h +++ b/eigenlib/Eigen/src/Core/MatrixBase.h @@ -237,7 +237,7 @@ template class MatrixBase // huuuge hack. make Eigen2's matrix.part() work in eigen3. Problem: Diagonal is now a class template instead // of an integer constant. Solution: overload the part() method template wrt template parameters list. - template class U> + template class U> const DiagonalWrapper part() const { return diagonal().asDiagonal(); } #endif // EIGEN2_SUPPORT diff --git a/eigenlib/Eigen/src/Core/StableNorm.h b/eigenlib/Eigen/src/Core/StableNorm.h index d8bf7db7..7499b195 100644 --- a/eigenlib/Eigen/src/Core/StableNorm.h +++ b/eigenlib/Eigen/src/Core/StableNorm.h @@ -131,7 +131,6 @@ MatrixBase::blueNorm() const abig = internal::sqrt(abig); if(abig > overfl) { - eigen_assert(false && "overflow"); return rbig; } if(amed > RealScalar(0)) diff --git a/eigenlib/Eigen/src/Core/TriangularMatrix.h b/eigenlib/Eigen/src/Core/TriangularMatrix.h index de954006..301b0ef2 100644 --- a/eigenlib/Eigen/src/Core/TriangularMatrix.h +++ b/eigenlib/Eigen/src/Core/TriangularMatrix.h @@ -511,6 +511,7 @@ template struct triangular_assignment_selector { typedef typename Derived1::Index Index; + typedef typename Derived1::Scalar Scalar; static inline void run(Derived1 &dst, const Derived2 &src) { for(Index j = 0; j < dst.cols(); ++j) @@ -520,7 +521,7 @@ struct triangular_assignment_selector1 ? alignedStart + ((size-alignedStart) & ~ResPacketAlignedMask) : 0; - const Index peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart; + const Index peeledSize = alignedSize - RhsPacketSize*peels - RhsPacketSize + 1; const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0; Index alignmentPattern = alignmentStep==0 ? AllAligned @@ -177,6 +176,8 @@ EIGEN_DONT_INLINE static void run( _EIGEN_ACCUMULATE_PACKETS(d,du,d); break; case FirstAligned: + { + Index j = alignedStart; if(peels>1) { LhsPacket A00, A01, A02, A03, A10, A11, A12, A13; @@ -186,7 +187,7 @@ EIGEN_DONT_INLINE static void run( A02 = pload(&lhs2[alignedStart-2]); A03 = pload(&lhs3[alignedStart-3]); - for (Index j = alignedStart; j(&lhs1[j-1+LhsPacketSize]); palign<1>(A01,A11); A12 = pload(&lhs2[j-2+LhsPacketSize]); palign<2>(A02,A12); @@ -210,9 +211,10 @@ EIGEN_DONT_INLINE static void run( pstore(&res[j+ResPacketSize],T1); } } - for (Index j = peeledSize; j1 ? alignedStart + ((depth-alignedStart) & ~RhsPacketAlignedMask) : 0; - const Index peeledSize = peels>1 ? alignedStart + ((alignedSize-alignedStart) & ~PeelAlignedMask) : alignedStart; + const Index peeledSize = alignedSize - RhsPacketSize*peels - RhsPacketSize + 1; const Index alignmentStep = LhsPacketSize>1 ? (LhsPacketSize - lhsStride % LhsPacketSize) & LhsPacketAlignedMask : 0; Index alignmentPattern = alignmentStep==0 ? AllAligned @@ -430,10 +431,12 @@ EIGEN_DONT_INLINE static void run( _EIGEN_ACCUMULATE_PACKETS(d,du,d); break; case FirstAligned: + { + Index j = alignedStart; if (peels>1) { /* Here we proccess 4 rows with with two peeled iterations to hide - * tghe overhead of unaligned loads. Moreover unaligned loads are handled + * the overhead of unaligned loads. Moreover unaligned loads are handled * using special shift/move operations between the two aligned packets * overlaping the desired unaligned packet. This is *much* more efficient * than basic unaligned loads. @@ -443,7 +446,7 @@ EIGEN_DONT_INLINE static void run( A02 = pload(&lhs2[alignedStart-2]); A03 = pload(&lhs3[alignedStart-3]); - for (Index j = alignedStart; j(&rhs[j]); A11 = pload(&lhs1[j-1+LhsPacketSize]); palign<1>(A01,A11); @@ -465,9 +468,10 @@ EIGEN_DONT_INLINE static void run( ptmp3 = pcj.pmadd(A13, b, ptmp3); } } - for (Index j = peeledSize; j { \ static inline void run(Index _rows, Index _cols, Index _depth, const Scalar* _lhs, Index lhsStride,\ - const Scalar* _rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha) { \ + const Scalar* _rhs, Index rhsStride, Scalar* res, Index resStride, Scalar alpha, level3_blocking& blocking) { \ product_triangular_matrix_matrix_trmm::run( \ - _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha); \ + _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ } \ }; @@ -96,7 +96,7 @@ struct product_triangular_matrix_matrix_trmm& blocking) \ { \ Index diagSize = (std::min)(_rows,_depth); \ Index rows = IsLower ? _rows : diagSize; \ @@ -115,16 +115,16 @@ struct product_triangular_matrix_matrix_trmm::run( \ - _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha); \ + _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ /*std::cout << "TRMM_L: A is not square! Go to Eigen TRMM implementation!\n";*/ \ } else { \ /* Make sense to call GEMM */ \ Map > lhsMap(_lhs,rows,depth,OuterStride<>(lhsStride)); \ MatrixLhs aa_tmp=lhsMap.template triangularView(); \ MKL_INT aStride = aa_tmp.outerStride(); \ - gemm_blocking_space blocking(_rows,_cols,_depth); \ + gemm_blocking_space gemm_blocking(_rows,_cols,_depth); \ general_matrix_matrix_product::run( \ - rows, cols, depth, aa_tmp.data(), aStride, _rhs, rhsStride, res, resStride, alpha, blocking, 0); \ + rows, cols, depth, aa_tmp.data(), aStride, _rhs, rhsStride, res, resStride, alpha, gemm_blocking, 0); \ \ /*std::cout << "TRMM_L: A is not square! Go to MKL GEMM implementation! " << nthr<<" \n";*/ \ } \ @@ -210,7 +210,7 @@ struct product_triangular_matrix_matrix_trmm& blocking) \ { \ Index diagSize = (std::min)(_cols,_depth); \ Index rows = _rows; \ @@ -229,16 +229,16 @@ struct product_triangular_matrix_matrix_trmm::run( \ - _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha); \ + _rows, _cols, _depth, _lhs, lhsStride, _rhs, rhsStride, res, resStride, alpha, blocking); \ /*std::cout << "TRMM_R: A is not square! Go to Eigen TRMM implementation!\n";*/ \ } else { \ /* Make sense to call GEMM */ \ Map > rhsMap(_rhs,depth,cols, OuterStride<>(rhsStride)); \ MatrixRhs aa_tmp=rhsMap.template triangularView(); \ MKL_INT aStride = aa_tmp.outerStride(); \ - gemm_blocking_space blocking(_rows,_cols,_depth); \ + gemm_blocking_space gemm_blocking(_rows,_cols,_depth); \ general_matrix_matrix_product::run( \ - rows, cols, depth, _lhs, lhsStride, aa_tmp.data(), aStride, res, resStride, alpha, blocking, 0); \ + rows, cols, depth, _lhs, lhsStride, aa_tmp.data(), aStride, res, resStride, alpha, gemm_blocking, 0); \ \ /*std::cout << "TRMM_R: A is not square! Go to MKL GEMM implementation! " << nthr<<" \n";*/ \ } \ diff --git a/eigenlib/Eigen/src/Core/products/TriangularMatrixVector_MKL.h b/eigenlib/Eigen/src/Core/products/TriangularMatrixVector_MKL.h index 3589b8c5..3c2c3049 100644 --- a/eigenlib/Eigen/src/Core/products/TriangularMatrixVector_MKL.h +++ b/eigenlib/Eigen/src/Core/products/TriangularMatrixVector_MKL.h @@ -82,11 +82,11 @@ struct triangular_matrix_vector_product_trmv& blocking) \ + const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \ { \ if (ConjLhs || IsZeroDiag) { \ triangular_matrix_vector_product::run( \ - _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha, blocking); \ + _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ return; \ }\ Index size = (std::min)(_rows,_cols); \ @@ -167,11 +167,11 @@ struct triangular_matrix_vector_product_trmv& blocking) \ + const EIGTYPE* _rhs, Index rhsIncr, EIGTYPE* _res, Index resIncr, EIGTYPE alpha) \ { \ if (IsZeroDiag) { \ triangular_matrix_vector_product::run( \ - _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha, blocking); \ + _rows, _cols, _lhs, lhsStride, _rhs, rhsIncr, _res, resIncr, alpha); \ return; \ }\ Index size = (std::min)(_rows,_cols); \ diff --git a/eigenlib/Eigen/src/Core/util/Macros.h b/eigenlib/Eigen/src/Core/util/Macros.h index d973a683..31436fe7 100644 --- a/eigenlib/Eigen/src/Core/util/Macros.h +++ b/eigenlib/Eigen/src/Core/util/Macros.h @@ -13,7 +13,7 @@ #define EIGEN_WORLD_VERSION 3 #define EIGEN_MAJOR_VERSION 1 -#define EIGEN_MINOR_VERSION 1 +#define EIGEN_MINOR_VERSION 2 #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/eigenlib/Eigen/src/Core/util/Memory.h b/eigenlib/Eigen/src/Core/util/Memory.h index 6e06ace4..69148703 100644 --- a/eigenlib/Eigen/src/Core/util/Memory.h +++ b/eigenlib/Eigen/src/Core/util/Memory.h @@ -204,7 +204,7 @@ inline void* aligned_malloc(size_t size) if(posix_memalign(&result, 16, size)) result = 0; #elif EIGEN_HAS_MM_MALLOC result = _mm_malloc(size, 16); - #elif (defined _MSC_VER) +#elif defined(_MSC_VER) && (!defined(_WIN32_WCE)) result = _aligned_malloc(size, 16); #else result = handmade_aligned_malloc(size); @@ -745,7 +745,7 @@ public: __asm__ __volatile__ ("cpuid": "=a" (abcd[0]), "=b" (abcd[1]), "=c" (abcd[2]), "=d" (abcd[3]) : "a" (func), "c" (id) ); # endif # elif defined(_MSC_VER) -# if (_MSC_VER > 1500) +# if (_MSC_VER > 1500) && ( defined(_M_IX86) || defined(_M_X64) ) # define EIGEN_CPUID(abcd,func,id) __cpuidex((int*)abcd,func,id) # endif # endif diff --git a/eigenlib/Eigen/src/Core/util/XprHelper.h b/eigenlib/Eigen/src/Core/util/XprHelper.h index 2a65c7cb..e6f8aaef 100644 --- a/eigenlib/Eigen/src/Core/util/XprHelper.h +++ b/eigenlib/Eigen/src/Core/util/XprHelper.h @@ -301,9 +301,9 @@ template::type> str // it's important that this value can still be squared without integer overflowing. DynamicAsInteger = 10000, ScalarReadCost = NumTraits::Scalar>::ReadCost, - ScalarReadCostAsInteger = ScalarReadCost == Dynamic ? DynamicAsInteger : ScalarReadCost, + ScalarReadCostAsInteger = ScalarReadCost == Dynamic ? int(DynamicAsInteger) : int(ScalarReadCost), CoeffReadCost = traits::CoeffReadCost, - CoeffReadCostAsInteger = CoeffReadCost == Dynamic ? DynamicAsInteger : CoeffReadCost, + CoeffReadCostAsInteger = CoeffReadCost == Dynamic ? int(DynamicAsInteger) : int(CoeffReadCost), NAsInteger = n == Dynamic ? int(DynamicAsInteger) : n, CostEvalAsInteger = (NAsInteger+1) * ScalarReadCostAsInteger + CoeffReadCostAsInteger, CostNoEvalAsInteger = NAsInteger * CoeffReadCostAsInteger diff --git a/eigenlib/Eigen/src/Eigenvalues/ComplexSchur.h b/eigenlib/Eigen/src/Eigenvalues/ComplexSchur.h index 16a9a03d..55aeedb9 100644 --- a/eigenlib/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/eigenlib/Eigen/src/Eigenvalues/ComplexSchur.h @@ -336,6 +336,7 @@ void ComplexSchur::reduceToTriangularForm(bool computeU) Index iu = m_matT.cols() - 1; Index il; Index iter = 0; // number of iterations we are working on the (iu,iu) element + Index totalIter = 0; // number of iterations for whole matrix while(true) { @@ -350,9 +351,10 @@ void ComplexSchur::reduceToTriangularForm(bool computeU) // if iu is zero then we are done; the whole matrix is triangularized if(iu==0) break; - // if we spent too many iterations on the current element, we give up + // if we spent too many iterations, we give up iter++; - if(iter > m_maxIterations * m_matT.cols()) break; + totalIter++; + if(totalIter > m_maxIterations * m_matT.cols()) break; // find il, the top row of the active submatrix il = iu-1; @@ -382,7 +384,7 @@ void ComplexSchur::reduceToTriangularForm(bool computeU) } } - if(iter <= m_maxIterations * m_matT.cols()) + if(totalIter <= m_maxIterations * m_matT.cols()) m_info = Success; else m_info = NoConvergence; diff --git a/eigenlib/Eigen/src/Eigenvalues/ComplexSchur_MKL.h b/eigenlib/Eigen/src/Eigenvalues/ComplexSchur_MKL.h index aa18e696..ada7a24e 100644 --- a/eigenlib/Eigen/src/Eigenvalues/ComplexSchur_MKL.h +++ b/eigenlib/Eigen/src/Eigenvalues/ComplexSchur_MKL.h @@ -40,7 +40,7 @@ namespace Eigen { /** \internal Specialization for the data types supported by MKL */ #define EIGEN_MKL_SCHUR_COMPLEX(EIGTYPE, MKLTYPE, MKLPREFIX, MKLPREFIX_U, EIGCOLROW, MKLCOLROW) \ -template<> inline\ +template<> inline \ ComplexSchur >& \ ComplexSchur >::compute(const Matrix& matrix, bool computeU) \ { \ diff --git a/eigenlib/Eigen/src/Eigenvalues/RealSchur.h b/eigenlib/Eigen/src/Eigenvalues/RealSchur.h index 781692ec..d1949b83 100644 --- a/eigenlib/Eigen/src/Eigenvalues/RealSchur.h +++ b/eigenlib/Eigen/src/Eigenvalues/RealSchur.h @@ -220,8 +220,9 @@ RealSchur& RealSchur::compute(const MatrixType& matrix, // Rows il,...,iu is the part we are working on (the active window). // Rows iu+1,...,end are already brought in triangular form. Index iu = m_matT.cols() - 1; - Index iter = 0; // iteration count - Scalar exshift(0); // sum of exceptional shifts + Index iter = 0; // iteration count for current eigenvalue + Index totalIter = 0; // iteration count for whole matrix + Scalar exshift(0); // sum of exceptional shifts Scalar norm = computeNormOfT(); if(norm!=0) @@ -251,14 +252,15 @@ RealSchur& RealSchur::compute(const MatrixType& matrix, Vector3s firstHouseholderVector(0,0,0), shiftInfo; computeShift(iu, iter, exshift, shiftInfo); iter = iter + 1; - if (iter > m_maxIterations * m_matT.cols()) break; + totalIter = totalIter + 1; + if (totalIter > m_maxIterations * matrix.cols()) break; Index im; initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector); performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace); } } } - if(iter <= m_maxIterations * m_matT.cols()) + if(totalIter <= m_maxIterations * matrix.cols()) m_info = Success; else m_info = NoConvergence; diff --git a/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index acc5576f..24c78b4b 100644 --- a/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -743,7 +743,16 @@ static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index sta // RealScalar e2 = abs2(subdiag[end-1]); // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); // This explain the following, somewhat more complicated, version: - RealScalar mu = diag[end] - (e / (td + (td>0 ? 1 : -1))) * (e / hypot(td,e)); + RealScalar mu = diag[end]; + if(td==0) + mu -= abs(e); + else + { + RealScalar e2 = abs2(subdiag[end-1]); + RealScalar h = hypot(td,e); + if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h); + else mu -= e2 / (td + (td>0 ? h : -h)); + } RealScalar x = diag[start] - mu; RealScalar z = subdiag[start]; diff --git a/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h b/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h index 9380956b..5de5f15d 100644 --- a/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h +++ b/eigenlib/Eigen/src/Eigenvalues/SelfAdjointEigenSolver_MKL.h @@ -40,7 +40,7 @@ namespace Eigen { /** \internal Specialization for the data types supported by MKL */ #define EIGEN_MKL_EIG_SELFADJ(EIGTYPE, MKLTYPE, MKLRTYPE, MKLNAME, EIGCOLROW, MKLCOLROW ) \ -template<> inline\ +template<> inline \ SelfAdjointEigenSolver >& \ SelfAdjointEigenSolver >::compute(const Matrix& matrix, int options) \ { \ diff --git a/eigenlib/Eigen/src/Geometry/AlignedBox.h b/eigenlib/Eigen/src/Geometry/AlignedBox.h index 5830fcd3..c0f97300 100644 --- a/eigenlib/Eigen/src/Geometry/AlignedBox.h +++ b/eigenlib/Eigen/src/Geometry/AlignedBox.h @@ -79,7 +79,7 @@ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_AmbientDim) ~AlignedBox() {} /** \returns the dimension in which the box holds */ - inline Index dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size()-1 : Index(AmbientDimAtCompileTime); } + inline Index dim() const { return AmbientDimAtCompileTime==Dynamic ? m_min.size() : Index(AmbientDimAtCompileTime); } /** \deprecated use isEmpty */ inline bool isNull() const { return isEmpty(); } diff --git a/eigenlib/Eigen/src/PardisoSupport/PardisoSupport.h b/eigenlib/Eigen/src/PardisoSupport/PardisoSupport.h index e6defc8c..d623bf51 100644 --- a/eigenlib/Eigen/src/PardisoSupport/PardisoSupport.h +++ b/eigenlib/Eigen/src/PardisoSupport/PardisoSupport.h @@ -108,6 +108,7 @@ class PardisoImpl typedef Matrix VectorType; typedef Matrix IntRowVectorType; typedef Matrix IntColVectorType; + typedef Array ParameterType; enum { ScalarIsComplex = NumTraits::IsComplex }; @@ -142,7 +143,7 @@ class PardisoImpl /** \warning for advanced usage only. * \returns a reference to the parameter array controlling PARDISO. * See the PARDISO manual to know how to use it. */ - Array& pardisoParameterArray() + ParameterType& pardisoParameterArray() { return m_iparm; } @@ -295,7 +296,7 @@ class PardisoImpl bool m_initialized, m_analysisIsOk, m_factorizationIsOk; Index m_type, m_msglvl; mutable void *m_pt[64]; - mutable Array m_iparm; + mutable ParameterType m_iparm; mutable IntColVectorType m_perm; Index m_size; diff --git a/eigenlib/Eigen/src/QR/ColPivHouseholderQR_MKL.h b/eigenlib/Eigen/src/QR/ColPivHouseholderQR_MKL.h index 745ecf8b..ebcafe7d 100644 --- a/eigenlib/Eigen/src/QR/ColPivHouseholderQR_MKL.h +++ b/eigenlib/Eigen/src/QR/ColPivHouseholderQR_MKL.h @@ -41,7 +41,7 @@ namespace Eigen { /** \internal Specialization for the data types supported by MKL */ #define EIGEN_MKL_QR_COLPIV(EIGTYPE, MKLTYPE, MKLPREFIX, EIGCOLROW, MKLCOLROW) \ -template<> inline\ +template<> inline \ ColPivHouseholderQR >& \ ColPivHouseholderQR >::compute( \ const Matrix& matrix) \ diff --git a/eigenlib/Eigen/src/SVD/JacobiSVD_MKL.h b/eigenlib/Eigen/src/SVD/JacobiSVD_MKL.h index 4d479f6b..decda754 100644 --- a/eigenlib/Eigen/src/SVD/JacobiSVD_MKL.h +++ b/eigenlib/Eigen/src/SVD/JacobiSVD_MKL.h @@ -40,7 +40,7 @@ namespace Eigen { /** \internal Specialization for the data types supported by MKL */ #define EIGEN_MKL_SVD(EIGTYPE, MKLTYPE, MKLRTYPE, MKLPREFIX, EIGCOLROW, MKLCOLROW) \ -template<> inline\ +template<> inline \ JacobiSVD, ColPivHouseholderQRPreconditioner>& \ JacobiSVD, ColPivHouseholderQRPreconditioner>::compute(const Matrix& matrix, unsigned int computationOptions) \ { \ diff --git a/eigenlib/Eigen/src/SparseCore/SparseDiagonalProduct.h b/eigenlib/Eigen/src/SparseCore/SparseDiagonalProduct.h index 095bf686..ccba0212 100644 --- a/eigenlib/Eigen/src/SparseCore/SparseDiagonalProduct.h +++ b/eigenlib/Eigen/src/SparseCore/SparseDiagonalProduct.h @@ -126,11 +126,15 @@ class sparse_diagonal_product_inner_iterator_selector SparseInnerVectorSet, typename Lhs::DiagonalVectorType>::InnerIterator Base; typedef typename Lhs::Index Index; + Index m_outer; public: inline sparse_diagonal_product_inner_iterator_selector( const SparseDiagonalProductType& expr, Index outer) - : Base(expr.rhs().innerVector(outer) .cwiseProduct(expr.lhs().diagonal()), 0) + : Base(expr.rhs().innerVector(outer) .cwiseProduct(expr.lhs().diagonal()), 0), m_outer(outer) {} + + inline Index outer() const { return m_outer; } + inline Index col() const { return m_outer; } }; template @@ -160,11 +164,15 @@ class sparse_diagonal_product_inner_iterator_selector SparseInnerVectorSet, Transpose >::InnerIterator Base; typedef typename Lhs::Index Index; + Index m_outer; public: inline sparse_diagonal_product_inner_iterator_selector( const SparseDiagonalProductType& expr, Index outer) - : Base(expr.lhs().innerVector(outer) .cwiseProduct(expr.rhs().diagonal().transpose()), 0) + : Base(expr.lhs().innerVector(outer) .cwiseProduct(expr.rhs().diagonal().transpose()), 0), m_outer(outer) {} + + inline Index outer() const { return m_outer; } + inline Index row() const { return m_outer; } }; } // end namespace internal diff --git a/eigenlib/Eigen/src/SparseCore/SparseMatrix.h b/eigenlib/Eigen/src/SparseCore/SparseMatrix.h index efb774f0..fc3749b5 100644 --- a/eigenlib/Eigen/src/SparseCore/SparseMatrix.h +++ b/eigenlib/Eigen/src/SparseCore/SparseMatrix.h @@ -572,6 +572,16 @@ class SparseMatrix *this = other.derived(); } + /** \brief Copy constructor with in-place evaluation */ + template + SparseMatrix(const ReturnByValue& other) + : Base(), m_outerSize(0), m_innerSize(0), m_outerIndex(0), m_innerNonZeros(0) + { + check_template_parameters(); + initAssignment(other); + other.evalTo(*this); + } + /** Swaps the content of two sparse matrices of the same type. * This is a fast operation that simply swaps the underlying pointers and parameters. */ inline void swap(SparseMatrix& other) @@ -613,7 +623,10 @@ class SparseMatrix template inline SparseMatrix& operator=(const ReturnByValue& other) - { return Base::operator=(other.derived()); } + { + initAssignment(other); + return Base::operator=(other.derived()); + } template inline SparseMatrix& operator=(const EigenBase& other) @@ -623,7 +636,6 @@ class SparseMatrix template EIGEN_DONT_INLINE SparseMatrix& operator=(const SparseMatrixBase& other) { - initAssignment(other.derived()); const bool needToTranspose = (Flags & RowMajorBit) != (OtherDerived::Flags & RowMajorBit); if (needToTranspose) { @@ -635,40 +647,45 @@ class SparseMatrix typedef typename internal::remove_all::type _OtherCopy; OtherCopy otherCopy(other.derived()); - Eigen::Map > (m_outerIndex,outerSize()).setZero(); + SparseMatrix dest(other.rows(),other.cols()); + Eigen::Map > (dest.m_outerIndex,dest.outerSize()).setZero(); + // pass 1 // FIXME the above copy could be merged with that pass for (Index j=0; jswap(dest); return *this; } else { + if(other.isRValue()) + initAssignment(other.derived()); // there is no special optimization return Base::operator=(other.derived()); } @@ -888,6 +905,7 @@ protected: m_data.value(p) = m_data.value(p-1); --p; } + eigen_assert((p<=startId || m_data.index(p-1)!=inner) && "you cannot insert an element that already exist, you must call coeffRef to this end"); m_innerNonZeros[outer]++; diff --git a/eigenlib/Eigen/src/SparseCore/SparseUtil.h b/eigenlib/Eigen/src/SparseCore/SparseUtil.h index 6062a086..a686e08d 100644 --- a/eigenlib/Eigen/src/SparseCore/SparseUtil.h +++ b/eigenlib/Eigen/src/SparseCore/SparseUtil.h @@ -113,9 +113,10 @@ template struct sparse_eval { template struct sparse_eval { typedef typename traits::Scalar _Scalar; - enum { _Flags = traits::Flags }; + typedef typename traits::Index _Index; + enum { _Options = ((traits::Flags&RowMajorBit)==RowMajorBit) ? RowMajor : ColMajor }; public: - typedef SparseMatrix<_Scalar, _Flags> type; + typedef SparseMatrix<_Scalar, _Options, _Index> type; }; template struct sparse_eval { diff --git a/eigenlib/Eigen/src/SuperLUSupport/SuperLUSupport.h b/eigenlib/Eigen/src/SuperLUSupport/SuperLUSupport.h index 11fb014d..d8a54e18 100644 --- a/eigenlib/Eigen/src/SuperLUSupport/SuperLUSupport.h +++ b/eigenlib/Eigen/src/SuperLUSupport/SuperLUSupport.h @@ -496,8 +496,8 @@ class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> > SuperLU(const MatrixType& matrix) : Base() { - Base::init(); - compute(matrix); + init(); + Base::compute(matrix); } ~SuperLU() @@ -833,7 +833,7 @@ class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> > SuperILU(const MatrixType& matrix) : Base() { init(); - compute(matrix); + Base::compute(matrix); } ~SuperILU() diff --git a/eigenlib/Eigen/src/plugins/ArrayCwiseBinaryOps.h b/eigenlib/Eigen/src/plugins/ArrayCwiseBinaryOps.h index 5b979ebf..1e751ad6 100644 --- a/eigenlib/Eigen/src/plugins/ArrayCwiseBinaryOps.h +++ b/eigenlib/Eigen/src/plugins/ArrayCwiseBinaryOps.h @@ -33,7 +33,8 @@ EIGEN_MAKE_CWISE_BINARY_OP(min,internal::scalar_min_op) * * \sa max() */ -EIGEN_STRONG_INLINE const CwiseBinaryOp, const Derived, const ConstantReturnType> +EIGEN_STRONG_INLINE const CwiseBinaryOp, const Derived, + const CwiseNullaryOp, PlainObject> > (min)(const Scalar &other) const { return (min)(Derived::PlainObject::Constant(rows(), cols(), other)); @@ -52,7 +53,8 @@ EIGEN_MAKE_CWISE_BINARY_OP(max,internal::scalar_max_op) * * \sa min() */ -EIGEN_STRONG_INLINE const CwiseBinaryOp, const Derived, const ConstantReturnType> +EIGEN_STRONG_INLINE const CwiseBinaryOp, const Derived, + const CwiseNullaryOp, PlainObject> > (max)(const Scalar &other) const { return (max)(Derived::PlainObject::Constant(rows(), cols(), other)); diff --git a/eigenlib/Eigen/src/plugins/ArrayCwiseUnaryOps.h b/eigenlib/Eigen/src/plugins/ArrayCwiseUnaryOps.h index 0dffaf41..a5963679 100644 --- a/eigenlib/Eigen/src/plugins/ArrayCwiseUnaryOps.h +++ b/eigenlib/Eigen/src/plugins/ArrayCwiseUnaryOps.h @@ -200,3 +200,4 @@ EIGEN_MAKE_SCALAR_CWISE_UNARY_OP(operator<=, std::less_equal) EIGEN_MAKE_SCALAR_CWISE_UNARY_OP(operator>, std::greater) EIGEN_MAKE_SCALAR_CWISE_UNARY_OP(operator>=, std::greater_equal) + diff --git a/eigenlib/howto.txt b/eigenlib/howto.txt index 4daadb8e..493e4bcb 100644 --- a/eigenlib/howto.txt +++ b/eigenlib/howto.txt @@ -1,4 +1,4 @@ -Current Eigen Version 3.1.1 (05.Oct.2012) +Current Eigen Version 3.1.2 (05.11.2012) To update the lib: - download Eigen - unzip it somewhere diff --git a/eigenlib/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h b/eigenlib/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h index 3a50514b..892d0c9a 100644 --- a/eigenlib/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h +++ b/eigenlib/unsupported/Eigen/src/MatrixFunctions/MatrixLogarithm.h @@ -195,24 +195,24 @@ template int MatrixLogarithmAtomic::getPadeDegree(long double normTminusI) { #if LDBL_MANT_DIG == 53 // double precision - const double maxNormForPade[] = { 1.6206284795015624e-2 /* degree = 3 */ , 5.3873532631381171e-2, - 1.1352802267628681e-1, 1.8662860613541288e-1, 2.642960831111435e-1 }; + const long double maxNormForPade[] = { 1.6206284795015624e-2L /* degree = 3 */ , 5.3873532631381171e-2L, + 1.1352802267628681e-1L, 1.8662860613541288e-1L, 2.642960831111435e-1L }; #elif LDBL_MANT_DIG <= 64 // extended precision - const double maxNormForPade[] = { 5.48256690357782863103e-3 /* degree = 3 */, 2.34559162387971167321e-2, - 5.84603923897347449857e-2, 1.08486423756725170223e-1, 1.68385767881294446649e-1, - 2.32777776523703892094e-1 }; + const long double maxNormForPade[] = { 5.48256690357782863103e-3L /* degree = 3 */, 2.34559162387971167321e-2L, + 5.84603923897347449857e-2L, 1.08486423756725170223e-1L, 1.68385767881294446649e-1L, + 2.32777776523703892094e-1L }; #elif LDBL_MANT_DIG <= 106 // double-double - const double maxNormForPade[] = { 8.58970550342939562202529664318890e-5 /* degree = 3 */, - 9.34074328446359654039446552677759e-4, 4.26117194647672175773064114582860e-3, - 1.21546224740281848743149666560464e-2, 2.61100544998339436713088248557444e-2, - 4.66170074627052749243018566390567e-2, 7.32585144444135027565872014932387e-2, - 1.05026503471351080481093652651105e-1 }; + const long double maxNormForPade[] = { 8.58970550342939562202529664318890e-5L /* degree = 3 */, + 9.34074328446359654039446552677759e-4L, 4.26117194647672175773064114582860e-3L, + 1.21546224740281848743149666560464e-2L, 2.61100544998339436713088248557444e-2L, + 4.66170074627052749243018566390567e-2L, 7.32585144444135027565872014932387e-2L, + 1.05026503471351080481093652651105e-1L }; #else // quadruple precision - const double maxNormForPade[] = { 4.7419931187193005048501568167858103e-5 /* degree = 3 */, - 5.8853168473544560470387769480192666e-4, 2.9216120366601315391789493628113520e-3, - 8.8415758124319434347116734705174308e-3, 1.9850836029449446668518049562565291e-2, - 3.6688019729653446926585242192447447e-2, 5.9290962294020186998954055264528393e-2, - 8.6998436081634343903250580992127677e-2, 1.1880960220216759245467951592883642e-1 }; + const long double maxNormForPade[] = { 4.7419931187193005048501568167858103e-5L /* degree = 3 */, + 5.8853168473544560470387769480192666e-4L, 2.9216120366601315391789493628113520e-3L, + 8.8415758124319434347116734705174308e-3L, 1.9850836029449446668518049562565291e-2L, + 3.6688019729653446926585242192447447e-2L, 5.9290962294020186998954055264528393e-2L, + 8.6998436081634343903250580992127677e-2L, 1.1880960220216759245467951592883642e-1L }; #endif for (int degree = 3; degree <= maxPadeDegree; ++degree) if (normTminusI <= maxNormForPade[degree - minPadeDegree]) @@ -423,8 +423,8 @@ void MatrixLogarithmAtomic::computePade11(MatrixType& result, const * This class holds the argument to the matrix function until it is * assigned or evaluated for some other reason (so the argument * should not be changed in the meantime). It is the return type of - * matrixBase::matrixLogarithm() and most of the time this is the - * only way it is used. + * matrixBase::log() and most of the time this is the only way it + * is used. */ template class MatrixLogarithmReturnValue : public ReturnByValue >