diff --git a/uppsrc/plugin/Eigen/Eigen/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/CMakeLists.txt deleted file mode 100644 index a92dd6f6c..000000000 --- a/uppsrc/plugin/Eigen/Eigen/CMakeLists.txt +++ /dev/null @@ -1,19 +0,0 @@ -include(RegexUtils) -test_escape_string_as_regex() - -file(GLOB Eigen_directory_files "*") - -escape_string_as_regex(ESCAPED_CMAKE_CURRENT_SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}") - -foreach(f ${Eigen_directory_files}) - if(NOT f MATCHES "\\.txt" AND NOT f MATCHES "${ESCAPED_CMAKE_CURRENT_SOURCE_DIR}/[.].+" AND NOT f MATCHES "${ESCAPED_CMAKE_CURRENT_SOURCE_DIR}/src") - list(APPEND Eigen_directory_files_to_install ${f}) - endif() -endforeach(f ${Eigen_directory_files}) - -install(FILES - ${Eigen_directory_files_to_install} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen COMPONENT Devel - ) - -add_subdirectory(src) diff --git a/uppsrc/plugin/Eigen/Eigen/Core b/uppsrc/plugin/Eigen/Eigen/Core index 9131cc3fc..c87f99df3 100644 --- a/uppsrc/plugin/Eigen/Eigen/Core +++ b/uppsrc/plugin/Eigen/Eigen/Core @@ -95,7 +95,7 @@ extern "C" { // In theory we should only include immintrin.h and not the other *mmintrin.h header files directly. // Doing so triggers some issues with ICC. However old gcc versions seems to not have this file, thus: - #ifdef __INTEL_COMPILER + #if defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1110 #include #else #include @@ -165,7 +165,7 @@ #endif // required for __cpuid, needs to be included after cmath -#if defined(_MSC_VER) && (defined(_M_IX86)||defined(_M_X64)) +#if defined(_MSC_VER) && (defined(_M_IX86)||defined(_M_X64)) && (!defined(_WIN32_WCE)) #include #endif diff --git a/uppsrc/plugin/Eigen/Eigen/src/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/CMakeLists.txt deleted file mode 100644 index c326f374d..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/CMakeLists.txt +++ /dev/null @@ -1,7 +0,0 @@ -file(GLOB Eigen_src_subdirectories "*") -escape_string_as_regex(ESCAPED_CMAKE_CURRENT_SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}") -foreach(f ${Eigen_src_subdirectories}) - if(NOT f MATCHES "\\.txt" AND NOT f MATCHES "${ESCAPED_CMAKE_CURRENT_SOURCE_DIR}/[.].+" ) - add_subdirectory(${f}) - endif() -endforeach() diff --git a/uppsrc/plugin/Eigen/Eigen/src/Cholesky/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Cholesky/CMakeLists.txt deleted file mode 100644 index d01488b41..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Cholesky/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Cholesky_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Cholesky_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Cholesky COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Cholesky/LDLT.h b/uppsrc/plugin/Eigen/Eigen/src/Cholesky/LDLT.h index d026418f8..02ab93880 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Cholesky/LDLT.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Cholesky/LDLT.h @@ -274,30 +274,13 @@ template<> struct ldlt_inplace return true; } - RealScalar cutoff(0), biggest_in_corner; - for (Index k = 0; k < size; ++k) { // Find largest diagonal element Index index_of_biggest_in_corner; - biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); + mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); index_of_biggest_in_corner += k; - if(k == 0) - { - // The biggest overall is the point of reference to which further diagonals - // are compared; if any diagonal is negligible compared - // to the largest overall, the algorithm bails. - cutoff = abs(NumTraits::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) { @@ -328,15 +311,20 @@ template<> struct ldlt_inplace if(k>0) { - temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint(); + temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint(); mat.coeffRef(k,k) -= (A10 * temp.head(k)).value(); if(rs>0) A21.noalias() -= A20 * temp.head(k); } - if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) - A21 /= mat.coeffRef(k,k); - + + // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot + // was smaller than the cutoff value. However, soince LDLT is not rank-revealing + // we should only make sure we do not introduce INF or NaN values. + // LAPACK also uses 0 as the cutoff value. RealScalar realAkk = numext::real(mat.coeffRef(k,k)); + if((rs>0) && (abs(realAkk) > RealScalar(0))) + A21 /= realAkk; + if (sign == PositiveSemiDef) { if (realAkk < 0) sign = Indefinite; } else if (sign == NegativeSemiDef) { @@ -454,6 +442,7 @@ LDLT& LDLT::compute(const MatrixType& a) m_transpositions.resize(size); m_isInitialized = false; m_temporary.resize(size); + m_sign = internal::ZeroSign; internal::ldlt_inplace::unblocked(m_matrix, m_transpositions, m_temporary, m_sign); @@ -514,16 +503,21 @@ struct solve_retval, Rhs> using std::abs; using std::max; typedef typename LDLTType::MatrixType MatrixType; - typedef typename LDLTType::Scalar Scalar; typedef typename LDLTType::RealScalar RealScalar; - const Diagonal vectorD = dec().vectorD(); - RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits::epsilon(), - RealScalar(1) / NumTraits::highest()); // motivated by LAPACK's xGELSS + const typename Diagonal::RealReturnType vectorD(dec().vectorD()); + // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon + // as motivated by LAPACK's xGELSS: + // RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() *NumTraits::epsilon(),RealScalar(1) / NumTraits::highest()); + // However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest + // diagonal element is not well justified and to numerical issues in some cases. + // Moreover, Lapack's xSYTRS routines use 0 for the tolerance. + RealScalar tolerance = RealScalar(1) / NumTraits::highest(); + for (Index i = 0; i < vectorD.size(); ++i) { if(abs(vectorD(i)) > tolerance) - dst.row(i) /= vectorD(i); + dst.row(i) /= vectorD(i); else - dst.row(i).setZero(); + dst.row(i).setZero(); } // dst = L^-T (D^-1 L^-1 P b) @@ -576,7 +570,7 @@ MatrixType LDLT::reconstructedMatrix() const // L^* P res = matrixU() * res; // D(L^*P) - res = vectorD().asDiagonal() * res; + res = vectorD().real().asDiagonal() * res; // L(DL^*P) res = matrixL() * res; // P^T (LDL^*P) diff --git a/uppsrc/plugin/Eigen/Eigen/src/CholmodSupport/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/CholmodSupport/CMakeLists.txt deleted file mode 100644 index 814dfa613..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/CholmodSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_CholmodSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_CholmodSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/CholmodSupport COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/ArrayWrapper.h b/uppsrc/plugin/Eigen/Eigen/src/Core/ArrayWrapper.h index a791bc358..b4641e2a0 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/ArrayWrapper.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/ArrayWrapper.h @@ -29,6 +29,11 @@ struct traits > : public traits::type > { typedef ArrayXpr XprKind; + // Let's remove NestByRefBit + enum { + Flags0 = traits::type >::Flags, + Flags = Flags0 & ~NestByRefBit + }; }; } @@ -149,6 +154,11 @@ struct traits > : public traits::type > { typedef MatrixXpr XprKind; + // Let's remove NestByRefBit + enum { + Flags0 = traits::type >::Flags, + Flags = Flags0 & ~NestByRefBit + }; }; } diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/Block.h b/uppsrc/plugin/Eigen/Eigen/src/Core/Block.h index 358b3188b..87bedfa46 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/Block.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/Block.h @@ -81,7 +81,7 @@ struct traits > : traits::Flags&LinearAccessBit))) ? LinearAccessBit : 0, FlagsLvalueBit = is_lvalue::value ? LvalueBit : 0, FlagsRowMajorBit = IsRowMajor ? RowMajorBit : 0, Flags0 = traits::Flags & ( (HereditaryBits & ~RowMajorBit) | diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Core/CMakeLists.txt deleted file mode 100644 index 2346fc2bb..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/CMakeLists.txt +++ /dev/null @@ -1,10 +0,0 @@ -FILE(GLOB Eigen_Core_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Core_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core COMPONENT Devel - ) - -ADD_SUBDIRECTORY(products) -ADD_SUBDIRECTORY(util) -ADD_SUBDIRECTORY(arch) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/CommaInitializer.h b/uppsrc/plugin/Eigen/Eigen/src/Core/CommaInitializer.h index a96867af4..a036d8c3b 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/CommaInitializer.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/CommaInitializer.h @@ -43,6 +43,17 @@ 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 + 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(o).m_row = m_xpr.rows(); + const_cast(o).m_col = m_xpr.cols(); + const_cast(o).m_currentBlockRows = 0; + } + /* inserts a scalar value in the target matrix */ CommaInitializer& operator,(const Scalar& s) { diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/DenseBase.h b/uppsrc/plugin/Eigen/Eigen/src/Core/DenseBase.h index c5800f6c8..04862f374 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/DenseBase.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/DenseBase.h @@ -462,8 +462,10 @@ template class DenseBase template RealScalar lpNorm() const; template - const Replicate replicate() const; - const Replicate replicate(Index rowFacor,Index colFactor) const; + inline const Replicate replicate() const; + + typedef Replicate ReplicateReturnType; + inline const ReplicateReturnType replicate(Index rowFacor,Index colFactor) const; typedef Reverse ReverseReturnType; typedef const Reverse ConstReverseReturnType; diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/DenseStorage.h b/uppsrc/plugin/Eigen/Eigen/src/Core/DenseStorage.h index 3e7f9c1b7..a72738e55 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/DenseStorage.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/DenseStorage.h @@ -24,6 +24,14 @@ namespace internal { struct constructor_without_unaligned_array_assert {}; +template void check_static_allocation_size() +{ + // if EIGEN_STACK_ALLOCATION_LIMIT is defined to 0, then no limit + #if EIGEN_STACK_ALLOCATION_LIMIT + EIGEN_STATIC_ASSERT(Size * sizeof(T) <= EIGEN_STACK_ALLOCATION_LIMIT, OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG); + #endif +} + /** \internal * Static array. If the MatrixOrArrayOptions require auto-alignment, the array will be automatically aligned: * to 16 bytes boundary if the total size is a multiple of 16 bytes. @@ -38,12 +46,12 @@ struct plain_array plain_array() { - EIGEN_STATIC_ASSERT(Size * sizeof(T) <= 128 * 128 * 8, OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG); + check_static_allocation_size(); } plain_array(constructor_without_unaligned_array_assert) { - EIGEN_STATIC_ASSERT(Size * sizeof(T) <= 128 * 128 * 8, OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG); + check_static_allocation_size(); } }; @@ -76,12 +84,12 @@ struct plain_array plain_array() { EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(0xf); - EIGEN_STATIC_ASSERT(Size * sizeof(T) <= 128 * 128 * 8, OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG); + check_static_allocation_size(); } plain_array(constructor_without_unaligned_array_assert) { - EIGEN_STATIC_ASSERT(Size * sizeof(T) <= 128 * 128 * 8, OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG); + check_static_allocation_size(); } }; diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/Diagonal.h b/uppsrc/plugin/Eigen/Eigen/src/Core/Diagonal.h index aab8007b3..68cf6d4b0 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/Diagonal.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/Diagonal.h @@ -190,18 +190,18 @@ MatrixBase::diagonal() const * * \sa MatrixBase::diagonal(), class Diagonal */ template -inline typename MatrixBase::template DiagonalIndexReturnType::Type +inline typename MatrixBase::DiagonalDynamicIndexReturnType MatrixBase::diagonal(Index index) { - return typename DiagonalIndexReturnType::Type(derived(), index); + return DiagonalDynamicIndexReturnType(derived(), index); } /** This is the const version of diagonal(Index). */ template -inline typename MatrixBase::template ConstDiagonalIndexReturnType::Type +inline typename MatrixBase::ConstDiagonalDynamicIndexReturnType MatrixBase::diagonal(Index index) const { - return typename ConstDiagonalIndexReturnType::Type(derived(), index); + return ConstDiagonalDynamicIndexReturnType(derived(), index); } /** \returns an expression of the \a DiagIndex-th sub or super diagonal of the matrix \c *this diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/Functors.h b/uppsrc/plugin/Eigen/Eigen/src/Core/Functors.h index 04fb21732..b08b967ff 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/Functors.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/Functors.h @@ -589,7 +589,7 @@ struct linspaced_op_impl template EIGEN_STRONG_INLINE const Packet packetOp(Index i) const - { return internal::padd(m_lowPacket, pmul(m_stepPacket, padd(pset1(i),m_interPacket))); } + { return internal::padd(m_lowPacket, pmul(m_stepPacket, padd(pset1(Scalar(i)),m_interPacket))); } const Scalar m_low; const Scalar m_step; @@ -609,7 +609,7 @@ template struct functor_traits< linspaced_o template struct linspaced_op { typedef typename packet_traits::type Packet; - linspaced_op(const Scalar& low, const Scalar& high, DenseIndex num_steps) : impl((num_steps==1 ? high : low), (num_steps==1 ? Scalar() : (high-low)/(num_steps-1))) {} + linspaced_op(const Scalar& low, const Scalar& high, DenseIndex num_steps) : impl((num_steps==1 ? high : low), (num_steps==1 ? Scalar() : (high-low)/Scalar(num_steps-1))) {} template EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return impl(i); } diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/GeneralProduct.h b/uppsrc/plugin/Eigen/Eigen/src/Core/GeneralProduct.h index 2a59d9464..9e805a80f 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/GeneralProduct.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/GeneralProduct.h @@ -232,7 +232,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest& // FIXME not very good if rhs is real and lhs complex while alpha is real too const Index cols = dest.cols(); for (Index j=0; j diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/MapBase.h b/uppsrc/plugin/Eigen/Eigen/src/Core/MapBase.h index 6876de588..b145a0b3d 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/MapBase.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/MapBase.h @@ -234,9 +234,15 @@ template class MapBase return derived(); } - using Base::Base::operator=; + // In theory MapBase should not make a using Base::operator=, + // and thus we should directly do: using Base::Base::operator=; + // However, this would confuse recent MSVC 2013 (bug 821), and since MapBase + // has operator= to make ICC 11 happy, we can also make MSVC 2013 happy as follow: + using Base::operator=; }; +#undef EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS + } // end namespace Eigen #endif // EIGEN_MAPBASE_H diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/MatrixBase.h b/uppsrc/plugin/Eigen/Eigen/src/Core/MatrixBase.h index 344b38f2f..cc3279746 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/MatrixBase.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/MatrixBase.h @@ -215,7 +215,7 @@ template class MatrixBase typedef Diagonal DiagonalReturnType; DiagonalReturnType diagonal(); - typedef typename internal::add_const >::type ConstDiagonalReturnType; + typedef typename internal::add_const >::type ConstDiagonalReturnType; ConstDiagonalReturnType diagonal() const; template struct DiagonalIndexReturnType { typedef Diagonal Type; }; @@ -223,16 +223,12 @@ template class MatrixBase template typename DiagonalIndexReturnType::Type diagonal(); template typename ConstDiagonalIndexReturnType::Type diagonal() const; + + typedef Diagonal DiagonalDynamicIndexReturnType; + typedef typename internal::add_const >::type ConstDiagonalDynamicIndexReturnType; - // Note: The "MatrixBase::" prefixes are added to help MSVC9 to match these declarations with the later implementations. - // On the other hand they confuse MSVC8... - #if (defined _MSC_VER) && (_MSC_VER >= 1500) // 2008 or later - typename MatrixBase::template DiagonalIndexReturnType::Type diagonal(Index index); - typename MatrixBase::template ConstDiagonalIndexReturnType::Type diagonal(Index index) const; - #else - typename DiagonalIndexReturnType::Type diagonal(Index index); - typename ConstDiagonalIndexReturnType::Type diagonal(Index index) const; - #endif + DiagonalDynamicIndexReturnType diagonal(Index index); + ConstDiagonalDynamicIndexReturnType diagonal(Index index) const; #ifdef EIGEN2_SUPPORT template typename internal::eigen2_part_return_type::type part(); diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/PermutationMatrix.h b/uppsrc/plugin/Eigen/Eigen/src/Core/PermutationMatrix.h index 1297b8413..f26f3e5cc 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/PermutationMatrix.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/PermutationMatrix.h @@ -555,7 +555,10 @@ struct permut_matrix_product_retval const Index n = Side==OnTheLeft ? rows() : cols(); // FIXME we need an is_same for expression that is not sensitive to constness. For instance // is_same_xpr, Block >::value should be true. - if(is_same::value && extract_data(dst) == extract_data(m_matrix)) + if( is_same::value + && blas_traits::HasUsableDirectAccess + && blas_traits::HasUsableDirectAccess + && extract_data(dst) == extract_data(m_matrix)) { // apply the permutation inplace Matrix mask(m_permutation.size()); diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/ProductBase.h b/uppsrc/plugin/Eigen/Eigen/src/Core/ProductBase.h index a494b5f87..cf74470a9 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/ProductBase.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/ProductBase.h @@ -85,7 +85,14 @@ class ProductBase : public MatrixBase public: +#ifndef EIGEN_NO_MALLOC + typedef typename Base::PlainObject BasePlainObject; + typedef Matrix DynPlainObject; + typedef typename internal::conditional<(BasePlainObject::SizeAtCompileTime==Dynamic) || (BasePlainObject::SizeAtCompileTime*int(sizeof(Scalar)) < int(EIGEN_STACK_ALLOCATION_LIMIT)), + BasePlainObject, DynPlainObject>::type PlainObject; +#else typedef typename Base::PlainObject PlainObject; +#endif ProductBase(const Lhs& a_lhs, const Rhs& a_rhs) : m_lhs(a_lhs), m_rhs(a_rhs) @@ -180,7 +187,12 @@ namespace internal { template struct nested, N, PlainObject> { - typedef PlainObject const& type; + typedef typename GeneralProduct::PlainObject const& type; +}; +template +struct nested, N, PlainObject> +{ + typedef typename GeneralProduct::PlainObject const& type; }; } diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/Ref.h b/uppsrc/plugin/Eigen/Eigen/src/Core/Ref.h index 00d9e6d2b..df87b1af4 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/Ref.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/Ref.h @@ -101,7 +101,7 @@ struct traits > template struct match { enum { HasDirectAccess = internal::has_direct_access::ret, - StorageOrderMatch = PlainObjectType::IsVectorAtCompileTime || ((PlainObjectType::Flags&RowMajorBit)==(Derived::Flags&RowMajorBit)), + StorageOrderMatch = PlainObjectType::IsVectorAtCompileTime || Derived::IsVectorAtCompileTime || ((PlainObjectType::Flags&RowMajorBit)==(Derived::Flags&RowMajorBit)), InnerStrideMatch = int(StrideType::InnerStrideAtCompileTime)==int(Dynamic) || int(StrideType::InnerStrideAtCompileTime)==int(Derived::InnerStrideAtCompileTime) || (int(StrideType::InnerStrideAtCompileTime)==0 && int(Derived::InnerStrideAtCompileTime)==1), @@ -172,8 +172,12 @@ protected: } else ::new (static_cast(this)) Base(expr.data(), expr.rows(), expr.cols()); - ::new (&m_stride) StrideBase(StrideType::OuterStrideAtCompileTime==0?0:expr.outerStride(), - StrideType::InnerStrideAtCompileTime==0?0:expr.innerStride()); + + if(Expression::IsVectorAtCompileTime && (!PlainObjectType::IsVectorAtCompileTime) && ((Expression::Flags&RowMajorBit)!=(PlainObjectType::Flags&RowMajorBit))) + ::new (&m_stride) StrideBase(expr.innerStride(), StrideType::InnerStrideAtCompileTime==0?0:1); + else + ::new (&m_stride) StrideBase(StrideType::OuterStrideAtCompileTime==0?0:expr.outerStride(), + StrideType::InnerStrideAtCompileTime==0?0:expr.innerStride()); } StrideBase m_stride; @@ -184,6 +188,8 @@ template class Ref : public RefBase > { typedef internal::traits Traits; + template + inline Ref(const PlainObjectBase& expr); public: typedef RefBase Base; @@ -192,20 +198,21 @@ template class Ref #ifndef EIGEN_PARSED_BY_DOXYGEN template - inline Ref(PlainObjectBase& expr, - typename internal::enable_if::MatchAtCompileTime),Derived>::type* = 0) + inline Ref(PlainObjectBase& expr) { - Base::construct(expr); + EIGEN_STATIC_ASSERT(static_cast(Traits::template match::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); + Base::construct(expr.derived()); } template - inline Ref(const DenseBase& expr, - typename internal::enable_if::value&&bool(Traits::template match::MatchAtCompileTime)),Derived>::type* = 0, - int = Derived::ThisConstantIsPrivateInPlainObjectBase) + inline Ref(const DenseBase& expr) #else template inline Ref(DenseBase& expr) #endif { + EIGEN_STATIC_ASSERT(static_cast(internal::is_lvalue::value), THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY); + EIGEN_STATIC_ASSERT(static_cast(Traits::template match::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH); + enum { THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY = Derived::ThisConstantIsPrivateInPlainObjectBase}; Base::construct(expr.const_cast_derived()); } diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/Replicate.h b/uppsrc/plugin/Eigen/Eigen/src/Core/Replicate.h index dde86a834..ac4537c14 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/Replicate.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/Replicate.h @@ -135,7 +135,7 @@ template class Replicate */ template template -inline const Replicate +const Replicate DenseBase::replicate() const { return Replicate(derived()); @@ -150,7 +150,7 @@ DenseBase::replicate() const * \sa VectorwiseOp::replicate(), DenseBase::replicate(), class Replicate */ template -inline const Replicate +const typename DenseBase::ReplicateReturnType DenseBase::replicate(Index rowFactor,Index colFactor) const { return Replicate(derived(),rowFactor,colFactor); diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/TriangularMatrix.h b/uppsrc/plugin/Eigen/Eigen/src/Core/TriangularMatrix.h index fba07365f..4d65392c6 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/TriangularMatrix.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/TriangularMatrix.h @@ -278,21 +278,21 @@ template class TriangularView /** Efficient triangular matrix times vector/matrix product */ template - TriangularProduct + TriangularProduct operator*(const MatrixBase& rhs) const { return TriangularProduct - + (m_matrix, rhs.derived()); } /** Efficient vector/matrix times triangular matrix product */ template friend - TriangularProduct + TriangularProduct operator*(const MatrixBase& lhs, const TriangularView& rhs) { return TriangularProduct - + (lhs.derived(),rhs.m_matrix); } @@ -380,19 +380,19 @@ template class TriangularView EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase& other) { setZero(); - return assignProduct(other,1); + return assignProduct(other.derived(),1); } template EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase& other) { - return assignProduct(other,1); + return assignProduct(other.derived(),1); } template EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase& other) { - return assignProduct(other,-1); + return assignProduct(other.derived(),-1); } @@ -400,25 +400,34 @@ template class TriangularView EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct& other) { setZero(); - return assignProduct(other,other.alpha()); + return assignProduct(other.derived(),other.alpha()); } template EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct& other) { - return assignProduct(other,other.alpha()); + return assignProduct(other.derived(),other.alpha()); } template EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct& other) { - return assignProduct(other,-other.alpha()); + return assignProduct(other.derived(),-other.alpha()); } protected: template EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase& prod, const Scalar& alpha); + + template + EIGEN_STRONG_INLINE TriangularView& assignProduct(const TriangularProduct& prod, const Scalar& alpha) + { + lazyAssign(alpha*prod.eval()); + return *this; + } MatrixTypeNested m_matrix; }; diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/AltiVec/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/AltiVec/CMakeLists.txt deleted file mode 100644 index 9f8d2e9c4..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/AltiVec/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Core_arch_AltiVec_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Core_arch_AltiVec_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/arch/AltiVec COMPONENT Devel -) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/CMakeLists.txt deleted file mode 100644 index 8456dec15..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/CMakeLists.txt +++ /dev/null @@ -1,4 +0,0 @@ -ADD_SUBDIRECTORY(SSE) -ADD_SUBDIRECTORY(AltiVec) -ADD_SUBDIRECTORY(NEON) -ADD_SUBDIRECTORY(Default) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/Default/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/Default/CMakeLists.txt deleted file mode 100644 index 339c091d1..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/Default/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Core_arch_Default_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Core_arch_Default_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/arch/Default COMPONENT Devel -) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/NEON/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/NEON/CMakeLists.txt deleted file mode 100644 index fd4d4af50..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/NEON/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Core_arch_NEON_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Core_arch_NEON_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/arch/NEON COMPONENT Devel -) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/NEON/PacketMath.h b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/NEON/PacketMath.h index 163bac215..a4dbdc1ea 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/NEON/PacketMath.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/NEON/PacketMath.h @@ -48,9 +48,18 @@ typedef uint32x4_t Packet4ui; #define EIGEN_INIT_NEON_PACKET2(X, Y) {X, Y} #define EIGEN_INIT_NEON_PACKET4(X, Y, Z, W) {X, Y, Z, W} #endif - -#ifndef __pld -#define __pld(x) asm volatile ( " pld [%[addr]]\n" :: [addr] "r" (x) : "cc" ); + +// arm64 does have the pld instruction. If available, let's trust the __builtin_prefetch built-in function +// which available on LLVM and GCC (at least) +#if EIGEN_HAS_BUILTIN(__builtin_prefetch) || defined(__GNUC__) + #define EIGEN_ARM_PREFETCH(ADDR) __builtin_prefetch(ADDR); +#elif defined __pld + #define EIGEN_ARM_PREFETCH(ADDR) __pld(ADDR) +#elif !defined(__aarch64__) + #define EIGEN_ARM_PREFETCH(ADDR) __asm__ __volatile__ ( " pld [%[addr]]\n" :: [addr] "r" (ADDR) : "cc" ); +#else + // by default no explicit prefetching + #define EIGEN_ARM_PREFETCH(ADDR) #endif template<> struct packet_traits : default_packet_traits diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/SSE/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/SSE/CMakeLists.txt deleted file mode 100644 index 46ea7cc62..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/SSE/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Core_arch_SSE_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Core_arch_SSE_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/arch/SSE COMPONENT Devel -) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/SSE/MathFunctions.h b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/SSE/MathFunctions.h index 99cbd0d95..d16f30bb0 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/arch/SSE/MathFunctions.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/arch/SSE/MathFunctions.h @@ -52,7 +52,7 @@ Packet4f plog(const Packet4f& _x) Packet4i emm0; - Packet4f invalid_mask = _mm_cmplt_ps(x, _mm_setzero_ps()); + Packet4f invalid_mask = _mm_cmpnge_ps(x, _mm_setzero_ps()); // not greater equal is true if x is NaN Packet4f iszero_mask = _mm_cmpeq_ps(x, _mm_setzero_ps()); x = pmax(x, p4f_min_norm_pos); /* cut off denormalized stuff */ @@ -166,7 +166,7 @@ Packet4f pexp(const Packet4f& _x) emm0 = _mm_cvttps_epi32(fx); emm0 = _mm_add_epi32(emm0, p4i_0x7f); emm0 = _mm_slli_epi32(emm0, 23); - return pmul(y, _mm_castsi128_ps(emm0)); + return pmax(pmul(y, Packet4f(_mm_castsi128_ps(emm0))), _x); } template<> EIGEN_DEFINE_FUNCTION_ALLOWING_MULTIPLE_DEFINITIONS EIGEN_UNUSED Packet2d pexp(const Packet2d& _x) @@ -239,7 +239,7 @@ Packet2d pexp(const Packet2d& _x) emm0 = _mm_add_epi32(emm0, p4i_1023_0); emm0 = _mm_slli_epi32(emm0, 20); emm0 = _mm_shuffle_epi32(emm0, _MM_SHUFFLE(1,2,0,3)); - return pmul(x, _mm_castsi128_pd(emm0)); + return pmax(pmul(x, Packet2d(_mm_castsi128_pd(emm0))), _x); } /* evaluation of 4 sines at onces, using SSE2 intrinsics. diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/products/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Core/products/CMakeLists.txt deleted file mode 100644 index 21fc94ae3..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/products/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Core_Product_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Core_Product_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/products COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/products/CoeffBasedProduct.h b/uppsrc/plugin/Eigen/Eigen/src/Core/products/CoeffBasedProduct.h index c06a0df1c..421f925e1 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/products/CoeffBasedProduct.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/products/CoeffBasedProduct.h @@ -90,6 +90,7 @@ struct traits > | (SameType && (CanVectorizeLhs || CanVectorizeRhs) ? PacketAccessBit : 0), CoeffReadCost = InnerSize == Dynamic ? Dynamic + : InnerSize == 0 ? 0 : InnerSize * (NumTraits::MulCost + LhsCoeffReadCost + RhsCoeffReadCost) + (InnerSize - 1) * NumTraits::AddCost, @@ -133,7 +134,7 @@ class CoeffBasedProduct }; typedef internal::product_coeff_impl ScalarCoeffImpl; typedef CoeffBasedProduct LazyCoeffBasedProductType; @@ -184,7 +185,7 @@ class CoeffBasedProduct { PacketScalar res; internal::product_packet_impl ::run(row, col, m_lhs, m_rhs, res); return res; @@ -262,10 +263,7 @@ struct product_coeff_impl typedef typename Lhs::Index Index; static EIGEN_STRONG_INLINE void run(Index row, Index col, const Lhs& lhs, const Rhs& rhs, RetScalar& res) { - eigen_assert(lhs.cols()>0 && "you are using a non initialized matrix"); - res = lhs.coeff(row, 0) * rhs.coeff(0, col); - for(Index i = 1; i < lhs.cols(); ++i) - res += lhs.coeff(row, i) * rhs.coeff(i, col); + res = (lhs.row(row).transpose().cwiseProduct( rhs.col(col) )).sum(); } }; diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/util/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Core/util/CMakeLists.txt deleted file mode 100644 index a1e2e521f..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/util/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Core_util_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Core_util_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Core/util COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/util/MKL_support.h b/uppsrc/plugin/Eigen/Eigen/src/Core/util/MKL_support.h index 1e6e355d6..8acca9c8c 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/util/MKL_support.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/util/MKL_support.h @@ -54,8 +54,25 @@ #endif #if defined EIGEN_USE_MKL +# include +/*Check IMKL version for compatibility: < 10.3 is not usable with Eigen*/ +# ifndef INTEL_MKL_VERSION +# undef EIGEN_USE_MKL /* INTEL_MKL_VERSION is not even defined on older versions */ +# elif INTEL_MKL_VERSION < 100305 /* the intel-mkl-103-release-notes say this was when the lapacke.h interface was added*/ +# undef EIGEN_USE_MKL +# endif +# ifndef EIGEN_USE_MKL + /*If the MKL version is too old, undef everything*/ +# undef EIGEN_USE_MKL_ALL +# undef EIGEN_USE_BLAS +# undef EIGEN_USE_LAPACKE +# undef EIGEN_USE_MKL_VML +# undef EIGEN_USE_LAPACKE_STRICT +# undef EIGEN_USE_LAPACKE +# endif +#endif -#include +#if defined EIGEN_USE_MKL #include #define EIGEN_MKL_VML_THRESHOLD 128 diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/util/Macros.h b/uppsrc/plugin/Eigen/Eigen/src/Core/util/Macros.h index 0088621ac..a4f6f22e7 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/util/Macros.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/util/Macros.h @@ -13,7 +13,7 @@ #define EIGEN_WORLD_VERSION 3 #define EIGEN_MAJOR_VERSION 2 -#define EIGEN_MINOR_VERSION 1 +#define EIGEN_MINOR_VERSION 3 #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ @@ -96,6 +96,13 @@ #define EIGEN_DEFAULT_DENSE_INDEX_TYPE std::ptrdiff_t #endif +// Cross compiler wrapper around LLVM's __has_builtin +#ifdef __has_builtin +# define EIGEN_HAS_BUILTIN(x) __has_builtin(x) +#else +# define EIGEN_HAS_BUILTIN(x) 0 +#endif + /** Allows to disable some optimizations which might affect the accuracy of the result. * Such optimization are enabled by default, and set EIGEN_FAST_MATH to 0 to disable them. * They currently include: @@ -247,7 +254,7 @@ namespace Eigen { #if !defined(EIGEN_ASM_COMMENT) #if (defined __GNUC__) && ( defined(__i386__) || defined(__x86_64__) ) - #define EIGEN_ASM_COMMENT(X) asm("#" X) + #define EIGEN_ASM_COMMENT(X) __asm__("#" X) #else #define EIGEN_ASM_COMMENT(X) #endif @@ -289,7 +296,8 @@ namespace Eigen { #endif #ifndef EIGEN_STACK_ALLOCATION_LIMIT -#define EIGEN_STACK_ALLOCATION_LIMIT 20000 +// 131072 == 128 KB +#define EIGEN_STACK_ALLOCATION_LIMIT 131072 #endif #ifndef EIGEN_DEFAULT_IO_FORMAT diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/util/Memory.h b/uppsrc/plugin/Eigen/Eigen/src/Core/util/Memory.h index cacbd02fd..07e54ef06 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/util/Memory.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/util/Memory.h @@ -63,7 +63,7 @@ // Currently, let's include it only on unix systems: #if defined(__unix__) || defined(__unix) #include - #if ((defined __QNXNTO__) || (defined _GNU_SOURCE) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0) + #if ((defined __QNXNTO__) || (defined _GNU_SOURCE) || (defined __PGI) || ((defined _XOPEN_SOURCE) && (_XOPEN_SOURCE >= 600))) && (defined _POSIX_ADVISORY_INFO) && (_POSIX_ADVISORY_INFO > 0) #define EIGEN_HAS_POSIX_MEMALIGN 1 #endif #endif @@ -272,12 +272,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); @@ -417,6 +417,8 @@ template inline T* conditional_aligned_realloc_new(T* pt template inline T* conditional_aligned_new_auto(size_t size) { + if(size==0) + return 0; // short-cut. Also fixes Bug 884 check_size_for_overflow(size); T *result = reinterpret_cast(conditional_aligned_malloc(sizeof(T)*size)); if(NumTraits::RequireInitialization) @@ -612,7 +614,6 @@ template class aligned_stack_memory_handler void* operator new(size_t size, const std::nothrow_t&) throw() { \ try { return Eigen::internal::conditional_aligned_malloc(size); } \ catch (...) { return 0; } \ - return 0; \ } #else #define EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \ @@ -777,9 +778,9 @@ namespace internal { #ifdef EIGEN_CPUID -inline bool cpuid_is_vendor(int abcd[4], const char* vendor) +inline bool cpuid_is_vendor(int abcd[4], const int vendor[3]) { - return abcd[1]==(reinterpret_cast(vendor))[0] && abcd[3]==(reinterpret_cast(vendor))[1] && abcd[2]==(reinterpret_cast(vendor))[2]; + return abcd[1]==vendor[0] && abcd[3]==vendor[1] && abcd[2]==vendor[2]; } inline void queryCacheSizes_intel_direct(int& l1, int& l2, int& l3) @@ -921,13 +922,16 @@ inline void queryCacheSizes(int& l1, int& l2, int& l3) { #ifdef EIGEN_CPUID int abcd[4]; + const int GenuineIntel[] = {0x756e6547, 0x49656e69, 0x6c65746e}; + const int AuthenticAMD[] = {0x68747541, 0x69746e65, 0x444d4163}; + const int AMDisbetter_[] = {0x69444d41, 0x74656273, 0x21726574}; // "AMDisbetter!" // identify the CPU vendor EIGEN_CPUID(abcd,0x0,0); int max_std_funcs = abcd[1]; - if(cpuid_is_vendor(abcd,"GenuineIntel")) + if(cpuid_is_vendor(abcd,GenuineIntel)) queryCacheSizes_intel(l1,l2,l3,max_std_funcs); - else if(cpuid_is_vendor(abcd,"AuthenticAMD") || cpuid_is_vendor(abcd,"AMDisbetter!")) + else if(cpuid_is_vendor(abcd,AuthenticAMD) || cpuid_is_vendor(abcd,AMDisbetter_)) queryCacheSizes_amd(l1,l2,l3); else // by default let's use Intel's API diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/util/StaticAssert.h b/uppsrc/plugin/Eigen/Eigen/src/Core/util/StaticAssert.h index 8872c5b64..bac5d9fe9 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/util/StaticAssert.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/util/StaticAssert.h @@ -90,7 +90,9 @@ YOU_PASSED_A_COLUMN_VECTOR_BUT_A_ROW_VECTOR_WAS_EXPECTED, THE_INDEX_TYPE_MUST_BE_A_SIGNED_TYPE, THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH, - OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG + OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG, + IMPLICIT_CONVERSION_TO_SCALAR_IS_FOR_INNER_PRODUCT_ONLY, + STORAGE_LAYOUT_DOES_NOT_MATCH }; }; diff --git a/uppsrc/plugin/Eigen/Eigen/src/Core/util/XprHelper.h b/uppsrc/plugin/Eigen/Eigen/src/Core/util/XprHelper.h index 3c4773054..781965d2c 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Core/util/XprHelper.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Core/util/XprHelper.h @@ -341,7 +341,7 @@ template::type> str }; template -T* const_cast_ptr(const T* ptr) +inline T* const_cast_ptr(const T* ptr) { return const_cast(ptr); } diff --git a/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/CMakeLists.txt deleted file mode 100644 index 7ae41b3cb..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -FILE(GLOB Eigen_Eigen2Support_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Eigen2Support_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support COMPONENT Devel - ) - -ADD_SUBDIRECTORY(Geometry) \ No newline at end of file diff --git a/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt deleted file mode 100644 index c347a8f26..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/Geometry/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Eigen2Support_Geometry_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Eigen2Support_Geometry_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigen2Support/Geometry - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/LeastSquares.h b/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/LeastSquares.h index 0e6fdb488..7992d4944 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/LeastSquares.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Eigen2Support/LeastSquares.h @@ -147,7 +147,6 @@ void fitHyperplane(int numPoints, // compute the covariance matrix CovMatrixType covMat = CovMatrixType::Zero(size, size); - VectorType remean = VectorType::Zero(size); for(int i = 0; i < numPoints; ++i) { VectorType diff = (*(points[i]) - mean).conjugate(); diff --git a/uppsrc/plugin/Eigen/Eigen/src/Eigenvalues/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Eigenvalues/CMakeLists.txt deleted file mode 100644 index 193e02685..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Eigenvalues/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_EIGENVALUES_SRCS "*.h") - -INSTALL(FILES - ${Eigen_EIGENVALUES_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigenvalues COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/uppsrc/plugin/Eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 3993046a8..be89de4a9 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -563,7 +563,6 @@ template struct direct_selfadjoint_eigenvalues::epsilon(); - safeNorm2 *= safeNorm2; if((eivals(2)-eivals(0))<=Eigen::NumTraits::epsilon()) { eivecs.setIdentity(); @@ -577,7 +576,7 @@ template struct direct_selfadjoint_eigenvalues d1 ? 2 : 0; - d0 = d0 > d1 ? d1 : d0; + d0 = d0 > d1 ? d0 : d1; tmp.diagonal().array () -= eivals(k); VectorType cross; @@ -585,19 +584,25 @@ template struct direct_selfadjoint_eigenvaluessafeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { // the input matrix and/or the eigenvaues probably contains some inf/NaN, @@ -617,12 +622,16 @@ template struct direct_selfadjoint_eigenvalues::epsilon()) + { eivecs.col(1) = eivecs.col(k).unitOrthogonal(); + } else { - n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm(); + n = (cross = eivecs.col(k).cross(tmp.row(0))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(1) = cross / sqrt(n); + } else { n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm(); @@ -636,13 +645,14 @@ template struct direct_selfadjoint_eigenvalues::epsilon()) + { + Matrix m; m << v0.transpose(), v1.transpose(); + JacobiSVD > svd(m, ComputeFullV); + result.normal() = svd.matrixV().col(2); + } + else + result.normal() /= norm; result.offset() = -p0.dot(result.normal()); return result; } diff --git a/uppsrc/plugin/Eigen/Eigen/src/Geometry/Quaternion.h b/uppsrc/plugin/Eigen/Eigen/src/Geometry/Quaternion.h index 9fee0c919..f06653f1c 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Geometry/Quaternion.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Geometry/Quaternion.h @@ -203,6 +203,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 */ @@ -344,7 +346,7 @@ class Map, _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. */ @@ -464,7 +466,7 @@ QuaternionBase::_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; @@ -584,7 +586,7 @@ inline Derived& QuaternionBase::setFromTwoVectors(const MatrixBase::dummy_precision()) { - c = max(c,-1); + c = (max)(c,Scalar(-1)); Matrix m; m << v0.transpose(), v1.transpose(); JacobiSVD > svd(m, ComputeFullV); Vector3 axis = svd.matrixV().col(2); @@ -667,10 +669,10 @@ QuaternionBase::angularDistance(const QuaternionBase& 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(2 * acos(d)); + return Scalar(2) * acos(d); } diff --git a/uppsrc/plugin/Eigen/Eigen/src/Geometry/Rotation2D.h b/uppsrc/plugin/Eigen/Eigen/src/Geometry/Rotation2D.h index 1cac343a5..776c36144 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Geometry/Rotation2D.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Geometry/Rotation2D.h @@ -59,7 +59,10 @@ protected: public: /** Construct a 2D counter clock wise rotation from the angle \a a in radian. */ - inline Rotation2D(const Scalar& a) : m_angle(a) {} + explicit inline Rotation2D(const Scalar& a) : m_angle(a) {} + + /** Default constructor wihtout initialization. The represented rotation is undefined. */ + Rotation2D() {} /** \returns the rotation angle */ inline Scalar angle() const { return m_angle; } @@ -81,10 +84,10 @@ public: /** Applies the rotation to a 2D vector */ Vector2 operator* (const Vector2& vec) const { return toRotationMatrix() * vec; } - + template Rotation2D& fromRotationMatrix(const MatrixBase& m); - Matrix2 toRotationMatrix(void) const; + Matrix2 toRotationMatrix() const; /** \returns the spherical interpolation between \c *this and \a other using * parameter \a t. It is in fact equivalent to a linear interpolation. diff --git a/uppsrc/plugin/Eigen/Eigen/src/Geometry/Transform.h b/uppsrc/plugin/Eigen/Eigen/src/Geometry/Transform.h index 498fea41a..e786e5356 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Geometry/Transform.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Geometry/Transform.h @@ -62,6 +62,8 @@ struct transform_construct_from_matrix; template struct transform_take_affine_part; +template struct transform_make_affine; + } // end namespace internal /** \geometry_module \ingroup Geometry_Module @@ -194,9 +196,9 @@ public: /** type of the matrix used to represent the linear part of the transformation */ typedef Matrix LinearMatrixType; /** type of read/write reference to the linear part of the transformation */ - typedef Block LinearPart; + typedef Block LinearPart; /** type of read reference to the linear part of the transformation */ - typedef const Block ConstLinearPart; + typedef const Block ConstLinearPart; /** type of read/write reference to the affine part of the transformation */ typedef typename internal::conditional::run(m_matrix); } inline Transform(const Transform& other) @@ -591,11 +592,7 @@ public: */ void makeAffine() { - if(int(Mode)!=int(AffineCompact)) - { - matrix().template block<1,Dim>(Dim,0).setZero(); - matrix().coeffRef(Dim,Dim) = Scalar(1); - } + internal::transform_make_affine::run(m_matrix); } /** \internal @@ -1079,6 +1076,24 @@ Transform::fromPositionOrientationScale(const MatrixBas namespace internal { +template +struct transform_make_affine +{ + template + static void run(MatrixType &mat) + { + static const int Dim = MatrixType::ColsAtCompileTime-1; + mat.template block<1,Dim>(Dim,0).setZero(); + mat.coeffRef(Dim,Dim) = typename MatrixType::Scalar(1); + } +}; + +template<> +struct transform_make_affine +{ + template static void run(MatrixType &) { } +}; + // selector needed to avoid taking the inverse of a 3x4 matrix template struct projective_transform_inverse diff --git a/uppsrc/plugin/Eigen/Eigen/src/Geometry/Umeyama.h b/uppsrc/plugin/Eigen/Eigen/src/Geometry/Umeyama.h index 345b47e0c..5e20662f8 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Geometry/Umeyama.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Geometry/Umeyama.h @@ -113,7 +113,7 @@ umeyama(const MatrixBase& src, const MatrixBase& dst, boo const Index n = src.cols(); // number of measurements // required for demeaning ... - const RealScalar one_over_n = 1 / static_cast(n); + const RealScalar one_over_n = RealScalar(1) / static_cast(n); // computation of mean const VectorType src_mean = src.rowwise().sum() * one_over_n; @@ -136,16 +136,16 @@ umeyama(const MatrixBase& src, const MatrixBase& dst, boo // Eq. (39) VectorType S = VectorType::Ones(m); - if (sigma.determinant()<0) S(m-1) = -1; + if (sigma.determinant() 0 ) { + if ( svd.matrixU().determinant() * svd.matrixV().determinant() > Scalar(0) ) { Rt.block(0,0,m,m).noalias() = svd.matrixU()*svd.matrixV().transpose(); } else { - const Scalar s = S(m-1); S(m-1) = -1; + const Scalar s = S(m-1); S(m-1) = Scalar(-1); Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose(); S(m-1) = s; } @@ -156,7 +156,7 @@ umeyama(const MatrixBase& src, const MatrixBase& dst, boo if (with_scaling) { // Eq. (42) - const Scalar c = 1/src_var * svd.singularValues().dot(S); + const Scalar c = Scalar(1)/src_var * svd.singularValues().dot(S); // Eq. (41) Rt.col(m).head(m) = dst_mean; diff --git a/uppsrc/plugin/Eigen/Eigen/src/Geometry/arch/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Geometry/arch/CMakeLists.txt deleted file mode 100644 index 1267a79c7..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Geometry/arch/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Geometry_arch_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Geometry_arch_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Geometry/arch COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Householder/BlockHouseholder.h b/uppsrc/plugin/Eigen/Eigen/src/Householder/BlockHouseholder.h index 1991c6527..60dbea5f5 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/Householder/BlockHouseholder.h +++ b/uppsrc/plugin/Eigen/Eigen/src/Householder/BlockHouseholder.h @@ -48,7 +48,7 @@ void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vec typedef typename MatrixType::Index Index; enum { TFactorSize = MatrixType::ColsAtCompileTime }; Index nbVecs = vectors.cols(); - Matrix T(nbVecs,nbVecs); + Matrix T(nbVecs,nbVecs); make_block_householder_triangular_factor(T, vectors, hCoeffs); const TriangularView& V(vectors); diff --git a/uppsrc/plugin/Eigen/Eigen/src/Householder/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Householder/CMakeLists.txt deleted file mode 100644 index ce4937db0..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Householder/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Householder_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Householder_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Householder COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/uppsrc/plugin/Eigen/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index 6fc6ab852..dd135c21f 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/uppsrc/plugin/Eigen/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -39,7 +39,6 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, int maxIters = iters; int n = mat.cols(); - x = precond.solve(x); VectorType r = rhs - mat * x; VectorType r0 = r; @@ -61,6 +60,7 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, VectorType s(n), t(n); RealScalar tol2 = tol*tol; + RealScalar eps2 = NumTraits::epsilon()*NumTraits::epsilon(); int i = 0; int restarts = 0; @@ -69,7 +69,7 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, Scalar rho_old = rho; rho = r0.dot(r); - if (internal::isMuchSmallerThan(rho,r0_sqnorm)) + if (abs(rho) < eps2*r0_sqnorm) { // The new residual vector became too orthogonal to the arbitrarily choosen direction r0 // Let's restart with a new r0: @@ -142,7 +142,7 @@ struct traits > * SparseMatrix A(n,n); * // fill A and b * BiCGSTAB > solver; - * solver(A); + * solver.compute(A); * x = solver.solve(b); * std::cout << "#iterations: " << solver.iterations() << std::endl; * std::cout << "estimated error: " << solver.error() << std::endl; diff --git a/uppsrc/plugin/Eigen/Eigen/src/IterativeLinearSolvers/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/IterativeLinearSolvers/CMakeLists.txt deleted file mode 100644 index 59ccc0072..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/IterativeLinearSolvers/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_IterativeLinearSolvers_SRCS "*.h") - -INSTALL(FILES - ${Eigen_IterativeLinearSolvers_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/IterativeLinearSolvers COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/Jacobi/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/Jacobi/CMakeLists.txt deleted file mode 100644 index 490dac626..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/Jacobi/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Jacobi_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Jacobi_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Jacobi COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/LU/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/LU/CMakeLists.txt deleted file mode 100644 index e0d8d78c1..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/LU/CMakeLists.txt +++ /dev/null @@ -1,8 +0,0 @@ -FILE(GLOB Eigen_LU_SRCS "*.h") - -INSTALL(FILES - ${Eigen_LU_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LU COMPONENT Devel - ) - -ADD_SUBDIRECTORY(arch) diff --git a/uppsrc/plugin/Eigen/Eigen/src/LU/FullPivLU.h b/uppsrc/plugin/Eigen/Eigen/src/LU/FullPivLU.h index dfe25f424..79ab6a8c8 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/LU/FullPivLU.h +++ b/uppsrc/plugin/Eigen/Eigen/src/LU/FullPivLU.h @@ -20,10 +20,11 @@ namespace Eigen { * * \param MatrixType the type of the matrix of which we are computing the LU decomposition * - * This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A - * is decomposed as A = PLUQ where L is unit-lower-triangular, U is upper-triangular, and P and Q - * are permutation matrices. This is a rank-revealing LU decomposition. The eigenvalues (diagonal - * coefficients) of U are sorted in such a way that any zeros are at the end. + * This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A is + * decomposed as \f$ A = P^{-1} L U Q^{-1} \f$ where L is unit-lower-triangular, U is + * upper-triangular, and P and Q are permutation matrices. This is a rank-revealing LU + * decomposition. The eigenvalues (diagonal coefficients) of U are sorted in such a way that any + * zeros are at the end. * * This decomposition provides the generic approach to solving systems of linear equations, computing * the rank, invertibility, inverse, kernel, and determinant. @@ -511,8 +512,8 @@ typename internal::traits::Scalar FullPivLU::determinant } /** \returns the matrix represented by the decomposition, - * i.e., it returns the product: P^{-1} L U Q^{-1}. - * This function is provided for debug purpose. */ + * i.e., it returns the product: \f$ P^{-1} L U Q^{-1} \f$. + * This function is provided for debug purposes. */ template MatrixType FullPivLU::reconstructedMatrix() const { diff --git a/uppsrc/plugin/Eigen/Eigen/src/LU/arch/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/LU/arch/CMakeLists.txt deleted file mode 100644 index f6b7ed9ec..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/LU/arch/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_LU_arch_SRCS "*.h") - -INSTALL(FILES - ${Eigen_LU_arch_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LU/arch COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/MetisSupport/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/MetisSupport/CMakeLists.txt deleted file mode 100644 index 2bad31416..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/MetisSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_MetisSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_MetisSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/MetisSupport COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/OrderingMethods/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/OrderingMethods/CMakeLists.txt deleted file mode 100644 index 9f4bb2758..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/OrderingMethods/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_OrderingMethods_SRCS "*.h") - -INSTALL(FILES - ${Eigen_OrderingMethods_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/OrderingMethods COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/OrderingMethods/Ordering.h b/uppsrc/plugin/Eigen/Eigen/src/OrderingMethods/Ordering.h index b4da6531a..f3c31f9cb 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/OrderingMethods/Ordering.h +++ b/uppsrc/plugin/Eigen/Eigen/src/OrderingMethods/Ordering.h @@ -109,7 +109,7 @@ class NaturalOrdering * \class COLAMDOrdering * * Functor computing the \em column \em approximate \em minimum \em degree ordering - * The matrix should be in column-major format + * The matrix should be in column-major and \b compressed format (see SparseMatrix::makeCompressed()). */ template class COLAMDOrdering @@ -118,10 +118,14 @@ class COLAMDOrdering typedef PermutationMatrix PermutationType; typedef Matrix IndexVector; - /** Compute the permutation vector form a sparse matrix */ + /** Compute the permutation vector \a perm form the sparse matrix \a mat + * \warning The input sparse matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). + */ template void operator() (const MatrixType& mat, PermutationType& perm) { + eigen_assert(mat.isCompressed() && "COLAMDOrdering requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to COLAMDOrdering"); + Index m = mat.rows(); Index n = mat.cols(); Index nnz = mat.nonZeros(); @@ -132,12 +136,12 @@ class COLAMDOrdering Index stats [COLAMD_STATS]; internal::colamd_set_defaults(knobs); - Index info; IndexVector p(n+1), A(Alen); for(Index i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; for(Index i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; // Call Colamd routine to compute the ordering - info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); + Index info = internal::colamd(m, n, Alen, A.data(), p.data(), knobs, stats); + EIGEN_UNUSED_VARIABLE(info); eigen_assert( info && "COLAMD failed " ); perm.resize(n); diff --git a/uppsrc/plugin/Eigen/Eigen/src/PaStiXSupport/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/PaStiXSupport/CMakeLists.txt deleted file mode 100644 index 28c657e9b..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/PaStiXSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_PastixSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_PastixSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/PaStiXSupport COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/PardisoSupport/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/PardisoSupport/CMakeLists.txt deleted file mode 100644 index a097ab401..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/PardisoSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_PardisoSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_PardisoSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/PardisoSupport COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/PardisoSupport/PardisoSupport.h b/uppsrc/plugin/Eigen/Eigen/src/PardisoSupport/PardisoSupport.h index 1c48f0df7..18cd7d88a 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/PardisoSupport/PardisoSupport.h +++ b/uppsrc/plugin/Eigen/Eigen/src/PardisoSupport/PardisoSupport.h @@ -219,7 +219,7 @@ class PardisoImpl void pardisoInit(int type) { m_type = type; - bool symmetric = abs(m_type) < 10; + bool symmetric = std::abs(m_type) < 10; m_iparm[0] = 1; // No solver default m_iparm[1] = 3; // use Metis for the ordering m_iparm[2] = 1; // Numbers of processors, value of OMP_NUM_THREADS diff --git a/uppsrc/plugin/Eigen/Eigen/src/QR/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/QR/CMakeLists.txt deleted file mode 100644 index 96f43d7f5..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/QR/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_QR_SRCS "*.h") - -INSTALL(FILES - ${Eigen_QR_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/QR COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/QR/ColPivHouseholderQR.h b/uppsrc/plugin/Eigen/Eigen/src/QR/ColPivHouseholderQR.h index bec85810c..773d1df8f 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/QR/ColPivHouseholderQR.h +++ b/uppsrc/plugin/Eigen/Eigen/src/QR/ColPivHouseholderQR.h @@ -76,7 +76,8 @@ template class ColPivHouseholderQR m_colsTranspositions(), m_temp(), m_colSqNorms(), - m_isInitialized(false) {} + m_isInitialized(false), + m_usePrescribedThreshold(false) {} /** \brief Default Constructor with memory preallocation * diff --git a/uppsrc/plugin/Eigen/Eigen/src/SPQRSupport/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/SPQRSupport/CMakeLists.txt deleted file mode 100644 index 4968beaf2..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/SPQRSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SPQRSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SPQRSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SPQRSupport/ COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/SVD/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/SVD/CMakeLists.txt deleted file mode 100644 index 55efc44b1..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/SVD/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SVD_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SVD_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SVD COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/SVD/JacobiSVD.h b/uppsrc/plugin/Eigen/Eigen/src/SVD/JacobiSVD.h index f44995cd3..c57f2974c 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SVD/JacobiSVD.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SVD/JacobiSVD.h @@ -375,17 +375,19 @@ struct svd_precondition_2x2_block_to_be_real Scalar z; JacobiRotation rot; RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); + if(n==0) { z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); work_matrix.row(p) *= z; if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); if(work_matrix.coeff(q,q)!=Scalar(0)) + { z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); - else - z = Scalar(0); - work_matrix.row(q) *= z; - if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); + work_matrix.row(q) *= z; + if(svd.computeU()) svd.m_matrixU.col(q) *= conj(z); + } + // otherwise the second row is already zero, so we have nothing to do. } else { @@ -415,6 +417,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, JacobiRotation *j_right) { using std::sqrt; + using std::abs; Matrix m; m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), numext::real(matrix.coeff(q,p)), numext::real(matrix.coeff(q,q)); @@ -428,9 +431,11 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q, } else { - RealScalar u = d / t; - rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); - rot1.s() = rot1.c() * u; + RealScalar t2d2 = numext::hypot(t,d); + rot1.c() = abs(t)/t2d2; + rot1.s() = d/t2d2; + if(tmakeJacobi(m,0,1); @@ -531,8 +536,9 @@ template class JacobiSVD JacobiSVD() : m_isInitialized(false), m_isAllocated(false), + m_usePrescribedThreshold(false), m_computationOptions(0), - m_rows(-1), m_cols(-1) + m_rows(-1), m_cols(-1), m_diagSize(0) {} @@ -545,6 +551,7 @@ template class JacobiSVD JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) : m_isInitialized(false), m_isAllocated(false), + m_usePrescribedThreshold(false), m_computationOptions(0), m_rows(-1), m_cols(-1) { @@ -564,6 +571,7 @@ template class JacobiSVD JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) : m_isInitialized(false), m_isAllocated(false), + m_usePrescribedThreshold(false), m_computationOptions(0), m_rows(-1), m_cols(-1) { @@ -665,6 +673,69 @@ template class JacobiSVD eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); return m_nonzeroSingularValues; } + + /** \returns the rank of the matrix of which \c *this is the SVD. + * + * \note This method has to determine which singular values should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline Index rank() const + { + using std::abs; + eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); + if(m_singularValues.size()==0) return 0; + RealScalar premultiplied_threshold = m_singularValues.coeff(0) * threshold(); + Index i = m_nonzeroSingularValues-1; + while(i>=0 && m_singularValues.coeff(i) < premultiplied_threshold) --i; + return i+1; + } + + /** Allows to prescribe a threshold to be used by certain methods, such as rank() and solve(), + * which need to determine when singular values are to be considered nonzero. + * This is not used for the SVD decomposition itself. + * + * When it needs to get the threshold value, Eigen calls threshold(). + * The default is \c NumTraits::epsilon() + * + * \param threshold The new value to use as the threshold. + * + * A singular value will be considered nonzero if its value is strictly greater than + * \f$ \vert singular value \vert \leqslant threshold \times \vert max singular value \vert \f$. + * + * If you want to come back to the default behavior, call setThreshold(Default_t) + */ + JacobiSVD& setThreshold(const RealScalar& threshold) + { + m_usePrescribedThreshold = true; + m_prescribedThreshold = threshold; + return *this; + } + + /** Allows to come back to the default behavior, letting Eigen use its default formula for + * determining the threshold. + * + * You should pass the special object Eigen::Default as parameter here. + * \code svd.setThreshold(Eigen::Default); \endcode + * + * See the documentation of setThreshold(const RealScalar&). + */ + JacobiSVD& setThreshold(Default_t) + { + m_usePrescribedThreshold = false; + return *this; + } + + /** Returns the threshold that will be used by certain methods such as rank(). + * + * See the documentation of setThreshold(const RealScalar&). + */ + RealScalar threshold() const + { + eigen_assert(m_isInitialized || m_usePrescribedThreshold); + return m_usePrescribedThreshold ? m_prescribedThreshold + : (std::max)(1,m_diagSize)*NumTraits::epsilon(); + } inline Index rows() const { return m_rows; } inline Index cols() const { return m_cols; } @@ -677,11 +748,12 @@ template class JacobiSVD MatrixVType m_matrixV; SingularValuesType m_singularValues; WorkMatrixType m_workMatrix; - bool m_isInitialized, m_isAllocated; + bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold; bool m_computeFullU, m_computeThinU; bool m_computeFullV, m_computeThinV; unsigned int m_computationOptions; Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize; + RealScalar m_prescribedThreshold; template friend struct internal::svd_precondition_2x2_block_to_be_real; @@ -690,6 +762,7 @@ template class JacobiSVD internal::qr_preconditioner_impl m_qr_precond_morecols; internal::qr_preconditioner_impl m_qr_precond_morerows; + MatrixType m_scaledMatrix; }; template @@ -736,8 +809,9 @@ void JacobiSVD::allocate(Index rows, Index cols, u : 0); m_workMatrix.resize(m_diagSize, m_diagSize); - if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); - if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this); + if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this); + if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this); + if(m_cols!=m_cols) m_scaledMatrix.resize(rows,cols); } template @@ -754,11 +828,21 @@ JacobiSVD::compute(const MatrixType& matrix, unsig // limit for very small denormal numbers to be considered zero in order to avoid infinite loops (see bug 286) const RealScalar considerAsZero = RealScalar(2) * std::numeric_limits::denorm_min(); + // Scaling factor to reduce over/under-flows + RealScalar scale = matrix.cwiseAbs().maxCoeff(); + if(scale==RealScalar(0)) scale = RealScalar(1); + /*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ - if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix)) + if(m_rows!=m_cols) { - m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize); + m_scaledMatrix = matrix / scale; + m_qr_precond_morecols.run(*this, m_scaledMatrix); + m_qr_precond_morerows.run(*this, m_scaledMatrix); + } + else + { + m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale; if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize); if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols); @@ -784,7 +868,8 @@ JacobiSVD::compute(const MatrixType& matrix, unsig using std::max; RealScalar threshold = (max)(considerAsZero, precision * (max)(abs(m_workMatrix.coeff(p,p)), abs(m_workMatrix.coeff(q,q)))); - if((max)(abs(m_workMatrix.coeff(p,q)),abs(m_workMatrix.coeff(q,p))) > threshold) + // We compare both values to threshold instead of calling max to be robust to NaN (See bug 791) + if(abs(m_workMatrix.coeff(p,q))>threshold || abs(m_workMatrix.coeff(q,p)) > threshold) { finished = false; @@ -833,6 +918,8 @@ JacobiSVD::compute(const MatrixType& matrix, unsig if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i)); } } + + m_singularValues *= scale; m_isInitialized = true; return *this; @@ -854,11 +941,11 @@ struct solve_retval, Rhs> // So A^{-1} = V S^{-1} U^* Matrix tmp; - Index nonzeroSingVals = dec().nonzeroSingularValues(); + Index rank = dec().rank(); - tmp.noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs(); - tmp = dec().singularValues().head(nonzeroSingVals).asDiagonal().inverse() * tmp; - dst = dec().matrixV().leftCols(nonzeroSingVals) * tmp; + tmp.noalias() = dec().matrixU().leftCols(rank).adjoint() * rhs(); + tmp = dec().singularValues().head(rank).asDiagonal().inverse() * tmp; + dst = dec().matrixV().leftCols(rank) * tmp; } }; } // end namespace internal diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCholesky/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/SparseCholesky/CMakeLists.txt deleted file mode 100644 index 375a59d7a..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCholesky/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SparseCholesky_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SparseCholesky_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseCholesky COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCholesky/SimplicialCholesky.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCholesky/SimplicialCholesky.h index f41d7e010..e1f96ba5a 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCholesky/SimplicialCholesky.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCholesky/SimplicialCholesky.h @@ -37,6 +37,7 @@ class SimplicialCholeskyBase : internal::noncopyable { public: typedef typename internal::traits::MatrixType MatrixType; + typedef typename internal::traits::OrderingType OrderingType; enum { UpLo = internal::traits::UpLo }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; @@ -240,15 +241,16 @@ class SimplicialCholeskyBase : internal::noncopyable RealScalar m_shiftScale; }; -template class SimplicialLLT; -template class SimplicialLDLT; -template class SimplicialCholesky; +template > class SimplicialLLT; +template > class SimplicialLDLT; +template > class SimplicialCholesky; namespace internal { -template struct traits > +template struct traits > { typedef _MatrixType MatrixType; + typedef _Ordering OrderingType; enum { UpLo = _UpLo }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; @@ -259,9 +261,10 @@ template struct traits struct traits > +template struct traits > { typedef _MatrixType MatrixType; + typedef _Ordering OrderingType; enum { UpLo = _UpLo }; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Index Index; @@ -272,9 +275,10 @@ template struct traits struct traits > +template struct traits > { typedef _MatrixType MatrixType; + typedef _Ordering OrderingType; enum { UpLo = _UpLo }; }; @@ -294,11 +298,12 @@ template struct traits * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower * or Upper. Default is Lower. + * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> * - * \sa class SimplicialLDLT + * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering */ -template - class SimplicialLLT : public SimplicialCholeskyBase > +template + class SimplicialLLT : public SimplicialCholeskyBase > { public: typedef _MatrixType MatrixType; @@ -382,11 +387,12 @@ public: * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower * or Upper. Default is Lower. + * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> * - * \sa class SimplicialLLT + * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering */ -template - class SimplicialLDLT : public SimplicialCholeskyBase > +template + class SimplicialLDLT : public SimplicialCholeskyBase > { public: typedef _MatrixType MatrixType; @@ -467,8 +473,8 @@ public: * * \sa class SimplicialLDLT, class SimplicialLLT */ -template - class SimplicialCholesky : public SimplicialCholeskyBase > +template + class SimplicialCholesky : public SimplicialCholeskyBase > { public: typedef _MatrixType MatrixType; @@ -612,15 +618,13 @@ void SimplicialCholeskyBase::ordering(const MatrixType& a, CholMatrixTy { eigen_assert(a.rows()==a.cols()); const Index size = a.rows(); - // TODO allows to configure the permutation // Note that amd compute the inverse permutation { CholMatrixType C; C = a.template selfadjointView(); - // remove diagonal entries: - // seems not to be needed - // C.prune(keep_diag()); - internal::minimum_degree_ordering(C, m_Pinv); + + OrderingType ordering; + ordering(C,m_Pinv); } if(m_Pinv.size()>0) diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/AmbiVector.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/AmbiVector.h index 17fff96a7..220c6451c 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/AmbiVector.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/AmbiVector.h @@ -69,7 +69,7 @@ class AmbiVector delete[] m_buffer; if (size<1000) { - Index allocSize = (size * sizeof(ListEl))/sizeof(Scalar); + Index allocSize = (size * sizeof(ListEl) + sizeof(Scalar) - 1)/sizeof(Scalar); m_allocatedElements = (allocSize*sizeof(Scalar))/sizeof(ListEl); m_buffer = new Scalar[allocSize]; } @@ -88,7 +88,7 @@ class AmbiVector Index copyElements = m_allocatedElements; m_allocatedElements = (std::min)(Index(m_allocatedElements*1.5),m_size); Index allocSize = m_allocatedElements * sizeof(ListEl); - allocSize = allocSize/sizeof(Scalar) + (allocSize%sizeof(Scalar)>0?1:0); + allocSize = (allocSize + sizeof(Scalar) - 1)/sizeof(Scalar); Scalar* newBuffer = new Scalar[allocSize]; memcpy(newBuffer, m_buffer, copyElements * sizeof(ListEl)); delete[] m_buffer; diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/CMakeLists.txt deleted file mode 100644 index d860452a6..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SparseCore_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SparseCore_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseCore COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/CompressedStorage.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/CompressedStorage.h index 3321fab4a..a667cb56e 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/CompressedStorage.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/CompressedStorage.h @@ -51,8 +51,8 @@ class CompressedStorage CompressedStorage& operator=(const CompressedStorage& other) { resize(other.size()); - memcpy(m_values, other.m_values, m_size * sizeof(Scalar)); - memcpy(m_indices, other.m_indices, m_size * sizeof(Index)); + internal::smart_copy(other.m_values, other.m_values + m_size, m_values); + internal::smart_copy(other.m_indices, other.m_indices + m_size, m_indices); return *this; } @@ -83,10 +83,10 @@ class CompressedStorage reallocate(m_size); } - void resize(size_t size, float reserveSizeFactor = 0) + void resize(size_t size, double reserveSizeFactor = 0) { if (m_allocatedSize m_outerSize; EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl) + private: + Index nonZeros() const; }; @@ -82,6 +84,7 @@ class BlockImpl,BlockRows,BlockCols,true typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; typedef typename internal::remove_all::type _MatrixTypeNested; typedef Block BlockType; + typedef Block ConstBlockType; public: enum { IsRowMajor = internal::traits::IsRowMajor }; EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) @@ -245,6 +248,93 @@ public: }; + +template +class BlockImpl,BlockRows,BlockCols,true,Sparse> + : public SparseMatrixBase,BlockRows,BlockCols,true> > +{ + typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType; + typedef typename internal::remove_all::type _MatrixTypeNested; + typedef Block BlockType; +public: + enum { IsRowMajor = internal::traits::IsRowMajor }; + EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType) +protected: + enum { OuterSize = IsRowMajor ? BlockRows : BlockCols }; +public: + + class InnerIterator: public SparseMatrixType::InnerIterator + { + public: + inline InnerIterator(const BlockType& xpr, Index outer) + : SparseMatrixType::InnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + class ReverseInnerIterator: public SparseMatrixType::ReverseInnerIterator + { + public: + inline ReverseInnerIterator(const BlockType& xpr, Index outer) + : SparseMatrixType::ReverseInnerIterator(xpr.m_matrix, xpr.m_outerStart + outer), m_outer(outer) + {} + inline Index row() const { return IsRowMajor ? m_outer : this->index(); } + inline Index col() const { return IsRowMajor ? this->index() : m_outer; } + protected: + Index m_outer; + }; + + inline BlockImpl(const SparseMatrixType& xpr, int i) + : m_matrix(xpr), m_outerStart(i), m_outerSize(OuterSize) + {} + + inline BlockImpl(const SparseMatrixType& xpr, int startRow, int startCol, int blockRows, int blockCols) + : m_matrix(xpr), m_outerStart(IsRowMajor ? startRow : startCol), m_outerSize(IsRowMajor ? blockRows : blockCols) + {} + + inline const Scalar* valuePtr() const + { return m_matrix.valuePtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + + inline const Index* innerIndexPtr() const + { return m_matrix.innerIndexPtr() + m_matrix.outerIndexPtr()[m_outerStart]; } + + inline const Index* outerIndexPtr() const + { return m_matrix.outerIndexPtr() + m_outerStart; } + + Index nonZeros() const + { + if(m_matrix.isCompressed()) + return std::size_t(m_matrix.outerIndexPtr()[m_outerStart+m_outerSize.value()]) + - std::size_t(m_matrix.outerIndexPtr()[m_outerStart]); + else if(m_outerSize.value()==0) + return 0; + else + return Map >(m_matrix.innerNonZeroPtr()+m_outerStart, m_outerSize.value()).sum(); + } + + const Scalar& lastCoeff() const + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(BlockImpl); + eigen_assert(nonZeros()>0); + if(m_matrix.isCompressed()) + return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart+1]-1]; + else + return m_matrix.valuePtr()[m_matrix.outerIndexPtr()[m_outerStart]+m_matrix.innerNonZeroPtr()[m_outerStart]-1]; + } + + EIGEN_STRONG_INLINE Index rows() const { return IsRowMajor ? m_outerSize.value() : m_matrix.rows(); } + EIGEN_STRONG_INLINE Index cols() const { return IsRowMajor ? m_matrix.cols() : m_outerSize.value(); } + + protected: + + typename SparseMatrixType::Nested m_matrix; + Index m_outerStart; + const internal::variable_if_dynamic m_outerSize; + +}; + //---------- /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h index ec86ca933..60ca7690c 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseCwiseBinaryOp.h @@ -73,7 +73,8 @@ class CwiseBinaryOpImpl::InnerIterator typedef internal::sparse_cwise_binary_op_inner_iterator_selector< BinaryOp,Lhs,Rhs, InnerIterator> Base; - EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOpImpl& binOp, Index outer) + // NOTE: we have to prefix Index by "typename Lhs::" to avoid an ICE with VC11 + EIGEN_STRONG_INLINE InnerIterator(const CwiseBinaryOpImpl& binOp, typename Lhs::Index outer) : Base(binOp.derived(),outer) {} }; diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseDenseProduct.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseDenseProduct.h index 54fd633a1..a9084192e 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseDenseProduct.h @@ -19,7 +19,10 @@ template struct SparseDenseProductRet template struct SparseDenseProductReturnType { - typedef SparseDenseOuterProduct Type; + typedef typename internal::conditional< + Lhs::IsRowMajor, + SparseDenseOuterProduct, + SparseDenseOuterProduct >::type Type; }; template struct DenseSparseProductReturnType @@ -29,7 +32,10 @@ template struct DenseSparseProductRet template struct DenseSparseProductReturnType { - typedef SparseDenseOuterProduct Type; + typedef typename internal::conditional< + Rhs::IsRowMajor, + SparseDenseOuterProduct, + SparseDenseOuterProduct >::type Type; }; namespace internal { @@ -114,17 +120,30 @@ class SparseDenseOuterProduct::InnerIterator : public _LhsNes typedef typename SparseDenseOuterProduct::Index Index; public: EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer) - : Base(prod.lhs(), 0), m_outer(outer), m_factor(prod.rhs().coeff(outer)) - { - } + : Base(prod.lhs(), 0), m_outer(outer), m_factor(get(prod.rhs(), outer, typename internal::traits::StorageKind() )) + { } inline Index outer() const { return m_outer; } - inline Index row() const { return Transpose ? Base::row() : m_outer; } - inline Index col() const { return Transpose ? m_outer : Base::row(); } + inline Index row() const { return Transpose ? m_outer : Base::index(); } + inline Index col() const { return Transpose ? Base::index() : m_outer; } inline Scalar value() const { return Base::value() * m_factor; } protected: + static Scalar get(const _RhsNested &rhs, Index outer, Dense = Dense()) + { + return rhs.coeff(outer); + } + + static Scalar get(const _RhsNested &rhs, Index outer, Sparse = Sparse()) + { + typename Traits::_RhsNested::InnerIterator it(rhs, outer); + if (it && it.index()==0) + return it.value(); + + return Scalar(0); + } + Index m_outer; Scalar m_factor; }; @@ -287,15 +306,6 @@ class DenseTimeSparseProduct DenseTimeSparseProduct& operator=(const DenseTimeSparseProduct&); }; -// sparse * dense -template -template -inline const typename SparseDenseProductReturnType::Type -SparseMatrixBase::operator*(const MatrixBase &other) const -{ - return typename SparseDenseProductReturnType::Type(derived(), other.derived()); -} - } // end namespace Eigen #endif // EIGEN_SPARSEDENSEPRODUCT_H diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseMatrix.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseMatrix.h index 01ce0dcfe..ba5e3a9b6 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseMatrix.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseMatrix.h @@ -940,7 +940,7 @@ void set_from_triplets(const InputIterator& begin, const InputIterator& end, Spa enum { IsRowMajor = SparseMatrixType::IsRowMajor }; typedef typename SparseMatrixType::Scalar Scalar; typedef typename SparseMatrixType::Index Index; - SparseMatrix trMat(mat.rows(),mat.cols()); + SparseMatrix trMat(mat.rows(),mat.cols()); if(begin!=end) { @@ -1178,7 +1178,7 @@ EIGEN_DONT_INLINE typename SparseMatrix<_Scalar,_Options,_Index>::Scalar& Sparse size_t p = m_outerIndex[outer+1]; ++m_outerIndex[outer+1]; - float reallocRatio = 1; + double reallocRatio = 1; if (m_data.allocatedSize()<=m_data.size()) { // if there is no preallocated memory, let's reserve a minimum of 32 elements @@ -1190,13 +1190,13 @@ EIGEN_DONT_INLINE typename SparseMatrix<_Scalar,_Options,_Index>::Scalar& Sparse { // we need to reallocate the data, to reduce multiple reallocations // we use a smart resize algorithm based on the current filling ratio - // in addition, we use float to avoid integers overflows - float nnzEstimate = float(m_outerIndex[outer])*float(m_outerSize)/float(outer+1); - reallocRatio = (nnzEstimate-float(m_data.size()))/float(m_data.size()); + // in addition, we use double to avoid integers overflows + double nnzEstimate = double(m_outerIndex[outer])*double(m_outerSize)/double(outer+1); + reallocRatio = (nnzEstimate-double(m_data.size()))/double(m_data.size()); // furthermore we bound the realloc ratio to: // 1) reduce multiple minor realloc when the matrix is almost filled // 2) avoid to allocate too much memory when the matrix is almost empty - reallocRatio = (std::min)((std::max)(reallocRatio,1.5f),8.f); + reallocRatio = (std::min)((std::max)(reallocRatio,1.5),8.); } } m_data.resize(m_data.size()+1,reallocRatio); diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseMatrixBase.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseMatrixBase.h index bbcf7fb1c..485e85bec 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseMatrixBase.h @@ -358,7 +358,8 @@ template class SparseMatrixBase : public EigenBase /** sparse * dense (returns a dense object unless it is an outer product) */ template const typename SparseDenseProductReturnType::Type - operator*(const MatrixBase &other) const; + operator*(const MatrixBase &other) const + { return typename SparseDenseProductReturnType::Type(derived(), other.derived()); } /** \returns an expression of P H P^-1 where H is the matrix represented by \c *this */ SparseSymmetricPermutationProduct twistedBy(const PermutationMatrix& perm) const diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparsePermutation.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparsePermutation.h index b85be93f6..75e210009 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparsePermutation.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparsePermutation.h @@ -61,7 +61,7 @@ struct permut_sparsematrix_product_retval for(Index j=0; j class TransposeImpl inline Index nonZeros() const { return derived().nestedExpression().nonZeros(); } }; -// NOTE: VC10 trigger an ICE if don't put typename TransposeImpl:: in front of Index, +// NOTE: VC10 and VC11 trigger an ICE if don't put typename TransposeImpl:: in front of Index, // a typedef typename TransposeImpl::Index Index; // does not fix the issue. // An alternative is to define the nested class in the parent class itself. @@ -40,8 +40,8 @@ template class TransposeImpl::InnerItera EIGEN_STRONG_INLINE InnerIterator(const TransposeImpl& trans, typename TransposeImpl::Index outer) : Base(trans.derived().nestedExpression(), outer) {} - Index row() const { return Base::col(); } - Index col() const { return Base::row(); } + typename TransposeImpl::Index row() const { return Base::col(); } + typename TransposeImpl::Index col() const { return Base::row(); } }; template class TransposeImpl::ReverseInnerIterator @@ -54,8 +54,8 @@ template class TransposeImpl::ReverseInn EIGEN_STRONG_INLINE ReverseInnerIterator(const TransposeImpl& xpr, typename TransposeImpl::Index outer) : Base(xpr.derived().nestedExpression(), outer) {} - Index row() const { return Base::col(); } - Index col() const { return Base::row(); } + typename TransposeImpl::Index row() const { return Base::col(); } + typename TransposeImpl::Index col() const { return Base::row(); } }; } // end namespace Eigen diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseUtil.h b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseUtil.h index 05023858b..0ba471320 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseUtil.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseCore/SparseUtil.h @@ -84,8 +84,10 @@ template class DenseTimeSparseProduct; template class SparseDenseOuterProduct; template struct SparseSparseProductReturnType; -template::ColsAtCompileTime> struct DenseSparseProductReturnType; -template::ColsAtCompileTime> struct SparseDenseProductReturnType; +template::ColsAtCompileTime,internal::traits::RowsAtCompileTime)> struct DenseSparseProductReturnType; +template::ColsAtCompileTime,internal::traits::RowsAtCompileTime)> struct SparseDenseProductReturnType; template class SparseSymmetricPermutationProduct; namespace internal { diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseLU/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/SparseLU/CMakeLists.txt deleted file mode 100644 index 69729ee89..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseLU/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SparseLU_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SparseLU_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseLU COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseLU/SparseLU.h b/uppsrc/plugin/Eigen/Eigen/src/SparseLU/SparseLU.h index 1d592f2c8..d6f5d59ca 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseLU/SparseLU.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseLU/SparseLU.h @@ -260,14 +260,13 @@ class SparseLU : public internal::SparseLUImplcols(); ++j) { for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) { - if(it.row() < j) continue; - if(it.row() == j) + if(it.index() == j) { det *= (std::abs)(it.value()); break; diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/uppsrc/plugin/Eigen/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h index ad6f2183f..b16afd6a4 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h @@ -189,8 +189,8 @@ class MappedSuperNodalMatrix::InnerIterator m_idval(mat.colIndexPtr()[outer]), m_startidval(m_idval), m_endidval(mat.colIndexPtr()[outer+1]), - m_idrow(mat.rowIndexPtr()[outer]), - m_endidrow(mat.rowIndexPtr()[outer+1]) + m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]), + m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1]) {} inline InnerIterator& operator++() { diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseQR/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/SparseQR/CMakeLists.txt deleted file mode 100644 index f9ddf2bdb..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseQR/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SparseQR_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SparseQR_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseQR/ COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/SparseQR/SparseQR.h b/uppsrc/plugin/Eigen/Eigen/src/SparseQR/SparseQR.h index afda43bfc..a00bd5db1 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/SparseQR/SparseQR.h +++ b/uppsrc/plugin/Eigen/Eigen/src/SparseQR/SparseQR.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2012-2013 Desire Nuentsa -// Copyright (C) 2012-2013 Gael Guennebaud +// Copyright (C) 2012-2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -58,6 +58,7 @@ namespace internal { * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module * OrderingMethods \endlink module for the list of built-in and external ordering methods. * + * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()). * */ template @@ -74,13 +75,26 @@ class SparseQR typedef Matrix ScalarVector; typedef PermutationMatrix PermutationType; public: - SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) + SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { } - SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) + /** Construct a QR factorization of the matrix \a mat. + * + * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). + * + * \sa compute() + */ + SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) { compute(mat); } + + /** Computes the QR factorization of the sparse matrix \a mat. + * + * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). + * + * \sa analyzePattern(), factorize() + */ void compute(const MatrixType& mat) { analyzePattern(mat); @@ -166,7 +180,7 @@ class SparseQR y.bottomRows(y.rows()-rank).setZero(); // Apply the column permutation - if (m_perm_c.size()) dest.topRows(cols()) = colsPermutation() * y.topRows(cols()); + if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols()); else dest = y.topRows(cols()); m_info = Success; @@ -206,7 +220,7 @@ class SparseQR /** \brief Reports whether previous computation was successful. * - * \returns \c Success if computation was succesful, + * \returns \c Success if computation was successful, * \c NumericalIssue if the QR factorization reports a numerical problem * \c InvalidInput if the input matrix is invalid * @@ -248,6 +262,7 @@ class SparseQR IndexVector m_etree; // Column elimination tree IndexVector m_firstRowElt; // First element in each row bool m_isQSorted; // whether Q is sorted or not + bool m_isEtreeOk; // whether the elimination tree match the initial input matrix template friend struct SparseQR_QProduct; template friend struct SparseQRMatrixQReturnType; @@ -255,20 +270,26 @@ class SparseQR }; /** \brief Preprocessing step of a QR factorization + * + * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). * * In this step, the fill-reducing permutation is computed and applied to the columns of A - * and the column elimination tree is computed as well. Only the sparcity pattern of \a mat is exploited. + * and the column elimination tree is computed as well. Only the sparsity pattern of \a mat is exploited. * * \note In this step it is assumed that there is no empty row in the matrix \a mat. */ template void SparseQR::analyzePattern(const MatrixType& mat) { + eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR"); + // Copy to a column major matrix if the input is rowmajor + typename internal::conditional::type matCpy(mat); // Compute the column fill reducing ordering OrderingType ord; - ord(mat, m_perm_c); + ord(matCpy, m_perm_c); Index n = mat.cols(); Index m = mat.rows(); + Index diagSize = (std::min)(m,n); if (!m_perm_c.size()) { @@ -278,22 +299,23 @@ void SparseQR::analyzePattern(const MatrixType& mat) // Compute the column elimination tree of the permuted matrix m_outputPerm_c = m_perm_c.inverse(); - internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); + internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); + m_isEtreeOk = true; - m_R.resize(n, n); - m_Q.resize(m, n); + m_R.resize(m, n); + m_Q.resize(m, diagSize); // Allocate space for nonzero elements : rough estimation m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree m_Q.reserve(2*mat.nonZeros()); - m_hcoeffs.resize(n); + m_hcoeffs.resize(diagSize); m_analysisIsok = true; } /** \brief Performs the numerical QR factorization of the input matrix * * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with - * a matrix having the same sparcity pattern than \a mat. + * a matrix having the same sparsity pattern than \a mat. * * \param mat The sparse column-major matrix */ @@ -306,23 +328,47 @@ void SparseQR::factorize(const MatrixType& mat) eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); Index m = mat.rows(); Index n = mat.cols(); - IndexVector mark(m); mark.setConstant(-1); // Record the visited nodes - IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q - Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q - ScalarVector tval(m); // The dense vector used to compute the current column - bool found_diag; - + Index diagSize = (std::min)(m,n); + IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes + IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q + Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q + ScalarVector tval(m); // The dense vector used to compute the current column + RealScalar pivotThreshold = m_threshold; + + m_R.setZero(); + m_Q.setZero(); m_pmat = mat; - m_pmat.uncompress(); // To have the innerNonZeroPtr allocated - // Apply the fill-in reducing permutation lazily: - for (int i = 0; i < n; i++) + if(!m_isEtreeOk) { - Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; - m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; - m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; + m_outputPerm_c = m_perm_c.inverse(); + internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); + m_isEtreeOk = true; + } + + m_pmat.uncompress(); // To have the innerNonZeroPtr allocated + + // Apply the fill-in reducing permutation lazily: + { + // If the input is row major, copy the original column indices, + // otherwise directly use the input matrix + // + IndexVector originalOuterIndicesCpy; + const Index *originalOuterIndices = mat.outerIndexPtr(); + if(MatrixType::IsRowMajor) + { + originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1); + originalOuterIndices = originalOuterIndicesCpy.data(); + } + + for (int i = 0; i < n; i++) + { + Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; + m_pmat.outerIndexPtr()[p] = originalOuterIndices[i]; + m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i]; + } } - /* Compute the default threshold, see : + /* Compute the default threshold as in MatLab, see: * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 */ @@ -330,33 +376,35 @@ void SparseQR::factorize(const MatrixType& mat) { RealScalar max2Norm = 0.0; for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm()); - m_threshold = 20 * (m + n) * max2Norm * NumTraits::epsilon(); + if(max2Norm==RealScalar(0)) + max2Norm = RealScalar(1); + pivotThreshold = 20 * (m + n) * max2Norm * NumTraits::epsilon(); } // Initialize the numerical permutation m_pivotperm.setIdentity(n); Index nonzeroCol = 0; // Record the number of valid pivots - + m_Q.startVec(0); + // Left looking rank-revealing QR factorization: compute a column of R and Q at a time - for (Index col = 0; col < (std::min)(n,m); ++col) + for (Index col = 0; col < n; ++col) { mark.setConstant(-1); m_R.startVec(col); - m_Q.startVec(col); mark(nonzeroCol) = col; Qidx(0) = nonzeroCol; nzcolR = 0; nzcolQ = 1; - found_diag = col>=m; + bool found_diag = nonzeroCol>=m; tval.setZero(); // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e., // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k. // Note: if the diagonal entry does not exist, then its contribution must be explicitly added, // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. - for (typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) + for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) { - Index curIdx = nonzeroCol ; + Index curIdx = nonzeroCol; if(itp) curIdx = itp.row(); if(curIdx == nonzeroCol) found_diag = true; @@ -398,7 +446,7 @@ void SparseQR::factorize(const MatrixType& mat) // Browse all the indexes of R(:,col) in reverse order for (Index i = nzcolR-1; i >= 0; i--) { - Index curIdx = m_pivotperm.indices()(Ridx(i)); + Index curIdx = Ridx(i); // Apply the curIdx-th householder vector to the current column (temporarily stored into tval) Scalar tdot(0); @@ -427,33 +475,36 @@ void SparseQR::factorize(const MatrixType& mat) } } } // End update current column - - // Compute the Householder reflection that eliminate the current column - // FIXME this step should call the Householder module. - Scalar tau; - RealScalar beta; - Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0); - // First, the squared norm of Q((col+1):m, col) - RealScalar sqrNorm = 0.; - for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); + Scalar tau = 0; + RealScalar beta = 0; - if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) + if(nonzeroCol < diagSize) { - tau = RealScalar(0); - beta = numext::real(c0); - tval(Qidx(0)) = 1; - } - else - { - beta = std::sqrt(numext::abs2(c0) + sqrNorm); - if(numext::real(c0) >= RealScalar(0)) - beta = -beta; - tval(Qidx(0)) = 1; - for (Index itq = 1; itq < nzcolQ; ++itq) - tval(Qidx(itq)) /= (c0 - beta); - tau = numext::conj((beta-c0) / beta); - + // Compute the Householder reflection that eliminate the current column + // FIXME this step should call the Householder module. + Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0); + + // First, the squared norm of Q((col+1):m, col) + RealScalar sqrNorm = 0.; + for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); + if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) + { + beta = numext::real(c0); + tval(Qidx(0)) = 1; + } + else + { + using std::sqrt; + beta = sqrt(numext::abs2(c0) + sqrNorm); + if(numext::real(c0) >= RealScalar(0)) + beta = -beta; + tval(Qidx(0)) = 1; + for (Index itq = 1; itq < nzcolQ; ++itq) + tval(Qidx(itq)) /= (c0 - beta); + tau = numext::conj((beta-c0) / beta); + + } } // Insert values in R @@ -467,45 +518,49 @@ void SparseQR::factorize(const MatrixType& mat) } } - if(abs(beta) >= m_threshold) + if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold) { m_R.insertBackByOuterInner(col, nonzeroCol) = beta; - nonzeroCol++; // The householder coefficient - m_hcoeffs(col) = tau; + m_hcoeffs(nonzeroCol) = tau; // Record the householder reflections for (Index itq = 0; itq < nzcolQ; ++itq) { Index iQ = Qidx(itq); - m_Q.insertBackByOuterInnerUnordered(col,iQ) = tval(iQ); + m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ); tval(iQ) = Scalar(0.); - } + } + nonzeroCol++; + if(nonzeroCol void evalTo(DesType& res) const { + Index m = m_qr.rows(); Index n = m_qr.cols(); + Index diagSize = (std::min)(m,n); res = m_other; if (m_transpose) { eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); //Compute res = Q' * other column by column for(Index j = 0; j < res.cols(); j++){ - for (Index k = 0; k < n; k++) + for (Index k = 0; k < diagSize; k++) { Scalar tau = Scalar(0); tau = m_qr.m_Q.col(k).dot(res.col(j)); @@ -581,10 +638,10 @@ struct SparseQR_QProduct : ReturnByValue=0; k--) + for (Index k = diagSize-1; k >=0; k--) { Scalar tau = Scalar(0); tau = m_qr.m_Q.col(k).dot(res.col(j)); @@ -618,7 +675,7 @@ struct SparseQRMatrixQReturnType : public EigenBase(m_qr); } inline Index rows() const { return m_qr.rows(); } - inline Index cols() const { return m_qr.cols(); } + inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); } // To use for operations with the transpose of Q SparseQRMatrixQTransposeReturnType transpose() const { diff --git a/uppsrc/plugin/Eigen/Eigen/src/StlSupport/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/StlSupport/CMakeLists.txt deleted file mode 100644 index 0f094f637..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/StlSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_StlSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_StlSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/StlSupport COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdDeque.h b/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdDeque.h index 4ee8e5c10..aaf66330b 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdDeque.h +++ b/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdDeque.h @@ -11,7 +11,7 @@ #ifndef EIGEN_STDDEQUE_H #define EIGEN_STDDEQUE_H -#include "Eigen/src/StlSupport/details.h" +#include "details.h" // Define the explicit instantiation (e.g. necessary for the Intel compiler) #if defined(__INTEL_COMPILER) || defined(__GNUC__) diff --git a/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdList.h b/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdList.h index 627381ece..3c742430c 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdList.h +++ b/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdList.h @@ -10,7 +10,7 @@ #ifndef EIGEN_STDLIST_H #define EIGEN_STDLIST_H -#include "Eigen/src/StlSupport/details.h" +#include "details.h" // Define the explicit instantiation (e.g. necessary for the Intel compiler) #if defined(__INTEL_COMPILER) || defined(__GNUC__) diff --git a/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdVector.h b/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdVector.h index 40a9abefa..611664a2e 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdVector.h +++ b/uppsrc/plugin/Eigen/Eigen/src/StlSupport/StdVector.h @@ -11,7 +11,7 @@ #ifndef EIGEN_STDVECTOR_H #define EIGEN_STDVECTOR_H -#include "Eigen/src/StlSupport/details.h" +#include "details.h" /** * This section contains a convenience MACRO which allows an easy specialization of diff --git a/uppsrc/plugin/Eigen/Eigen/src/SuperLUSupport/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/SuperLUSupport/CMakeLists.txt deleted file mode 100644 index b28ebe583..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/SuperLUSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SuperLUSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SuperLUSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SuperLUSupport COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/UmfPackSupport/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/UmfPackSupport/CMakeLists.txt deleted file mode 100644 index a57de0020..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/UmfPackSupport/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_UmfPackSupport_SRCS "*.h") - -INSTALL(FILES - ${Eigen_UmfPackSupport_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/UmfPackSupport COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/UmfPackSupport/UmfPackSupport.h b/uppsrc/plugin/Eigen/Eigen/src/UmfPackSupport/UmfPackSupport.h index 3a48cecf7..29c60c378 100644 --- a/uppsrc/plugin/Eigen/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/uppsrc/plugin/Eigen/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -107,6 +107,16 @@ inline int umfpack_get_determinant(std::complex *Mx, double *Ex, void *N return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info); } +namespace internal { + template struct umfpack_helper_is_sparse_plain : false_type {}; + template + struct umfpack_helper_is_sparse_plain > + : true_type {}; + template + struct umfpack_helper_is_sparse_plain > + : true_type {}; +} + /** \ingroup UmfPackSupport_Module * \brief A sparse LU factorization and solver based on UmfPack * @@ -192,10 +202,14 @@ class UmfPackLU : internal::noncopyable * Note that the matrix should be column-major, and in compressed format for best performance. * \sa SparseMatrix::makeCompressed(). */ - void compute(const MatrixType& matrix) + template + void compute(const InputMatrixType& matrix) { - analyzePattern(matrix); - factorize(matrix); + if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); + if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); + grapInput(matrix.derived()); + analyzePattern_impl(); + factorize_impl(); } /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. @@ -230,23 +244,15 @@ class UmfPackLU : internal::noncopyable * * \sa factorize(), compute() */ - void analyzePattern(const MatrixType& matrix) + template + void analyzePattern(const InputMatrixType& matrix) { - if(m_symbolic) - umfpack_free_symbolic(&m_symbolic,Scalar()); - if(m_numeric) - umfpack_free_numeric(&m_numeric,Scalar()); + if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); + if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); - grapInput(matrix); + grapInput(matrix.derived()); - int errorCode = 0; - errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, - &m_symbolic, 0, 0); - - m_isInitialized = true; - m_info = errorCode ? InvalidInput : Success; - m_analysisIsOk = true; - m_factorizationIsOk = false; + analyzePattern_impl(); } /** Performs a numeric decomposition of \a matrix @@ -255,20 +261,16 @@ class UmfPackLU : internal::noncopyable * * \sa analyzePattern(), compute() */ - void factorize(const MatrixType& matrix) + template + void factorize(const InputMatrixType& matrix) { eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()"); if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); - grapInput(matrix); - - int errorCode; - errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, - m_symbolic, &m_numeric, 0, 0); - - m_info = errorCode ? NumericalIssue : Success; - m_factorizationIsOk = true; + grapInput(matrix.derived()); + + factorize_impl(); } #ifndef EIGEN_PARSED_BY_DOXYGEN @@ -283,19 +285,20 @@ class UmfPackLU : internal::noncopyable protected: - void init() { - m_info = InvalidInput; - m_isInitialized = false; - m_numeric = 0; - m_symbolic = 0; - m_outerIndexPtr = 0; - m_innerIndexPtr = 0; - m_valuePtr = 0; + m_info = InvalidInput; + m_isInitialized = false; + m_numeric = 0; + m_symbolic = 0; + m_outerIndexPtr = 0; + m_innerIndexPtr = 0; + m_valuePtr = 0; + m_extractedDataAreDirty = true; } - void grapInput(const MatrixType& mat) + template + void grapInput_impl(const InputMatrixType& mat, internal::true_type) { m_copyMatrix.resize(mat.rows(), mat.cols()); if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() ) @@ -313,6 +316,45 @@ class UmfPackLU : internal::noncopyable m_valuePtr = mat.valuePtr(); } } + + template + void grapInput_impl(const InputMatrixType& mat, internal::false_type) + { + m_copyMatrix = mat; + m_outerIndexPtr = m_copyMatrix.outerIndexPtr(); + m_innerIndexPtr = m_copyMatrix.innerIndexPtr(); + m_valuePtr = m_copyMatrix.valuePtr(); + } + + template + void grapInput(const InputMatrixType& mat) + { + grapInput_impl(mat, internal::umfpack_helper_is_sparse_plain()); + } + + void analyzePattern_impl() + { + int errorCode = 0; + errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, + &m_symbolic, 0, 0); + + m_isInitialized = true; + m_info = errorCode ? InvalidInput : Success; + m_analysisIsOk = true; + m_factorizationIsOk = false; + m_extractedDataAreDirty = true; + } + + void factorize_impl() + { + int errorCode; + errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, + m_symbolic, &m_numeric, 0, 0); + + m_info = errorCode ? NumericalIssue : Success; + m_factorizationIsOk = true; + m_extractedDataAreDirty = true; + } // cached data to reduce reallocation, etc. mutable LUMatrixType m_l; diff --git a/uppsrc/plugin/Eigen/Eigen/src/misc/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/misc/CMakeLists.txt deleted file mode 100644 index a58ffb745..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/misc/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_misc_SRCS "*.h") - -INSTALL(FILES - ${Eigen_misc_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/misc COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/Eigen/src/plugins/CMakeLists.txt b/uppsrc/plugin/Eigen/Eigen/src/plugins/CMakeLists.txt deleted file mode 100644 index 1a1d3ffbd..000000000 --- a/uppsrc/plugin/Eigen/Eigen/src/plugins/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_plugins_SRCS "*.h") - -INSTALL(FILES - ${Eigen_plugins_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/plugins COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/CMakeLists.txt deleted file mode 100644 index e06f1238b..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/CMakeLists.txt +++ /dev/null @@ -1,11 +0,0 @@ -set(Eigen_HEADERS AdolcForward BVH IterativeSolvers MatrixFunctions MoreVectorization AutoDiff AlignedVector3 Polynomials - FFT NonLinearOptimization SparseExtra IterativeSolvers - NumericalDiff Skyline MPRealSupport OpenGLSupport KroneckerProduct Splines LevenbergMarquardt - ) - -install(FILES - ${Eigen_HEADERS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen COMPONENT Devel - ) - -add_subdirectory(src) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/OpenGLSupport b/uppsrc/plugin/Eigen/unsupported/Eigen/OpenGLSupport index c4090ab11..e2769449c 100644 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/OpenGLSupport +++ b/uppsrc/plugin/Eigen/unsupported/Eigen/OpenGLSupport @@ -178,11 +178,11 @@ template void glLoadMatrix(const Transform& t) template void glLoadMatrix(const Transform& t) { glLoadMatrix(t.matrix()); } template void glLoadMatrix(const Transform& t) { glLoadMatrix(Transform(t).matrix()); } -static void glRotate(const Rotation2D& rot) +inline void glRotate(const Rotation2D& rot) { glRotatef(rot.angle()*180.f/float(M_PI), 0.f, 0.f, 1.f); } -static void glRotate(const Rotation2D& rot) +inline void glRotate(const Rotation2D& rot) { glRotated(rot.angle()*180.0/M_PI, 0.0, 0.0, 1.0); } @@ -246,18 +246,18 @@ EIGEN_GL_FUNC1_SPECIALIZATION_MAT(glGet,GLenum,_,double, 4,4,Doublev) #ifdef GL_VERSION_2_0 -static void glUniform2fv_ei (GLint loc, const float* v) { glUniform2fv(loc,1,v); } -static void glUniform2iv_ei (GLint loc, const int* v) { glUniform2iv(loc,1,v); } +inline void glUniform2fv_ei (GLint loc, const float* v) { glUniform2fv(loc,1,v); } +inline void glUniform2iv_ei (GLint loc, const int* v) { glUniform2iv(loc,1,v); } -static void glUniform3fv_ei (GLint loc, const float* v) { glUniform3fv(loc,1,v); } -static void glUniform3iv_ei (GLint loc, const int* v) { glUniform3iv(loc,1,v); } +inline void glUniform3fv_ei (GLint loc, const float* v) { glUniform3fv(loc,1,v); } +inline void glUniform3iv_ei (GLint loc, const int* v) { glUniform3iv(loc,1,v); } -static void glUniform4fv_ei (GLint loc, const float* v) { glUniform4fv(loc,1,v); } -static void glUniform4iv_ei (GLint loc, const int* v) { glUniform4iv(loc,1,v); } +inline void glUniform4fv_ei (GLint loc, const float* v) { glUniform4fv(loc,1,v); } +inline void glUniform4iv_ei (GLint loc, const int* v) { glUniform4iv(loc,1,v); } -static void glUniformMatrix2fv_ei (GLint loc, const float* v) { glUniformMatrix2fv(loc,1,false,v); } -static void glUniformMatrix3fv_ei (GLint loc, const float* v) { glUniformMatrix3fv(loc,1,false,v); } -static void glUniformMatrix4fv_ei (GLint loc, const float* v) { glUniformMatrix4fv(loc,1,false,v); } +inline void glUniformMatrix2fv_ei (GLint loc, const float* v) { glUniformMatrix2fv(loc,1,false,v); } +inline void glUniformMatrix3fv_ei (GLint loc, const float* v) { glUniformMatrix3fv(loc,1,false,v); } +inline void glUniformMatrix4fv_ei (GLint loc, const float* v) { glUniformMatrix4fv(loc,1,false,v); } EIGEN_GL_FUNC1_DECLARATION (glUniform,GLint,const) @@ -294,9 +294,9 @@ EIGEN_GL_FUNC1_SPECIALIZATION_MAT(glUniform,GLint,const,float, 4,3,Matrix #ifdef GL_VERSION_3_0 -static void glUniform2uiv_ei (GLint loc, const unsigned int* v) { glUniform2uiv(loc,1,v); } -static void glUniform3uiv_ei (GLint loc, const unsigned int* v) { glUniform3uiv(loc,1,v); } -static void glUniform4uiv_ei (GLint loc, const unsigned int* v) { glUniform4uiv(loc,1,v); } +inline void glUniform2uiv_ei (GLint loc, const unsigned int* v) { glUniform2uiv(loc,1,v); } +inline void glUniform3uiv_ei (GLint loc, const unsigned int* v) { glUniform3uiv(loc,1,v); } +inline void glUniform4uiv_ei (GLint loc, const unsigned int* v) { glUniform4uiv(loc,1,v); } EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 2,2uiv_ei) EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 3,3uiv_ei) @@ -305,9 +305,9 @@ EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,unsigned int, 4,4uiv_ei) #endif #ifdef GL_ARB_gpu_shader_fp64 -static void glUniform2dv_ei (GLint loc, const double* v) { glUniform2dv(loc,1,v); } -static void glUniform3dv_ei (GLint loc, const double* v) { glUniform3dv(loc,1,v); } -static void glUniform4dv_ei (GLint loc, const double* v) { glUniform4dv(loc,1,v); } +inline void glUniform2dv_ei (GLint loc, const double* v) { glUniform2dv(loc,1,v); } +inline void glUniform3dv_ei (GLint loc, const double* v) { glUniform3dv(loc,1,v); } +inline void glUniform4dv_ei (GLint loc, const double* v) { glUniform4dv(loc,1,v); } EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,double, 2,2dv_ei) EIGEN_GL_FUNC1_SPECIALIZATION_VEC(glUniform,GLint,const,double, 3,3dv_ei) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/AutoDiff/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/AutoDiff/CMakeLists.txt deleted file mode 100644 index ad91fd9c4..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/AutoDiff/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_AutoDiff_SRCS "*.h") - -INSTALL(FILES - ${Eigen_AutoDiff_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/AutoDiff COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/BVH/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/BVH/CMakeLists.txt deleted file mode 100644 index b377d865c..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/BVH/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_BVH_SRCS "*.h") - -INSTALL(FILES - ${Eigen_BVH_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/BVH COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/CMakeLists.txt deleted file mode 100644 index f3180b52b..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/CMakeLists.txt +++ /dev/null @@ -1,13 +0,0 @@ -ADD_SUBDIRECTORY(AutoDiff) -ADD_SUBDIRECTORY(BVH) -ADD_SUBDIRECTORY(FFT) -ADD_SUBDIRECTORY(IterativeSolvers) -ADD_SUBDIRECTORY(MatrixFunctions) -ADD_SUBDIRECTORY(MoreVectorization) -ADD_SUBDIRECTORY(NonLinearOptimization) -ADD_SUBDIRECTORY(NumericalDiff) -ADD_SUBDIRECTORY(Polynomials) -ADD_SUBDIRECTORY(Skyline) -ADD_SUBDIRECTORY(SparseExtra) -ADD_SUBDIRECTORY(KroneckerProduct) -ADD_SUBDIRECTORY(Splines) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/FFT/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/FFT/CMakeLists.txt deleted file mode 100644 index edcffcb18..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/FFT/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_FFT_SRCS "*.h") - -INSTALL(FILES - ${Eigen_FFT_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/FFT COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/IterativeSolvers/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/IterativeSolvers/CMakeLists.txt deleted file mode 100644 index 7986afc5e..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/IterativeSolvers/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_IterativeSolvers_SRCS "*.h") - -INSTALL(FILES - ${Eigen_IterativeSolvers_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/IterativeSolvers COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h b/uppsrc/plugin/Eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h index 073367506..c8c84069e 100644 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h +++ b/uppsrc/plugin/Eigen/unsupported/Eigen/src/IterativeSolvers/GMRES.h @@ -2,7 +2,7 @@ // for linear algebra. // // Copyright (C) 2011 Gael Guennebaud -// Copyright (C) 2012 Kolja Brix +// Copyright (C) 2012, 2014 Kolja Brix // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed @@ -72,16 +72,20 @@ bool gmres(const MatrixType & mat, const Rhs & rhs, Dest & x, const Precondition VectorType p0 = rhs - mat*x; VectorType r0 = precond.solve(p0); -// RealScalar r0_sqnorm = r0.squaredNorm(); + + // is initial guess already good enough? + if(abs(r0.norm()) < tol) { + return true; + } VectorType w = VectorType::Zero(restart + 1); - FMatrixType H = FMatrixType::Zero(m, restart + 1); + FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix VectorType tau = VectorType::Zero(restart + 1); std::vector < JacobiRotation < Scalar > > G(restart); // generate first Householder vector - VectorType e; + VectorType e(m-1); RealScalar beta; r0.makeHouseholder(e, tau.coeffRef(0), beta); w(0)=(Scalar) beta; diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/KroneckerProduct/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/KroneckerProduct/CMakeLists.txt deleted file mode 100644 index 4daefebee..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/KroneckerProduct/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_KroneckerProduct_SRCS "*.h") - -INSTALL(FILES - ${Eigen_KroneckerProduct_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/KroneckerProduct COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt deleted file mode 100644 index 8513803ce..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/LevenbergMarquardt/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_LevenbergMarquardt_SRCS "*.h") - -INSTALL(FILES - ${Eigen_LevenbergMarquardt_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/LevenbergMarquardt COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/MatrixFunctions/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/MatrixFunctions/CMakeLists.txt deleted file mode 100644 index cdde64d2c..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/MatrixFunctions/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_MatrixFunctions_SRCS "*.h") - -INSTALL(FILES - ${Eigen_MatrixFunctions_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/MatrixFunctions COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h b/uppsrc/plugin/Eigen/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h index c32437281..78a307e96 100644 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h +++ b/uppsrc/plugin/Eigen/unsupported/Eigen/src/MatrixFunctions/MatrixPower.h @@ -110,7 +110,6 @@ void MatrixPowerAtomic::compute2x2(MatrixType& res, RealScalar p) co using std::abs; using std::pow; - ArrayType logTdiag = m_A.diagonal().array().log(); res.coeffRef(0,0) = pow(m_A.coeff(0,0), p); for (Index i=1; i < m_A.cols(); ++i) { diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/MoreVectorization/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/MoreVectorization/CMakeLists.txt deleted file mode 100644 index 1b887cc8e..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/MoreVectorization/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_MoreVectorization_SRCS "*.h") - -INSTALL(FILES - ${Eigen_MoreVectorization_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/MoreVectorization COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/NonLinearOptimization/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/NonLinearOptimization/CMakeLists.txt deleted file mode 100644 index 9322ddadf..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/NonLinearOptimization/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_NonLinearOptimization_SRCS "*.h") - -INSTALL(FILES - ${Eigen_NonLinearOptimization_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/NonLinearOptimization COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/NumericalDiff/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/NumericalDiff/CMakeLists.txt deleted file mode 100644 index 1199aca2f..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/NumericalDiff/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_NumericalDiff_SRCS "*.h") - -INSTALL(FILES - ${Eigen_NumericalDiff_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/NumericalDiff COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/Polynomials/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/Polynomials/CMakeLists.txt deleted file mode 100644 index 51f13f3cb..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/Polynomials/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Polynomials_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Polynomials_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Polynomials COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/SVD/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/SVD/CMakeLists.txt deleted file mode 100644 index b40baf092..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/SVD/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SVD_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SVD_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}unsupported/Eigen/src/SVD COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/Skyline/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/Skyline/CMakeLists.txt deleted file mode 100644 index 3bf1b0dd4..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/Skyline/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Skyline_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Skyline_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Skyline COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/SparseExtra/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/SparseExtra/CMakeLists.txt deleted file mode 100644 index 7ea32ca5e..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/SparseExtra/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_SparseExtra_SRCS "*.h") - -INSTALL(FILES - ${Eigen_SparseExtra_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/SparseExtra COMPONENT Devel - ) diff --git a/uppsrc/plugin/Eigen/unsupported/Eigen/src/Splines/CMakeLists.txt b/uppsrc/plugin/Eigen/unsupported/Eigen/src/Splines/CMakeLists.txt deleted file mode 100644 index 55c6271e9..000000000 --- a/uppsrc/plugin/Eigen/unsupported/Eigen/src/Splines/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ -FILE(GLOB Eigen_Splines_SRCS "*.h") - -INSTALL(FILES - ${Eigen_Splines_SRCS} - DESTINATION ${INCLUDE_INSTALL_DIR}/unsupported/Eigen/src/Splines COMPONENT Devel - )