manolas
/
vcglib
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Updated Eigen library to the 3.2.2 stable version.

This commit is contained in:
Paolo Cignoni 2014-10-21 16:51:37 +00:00
parent af0f42fedf
commit 54875ffa0a
35 changed files with 407 additions and 219 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
eigenlib/Eigen/Core
View File

@ -95,7 +95,7 @@
extern "C" { extern "C" {
// In theory we should only include immintrin.h and not the other *mmintrin.h header files directly. // 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: // 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 <immintrin.h> #include <immintrin.h>
#else #else
#include <emmintrin.h> #include <emmintrin.h>
@ -165,7 +165,7 @@
#endif #endif
// required for __cpuid, needs to be included after cmath // 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 <intrin.h> #include <intrin.h>
#endif #endif

50
eigenlib/Eigen/src/Cholesky/LDLT.h
View File

@ -274,30 +274,13 @@ template<> struct ldlt_inplace<Lower>
return true; return true;
} }
RealScalar cutoff(0), biggest_in_corner;
for (Index k = 0; k < size; ++k) for (Index k = 0; k < size; ++k)
{ {
// Find largest diagonal element // Find largest diagonal element
Index index_of_biggest_in_corner; 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; 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<Scalar>::epsilon() * biggest_in_corner);
}
// Finish early if the matrix is not full rank.
if(biggest_in_corner < cutoff)
{
for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
break;
}
transpositions.coeffRef(k) = index_of_biggest_in_corner; transpositions.coeffRef(k) = index_of_biggest_in_corner;
if(k != index_of_biggest_in_corner) if(k != index_of_biggest_in_corner)
{ {
@ -328,15 +311,20 @@ template<> struct ldlt_inplace<Lower>
if(k>0) 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(); mat.coeffRef(k,k) -= (A10 * temp.head(k)).value();
if(rs>0) if(rs>0)
A21.noalias() -= A20 * temp.head(k); 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)); RealScalar realAkk = numext::real(mat.coeffRef(k,k));
if((rs>0) && (abs(realAkk) > RealScalar(0)))
A21 /= realAkk;
if (sign == PositiveSemiDef) { if (sign == PositiveSemiDef) {
if (realAkk < 0) sign = Indefinite; if (realAkk < 0) sign = Indefinite;
} else if (sign == NegativeSemiDef) { } else if (sign == NegativeSemiDef) {
@ -516,14 +504,20 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
typedef typename LDLTType::MatrixType MatrixType; typedef typename LDLTType::MatrixType MatrixType;
typedef typename LDLTType::Scalar Scalar; typedef typename LDLTType::Scalar Scalar;
typedef typename LDLTType::RealScalar RealScalar; typedef typename LDLTType::RealScalar RealScalar;
const Diagonal<const MatrixType> vectorD = dec().vectorD(); const typename Diagonal<const MatrixType>::RealReturnType vectorD(dec().vectorD());
RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(), // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon
RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS // as motivated by LAPACK's xGELSS:
// RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() *NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::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<RealScalar>::highest();
for (Index i = 0; i < vectorD.size(); ++i) { for (Index i = 0; i < vectorD.size(); ++i) {
if(abs(vectorD(i)) > tolerance) if(abs(vectorD(i)) > tolerance)
dst.row(i) /= vectorD(i); dst.row(i) /= vectorD(i);
else else
dst.row(i).setZero(); dst.row(i).setZero();
} }
// dst = L^-T (D^-1 L^-1 P b) // dst = L^-T (D^-1 L^-1 P b)
@ -576,7 +570,7 @@ MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const
// L^* P // L^* P
res = matrixU() * res; res = matrixU() * res;
// D(L^*P) // D(L^*P)
res = vectorD().asDiagonal() * res; res = vectorD().real().asDiagonal() * res;
// L(DL^*P) // L(DL^*P)
res = matrixL() * res; res = matrixL() * res;
// P^T (LDL^*P) // P^T (LDL^*P)

2
eigenlib/Eigen/src/Core/Block.h
View File

@ -81,7 +81,7 @@ struct traits<Block<XprType, BlockRows, BlockCols, InnerPanel> > : traits<XprTyp
&& (InnerStrideAtCompileTime == 1) && (InnerStrideAtCompileTime == 1)
? PacketAccessBit : 0, ? PacketAccessBit : 0,
MaskAlignedBit = (InnerPanel && (OuterStrideAtCompileTime!=Dynamic) && (((OuterStrideAtCompileTime * int(sizeof(Scalar))) % 16) == 0)) ? AlignedBit : 0, MaskAlignedBit = (InnerPanel && (OuterStrideAtCompileTime!=Dynamic) && (((OuterStrideAtCompileTime * int(sizeof(Scalar))) % 16) == 0)) ? AlignedBit : 0,
FlagsLinearAccessBit = (RowsAtCompileTime == 1 || ColsAtCompileTime == 1) ? LinearAccessBit : 0, FlagsLinearAccessBit = (RowsAtCompileTime == 1 || ColsAtCompileTime == 1 || (InnerPanel && (traits<XprType>::Flags&LinearAccessBit))) ? LinearAccessBit : 0,
FlagsLvalueBit = is_lvalue<XprType>::value ? LvalueBit : 0, FlagsLvalueBit = is_lvalue<XprType>::value ? LvalueBit : 0,
FlagsRowMajorBit = IsRowMajor ? RowMajorBit : 0, FlagsRowMajorBit = IsRowMajor ? RowMajorBit : 0,
Flags0 = traits<XprType>::Flags & ( (HereditaryBits & ~RowMajorBit) | Flags0 = traits<XprType>::Flags & ( (HereditaryBits & ~RowMajorBit) |

11
eigenlib/Eigen/src/Core/CommaInitializer.h
View File

@ -43,6 +43,17 @@ struct CommaInitializer
m_xpr.block(0, 0, other.rows(), other.cols()) = other; 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<CommaInitializer&>(o).m_row = m_xpr.rows();
const_cast<CommaInitializer&>(o).m_col = m_xpr.cols();
const_cast<CommaInitializer&>(o).m_currentBlockRows = 0;
}
/* inserts a scalar value in the target matrix */ /* inserts a scalar value in the target matrix */
CommaInitializer& operator,(const Scalar& s) CommaInitializer& operator,(const Scalar& s)
{ {

16
eigenlib/Eigen/src/Core/DenseStorage.h
View File

@ -24,6 +24,14 @@ namespace internal {
struct constructor_without_unaligned_array_assert {}; struct constructor_without_unaligned_array_assert {};
template<typename T, int Size> 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 /** \internal
* Static array. If the MatrixOrArrayOptions require auto-alignment, the array will be automatically aligned: * 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. * to 16 bytes boundary if the total size is a multiple of 16 bytes.
@ -38,12 +46,12 @@ struct plain_array
plain_array() plain_array()
{ {
EIGEN_STATIC_ASSERT(Size * sizeof(T) <= 128 * 128 * 8, OBJECT_ALLOCATED_ON_STACK_IS_TOO_BIG); check_static_allocation_size<T,Size>();
} }
plain_array(constructor_without_unaligned_array_assert) 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<T,Size>();
} }
}; };
@ -76,12 +84,12 @@ struct plain_array<T, Size, MatrixOrArrayOptions, 16>
plain_array() plain_array()
{ {
EIGEN_MAKE_UNALIGNED_ARRAY_ASSERT(0xf); 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<T,Size>();
} }
plain_array(constructor_without_unaligned_array_assert) 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<T,Size>();
} }
}; };

4
eigenlib/Eigen/src/Core/Functors.h
View File

@ -589,7 +589,7 @@ struct linspaced_op_impl<Scalar,true>
template<typename Index> template<typename Index>
EIGEN_STRONG_INLINE const Packet packetOp(Index i) const EIGEN_STRONG_INLINE const Packet packetOp(Index i) const
{ return internal::padd(m_lowPacket, pmul(m_stepPacket, padd(pset1<Packet>(i),m_interPacket))); } { return internal::padd(m_lowPacket, pmul(m_stepPacket, padd(pset1<Packet>(Scalar(i)),m_interPacket))); }
const Scalar m_low; const Scalar m_low;
const Scalar m_step; const Scalar m_step;
@ -609,7 +609,7 @@ template <typename Scalar, bool RandomAccess> struct functor_traits< linspaced_o
template <typename Scalar, bool RandomAccess> struct linspaced_op template <typename Scalar, bool RandomAccess> struct linspaced_op
{ {
typedef typename packet_traits<Scalar>::type Packet; typedef typename packet_traits<Scalar>::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<typename Index> template<typename Index>
EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return impl(i); } EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return impl(i); }

2
eigenlib/Eigen/src/Core/MapBase.h
View File

@ -237,6 +237,8 @@ template<typename Derived> class MapBase<Derived, WriteAccessors>
using Base::Base::operator=; using Base::Base::operator=;
}; };
#undef EIGEN_STATIC_ASSERT_INDEX_BASED_ACCESS
} // end namespace Eigen } // end namespace Eigen
#endif // EIGEN_MAPBASE_H #endif // EIGEN_MAPBASE_H

10
eigenlib/Eigen/src/Core/Ref.h
View File

@ -101,7 +101,7 @@ struct traits<Ref<_PlainObjectType, _Options, _StrideType> >
template<typename Derived> struct match { template<typename Derived> struct match {
enum { enum {
HasDirectAccess = internal::has_direct_access<Derived>::ret, HasDirectAccess = internal::has_direct_access<Derived>::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) InnerStrideMatch = int(StrideType::InnerStrideAtCompileTime)==int(Dynamic)
|| int(StrideType::InnerStrideAtCompileTime)==int(Derived::InnerStrideAtCompileTime) || int(StrideType::InnerStrideAtCompileTime)==int(Derived::InnerStrideAtCompileTime)
|| (int(StrideType::InnerStrideAtCompileTime)==0 && int(Derived::InnerStrideAtCompileTime)==1), || (int(StrideType::InnerStrideAtCompileTime)==0 && int(Derived::InnerStrideAtCompileTime)==1),
@ -172,8 +172,12 @@ protected:
} }
else else
::new (static_cast<Base*>(this)) Base(expr.data(), expr.rows(), expr.cols()); ::new (static_cast<Base*>(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; StrideBase m_stride;

8
eigenlib/Eigen/src/Core/TriangularMatrix.h
View File

@ -278,21 +278,21 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
/** Efficient triangular matrix times vector/matrix product */ /** Efficient triangular matrix times vector/matrix product */
template<typename OtherDerived> template<typename OtherDerived>
TriangularProduct<Mode,true,MatrixType,false,OtherDerived, OtherDerived::IsVectorAtCompileTime> TriangularProduct<Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
operator*(const MatrixBase<OtherDerived>& rhs) const operator*(const MatrixBase<OtherDerived>& rhs) const
{ {
return TriangularProduct return TriangularProduct
<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime> <Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
(m_matrix, rhs.derived()); (m_matrix, rhs.derived());
} }
/** Efficient vector/matrix times triangular matrix product */ /** Efficient vector/matrix times triangular matrix product */
template<typename OtherDerived> friend template<typename OtherDerived> friend
TriangularProduct<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false> TriangularProduct<Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs) operator*(const MatrixBase<OtherDerived>& lhs, const TriangularView& rhs)
{ {
return TriangularProduct return TriangularProduct
<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false> <Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
(lhs.derived(),rhs.m_matrix); (lhs.derived(),rhs.m_matrix);
} }

19
eigenlib/Eigen/src/Core/util/MKL_support.h
View File

@ -54,8 +54,25 @@
#endif #endif
#if defined EIGEN_USE_MKL #if defined EIGEN_USE_MKL
# include <mkl.h>
/*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 <mkl.h> #if defined EIGEN_USE_MKL
#include <mkl_lapacke.h> #include <mkl_lapacke.h>
#define EIGEN_MKL_VML_THRESHOLD 128 #define EIGEN_MKL_VML_THRESHOLD 128

5
eigenlib/Eigen/src/Core/util/Macros.h
View File

@ -13,7 +13,7 @@
#define EIGEN_WORLD_VERSION 3 #define EIGEN_WORLD_VERSION 3
#define EIGEN_MAJOR_VERSION 2 #define EIGEN_MAJOR_VERSION 2
#define EIGEN_MINOR_VERSION 1 #define EIGEN_MINOR_VERSION 2
#define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \ #define EIGEN_VERSION_AT_LEAST(x,y,z) (EIGEN_WORLD_VERSION>x || (EIGEN_WORLD_VERSION>=x && \
(EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \ (EIGEN_MAJOR_VERSION>y || (EIGEN_MAJOR_VERSION>=y && \
@ -289,7 +289,8 @@ namespace Eigen {
#endif #endif
#ifndef EIGEN_STACK_ALLOCATION_LIMIT #ifndef EIGEN_STACK_ALLOCATION_LIMIT
#define EIGEN_STACK_ALLOCATION_LIMIT 20000 // 131072 == 128 KB
#define EIGEN_STACK_ALLOCATION_LIMIT 131072
#endif #endif
#ifndef EIGEN_DEFAULT_IO_FORMAT #ifndef EIGEN_DEFAULT_IO_FORMAT

15
eigenlib/Eigen/src/Core/util/Memory.h
View File

@ -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 // The defined(_mm_free) is just here to verify that this MSVC version
// implements _mm_malloc/_mm_free based on the corresponding _aligned_ // 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. // 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); result = _aligned_realloc(ptr,new_size,16);
#else #else
result = generic_aligned_realloc(ptr,new_size,old_size); result = generic_aligned_realloc(ptr,new_size,old_size);
#endif #endif
#elif defined(_MSC_VER) #elif defined(_MSC_VER) && (!defined(_WIN32_WCE))
result = _aligned_realloc(ptr,new_size,16); result = _aligned_realloc(ptr,new_size,16);
#else #else
result = handmade_aligned_realloc(ptr,new_size,old_size); result = handmade_aligned_realloc(ptr,new_size,old_size);
@ -777,9 +777,9 @@ namespace internal {
#ifdef EIGEN_CPUID #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<const int*>(vendor))[0] && abcd[3]==(reinterpret_cast<const int*>(vendor))[1] && abcd[2]==(reinterpret_cast<const int*>(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) inline void queryCacheSizes_intel_direct(int& l1, int& l2, int& l3)
@ -921,13 +921,16 @@ inline void queryCacheSizes(int& l1, int& l2, int& l3)
{ {
#ifdef EIGEN_CPUID #ifdef EIGEN_CPUID
int abcd[4]; 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 // identify the CPU vendor
EIGEN_CPUID(abcd,0x0,0); EIGEN_CPUID(abcd,0x0,0);
int max_std_funcs = abcd[1]; 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); 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); queryCacheSizes_amd(l1,l2,l3);
else else
// by default let's use Intel's API // by default let's use Intel's API

14
eigenlib/Eigen/src/Geometry/Quaternion.h
View File

@ -203,6 +203,8 @@ public:
* \li \c Quaternionf for \c float * \li \c Quaternionf for \c float
* \li \c Quaterniond for \c double * \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 * \sa class AngleAxis, class Transform
*/ */
@ -344,7 +346,7 @@ class Map<const Quaternion<_Scalar>, _Options >
/** Constructs a Mapped Quaternion object from the pointer \a coeffs /** 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 * \code *coeffs == {x, y, z, w} \endcode
* *
* If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */ * If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */
@ -464,7 +466,7 @@ QuaternionBase<Derived>::_transformVector(Vector3 v) const
// Note that this algorithm comes from the optimization by hand // Note that this algorithm comes from the optimization by hand
// of the conversion to a Matrix followed by a Matrix/Vector product. // of the conversion to a Matrix followed by a Matrix/Vector product.
// It appears to be much faster than the common algorithm found // 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 as temporaries.
Vector3 uv = this->vec().cross(v); Vector3 uv = this->vec().cross(v);
uv += uv; uv += uv;
@ -584,7 +586,7 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri
// which yields a singular value problem // which yields a singular value problem
if (c < Scalar(-1)+NumTraits<Scalar>::dummy_precision()) if (c < Scalar(-1)+NumTraits<Scalar>::dummy_precision())
{ {
c = max<Scalar>(c,-1); c = (max)(c,Scalar(-1));
Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose(); Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV); JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
Vector3 axis = svd.matrixV().col(2); Vector3 axis = svd.matrixV().col(2);
@ -667,10 +669,10 @@ QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& oth
{ {
using std::acos; using std::acos;
using std::abs; using std::abs;
double d = abs(this->dot(other)); Scalar d = abs(this->dot(other));
if (d>=1.0) if (d>=Scalar(1))
return Scalar(0); return Scalar(0);
return static_cast<Scalar>(2 * acos(d)); return Scalar(2) * acos(d);
} }

4
eigenlib/Eigen/src/Geometry/Transform.h
View File

@ -194,9 +194,9 @@ public:
/** type of the matrix used to represent the linear part of the transformation */ /** type of the matrix used to represent the linear part of the transformation */
typedef Matrix<Scalar,Dim,Dim,Options> LinearMatrixType; typedef Matrix<Scalar,Dim,Dim,Options> LinearMatrixType;
/** type of read/write reference to the linear part of the transformation */ /** type of read/write reference to the linear part of the transformation */
typedef Block<MatrixType,Dim,Dim,int(Mode)==(AffineCompact)> LinearPart; typedef Block<MatrixType,Dim,Dim,int(Mode)==(AffineCompact) && (Options&RowMajor)==0> LinearPart;
/** type of read reference to the linear part of the transformation */ /** type of read reference to the linear part of the transformation */
typedef const Block<ConstMatrixType,Dim,Dim,int(Mode)==(AffineCompact)> ConstLinearPart; typedef const Block<ConstMatrixType,Dim,Dim,int(Mode)==(AffineCompact) && (Options&RowMajor)==0> ConstLinearPart;
/** type of read/write reference to the affine part of the transformation */ /** type of read/write reference to the affine part of the transformation */
typedef typename internal::conditional<int(Mode)==int(AffineCompact), typedef typename internal::conditional<int(Mode)==int(AffineCompact),
MatrixType&, MatrixType&,

10
eigenlib/Eigen/src/Geometry/Umeyama.h
View File

@ -113,7 +113,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
const Index n = src.cols(); // number of measurements const Index n = src.cols(); // number of measurements
// required for demeaning ... // required for demeaning ...
const RealScalar one_over_n = 1 / static_cast<RealScalar>(n); const RealScalar one_over_n = RealScalar(1) / static_cast<RealScalar>(n);
// computation of mean // computation of mean
const VectorType src_mean = src.rowwise().sum() * one_over_n; const VectorType src_mean = src.rowwise().sum() * one_over_n;
@ -136,16 +136,16 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
// Eq. (39) // Eq. (39)
VectorType S = VectorType::Ones(m); VectorType S = VectorType::Ones(m);
if (sigma.determinant()<0) S(m-1) = -1; if (sigma.determinant()<Scalar(0)) S(m-1) = Scalar(-1);
// Eq. (40) and (43) // Eq. (40) and (43)
const VectorType& d = svd.singularValues(); const VectorType& d = svd.singularValues();
Index rank = 0; for (Index i=0; i<m; ++i) if (!internal::isMuchSmallerThan(d.coeff(i),d.coeff(0))) ++rank; Index rank = 0; for (Index i=0; i<m; ++i) if (!internal::isMuchSmallerThan(d.coeff(i),d.coeff(0))) ++rank;
if (rank == m-1) { if (rank == m-1) {
if ( svd.matrixU().determinant() * svd.matrixV().determinant() > 0 ) { if ( svd.matrixU().determinant() * svd.matrixV().determinant() > Scalar(0) ) {
Rt.block(0,0,m,m).noalias() = svd.matrixU()*svd.matrixV().transpose(); Rt.block(0,0,m,m).noalias() = svd.matrixU()*svd.matrixV().transpose();
} else { } 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(); Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose();
S(m-1) = s; S(m-1) = s;
} }
@ -156,7 +156,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
if (with_scaling) if (with_scaling)
{ {
// Eq. (42) // 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) // Eq. (41)
Rt.col(m).head(m) = dst_mean; Rt.col(m).head(m) = dst_mean;

2
eigenlib/Eigen/src/Householder/BlockHouseholder.h
View File

@ -48,7 +48,7 @@ void apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vec
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
enum { TFactorSize = MatrixType::ColsAtCompileTime }; enum { TFactorSize = MatrixType::ColsAtCompileTime };
Index nbVecs = vectors.cols(); Index nbVecs = vectors.cols();
Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize> T(nbVecs,nbVecs); Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize, ColMajor> T(nbVecs,nbVecs);
make_block_householder_triangular_factor(T, vectors, hCoeffs); make_block_householder_triangular_factor(T, vectors, hCoeffs);
const TriangularView<const VectorsType, UnitLower>& V(vectors); const TriangularView<const VectorsType, UnitLower>& V(vectors);

3
eigenlib/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
View File

@ -61,6 +61,7 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
VectorType s(n), t(n); VectorType s(n), t(n);
RealScalar tol2 = tol*tol; RealScalar tol2 = tol*tol;
RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon();
int i = 0; int i = 0;
int restarts = 0; int restarts = 0;
@ -69,7 +70,7 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
Scalar rho_old = rho; Scalar rho_old = rho;
rho = r0.dot(r); 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 // The new residual vector became too orthogonal to the arbitrarily choosen direction r0
// Let's restart with a new r0: // Let's restart with a new r0:

13
eigenlib/Eigen/src/LU/FullPivLU.h
View File

@ -20,10 +20,11 @@ namespace Eigen {
* *
* \param MatrixType the type of the matrix of which we are computing the LU decomposition * \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 * This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A is
* is decomposed as A = PLUQ where L is unit-lower-triangular, U is upper-triangular, and P and Q * decomposed as \f$ A = P^{-1} L U Q^{-1} \f$ where L is unit-lower-triangular, U is
* are permutation matrices. This is a rank-revealing LU decomposition. The eigenvalues (diagonal * upper-triangular, and P and Q are permutation matrices. This is a rank-revealing LU
* coefficients) of U are sorted in such a way that any zeros are at the end. * 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 * This decomposition provides the generic approach to solving systems of linear equations, computing
* the rank, invertibility, inverse, kernel, and determinant. * the rank, invertibility, inverse, kernel, and determinant.
@ -511,8 +512,8 @@ typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant
} }
/** \returns the matrix represented by the decomposition, /** \returns the matrix represented by the decomposition,
* i.e., it returns the product: P^{-1} L U Q^{-1}. * i.e., it returns the product: \f$ P^{-1} L U Q^{-1} \f$.
* This function is provided for debug purpose. */ * This function is provided for debug purposes. */
template<typename MatrixType> template<typename MatrixType>
MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
{ {

12
eigenlib/Eigen/src/OrderingMethods/Ordering.h
View File

@ -109,7 +109,7 @@ class NaturalOrdering
* \class COLAMDOrdering * \class COLAMDOrdering
* *
* Functor computing the \em column \em approximate \em minimum \em degree ordering * 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<typename Index> template<typename Index>
class COLAMDOrdering class COLAMDOrdering
@ -118,10 +118,14 @@ class COLAMDOrdering
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
typedef Matrix<Index, Dynamic, 1> IndexVector; typedef Matrix<Index, Dynamic, 1> 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 <typename MatrixType> template <typename MatrixType>
void operator() (const MatrixType& mat, PermutationType& perm) 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 m = mat.rows();
Index n = mat.cols(); Index n = mat.cols();
Index nnz = mat.nonZeros(); Index nnz = mat.nonZeros();
@ -132,12 +136,12 @@ class COLAMDOrdering
Index stats [COLAMD_STATS]; Index stats [COLAMD_STATS];
internal::colamd_set_defaults(knobs); internal::colamd_set_defaults(knobs);
Index info;
IndexVector p(n+1), A(Alen); IndexVector p(n+1), A(Alen);
for(Index i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i]; for(Index i=0; i <= n; i++) p(i) = mat.outerIndexPtr()[i];
for(Index i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i]; for(Index i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()[i];
// Call Colamd routine to compute the ordering // 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 " ); eigen_assert( info && "COLAMD failed " );
perm.resize(n); perm.resize(n);

3
eigenlib/Eigen/src/QR/ColPivHouseholderQR.h
View File

@ -76,7 +76,8 @@ template<typename _MatrixType> class ColPivHouseholderQR
m_colsTranspositions(), m_colsTranspositions(),
m_temp(), m_temp(),
m_colSqNorms(), m_colSqNorms(),
m_isInitialized(false) {} m_isInitialized(false),
m_usePrescribedThreshold(false) {}
/** \brief Default Constructor with memory preallocation /** \brief Default Constructor with memory preallocation
* *

105
eigenlib/Eigen/src/SVD/JacobiSVD.h
View File

@ -375,17 +375,19 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
Scalar z; Scalar z;
JacobiRotation<Scalar> rot; JacobiRotation<Scalar> rot;
RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p))); RealScalar n = sqrt(numext::abs2(work_matrix.coeff(p,p)) + numext::abs2(work_matrix.coeff(q,p)));
if(n==0) if(n==0)
{ {
z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); z = abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q);
work_matrix.row(p) *= z; work_matrix.row(p) *= z;
if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z); if(svd.computeU()) svd.m_matrixU.col(p) *= conj(z);
if(work_matrix.coeff(q,q)!=Scalar(0)) if(work_matrix.coeff(q,q)!=Scalar(0))
{
z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); z = abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q);
else work_matrix.row(q) *= z;
z = Scalar(0); 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 else
{ {
@ -415,6 +417,7 @@ void real_2x2_jacobi_svd(const MatrixType& matrix, Index p, Index q,
JacobiRotation<RealScalar> *j_right) JacobiRotation<RealScalar> *j_right)
{ {
using std::sqrt; using std::sqrt;
using std::abs;
Matrix<RealScalar,2,2> m; Matrix<RealScalar,2,2> m;
m << numext::real(matrix.coeff(p,p)), numext::real(matrix.coeff(p,q)), 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)); 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 else
{ {
RealScalar u = d / t; RealScalar t2d2 = numext::hypot(t,d);
rot1.c() = RealScalar(1) / sqrt(RealScalar(1) + numext::abs2(u)); rot1.c() = abs(t)/t2d2;
rot1.s() = rot1.c() * u; rot1.s() = d/t2d2;
if(t<RealScalar(0))
rot1.s() = -rot1.s();
} }
m.applyOnTheLeft(0,1,rot1); m.applyOnTheLeft(0,1,rot1);
j_right->makeJacobi(m,0,1); j_right->makeJacobi(m,0,1);
@ -531,8 +536,9 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
JacobiSVD() JacobiSVD()
: m_isInitialized(false), : m_isInitialized(false),
m_isAllocated(false), m_isAllocated(false),
m_usePrescribedThreshold(false),
m_computationOptions(0), m_computationOptions(0),
m_rows(-1), m_cols(-1) m_rows(-1), m_cols(-1), m_diagSize(0)
{} {}
@ -545,6 +551,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0)
: m_isInitialized(false), : m_isInitialized(false),
m_isAllocated(false), m_isAllocated(false),
m_usePrescribedThreshold(false),
m_computationOptions(0), m_computationOptions(0),
m_rows(-1), m_cols(-1) m_rows(-1), m_cols(-1)
{ {
@ -564,6 +571,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0)
: m_isInitialized(false), : m_isInitialized(false),
m_isAllocated(false), m_isAllocated(false),
m_usePrescribedThreshold(false),
m_computationOptions(0), m_computationOptions(0),
m_rows(-1), m_cols(-1) m_rows(-1), m_cols(-1)
{ {
@ -665,6 +673,69 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
eigen_assert(m_isInitialized && "JacobiSVD is not initialized."); eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
return m_nonzeroSingularValues; 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<Scalar>::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<Index>)(1,m_diagSize)*NumTraits<Scalar>::epsilon();
}
inline Index rows() const { return m_rows; } inline Index rows() const { return m_rows; }
inline Index cols() const { return m_cols; } inline Index cols() const { return m_cols; }
@ -677,11 +748,12 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
MatrixVType m_matrixV; MatrixVType m_matrixV;
SingularValuesType m_singularValues; SingularValuesType m_singularValues;
WorkMatrixType m_workMatrix; WorkMatrixType m_workMatrix;
bool m_isInitialized, m_isAllocated; bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
bool m_computeFullU, m_computeThinU; bool m_computeFullU, m_computeThinU;
bool m_computeFullV, m_computeThinV; bool m_computeFullV, m_computeThinV;
unsigned int m_computationOptions; unsigned int m_computationOptions;
Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize; Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
RealScalar m_prescribedThreshold;
template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex> template<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
friend struct internal::svd_precondition_2x2_block_to_be_real; friend struct internal::svd_precondition_2x2_block_to_be_real;
@ -764,6 +836,11 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols); if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize); if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
} }
// Scaling factor to reduce over/under-flows
RealScalar scale = m_workMatrix.cwiseAbs().maxCoeff();
if(scale==RealScalar(0)) scale = RealScalar(1);
m_workMatrix /= scale;
/*** step 2. The main Jacobi SVD iteration. ***/ /*** step 2. The main Jacobi SVD iteration. ***/
@ -833,6 +910,8 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i)); if(computeV()) m_matrixV.col(pos).swap(m_matrixV.col(i));
} }
} }
m_singularValues *= scale;
m_isInitialized = true; m_isInitialized = true;
return *this; return *this;
@ -854,11 +933,11 @@ struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
// So A^{-1} = V S^{-1} U^* // So A^{-1} = V S^{-1} U^*
Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp; Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp;
Index nonzeroSingVals = dec().nonzeroSingularValues(); Index rank = dec().rank();
tmp.noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs(); tmp.noalias() = dec().matrixU().leftCols(rank).adjoint() * rhs();
tmp = dec().singularValues().head(nonzeroSingVals).asDiagonal().inverse() * tmp; tmp = dec().singularValues().head(rank).asDiagonal().inverse() * tmp;
dst = dec().matrixV().leftCols(nonzeroSingVals) * tmp; dst = dec().matrixV().leftCols(rank) * tmp;
} }
}; };
} // end namespace internal } // end namespace internal

42
eigenlib/Eigen/src/SparseCholesky/SimplicialCholesky.h
View File

@ -37,6 +37,7 @@ class SimplicialCholeskyBase : internal::noncopyable
{ {
public: public:
typedef typename internal::traits<Derived>::MatrixType MatrixType; typedef typename internal::traits<Derived>::MatrixType MatrixType;
typedef typename internal::traits<Derived>::OrderingType OrderingType;
enum { UpLo = internal::traits<Derived>::UpLo }; enum { UpLo = internal::traits<Derived>::UpLo };
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar; typedef typename MatrixType::RealScalar RealScalar;
@ -240,15 +241,16 @@ class SimplicialCholeskyBase : internal::noncopyable
RealScalar m_shiftScale; RealScalar m_shiftScale;
}; };
template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT; template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT;
template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT; template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT;
template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky; template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialCholesky;
namespace internal { namespace internal {
template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> > template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
{ {
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef _Ordering OrderingType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
@ -259,9 +261,10 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixTyp
static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
}; };
template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> > template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
{ {
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef _Ordering OrderingType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::Index Index; typedef typename MatrixType::Index Index;
@ -272,9 +275,10 @@ template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixTyp
static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
}; };
template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> > template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
{ {
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
typedef _Ordering OrderingType;
enum { UpLo = _UpLo }; enum { UpLo = _UpLo };
}; };
@ -294,11 +298,12 @@ template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_Matr
* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * \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 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is 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<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo, typename _Ordering>
class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> > class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
@ -382,11 +387,12 @@ public:
* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> * \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 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is 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<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo, typename _Ordering>
class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> > class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
@ -467,8 +473,8 @@ public:
* *
* \sa class SimplicialLDLT, class SimplicialLLT * \sa class SimplicialLDLT, class SimplicialLLT
*/ */
template<typename _MatrixType, int _UpLo> template<typename _MatrixType, int _UpLo, typename _Ordering>
class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> > class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
{ {
public: public:
typedef _MatrixType MatrixType; typedef _MatrixType MatrixType;
@ -612,15 +618,13 @@ void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixTy
{ {
eigen_assert(a.rows()==a.cols()); eigen_assert(a.rows()==a.cols());
const Index size = a.rows(); const Index size = a.rows();
// TODO allows to configure the permutation
// Note that amd compute the inverse permutation // Note that amd compute the inverse permutation
{ {
CholMatrixType C; CholMatrixType C;
C = a.template selfadjointView<UpLo>(); C = a.template selfadjointView<UpLo>();
// remove diagonal entries:
// seems not to be needed OrderingType ordering;
// C.prune(keep_diag()); ordering(C,m_Pinv);
internal::minimum_degree_ordering(C, m_Pinv);
} }
if(m_Pinv.size()>0) if(m_Pinv.size()>0)

8
eigenlib/Eigen/src/SparseCore/CompressedStorage.h
View File

@ -51,8 +51,8 @@ class CompressedStorage
CompressedStorage& operator=(const CompressedStorage& other) CompressedStorage& operator=(const CompressedStorage& other)
{ {
resize(other.size()); resize(other.size());
memcpy(m_values, other.m_values, m_size * sizeof(Scalar)); internal::smart_copy(other.m_values, other.m_values + m_size, m_values);
memcpy(m_indices, other.m_indices, m_size * sizeof(Index)); internal::smart_copy(other.m_indices, other.m_indices + m_size, m_indices);
return *this; return *this;
} }
@ -83,10 +83,10 @@ class CompressedStorage
reallocate(m_size); reallocate(m_size);
} }
void resize(size_t size, float reserveSizeFactor = 0) void resize(size_t size, double reserveSizeFactor = 0)
{ {
if (m_allocatedSize<size) if (m_allocatedSize<size)
reallocate(size + size_t(reserveSizeFactor*size)); reallocate(size + size_t(reserveSizeFactor*double(size)));
m_size = size; m_size = size;
} }

3
eigenlib/Eigen/src/SparseCore/SparseCwiseBinaryOp.h
View File

@ -73,7 +73,8 @@ class CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::InnerIterator
typedef internal::sparse_cwise_binary_op_inner_iterator_selector< typedef internal::sparse_cwise_binary_op_inner_iterator_selector<
BinaryOp,Lhs,Rhs, InnerIterator> Base; 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) : Base(binOp.derived(),outer)
{} {}
}; };

33
eigenlib/Eigen/src/SparseCore/SparseDenseProduct.h
View File

@ -19,7 +19,10 @@ template<typename Lhs, typename Rhs, int InnerSize> struct SparseDenseProductRet
template<typename Lhs, typename Rhs> struct SparseDenseProductReturnType<Lhs,Rhs,1> template<typename Lhs, typename Rhs> struct SparseDenseProductReturnType<Lhs,Rhs,1>
{ {
typedef SparseDenseOuterProduct<Lhs,Rhs,false> Type; typedef typename internal::conditional<
Lhs::IsRowMajor,
SparseDenseOuterProduct<Rhs,Lhs,true>,
SparseDenseOuterProduct<Lhs,Rhs,false> >::type Type;
}; };
template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductReturnType template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductReturnType
@ -29,7 +32,10 @@ template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductRet
template<typename Lhs, typename Rhs> struct DenseSparseProductReturnType<Lhs,Rhs,1> template<typename Lhs, typename Rhs> struct DenseSparseProductReturnType<Lhs,Rhs,1>
{ {
typedef SparseDenseOuterProduct<Rhs,Lhs,true> Type; typedef typename internal::conditional<
Rhs::IsRowMajor,
SparseDenseOuterProduct<Rhs,Lhs,true>,
SparseDenseOuterProduct<Lhs,Rhs,false> >::type Type;
}; };
namespace internal { namespace internal {
@ -114,17 +120,30 @@ class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::InnerIterator : public _LhsNes
typedef typename SparseDenseOuterProduct::Index Index; typedef typename SparseDenseOuterProduct::Index Index;
public: public:
EIGEN_STRONG_INLINE InnerIterator(const SparseDenseOuterProduct& prod, Index outer) 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<Rhs>::StorageKind() ))
{ { }
}
inline Index outer() const { return m_outer; } inline Index outer() const { return m_outer; }
inline Index row() const { return Transpose ? Base::row() : m_outer; } inline Index row() const { return Transpose ? m_outer : Base::index(); }
inline Index col() const { return Transpose ? m_outer : Base::row(); } inline Index col() const { return Transpose ? Base::index() : m_outer; }
inline Scalar value() const { return Base::value() * m_factor; } inline Scalar value() const { return Base::value() * m_factor; }
protected: 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; Index m_outer;
Scalar m_factor; Scalar m_factor;
}; };

12
eigenlib/Eigen/src/SparseCore/SparseMatrix.h
View File

@ -940,7 +940,7 @@ void set_from_triplets(const InputIterator& begin, const InputIterator& end, Spa
enum { IsRowMajor = SparseMatrixType::IsRowMajor }; enum { IsRowMajor = SparseMatrixType::IsRowMajor };
typedef typename SparseMatrixType::Scalar Scalar; typedef typename SparseMatrixType::Scalar Scalar;
typedef typename SparseMatrixType::Index Index; typedef typename SparseMatrixType::Index Index;
SparseMatrix<Scalar,IsRowMajor?ColMajor:RowMajor> trMat(mat.rows(),mat.cols()); SparseMatrix<Scalar,IsRowMajor?ColMajor:RowMajor,Index> trMat(mat.rows(),mat.cols());
if(begin!=end) if(begin!=end)
{ {
@ -1178,7 +1178,7 @@ EIGEN_DONT_INLINE typename SparseMatrix<_Scalar,_Options,_Index>::Scalar& Sparse
size_t p = m_outerIndex[outer+1]; size_t p = m_outerIndex[outer+1];
++m_outerIndex[outer+1]; ++m_outerIndex[outer+1];
float reallocRatio = 1; double reallocRatio = 1;
if (m_data.allocatedSize()<=m_data.size()) if (m_data.allocatedSize()<=m_data.size())
{ {
// if there is no preallocated memory, let's reserve a minimum of 32 elements // 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 need to reallocate the data, to reduce multiple reallocations
// we use a smart resize algorithm based on the current filling ratio // we use a smart resize algorithm based on the current filling ratio
// in addition, we use float to avoid integers overflows // in addition, we use double to avoid integers overflows
float nnzEstimate = float(m_outerIndex[outer])*float(m_outerSize)/float(outer+1); double nnzEstimate = double(m_outerIndex[outer])*double(m_outerSize)/double(outer+1);
reallocRatio = (nnzEstimate-float(m_data.size()))/float(m_data.size()); reallocRatio = (nnzEstimate-double(m_data.size()))/double(m_data.size());
// furthermore we bound the realloc ratio to: // furthermore we bound the realloc ratio to:
// 1) reduce multiple minor realloc when the matrix is almost filled // 1) reduce multiple minor realloc when the matrix is almost filled
// 2) avoid to allocate too much memory when the matrix is almost empty // 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); m_data.resize(m_data.size()+1,reallocRatio);

10
eigenlib/Eigen/src/SparseCore/SparseTranspose.h
View File

@ -26,7 +26,7 @@ template<typename MatrixType> class TransposeImpl<MatrixType,Sparse>
inline Index nonZeros() const { return derived().nestedExpression().nonZeros(); } inline Index nonZeros() const { return derived().nestedExpression().nonZeros(); }
}; };
// NOTE: VC10 trigger an ICE if don't put typename TransposeImpl<MatrixType,Sparse>:: in front of Index, // NOTE: VC10 and VC11 trigger an ICE if don't put typename TransposeImpl<MatrixType,Sparse>:: in front of Index,
// a typedef typename TransposeImpl<MatrixType,Sparse>::Index Index; // a typedef typename TransposeImpl<MatrixType,Sparse>::Index Index;
// does not fix the issue. // does not fix the issue.
// An alternative is to define the nested class in the parent class itself. // An alternative is to define the nested class in the parent class itself.
@ -40,8 +40,8 @@ template<typename MatrixType> class TransposeImpl<MatrixType,Sparse>::InnerItera
EIGEN_STRONG_INLINE InnerIterator(const TransposeImpl& trans, typename TransposeImpl<MatrixType,Sparse>::Index outer) EIGEN_STRONG_INLINE InnerIterator(const TransposeImpl& trans, typename TransposeImpl<MatrixType,Sparse>::Index outer)
: Base(trans.derived().nestedExpression(), outer) : Base(trans.derived().nestedExpression(), outer)
{} {}
Index row() const { return Base::col(); } typename TransposeImpl<MatrixType,Sparse>::Index row() const { return Base::col(); }
Index col() const { return Base::row(); } typename TransposeImpl<MatrixType,Sparse>::Index col() const { return Base::row(); }
}; };
template<typename MatrixType> class TransposeImpl<MatrixType,Sparse>::ReverseInnerIterator template<typename MatrixType> class TransposeImpl<MatrixType,Sparse>::ReverseInnerIterator
@ -54,8 +54,8 @@ template<typename MatrixType> class TransposeImpl<MatrixType,Sparse>::ReverseInn
EIGEN_STRONG_INLINE ReverseInnerIterator(const TransposeImpl& xpr, typename TransposeImpl<MatrixType,Sparse>::Index outer) EIGEN_STRONG_INLINE ReverseInnerIterator(const TransposeImpl& xpr, typename TransposeImpl<MatrixType,Sparse>::Index outer)
: Base(xpr.derived().nestedExpression(), outer) : Base(xpr.derived().nestedExpression(), outer)
{} {}
Index row() const { return Base::col(); } typename TransposeImpl<MatrixType,Sparse>::Index row() const { return Base::col(); }
Index col() const { return Base::row(); } typename TransposeImpl<MatrixType,Sparse>::Index col() const { return Base::row(); }
}; };
} // end namespace Eigen } // end namespace Eigen

6
eigenlib/Eigen/src/SparseCore/SparseUtil.h
View File

@ -84,8 +84,10 @@ template<typename Lhs, typename Rhs> class DenseTimeSparseProduct;
template<typename Lhs, typename Rhs, bool Transpose> class SparseDenseOuterProduct; template<typename Lhs, typename Rhs, bool Transpose> class SparseDenseOuterProduct;
template<typename Lhs, typename Rhs> struct SparseSparseProductReturnType; template<typename Lhs, typename Rhs> struct SparseSparseProductReturnType;
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct DenseSparseProductReturnType; template<typename Lhs, typename Rhs,
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct SparseDenseProductReturnType; int InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(internal::traits<Lhs>::ColsAtCompileTime,internal::traits<Rhs>::RowsAtCompileTime)> struct DenseSparseProductReturnType;
template<typename Lhs, typename Rhs,
int InnerSize = EIGEN_SIZE_MIN_PREFER_FIXED(internal::traits<Lhs>::ColsAtCompileTime,internal::traits<Rhs>::RowsAtCompileTime)> struct SparseDenseProductReturnType;
template<typename MatrixType,int UpLo> class SparseSymmetricPermutationProduct; template<typename MatrixType,int UpLo> class SparseSymmetricPermutationProduct;
namespace internal { namespace internal {

136
eigenlib/Eigen/src/SparseQR/SparseQR.h
View File

@ -2,7 +2,7 @@
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> // Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
// Copyright (C) 2012-2013 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr>
// //
// This Source Code Form is subject to the terms of the Mozilla // 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 // 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 * \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. * 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<typename _MatrixType, typename _OrderingType> template<typename _MatrixType, typename _OrderingType>
@ -77,10 +78,23 @@ class SparseQR
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)
{ } { }
/** 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) SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false)
{ {
compute(mat); 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) void compute(const MatrixType& mat)
{ {
analyzePattern(mat); analyzePattern(mat);
@ -166,7 +180,7 @@ class SparseQR
y.bottomRows(y.rows()-rank).setZero(); y.bottomRows(y.rows()-rank).setZero();
// Apply the column permutation // 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()); else dest = y.topRows(cols());
m_info = Success; m_info = Success;
@ -206,7 +220,7 @@ class SparseQR
/** \brief Reports whether previous computation was successful. /** \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 NumericalIssue if the QR factorization reports a numerical problem
* \c InvalidInput if the input matrix is invalid * \c InvalidInput if the input matrix is invalid
* *
@ -255,20 +269,24 @@ class SparseQR
}; };
/** \brief Preprocessing step of a QR factorization /** \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 * 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. * \note In this step it is assumed that there is no empty row in the matrix \a mat.
*/ */
template <typename MatrixType, typename OrderingType> template <typename MatrixType, typename OrderingType>
void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
{ {
eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR");
// Compute the column fill reducing ordering // Compute the column fill reducing ordering
OrderingType ord; OrderingType ord;
ord(mat, m_perm_c); ord(mat, m_perm_c);
Index n = mat.cols(); Index n = mat.cols();
Index m = mat.rows(); Index m = mat.rows();
Index diagSize = (std::min)(m,n);
if (!m_perm_c.size()) if (!m_perm_c.size())
{ {
@ -280,20 +298,20 @@ void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat)
m_outputPerm_c = m_perm_c.inverse(); m_outputPerm_c = m_perm_c.inverse();
internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data());
m_R.resize(n, n); m_R.resize(m, n);
m_Q.resize(m, n); m_Q.resize(m, diagSize);
// Allocate space for nonzero elements : rough estimation // 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_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree
m_Q.reserve(2*mat.nonZeros()); m_Q.reserve(2*mat.nonZeros());
m_hcoeffs.resize(n); m_hcoeffs.resize(diagSize);
m_analysisIsok = true; m_analysisIsok = true;
} }
/** \brief Performs the numerical QR factorization of the input matrix /** \brief Performs the numerical QR factorization of the input matrix
* *
* The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with * 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 * \param mat The sparse column-major matrix
*/ */
@ -306,11 +324,12 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step");
Index m = mat.rows(); Index m = mat.rows();
Index n = mat.cols(); Index n = mat.cols();
IndexVector mark(m); mark.setConstant(-1); // Record the visited nodes Index diagSize = (std::min)(m,n);
IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes
Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q
ScalarVector tval(m); // The dense vector used to compute the current column Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q
bool found_diag; ScalarVector tval(m); // The dense vector used to compute the current column
RealScalar pivotThreshold = m_threshold;
m_pmat = mat; m_pmat = mat;
m_pmat.uncompress(); // To have the innerNonZeroPtr allocated m_pmat.uncompress(); // To have the innerNonZeroPtr allocated
@ -322,7 +341,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i];
} }
/* Compute the default threshold, see : /* Compute the default threshold as in MatLab, see:
* Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing
* Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3
*/ */
@ -330,24 +349,24 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
{ {
RealScalar max2Norm = 0.0; RealScalar max2Norm = 0.0;
for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm()); for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm());
m_threshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon(); pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon();
} }
// Initialize the numerical permutation // Initialize the numerical permutation
m_pivotperm.setIdentity(n); m_pivotperm.setIdentity(n);
Index nonzeroCol = 0; // Record the number of valid pivots 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 // 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); mark.setConstant(-1);
m_R.startVec(col); m_R.startVec(col);
m_Q.startVec(col);
mark(nonzeroCol) = col; mark(nonzeroCol) = col;
Qidx(0) = nonzeroCol; Qidx(0) = nonzeroCol;
nzcolR = 0; nzcolQ = 1; nzcolR = 0; nzcolQ = 1;
found_diag = col>=m; bool found_diag = nonzeroCol>=m;
tval.setZero(); tval.setZero();
// Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e., // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e.,
@ -356,7 +375,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
// thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. // 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 MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp)
{ {
Index curIdx = nonzeroCol ; Index curIdx = nonzeroCol;
if(itp) curIdx = itp.row(); if(itp) curIdx = itp.row();
if(curIdx == nonzeroCol) found_diag = true; if(curIdx == nonzeroCol) found_diag = true;
@ -398,7 +417,7 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
// Browse all the indexes of R(:,col) in reverse order // Browse all the indexes of R(:,col) in reverse order
for (Index i = nzcolR-1; i >= 0; i--) 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) // Apply the curIdx-th householder vector to the current column (temporarily stored into tval)
Scalar tdot(0); Scalar tdot(0);
@ -427,33 +446,37 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
} }
} }
} // End update current column } // End update current column
// Compute the Householder reflection that eliminate the current column
// FIXME this step should call the Householder module.
Scalar tau; Scalar tau;
RealScalar beta; RealScalar beta = 0;
Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
// First, the squared norm of Q((col+1):m, col) if(nonzeroCol < diagSize)
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))
{ {
tau = RealScalar(0); // Compute the Householder reflection that eliminate the current column
beta = numext::real(c0); // FIXME this step should call the Householder module.
tval(Qidx(0)) = 1; Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0);
}
else // First, the squared norm of Q((col+1):m, col)
{ RealScalar sqrNorm = 0.;
beta = std::sqrt(numext::abs2(c0) + sqrNorm); for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq)));
if(numext::real(c0) >= RealScalar(0)) if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0))
beta = -beta; {
tval(Qidx(0)) = 1; tau = RealScalar(0);
for (Index itq = 1; itq < nzcolQ; ++itq) beta = numext::real(c0);
tval(Qidx(itq)) /= (c0 - beta); tval(Qidx(0)) = 1;
tau = numext::conj((beta-c0) / beta); }
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 // Insert values in R
@ -467,24 +490,25 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
} }
} }
if(abs(beta) >= m_threshold) if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold)
{ {
m_R.insertBackByOuterInner(col, nonzeroCol) = beta; m_R.insertBackByOuterInner(col, nonzeroCol) = beta;
nonzeroCol++;
// The householder coefficient // The householder coefficient
m_hcoeffs(col) = tau; m_hcoeffs(nonzeroCol) = tau;
// Record the householder reflections // Record the householder reflections
for (Index itq = 0; itq < nzcolQ; ++itq) for (Index itq = 0; itq < nzcolQ; ++itq)
{ {
Index iQ = Qidx(itq); Index iQ = Qidx(itq);
m_Q.insertBackByOuterInnerUnordered(col,iQ) = tval(iQ); m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ);
tval(iQ) = Scalar(0.); tval(iQ) = Scalar(0.);
} }
nonzeroCol++;
if(nonzeroCol<diagSize)
m_Q.startVec(nonzeroCol);
} }
else else
{ {
// Zero pivot found: move implicitly this column to the end // Zero pivot found: move implicitly this column to the end
m_hcoeffs(col) = Scalar(0);
for (Index j = nonzeroCol; j < n-1; j++) for (Index j = nonzeroCol; j < n-1; j++)
std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]); std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
@ -493,6 +517,8 @@ void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat)
} }
} }
m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
// Finalize the column pointers of the sparse matrices R and Q // Finalize the column pointers of the sparse matrices R and Q
m_Q.finalize(); m_Q.finalize();
m_Q.makeCompressed(); m_Q.makeCompressed();
@ -561,14 +587,16 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
template<typename DesType> template<typename DesType>
void evalTo(DesType& res) const void evalTo(DesType& res) const
{ {
Index m = m_qr.rows();
Index n = m_qr.cols(); Index n = m_qr.cols();
Index diagSize = (std::min)(m,n);
res = m_other; res = m_other;
if (m_transpose) if (m_transpose)
{ {
eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
//Compute res = Q' * other column by column //Compute res = Q' * other column by column
for(Index j = 0; j < res.cols(); j++){ 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); Scalar tau = Scalar(0);
tau = m_qr.m_Q.col(k).dot(res.col(j)); tau = m_qr.m_Q.col(k).dot(res.col(j));
@ -581,10 +609,10 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
else else
{ {
eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes");
// Compute res = Q' * other column by column // Compute res = Q * other column by column
for(Index j = 0; j < res.cols(); j++) for(Index j = 0; j < res.cols(); j++)
{ {
for (Index k = n-1; k >=0; k--) for (Index k = diagSize-1; k >=0; k--)
{ {
Scalar tau = Scalar(0); Scalar tau = Scalar(0);
tau = m_qr.m_Q.col(k).dot(res.col(j)); tau = m_qr.m_Q.col(k).dot(res.col(j));
@ -618,7 +646,7 @@ struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<Sp
return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr);
} }
inline Index rows() const { return m_qr.rows(); } 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 // To use for operations with the transpose of Q
SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const
{ {

2
eigenlib/Eigen/src/StlSupport/StdDeque.h
View File

@ -11,7 +11,7 @@
#ifndef EIGEN_STDDEQUE_H #ifndef EIGEN_STDDEQUE_H
#define 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) // Define the explicit instantiation (e.g. necessary for the Intel compiler)
#if defined(__INTEL_COMPILER) || defined(__GNUC__) #if defined(__INTEL_COMPILER) || defined(__GNUC__)

2
eigenlib/Eigen/src/StlSupport/StdList.h
View File

@ -10,7 +10,7 @@
#ifndef EIGEN_STDLIST_H #ifndef EIGEN_STDLIST_H
#define 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) // Define the explicit instantiation (e.g. necessary for the Intel compiler)
#if defined(__INTEL_COMPILER) || defined(__GNUC__) #if defined(__INTEL_COMPILER) || defined(__GNUC__)

2
eigenlib/Eigen/src/StlSupport/StdVector.h
View File

@ -11,7 +11,7 @@
#ifndef EIGEN_STDVECTOR_H #ifndef EIGEN_STDVECTOR_H
#define 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 * This section contains a convenience MACRO which allows an easy specialization of

1
eigenlib/howto.txt
View File

@ -1,5 +1,6 @@
Current Eigen Version 3.1.2 (05.11.2012) Current Eigen Version 3.1.2 (05.11.2012)
Current Eigen Version 3.2.1 (26.02.2014) updated on 14/05/2014 Current Eigen Version 3.2.1 (26.02.2014) updated on 14/05/2014
Current Eigen Version 3.2.2 (04.08.2014) updated on 21/10/2014
To update the lib: To update the lib:
- download Eigen - download Eigen

12
eigenlib/unsupported/Eigen/src/IterativeSolvers/GMRES.h
View File

@ -2,7 +2,7 @@
// for linear algebra. // for linear algebra.
// //
// Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr> // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2012 Kolja Brix <brix@igpm.rwth-aaachen.de> // Copyright (C) 2012, 2014 Kolja Brix <brix@igpm.rwth-aaachen.de>
// //
// This Source Code Form is subject to the terms of the Mozilla // 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 // 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 p0 = rhs - mat*x;
VectorType r0 = precond.solve(p0); 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); 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); VectorType tau = VectorType::Zero(restart + 1);
std::vector < JacobiRotation < Scalar > > G(restart); std::vector < JacobiRotation < Scalar > > G(restart);
// generate first Householder vector // generate first Householder vector
VectorType e; VectorType e(m-1);
RealScalar beta; RealScalar beta;
r0.makeHouseholder(e, tau.coeffRef(0), beta); r0.makeHouseholder(e, tau.coeffRef(0), beta);
w(0)=(Scalar) beta; w(0)=(Scalar) beta;

45
eigenlib/unsupported/test/autodiff.cpp
View File

@ -127,46 +127,47 @@ template<typename Func> void forward_jacobian(const Func& f)
VERIFY_IS_APPROX(j, jref); VERIFY_IS_APPROX(j, jref);
} }
// TODO also check actual derivatives!
void test_autodiff_scalar() void test_autodiff_scalar()
{ {
std::cerr << foo<float>(1,2) << "\n"; Vector2f p = Vector2f::Random();
typedef AutoDiffScalar<Vector2f> AD; typedef AutoDiffScalar<Vector2f> AD;
AD ax(1,Vector2f::UnitX()); AD ax(p.x(),Vector2f::UnitX());
AD ay(2,Vector2f::UnitY()); AD ay(p.y(),Vector2f::UnitY());
AD res = foo<AD>(ax,ay); AD res = foo<AD>(ax,ay);
std::cerr << res.value() << " <> " VERIFY_IS_APPROX(res.value(), foo(p.x(),p.y()));
<< res.derivatives().transpose() << "\n\n";
} }
// TODO also check actual derivatives!
void test_autodiff_vector() void test_autodiff_vector()
{ {
std::cerr << foo<Vector2f>(Vector2f(1,2)) << "\n"; Vector2f p = Vector2f::Random();
typedef AutoDiffScalar<Vector2f> AD; typedef AutoDiffScalar<Vector2f> AD;
typedef Matrix<AD,2,1> VectorAD; typedef Matrix<AD,2,1> VectorAD;
VectorAD p(AD(1),AD(-1)); VectorAD ap = p.cast<AD>();
p.x().derivatives() = Vector2f::UnitX(); ap.x().derivatives() = Vector2f::UnitX();
p.y().derivatives() = Vector2f::UnitY(); ap.y().derivatives() = Vector2f::UnitY();
AD res = foo<VectorAD>(p); AD res = foo<VectorAD>(ap);
std::cerr << res.value() << " <> " VERIFY_IS_APPROX(res.value(), foo(p));
<< res.derivatives().transpose() << "\n\n";
} }
void test_autodiff_jacobian() void test_autodiff_jacobian()
{ {
for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,2>()) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,2>()) )); CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,3>()) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double,2,3>()) )); CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,2>()) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,2>()) )); CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,3>()) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double,3,3>()) )); CALL_SUBTEST(( forward_jacobian(TestFunc1<double>(3,3)) ));
CALL_SUBTEST(( forward_jacobian(TestFunc1<double>(3,3)) ));
}
} }
void test_autodiff() void test_autodiff()
{ {
test_autodiff_scalar(); for(int i = 0; i < g_repeat; i++) {
test_autodiff_vector(); CALL_SUBTEST_1( test_autodiff_scalar() );
// test_autodiff_jacobian(); CALL_SUBTEST_2( test_autodiff_vector() );
CALL_SUBTEST_3( test_autodiff_jacobian() );
}
} }