mirror of
https://github.com/ultimatepp/ultimatepp.git
synced 2026-07-11 22:03:01 -06:00
*Eigen: Updated
git-svn-id: svn://ultimatepp.org/upp/trunk@8049 f0d560ea-af0d-0410-9eb7-867de7ffcac7
This commit is contained in:
koldo
parent
2ab4268f05
commit
5181211520
115 changed files with 797 additions and 724 deletions
|
|
@ -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)
|
||||
|
|
@ -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 <immintrin.h>
|
||||
#else
|
||||
#include <emmintrin.h>
|
||||
|
|
@ -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 <intrin.h>
|
||||
#endif
|
||||
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_Cholesky_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_Cholesky_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Cholesky COMPONENT Devel
|
||||
)
|
||||
|
|
@ -274,30 +274,13 @@ template<> struct ldlt_inplace<Lower>
|
|||
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<Scalar>::epsilon() * biggest_in_corner);
|
||||
}
|
||||
|
||||
// Finish early if the matrix is not full rank.
|
||||
if(biggest_in_corner < cutoff)
|
||||
{
|
||||
for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i;
|
||||
break;
|
||||
}
|
||||
|
||||
transpositions.coeffRef(k) = index_of_biggest_in_corner;
|
||||
if(k != index_of_biggest_in_corner)
|
||||
{
|
||||
|
|
@ -328,15 +311,20 @@ template<> struct ldlt_inplace<Lower>
|
|||
|
||||
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<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a)
|
|||
m_transpositions.resize(size);
|
||||
m_isInitialized = false;
|
||||
m_temporary.resize(size);
|
||||
m_sign = internal::ZeroSign;
|
||||
|
||||
internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign);
|
||||
|
||||
|
|
@ -514,16 +503,21 @@ struct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs>
|
|||
using std::abs;
|
||||
using std::max;
|
||||
typedef typename LDLTType::MatrixType MatrixType;
|
||||
typedef typename LDLTType::Scalar Scalar;
|
||||
typedef typename LDLTType::RealScalar RealScalar;
|
||||
const Diagonal<const MatrixType> vectorD = dec().vectorD();
|
||||
RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() * NumTraits<Scalar>::epsilon(),
|
||||
RealScalar(1) / NumTraits<RealScalar>::highest()); // motivated by LAPACK's xGELSS
|
||||
const typename Diagonal<const MatrixType>::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<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) {
|
||||
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<MatrixType,_UpLo>::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)
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_CholmodSupport_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_CholmodSupport_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/CholmodSupport COMPONENT Devel
|
||||
)
|
||||
|
|
@ -29,6 +29,11 @@ struct traits<ArrayWrapper<ExpressionType> >
|
|||
: public traits<typename remove_all<typename ExpressionType::Nested>::type >
|
||||
{
|
||||
typedef ArrayXpr XprKind;
|
||||
// Let's remove NestByRefBit
|
||||
enum {
|
||||
Flags0 = traits<typename remove_all<typename ExpressionType::Nested>::type >::Flags,
|
||||
Flags = Flags0 & ~NestByRefBit
|
||||
};
|
||||
};
|
||||
}
|
||||
|
||||
|
|
@ -149,6 +154,11 @@ struct traits<MatrixWrapper<ExpressionType> >
|
|||
: public traits<typename remove_all<typename ExpressionType::Nested>::type >
|
||||
{
|
||||
typedef MatrixXpr XprKind;
|
||||
// Let's remove NestByRefBit
|
||||
enum {
|
||||
Flags0 = traits<typename remove_all<typename ExpressionType::Nested>::type >::Flags,
|
||||
Flags = Flags0 & ~NestByRefBit
|
||||
};
|
||||
};
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -81,7 +81,7 @@ struct traits<Block<XprType, BlockRows, BlockCols, InnerPanel> > : traits<XprTyp
|
|||
&& (InnerStrideAtCompileTime == 1)
|
||||
? PacketAccessBit : 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,
|
||||
FlagsRowMajorBit = IsRowMajor ? RowMajorBit : 0,
|
||||
Flags0 = traits<XprType>::Flags & ( (HereditaryBits & ~RowMajorBit) |
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
@ -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<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 */
|
||||
CommaInitializer& operator,(const Scalar& s)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -462,8 +462,10 @@ template<typename Derived> class DenseBase
|
|||
template<int p> RealScalar lpNorm() const;
|
||||
|
||||
template<int RowFactor, int ColFactor>
|
||||
const Replicate<Derived,RowFactor,ColFactor> replicate() const;
|
||||
const Replicate<Derived,Dynamic,Dynamic> replicate(Index rowFacor,Index colFactor) const;
|
||||
inline const Replicate<Derived,RowFactor,ColFactor> replicate() const;
|
||||
|
||||
typedef Replicate<Derived,Dynamic,Dynamic> ReplicateReturnType;
|
||||
inline const ReplicateReturnType replicate(Index rowFacor,Index colFactor) const;
|
||||
|
||||
typedef Reverse<Derived, BothDirections> ReverseReturnType;
|
||||
typedef const Reverse<const Derived, BothDirections> ConstReverseReturnType;
|
||||
|
|
|
|||
|
|
@ -24,6 +24,14 @@ namespace internal {
|
|||
|
||||
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
|
||||
* 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<T,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<T,Size>();
|
||||
}
|
||||
};
|
||||
|
||||
|
|
@ -76,12 +84,12 @@ struct plain_array<T, Size, MatrixOrArrayOptions, 16>
|
|||
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<T,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<T,Size>();
|
||||
}
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -190,18 +190,18 @@ MatrixBase<Derived>::diagonal() const
|
|||
*
|
||||
* \sa MatrixBase::diagonal(), class Diagonal */
|
||||
template<typename Derived>
|
||||
inline typename MatrixBase<Derived>::template DiagonalIndexReturnType<DynamicIndex>::Type
|
||||
inline typename MatrixBase<Derived>::DiagonalDynamicIndexReturnType
|
||||
MatrixBase<Derived>::diagonal(Index index)
|
||||
{
|
||||
return typename DiagonalIndexReturnType<DynamicIndex>::Type(derived(), index);
|
||||
return DiagonalDynamicIndexReturnType(derived(), index);
|
||||
}
|
||||
|
||||
/** This is the const version of diagonal(Index). */
|
||||
template<typename Derived>
|
||||
inline typename MatrixBase<Derived>::template ConstDiagonalIndexReturnType<DynamicIndex>::Type
|
||||
inline typename MatrixBase<Derived>::ConstDiagonalDynamicIndexReturnType
|
||||
MatrixBase<Derived>::diagonal(Index index) const
|
||||
{
|
||||
return typename ConstDiagonalIndexReturnType<DynamicIndex>::Type(derived(), index);
|
||||
return ConstDiagonalDynamicIndexReturnType(derived(), index);
|
||||
}
|
||||
|
||||
/** \returns an expression of the \a DiagIndex-th sub or super diagonal of the matrix \c *this
|
||||
|
|
|
|||
|
|
@ -589,7 +589,7 @@ struct linspaced_op_impl<Scalar,true>
|
|||
|
||||
template<typename Index>
|
||||
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_step;
|
||||
|
|
@ -609,7 +609,7 @@ template <typename Scalar, bool RandomAccess> struct functor_traits< linspaced_o
|
|||
template <typename Scalar, bool RandomAccess> struct linspaced_op
|
||||
{
|
||||
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>
|
||||
EIGEN_STRONG_INLINE const Scalar operator() (Index i) const { return impl(i); }
|
||||
|
|
|
|||
|
|
@ -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<cols; ++j)
|
||||
func(dest.col(j), prod.rhs().coeff(j) * prod.lhs());
|
||||
func(dest.col(j), prod.rhs().coeff(0,j) * prod.lhs());
|
||||
}
|
||||
|
||||
// Row major
|
||||
|
|
@ -243,7 +243,7 @@ EIGEN_DONT_INLINE void outer_product_selector_run(const ProductType& prod, Dest&
|
|||
// FIXME not very good if lhs is real and rhs complex while alpha is real too
|
||||
const Index rows = dest.rows();
|
||||
for (Index i=0; i<rows; ++i)
|
||||
func(dest.row(i), prod.lhs().coeff(i) * prod.rhs());
|
||||
func(dest.row(i), prod.lhs().coeff(i,0) * prod.rhs());
|
||||
}
|
||||
|
||||
template<typename Lhs, typename Rhs>
|
||||
|
|
|
|||
|
|
@ -234,9 +234,15 @@ template<typename Derived> class MapBase<Derived, WriteAccessors>
|
|||
return derived();
|
||||
}
|
||||
|
||||
using Base::Base::operator=;
|
||||
// In theory MapBase<Derived, ReadOnlyAccessors> 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<Derived, ReadOnlyAccessors>
|
||||
// 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
|
||||
|
|
|
|||
|
|
@ -215,7 +215,7 @@ template<typename Derived> class MatrixBase
|
|||
|
||||
typedef Diagonal<Derived> DiagonalReturnType;
|
||||
DiagonalReturnType diagonal();
|
||||
typedef typename internal::add_const<Diagonal<const Derived> >::type ConstDiagonalReturnType;
|
||||
typedef typename internal::add_const<Diagonal<const Derived> >::type ConstDiagonalReturnType;
|
||||
ConstDiagonalReturnType diagonal() const;
|
||||
|
||||
template<int Index> struct DiagonalIndexReturnType { typedef Diagonal<Derived,Index> Type; };
|
||||
|
|
@ -223,16 +223,12 @@ template<typename Derived> class MatrixBase
|
|||
|
||||
template<int Index> typename DiagonalIndexReturnType<Index>::Type diagonal();
|
||||
template<int Index> typename ConstDiagonalIndexReturnType<Index>::Type diagonal() const;
|
||||
|
||||
typedef Diagonal<Derived,DynamicIndex> DiagonalDynamicIndexReturnType;
|
||||
typedef typename internal::add_const<Diagonal<const Derived,DynamicIndex> >::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<DynamicIndex>::Type diagonal(Index index);
|
||||
typename MatrixBase::template ConstDiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index) const;
|
||||
#else
|
||||
typename DiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index);
|
||||
typename ConstDiagonalIndexReturnType<DynamicIndex>::Type diagonal(Index index) const;
|
||||
#endif
|
||||
DiagonalDynamicIndexReturnType diagonal(Index index);
|
||||
ConstDiagonalDynamicIndexReturnType diagonal(Index index) const;
|
||||
|
||||
#ifdef EIGEN2_SUPPORT
|
||||
template<unsigned int Mode> typename internal::eigen2_part_return_type<Derived, Mode>::type part();
|
||||
|
|
|
|||
|
|
@ -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<const Matrix>, Block<Matrix> >::value should be true.
|
||||
if(is_same<MatrixTypeNestedCleaned,Dest>::value && extract_data(dst) == extract_data(m_matrix))
|
||||
if( is_same<MatrixTypeNestedCleaned,Dest>::value
|
||||
&& blas_traits<MatrixTypeNestedCleaned>::HasUsableDirectAccess
|
||||
&& blas_traits<Dest>::HasUsableDirectAccess
|
||||
&& extract_data(dst) == extract_data(m_matrix))
|
||||
{
|
||||
// apply the permutation inplace
|
||||
Matrix<bool,PermutationType::RowsAtCompileTime,1,0,PermutationType::MaxRowsAtCompileTime> mask(m_permutation.size());
|
||||
|
|
|
|||
|
|
@ -85,7 +85,14 @@ class ProductBase : public MatrixBase<Derived>
|
|||
|
||||
public:
|
||||
|
||||
#ifndef EIGEN_NO_MALLOC
|
||||
typedef typename Base::PlainObject BasePlainObject;
|
||||
typedef Matrix<Scalar,RowsAtCompileTime==1?1:Dynamic,ColsAtCompileTime==1?1:Dynamic,BasePlainObject::Options> 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<typename Lhs, typename Rhs, int Mode, int N, typename PlainObject>
|
||||
struct nested<GeneralProduct<Lhs,Rhs,Mode>, N, PlainObject>
|
||||
{
|
||||
typedef PlainObject const& type;
|
||||
typedef typename GeneralProduct<Lhs,Rhs,Mode>::PlainObject const& type;
|
||||
};
|
||||
template<typename Lhs, typename Rhs, int Mode, int N, typename PlainObject>
|
||||
struct nested<const GeneralProduct<Lhs,Rhs,Mode>, N, PlainObject>
|
||||
{
|
||||
typedef typename GeneralProduct<Lhs,Rhs,Mode>::PlainObject const& type;
|
||||
};
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -101,7 +101,7 @@ struct traits<Ref<_PlainObjectType, _Options, _StrideType> >
|
|||
template<typename Derived> struct match {
|
||||
enum {
|
||||
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)
|
||||
|| int(StrideType::InnerStrideAtCompileTime)==int(Derived::InnerStrideAtCompileTime)
|
||||
|| (int(StrideType::InnerStrideAtCompileTime)==0 && int(Derived::InnerStrideAtCompileTime)==1),
|
||||
|
|
@ -172,8 +172,12 @@ protected:
|
|||
}
|
||||
else
|
||||
::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;
|
||||
|
|
@ -184,6 +188,8 @@ template<typename PlainObjectType, int Options, typename StrideType> class Ref
|
|||
: public RefBase<Ref<PlainObjectType, Options, StrideType> >
|
||||
{
|
||||
typedef internal::traits<Ref> Traits;
|
||||
template<typename Derived>
|
||||
inline Ref(const PlainObjectBase<Derived>& expr);
|
||||
public:
|
||||
|
||||
typedef RefBase<Ref> Base;
|
||||
|
|
@ -192,20 +198,21 @@ template<typename PlainObjectType, int Options, typename StrideType> class Ref
|
|||
|
||||
#ifndef EIGEN_PARSED_BY_DOXYGEN
|
||||
template<typename Derived>
|
||||
inline Ref(PlainObjectBase<Derived>& expr,
|
||||
typename internal::enable_if<bool(Traits::template match<Derived>::MatchAtCompileTime),Derived>::type* = 0)
|
||||
inline Ref(PlainObjectBase<Derived>& expr)
|
||||
{
|
||||
Base::construct(expr);
|
||||
EIGEN_STATIC_ASSERT(static_cast<bool>(Traits::template match<Derived>::MatchAtCompileTime), STORAGE_LAYOUT_DOES_NOT_MATCH);
|
||||
Base::construct(expr.derived());
|
||||
}
|
||||
template<typename Derived>
|
||||
inline Ref(const DenseBase<Derived>& expr,
|
||||
typename internal::enable_if<bool(internal::is_lvalue<Derived>::value&&bool(Traits::template match<Derived>::MatchAtCompileTime)),Derived>::type* = 0,
|
||||
int = Derived::ThisConstantIsPrivateInPlainObjectBase)
|
||||
inline Ref(const DenseBase<Derived>& expr)
|
||||
#else
|
||||
template<typename Derived>
|
||||
inline Ref(DenseBase<Derived>& expr)
|
||||
#endif
|
||||
{
|
||||
EIGEN_STATIC_ASSERT(static_cast<bool>(internal::is_lvalue<Derived>::value), THIS_EXPRESSION_IS_NOT_A_LVALUE__IT_IS_READ_ONLY);
|
||||
EIGEN_STATIC_ASSERT(static_cast<bool>(Traits::template match<Derived>::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());
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -135,7 +135,7 @@ template<typename MatrixType,int RowFactor,int ColFactor> class Replicate
|
|||
*/
|
||||
template<typename Derived>
|
||||
template<int RowFactor, int ColFactor>
|
||||
inline const Replicate<Derived,RowFactor,ColFactor>
|
||||
const Replicate<Derived,RowFactor,ColFactor>
|
||||
DenseBase<Derived>::replicate() const
|
||||
{
|
||||
return Replicate<Derived,RowFactor,ColFactor>(derived());
|
||||
|
|
@ -150,7 +150,7 @@ DenseBase<Derived>::replicate() const
|
|||
* \sa VectorwiseOp::replicate(), DenseBase::replicate<int,int>(), class Replicate
|
||||
*/
|
||||
template<typename Derived>
|
||||
inline const Replicate<Derived,Dynamic,Dynamic>
|
||||
const typename DenseBase<Derived>::ReplicateReturnType
|
||||
DenseBase<Derived>::replicate(Index rowFactor,Index colFactor) const
|
||||
{
|
||||
return Replicate<Derived,Dynamic,Dynamic>(derived(),rowFactor,colFactor);
|
||||
|
|
|
|||
|
|
@ -278,21 +278,21 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
|
|||
|
||||
/** Efficient triangular matrix times vector/matrix product */
|
||||
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
|
||||
{
|
||||
return TriangularProduct
|
||||
<Mode,true,MatrixType,false,OtherDerived,OtherDerived::IsVectorAtCompileTime>
|
||||
<Mode, true, MatrixType, false, OtherDerived, OtherDerived::ColsAtCompileTime==1>
|
||||
(m_matrix, rhs.derived());
|
||||
}
|
||||
|
||||
/** Efficient vector/matrix times triangular matrix product */
|
||||
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)
|
||||
{
|
||||
return TriangularProduct
|
||||
<Mode,false,OtherDerived,OtherDerived::IsVectorAtCompileTime,MatrixType,false>
|
||||
<Mode, false, OtherDerived, OtherDerived::RowsAtCompileTime==1, MatrixType, false>
|
||||
(lhs.derived(),rhs.m_matrix);
|
||||
}
|
||||
|
||||
|
|
@ -380,19 +380,19 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
|
|||
EIGEN_STRONG_INLINE TriangularView& operator=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
|
||||
{
|
||||
setZero();
|
||||
return assignProduct(other,1);
|
||||
return assignProduct(other.derived(),1);
|
||||
}
|
||||
|
||||
template<typename ProductDerived, typename Lhs, typename Rhs>
|
||||
EIGEN_STRONG_INLINE TriangularView& operator+=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
|
||||
{
|
||||
return assignProduct(other,1);
|
||||
return assignProduct(other.derived(),1);
|
||||
}
|
||||
|
||||
template<typename ProductDerived, typename Lhs, typename Rhs>
|
||||
EIGEN_STRONG_INLINE TriangularView& operator-=(const ProductBase<ProductDerived, Lhs,Rhs>& other)
|
||||
{
|
||||
return assignProduct(other,-1);
|
||||
return assignProduct(other.derived(),-1);
|
||||
}
|
||||
|
||||
|
||||
|
|
@ -400,25 +400,34 @@ template<typename _MatrixType, unsigned int _Mode> class TriangularView
|
|||
EIGEN_STRONG_INLINE TriangularView& operator=(const ScaledProduct<ProductDerived>& other)
|
||||
{
|
||||
setZero();
|
||||
return assignProduct(other,other.alpha());
|
||||
return assignProduct(other.derived(),other.alpha());
|
||||
}
|
||||
|
||||
template<typename ProductDerived>
|
||||
EIGEN_STRONG_INLINE TriangularView& operator+=(const ScaledProduct<ProductDerived>& other)
|
||||
{
|
||||
return assignProduct(other,other.alpha());
|
||||
return assignProduct(other.derived(),other.alpha());
|
||||
}
|
||||
|
||||
template<typename ProductDerived>
|
||||
EIGEN_STRONG_INLINE TriangularView& operator-=(const ScaledProduct<ProductDerived>& other)
|
||||
{
|
||||
return assignProduct(other,-other.alpha());
|
||||
return assignProduct(other.derived(),-other.alpha());
|
||||
}
|
||||
|
||||
protected:
|
||||
|
||||
template<typename ProductDerived, typename Lhs, typename Rhs>
|
||||
EIGEN_STRONG_INLINE TriangularView& assignProduct(const ProductBase<ProductDerived, Lhs,Rhs>& prod, const Scalar& alpha);
|
||||
|
||||
template<int Mode, bool LhsIsTriangular,
|
||||
typename Lhs, bool LhsIsVector,
|
||||
typename Rhs, bool RhsIsVector>
|
||||
EIGEN_STRONG_INLINE TriangularView& assignProduct(const TriangularProduct<Mode, LhsIsTriangular, Lhs, LhsIsVector, Rhs, RhsIsVector>& prod, const Scalar& alpha)
|
||||
{
|
||||
lazyAssign(alpha*prod.eval());
|
||||
return *this;
|
||||
}
|
||||
|
||||
MatrixTypeNested m_matrix;
|
||||
};
|
||||
|
|
|
|||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -1,4 +0,0 @@
|
|||
ADD_SUBDIRECTORY(SSE)
|
||||
ADD_SUBDIRECTORY(AltiVec)
|
||||
ADD_SUBDIRECTORY(NEON)
|
||||
ADD_SUBDIRECTORY(Default)
|
||||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -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<float> : default_packet_traits
|
||||
|
|
|
|||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -52,7 +52,7 @@ Packet4f plog<Packet4f>(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<Packet4f>(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<Packet2d>(const Packet2d& _x)
|
||||
|
|
@ -239,7 +239,7 @@ Packet2d pexp<Packet2d>(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.
|
||||
|
|
|
|||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -90,6 +90,7 @@ struct traits<CoeffBasedProduct<LhsNested,RhsNested,NestingFlags> >
|
|||
| (SameType && (CanVectorizeLhs || CanVectorizeRhs) ? PacketAccessBit : 0),
|
||||
|
||||
CoeffReadCost = InnerSize == Dynamic ? Dynamic
|
||||
: InnerSize == 0 ? 0
|
||||
: InnerSize * (NumTraits<Scalar>::MulCost + LhsCoeffReadCost + RhsCoeffReadCost)
|
||||
+ (InnerSize - 1) * NumTraits<Scalar>::AddCost,
|
||||
|
||||
|
|
@ -133,7 +134,7 @@ class CoeffBasedProduct
|
|||
};
|
||||
|
||||
typedef internal::product_coeff_impl<CanVectorizeInner ? InnerVectorizedTraversal : DefaultTraversal,
|
||||
Unroll ? InnerSize-1 : Dynamic,
|
||||
Unroll ? (InnerSize==0 ? 0 : InnerSize-1) : Dynamic,
|
||||
_LhsNested, _RhsNested, Scalar> ScalarCoeffImpl;
|
||||
|
||||
typedef CoeffBasedProduct<LhsNested,RhsNested,NestByRefBit> LazyCoeffBasedProductType;
|
||||
|
|
@ -184,7 +185,7 @@ class CoeffBasedProduct
|
|||
{
|
||||
PacketScalar res;
|
||||
internal::product_packet_impl<Flags&RowMajorBit ? RowMajor : ColMajor,
|
||||
Unroll ? InnerSize-1 : Dynamic,
|
||||
Unroll ? (InnerSize==0 ? 0 : InnerSize-1) : Dynamic,
|
||||
_LhsNested, _RhsNested, PacketScalar, LoadMode>
|
||||
::run(row, col, m_lhs, m_rhs, res);
|
||||
return res;
|
||||
|
|
@ -262,10 +263,7 @@ struct product_coeff_impl<DefaultTraversal, Dynamic, Lhs, Rhs, RetScalar>
|
|||
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();
|
||||
}
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -54,8 +54,25 @@
|
|||
#endif
|
||||
|
||||
#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>
|
||||
#define EIGEN_MKL_VML_THRESHOLD 128
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -63,7 +63,7 @@
|
|||
// Currently, let's include it only on unix systems:
|
||||
#if defined(__unix__) || defined(__unix)
|
||||
#include <unistd.h>
|
||||
#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<typename T, bool Align> inline T* conditional_aligned_realloc_new(T* pt
|
|||
|
||||
template<typename T, bool Align> inline T* conditional_aligned_new_auto(size_t size)
|
||||
{
|
||||
if(size==0)
|
||||
return 0; // short-cut. Also fixes Bug 884
|
||||
check_size_for_overflow<T>(size);
|
||||
T *result = reinterpret_cast<T*>(conditional_aligned_malloc<Align>(sizeof(T)*size));
|
||||
if(NumTraits<T>::RequireInitialization)
|
||||
|
|
@ -612,7 +614,6 @@ template<typename T> class aligned_stack_memory_handler
|
|||
void* operator new(size_t size, const std::nothrow_t&) throw() { \
|
||||
try { return Eigen::internal::conditional_aligned_malloc<NeedsToAlign>(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<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)
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
};
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -341,7 +341,7 @@ template<typename T, int n=1, typename PlainObject = typename eval<T>::type> str
|
|||
};
|
||||
|
||||
template<typename T>
|
||||
T* const_cast_ptr(const T* ptr)
|
||||
inline T* const_cast_ptr(const T* ptr)
|
||||
{
|
||||
return const_cast<T*>(ptr);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -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();
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_EIGENVALUES_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_EIGENVALUES_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigenvalues COMPONENT Devel
|
||||
)
|
||||
|
|
@ -563,7 +563,6 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
|
|||
if(computeEigenvectors)
|
||||
{
|
||||
Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon();
|
||||
safeNorm2 *= safeNorm2;
|
||||
if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
|
||||
{
|
||||
eivecs.setIdentity();
|
||||
|
|
@ -577,7 +576,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
|
|||
Scalar d0 = eivals(2) - eivals(1);
|
||||
Scalar d1 = eivals(1) - eivals(0);
|
||||
int k = d0 > 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<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
|
|||
n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm();
|
||||
|
||||
if(n>safeNorm2)
|
||||
{
|
||||
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<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
|
|||
tmp.diagonal().array() -= eivals(1);
|
||||
|
||||
if(d0<=Eigen::NumTraits<Scalar>::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<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
|
|||
else
|
||||
{
|
||||
// we should never reach this point,
|
||||
// if so the last two eigenvalues are likely to ve very closed to each other
|
||||
// if so the last two eigenvalues are likely to be very close to each other
|
||||
eivecs.col(1) = eivecs.col(k).unitOrthogonal();
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// make sure that eivecs[1] is orthogonal to eivecs[2]
|
||||
// FIXME: this step should not be needed
|
||||
Scalar d = eivecs.col(1).dot(eivecs.col(k));
|
||||
eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized();
|
||||
}
|
||||
|
|
|
|||
|
|
@ -1,8 +0,0 @@
|
|||
FILE(GLOB Eigen_Geometry_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_Geometry_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Geometry COMPONENT Devel
|
||||
)
|
||||
|
||||
ADD_SUBDIRECTORY(arch)
|
||||
|
|
@ -100,7 +100,17 @@ public:
|
|||
{
|
||||
EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(VectorType, 3)
|
||||
Hyperplane result(p0.size());
|
||||
result.normal() = (p2 - p0).cross(p1 - p0).normalized();
|
||||
VectorType v0(p2 - p0), v1(p1 - p0);
|
||||
result.normal() = v0.cross(v1);
|
||||
RealScalar norm = result.normal().norm();
|
||||
if(norm <= v0.norm() * v1.norm() * NumTraits<RealScalar>::epsilon())
|
||||
{
|
||||
Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
|
||||
JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
|
||||
result.normal() = svd.matrixV().col(2);
|
||||
}
|
||||
else
|
||||
result.normal() /= norm;
|
||||
result.offset() = -p0.dot(result.normal());
|
||||
return result;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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<const Quaternion<_Scalar>, _Options >
|
|||
|
||||
/** Constructs a Mapped Quaternion object from the pointer \a coeffs
|
||||
*
|
||||
* The pointer \a coeffs must reference the four coeffecients of Quaternion in the following order:
|
||||
* The pointer \a coeffs must reference the four coefficients of Quaternion in the following order:
|
||||
* \code *coeffs == {x, y, z, w} \endcode
|
||||
*
|
||||
* If the template parameter _Options is set to #Aligned, then the pointer coeffs must be aligned. */
|
||||
|
|
@ -464,7 +466,7 @@ QuaternionBase<Derived>::_transformVector(Vector3 v) const
|
|||
// Note that this algorithm comes from the optimization by hand
|
||||
// of the conversion to a Matrix followed by a Matrix/Vector product.
|
||||
// It appears to be much faster than the common algorithm found
|
||||
// in the litterature (30 versus 39 flops). It also requires two
|
||||
// in the literature (30 versus 39 flops). It also requires two
|
||||
// Vector3 as temporaries.
|
||||
Vector3 uv = this->vec().cross(v);
|
||||
uv += uv;
|
||||
|
|
@ -584,7 +586,7 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri
|
|||
// which yields a singular value problem
|
||||
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();
|
||||
JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
|
||||
Vector3 axis = svd.matrixV().col(2);
|
||||
|
|
@ -667,10 +669,10 @@ QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& oth
|
|||
{
|
||||
using std::acos;
|
||||
using std::abs;
|
||||
double d = abs(this->dot(other));
|
||||
if (d>=1.0)
|
||||
Scalar d = abs(this->dot(other));
|
||||
if (d>=Scalar(1))
|
||||
return Scalar(0);
|
||||
return static_cast<Scalar>(2 * acos(d));
|
||||
return Scalar(2) * acos(d);
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -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<typename Derived>
|
||||
Rotation2D& fromRotationMatrix(const MatrixBase<Derived>& 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.
|
||||
|
|
|
|||
|
|
@ -62,6 +62,8 @@ struct transform_construct_from_matrix;
|
|||
|
||||
template<typename TransformType> struct transform_take_affine_part;
|
||||
|
||||
template<int Mode> 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<Scalar,Dim,Dim,Options> LinearMatrixType;
|
||||
/** 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 */
|
||||
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 */
|
||||
typedef typename internal::conditional<int(Mode)==int(AffineCompact),
|
||||
MatrixType&,
|
||||
|
|
@ -230,8 +232,7 @@ public:
|
|||
inline Transform()
|
||||
{
|
||||
check_template_params();
|
||||
if (int(Mode)==Affine)
|
||||
makeAffine();
|
||||
internal::transform_make_affine<(int(Mode)==Affine) ? Affine : AffineCompact>::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<int(Mode)>::run(m_matrix);
|
||||
}
|
||||
|
||||
/** \internal
|
||||
|
|
@ -1079,6 +1076,24 @@ Transform<Scalar,Dim,Mode,Options>::fromPositionOrientationScale(const MatrixBas
|
|||
|
||||
namespace internal {
|
||||
|
||||
template<int Mode>
|
||||
struct transform_make_affine
|
||||
{
|
||||
template<typename MatrixType>
|
||||
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<AffineCompact>
|
||||
{
|
||||
template<typename MatrixType> static void run(MatrixType &) { }
|
||||
};
|
||||
|
||||
// selector needed to avoid taking the inverse of a 3x4 matrix
|
||||
template<typename TransformType, int Mode=TransformType::Mode>
|
||||
struct projective_transform_inverse
|
||||
|
|
|
|||
|
|
@ -113,7 +113,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
|
|||
const Index n = src.cols(); // number of measurements
|
||||
|
||||
// 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
|
||||
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)
|
||||
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)
|
||||
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;
|
||||
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();
|
||||
} 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<Derived>& src, const MatrixBase<OtherDerived>& 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;
|
||||
|
|
|
|||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -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<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);
|
||||
|
||||
const TriangularView<const VectorsType, UnitLower>& V(vectors);
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_Householder_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_Householder_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Householder COMPONENT Devel
|
||||
)
|
||||
|
|
@ -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<Scalar>::epsilon()*NumTraits<Scalar>::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<BiCGSTAB<_MatrixType,_Preconditioner> >
|
|||
* SparseMatrix<double> A(n,n);
|
||||
* // fill A and b
|
||||
* BiCGSTAB<SparseMatrix<double> > 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;
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_IterativeLinearSolvers_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_IterativeLinearSolvers_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/IterativeLinearSolvers COMPONENT Devel
|
||||
)
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_Jacobi_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_Jacobi_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Jacobi COMPONENT Devel
|
||||
)
|
||||
|
|
@ -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)
|
||||
|
|
@ -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<MatrixType>::Scalar FullPivLU<MatrixType>::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<typename MatrixType>
|
||||
MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
|
||||
{
|
||||
|
|
|
|||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_MetisSupport_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_MetisSupport_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/MetisSupport COMPONENT Devel
|
||||
)
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_OrderingMethods_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_OrderingMethods_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/OrderingMethods COMPONENT Devel
|
||||
)
|
||||
|
|
@ -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<typename Index>
|
||||
class COLAMDOrdering
|
||||
|
|
@ -118,10 +118,14 @@ class COLAMDOrdering
|
|||
typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
|
||||
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>
|
||||
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);
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_PastixSupport_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_PastixSupport_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/PaStiXSupport COMPONENT Devel
|
||||
)
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_PardisoSupport_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_PardisoSupport_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/PardisoSupport COMPONENT Devel
|
||||
)
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_QR_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_QR_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/QR COMPONENT Devel
|
||||
)
|
||||
|
|
@ -76,7 +76,8 @@ template<typename _MatrixType> class ColPivHouseholderQR
|
|||
m_colsTranspositions(),
|
||||
m_temp(),
|
||||
m_colSqNorms(),
|
||||
m_isInitialized(false) {}
|
||||
m_isInitialized(false),
|
||||
m_usePrescribedThreshold(false) {}
|
||||
|
||||
/** \brief Default Constructor with memory preallocation
|
||||
*
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_SPQRSupport_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_SPQRSupport_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SPQRSupport/ COMPONENT Devel
|
||||
)
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_SVD_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_SVD_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SVD COMPONENT Devel
|
||||
)
|
||||
|
|
@ -375,17 +375,19 @@ struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
|
|||
Scalar z;
|
||||
JacobiRotation<Scalar> 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<RealScalar> *j_right)
|
||||
{
|
||||
using std::sqrt;
|
||||
using std::abs;
|
||||
Matrix<RealScalar,2,2> 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(t<RealScalar(0))
|
||||
rot1.s() = -rot1.s();
|
||||
}
|
||||
m.applyOnTheLeft(0,1,rot1);
|
||||
j_right->makeJacobi(m,0,1);
|
||||
|
|
@ -531,8 +536,9 @@ template<typename _MatrixType, int QRPreconditioner> 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<typename _MatrixType, int QRPreconditioner> 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<typename _MatrixType, int QRPreconditioner> 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<typename _MatrixType, int QRPreconditioner> 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<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 cols() const { return m_cols; }
|
||||
|
|
@ -677,11 +748,12 @@ template<typename _MatrixType, int QRPreconditioner> 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<typename __MatrixType, int _QRPreconditioner, bool _IsComplex>
|
||||
friend struct internal::svd_precondition_2x2_block_to_be_real;
|
||||
|
|
@ -690,6 +762,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
|
|||
|
||||
internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
|
||||
internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
|
||||
MatrixType m_scaledMatrix;
|
||||
};
|
||||
|
||||
template<typename MatrixType, int QRPreconditioner>
|
||||
|
|
@ -736,8 +809,9 @@ void JacobiSVD<MatrixType, QRPreconditioner>::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<typename MatrixType, int QRPreconditioner>
|
||||
|
|
@ -754,11 +828,21 @@ JacobiSVD<MatrixType, QRPreconditioner>::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<RealScalar>::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<MatrixType, QRPreconditioner>::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<MatrixType, QRPreconditioner>::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<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
|
|||
// So A^{-1} = V S^{-1} U^*
|
||||
|
||||
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 = 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
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_SparseCholesky_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_SparseCholesky_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseCholesky COMPONENT Devel
|
||||
)
|
||||
|
|
@ -37,6 +37,7 @@ class SimplicialCholeskyBase : internal::noncopyable
|
|||
{
|
||||
public:
|
||||
typedef typename internal::traits<Derived>::MatrixType MatrixType;
|
||||
typedef typename internal::traits<Derived>::OrderingType OrderingType;
|
||||
enum { UpLo = internal::traits<Derived>::UpLo };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
typedef typename MatrixType::RealScalar RealScalar;
|
||||
|
|
@ -240,15 +241,16 @@ class SimplicialCholeskyBase : internal::noncopyable
|
|||
RealScalar m_shiftScale;
|
||||
};
|
||||
|
||||
template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
|
||||
template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
|
||||
template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
|
||||
template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT;
|
||||
template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT;
|
||||
template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialCholesky;
|
||||
|
||||
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 _Ordering OrderingType;
|
||||
enum { UpLo = _UpLo };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
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(); }
|
||||
};
|
||||
|
||||
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 _Ordering OrderingType;
|
||||
enum { UpLo = _UpLo };
|
||||
typedef typename MatrixType::Scalar Scalar;
|
||||
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(); }
|
||||
};
|
||||
|
||||
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 _Ordering OrderingType;
|
||||
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 _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<typename _MatrixType, int _UpLo>
|
||||
class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
|
||||
template<typename _MatrixType, int _UpLo, typename _Ordering>
|
||||
class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
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<typename _MatrixType, int _UpLo>
|
||||
class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
|
||||
template<typename _MatrixType, int _UpLo, typename _Ordering>
|
||||
class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
public:
|
||||
typedef _MatrixType MatrixType;
|
||||
|
|
@ -467,8 +473,8 @@ public:
|
|||
*
|
||||
* \sa class SimplicialLDLT, class SimplicialLLT
|
||||
*/
|
||||
template<typename _MatrixType, int _UpLo>
|
||||
class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
|
||||
template<typename _MatrixType, int _UpLo, typename _Ordering>
|
||||
class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
|
||||
{
|
||||
public:
|
||||
typedef _MatrixType MatrixType;
|
||||
|
|
@ -612,15 +618,13 @@ void SimplicialCholeskyBase<Derived>::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<UpLo>();
|
||||
// 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)
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_SparseCore_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_SparseCore_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseCore COMPONENT Devel
|
||||
)
|
||||
|
|
@ -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<size)
|
||||
reallocate(size + size_t(reserveSizeFactor*size));
|
||||
reallocate(size + size_t(reserveSizeFactor*double(size)));
|
||||
m_size = size;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -68,6 +68,8 @@ public:
|
|||
const internal::variable_if_dynamic<Index, OuterSize> m_outerSize;
|
||||
|
||||
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl)
|
||||
private:
|
||||
Index nonZeros() const;
|
||||
};
|
||||
|
||||
|
||||
|
|
@ -82,6 +84,7 @@ class BlockImpl<SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true
|
|||
typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType;
|
||||
typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested;
|
||||
typedef Block<SparseMatrixType, BlockRows, BlockCols, true> BlockType;
|
||||
typedef Block<const SparseMatrixType, BlockRows, BlockCols, true> ConstBlockType;
|
||||
public:
|
||||
enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
|
||||
EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType)
|
||||
|
|
@ -245,6 +248,93 @@ public:
|
|||
|
||||
};
|
||||
|
||||
|
||||
template<typename _Scalar, int _Options, typename _Index, int BlockRows, int BlockCols>
|
||||
class BlockImpl<const SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true,Sparse>
|
||||
: public SparseMatrixBase<Block<const SparseMatrix<_Scalar, _Options, _Index>,BlockRows,BlockCols,true> >
|
||||
{
|
||||
typedef SparseMatrix<_Scalar, _Options, _Index> SparseMatrixType;
|
||||
typedef typename internal::remove_all<typename SparseMatrixType::Nested>::type _MatrixTypeNested;
|
||||
typedef Block<const SparseMatrixType, BlockRows, BlockCols, true> BlockType;
|
||||
public:
|
||||
enum { IsRowMajor = internal::traits<BlockType>::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<const Matrix<Index,OuterSize,1> >(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<Index, OuterSize> m_outerSize;
|
||||
|
||||
};
|
||||
|
||||
//----------
|
||||
|
||||
/** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
|
||||
|
|
|
|||
|
|
@ -73,7 +73,8 @@ class CwiseBinaryOpImpl<BinaryOp,Lhs,Rhs,Sparse>::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)
|
||||
{}
|
||||
};
|
||||
|
|
|
|||
|
|
@ -19,7 +19,10 @@ template<typename Lhs, typename Rhs, int InnerSize> struct SparseDenseProductRet
|
|||
|
||||
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
|
||||
|
|
@ -29,7 +32,10 @@ template<typename Lhs, typename Rhs, int InnerSize> struct DenseSparseProductRet
|
|||
|
||||
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 {
|
||||
|
|
@ -114,17 +120,30 @@ class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::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<Rhs>::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<typename Derived>
|
||||
template<typename OtherDerived>
|
||||
inline const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
|
||||
SparseMatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
|
||||
{
|
||||
return typename SparseDenseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
|
||||
}
|
||||
|
||||
} // end namespace Eigen
|
||||
|
||||
#endif // EIGEN_SPARSEDENSEPRODUCT_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<Scalar,IsRowMajor?ColMajor:RowMajor> trMat(mat.rows(),mat.cols());
|
||||
SparseMatrix<Scalar,IsRowMajor?ColMajor:RowMajor,Index> 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);
|
||||
|
|
|
|||
|
|
@ -358,7 +358,8 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived>
|
|||
/** sparse * dense (returns a dense object unless it is an outer product) */
|
||||
template<typename OtherDerived>
|
||||
const typename SparseDenseProductReturnType<Derived,OtherDerived>::Type
|
||||
operator*(const MatrixBase<OtherDerived> &other) const;
|
||||
operator*(const MatrixBase<OtherDerived> &other) const
|
||||
{ return typename SparseDenseProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); }
|
||||
|
||||
/** \returns an expression of P H P^-1 where H is the matrix represented by \c *this */
|
||||
SparseSymmetricPermutationProduct<Derived,Upper|Lower> twistedBy(const PermutationMatrix<Dynamic,Dynamic,Index>& perm) const
|
||||
|
|
|
|||
|
|
@ -61,7 +61,7 @@ struct permut_sparsematrix_product_retval
|
|||
for(Index j=0; j<m_matrix.outerSize(); ++j)
|
||||
{
|
||||
Index jp = m_permutation.indices().coeff(j);
|
||||
sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).size();
|
||||
sizes[((Side==OnTheLeft) ^ Transposed) ? jp : j] = m_matrix.innerVector(((Side==OnTheRight) ^ Transposed) ? jp : j).nonZeros();
|
||||
}
|
||||
tmp.reserve(sizes);
|
||||
for(Index j=0; j<m_matrix.outerSize(); ++j)
|
||||
|
|
|
|||
|
|
@ -26,7 +26,7 @@ template<typename MatrixType> class TransposeImpl<MatrixType,Sparse>
|
|||
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;
|
||||
// does not fix the issue.
|
||||
// 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)
|
||||
: Base(trans.derived().nestedExpression(), outer)
|
||||
{}
|
||||
Index row() const { return Base::col(); }
|
||||
Index col() const { return Base::row(); }
|
||||
typename TransposeImpl<MatrixType,Sparse>::Index row() const { return Base::col(); }
|
||||
typename TransposeImpl<MatrixType,Sparse>::Index col() const { return Base::row(); }
|
||||
};
|
||||
|
||||
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)
|
||||
: Base(xpr.derived().nestedExpression(), outer)
|
||||
{}
|
||||
Index row() const { return Base::col(); }
|
||||
Index col() const { return Base::row(); }
|
||||
typename TransposeImpl<MatrixType,Sparse>::Index row() const { return Base::col(); }
|
||||
typename TransposeImpl<MatrixType,Sparse>::Index col() const { return Base::row(); }
|
||||
};
|
||||
|
||||
} // end namespace Eigen
|
||||
|
|
|
|||
|
|
@ -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> struct SparseSparseProductReturnType;
|
||||
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct DenseSparseProductReturnType;
|
||||
template<typename Lhs, typename Rhs, int InnerSize = internal::traits<Lhs>::ColsAtCompileTime> struct SparseDenseProductReturnType;
|
||||
template<typename Lhs, typename Rhs,
|
||||
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;
|
||||
|
||||
namespace internal {
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_SparseLU_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_SparseLU_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseLU COMPONENT Devel
|
||||
)
|
||||
|
|
@ -260,14 +260,13 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ
|
|||
eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
|
||||
// Initialize with the determinant of the row matrix
|
||||
Scalar det = Scalar(1.);
|
||||
//Note that the diagonal blocks of U are stored in supernodes,
|
||||
// Note that the diagonal blocks of U are stored in supernodes,
|
||||
// which are available in the L part :)
|
||||
for (Index j = 0; j < this->cols(); ++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;
|
||||
|
|
|
|||
|
|
@ -189,8 +189,8 @@ class MappedSuperNodalMatrix<Scalar,Index>::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++()
|
||||
{
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_SparseQR_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_SparseQR_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SparseQR/ COMPONENT Devel
|
||||
)
|
||||
|
|
@ -2,7 +2,7 @@
|
|||
// for linear algebra.
|
||||
//
|
||||
// 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
|
||||
// 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<typename _MatrixType, typename _OrderingType>
|
||||
|
|
@ -74,13 +75,26 @@ class SparseQR
|
|||
typedef Matrix<Scalar, Dynamic, 1> ScalarVector;
|
||||
typedef PermutationMatrix<Dynamic, Dynamic, Index> 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 <typename, typename > friend struct SparseQR_QProduct;
|
||||
template <typename > 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 <typename MatrixType, typename OrderingType>
|
||||
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");
|
||||
// Copy to a column major matrix if the input is rowmajor
|
||||
typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::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<MatrixType,OrderingType>::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<MatrixType,OrderingType>::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<MatrixType,OrderingType>::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<RealScalar>::epsilon();
|
||||
if(max2Norm==RealScalar(0))
|
||||
max2Norm = RealScalar(1);
|
||||
pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::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<MatrixType,OrderingType>::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<MatrixType,OrderingType>::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<MatrixType,OrderingType>::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<diagSize)
|
||||
m_Q.startVec(nonzeroCol);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Zero pivot found: move implicitly this column to the end
|
||||
m_hcoeffs(col) = Scalar(0);
|
||||
for (Index j = nonzeroCol; j < n-1; j++)
|
||||
std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]);
|
||||
|
||||
// Recompute the column elimination tree
|
||||
internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data());
|
||||
m_isEtreeOk = false;
|
||||
}
|
||||
}
|
||||
|
||||
m_hcoeffs.tail(diagSize-nonzeroCol).setZero();
|
||||
|
||||
// Finalize the column pointers of the sparse matrices R and Q
|
||||
m_Q.finalize();
|
||||
m_Q.makeCompressed();
|
||||
m_R.finalize();
|
||||
m_R.makeCompressed();
|
||||
m_isQSorted = false;
|
||||
|
||||
|
||||
m_nonzeropivots = nonzeroCol;
|
||||
|
||||
if(nonzeroCol<n)
|
||||
{
|
||||
// Permute the triangular factor to put the 'dead' columns to the end
|
||||
MatrixType tempR(m_R);
|
||||
QRMatrixType tempR(m_R);
|
||||
m_R = tempR * m_pivotperm;
|
||||
|
||||
// Update the column permutation
|
||||
|
|
@ -561,14 +616,16 @@ struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived
|
|||
template<typename DesType>
|
||||
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<SparseQR_QProduct<SparseQRType, Derived
|
|||
else
|
||||
{
|
||||
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 k = n-1; k >=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<SparseQRMatrixQReturnType<Sp
|
|||
return SparseQRMatrixQTransposeReturnType<SparseQRType>(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<SparseQRType> transpose() const
|
||||
{
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_StlSupport_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_StlSupport_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/StlSupport COMPONENT Devel
|
||||
)
|
||||
|
|
@ -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__)
|
||||
|
|
|
|||
|
|
@ -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__)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_SuperLUSupport_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_SuperLUSupport_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/SuperLUSupport COMPONENT Devel
|
||||
)
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_UmfPackSupport_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_UmfPackSupport_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/UmfPackSupport COMPONENT Devel
|
||||
)
|
||||
|
|
@ -107,6 +107,16 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N
|
|||
return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
|
||||
}
|
||||
|
||||
namespace internal {
|
||||
template<typename T> struct umfpack_helper_is_sparse_plain : false_type {};
|
||||
template<typename Scalar, int Options, typename StorageIndex>
|
||||
struct umfpack_helper_is_sparse_plain<SparseMatrix<Scalar,Options,StorageIndex> >
|
||||
: true_type {};
|
||||
template<typename Scalar, int Options, typename StorageIndex>
|
||||
struct umfpack_helper_is_sparse_plain<MappedSparseMatrix<Scalar,Options,StorageIndex> >
|
||||
: 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<typename InputMatrixType>
|
||||
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<typename InputMatrixType>
|
||||
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<typename InputMatrixType>
|
||||
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<typename InputMatrixType>
|
||||
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<typename InputMatrixType>
|
||||
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<typename InputMatrixType>
|
||||
void grapInput(const InputMatrixType& mat)
|
||||
{
|
||||
grapInput_impl(mat, internal::umfpack_helper_is_sparse_plain<InputMatrixType>());
|
||||
}
|
||||
|
||||
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;
|
||||
|
|
|
|||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_misc_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_misc_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/misc COMPONENT Devel
|
||||
)
|
||||
|
|
@ -1,6 +0,0 @@
|
|||
FILE(GLOB Eigen_plugins_SRCS "*.h")
|
||||
|
||||
INSTALL(FILES
|
||||
${Eigen_plugins_SRCS}
|
||||
DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/plugins COMPONENT Devel
|
||||
)
|
||||
|
|
@ -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)
|
||||
|
|
@ -178,11 +178,11 @@ template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,Affine>& t)
|
|||
template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,Projective>& t) { glLoadMatrix(t.matrix()); }
|
||||
template<typename Scalar> void glLoadMatrix(const Transform<Scalar,3,AffineCompact>& t) { glLoadMatrix(Transform<Scalar,3,Affine>(t).matrix()); }
|
||||
|
||||
static void glRotate(const Rotation2D<float>& rot)
|
||||
inline void glRotate(const Rotation2D<float>& rot)
|
||||
{
|
||||
glRotatef(rot.angle()*180.f/float(M_PI), 0.f, 0.f, 1.f);
|
||||
}
|
||||
static void glRotate(const Rotation2D<double>& rot)
|
||||
inline void glRotate(const Rotation2D<double>& 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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -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
|
||||
)
|
||||
|
|
@ -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)
|
||||
Some files were not shown because too many files have changed in this diff Show more
Loading…
Add table
Add a link
Reference in a new issue