diff options
author | Gael Guennebaud <g.gael@free.fr> | 2014-02-18 15:45:39 +0100 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2014-02-18 15:45:39 +0100 |
commit | 463554c25456613e55dbac3e5f734b59f1708dca (patch) | |
tree | 47395219cf1d2301d1b3aac8f0bc1a01857919a8 /Eigen | |
parent | 82c066b3c4a0d64962ad0bf268bf7a4a28b2704b (diff) | |
parent | 37a1d736bfd6fdaa0e1a93aee585d1861ab4004a (diff) |
Merge with default branch
Diffstat (limited to 'Eigen')
28 files changed, 543 insertions, 320 deletions
@@ -15,7 +15,9 @@ * * This module provides various QR decompositions * This module also provides some MatrixBase methods, including: - * - MatrixBase::qr(), + * - MatrixBase::householderQr() + * - MatrixBase::colPivHouseholderQr() + * - MatrixBase::fullPivHouseholderQr() * * \code * #include <Eigen/QR> diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 4faedd257..f34d26465 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -16,7 +16,10 @@ namespace Eigen { namespace internal { -template<typename MatrixType, int UpLo> struct LDLT_Traits; + template<typename MatrixType, int UpLo> struct LDLT_Traits; + + // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef + enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite }; } /** \ingroup Cholesky_Module @@ -69,7 +72,12 @@ template<typename _MatrixType, int _UpLo> class LDLT * The default constructor is useful in cases in which the user intends to * perform decompositions via LDLT::compute(const MatrixType&). */ - LDLT() : m_matrix(), m_transpositions(), m_isInitialized(false) {} + LDLT() + : m_matrix(), + m_transpositions(), + m_sign(internal::ZeroSign), + m_isInitialized(false) + {} /** \brief Default Constructor with memory preallocation * @@ -81,6 +89,7 @@ template<typename _MatrixType, int _UpLo> class LDLT : m_matrix(size, size), m_transpositions(size), m_temporary(size), + m_sign(internal::ZeroSign), m_isInitialized(false) {} @@ -93,6 +102,7 @@ template<typename _MatrixType, int _UpLo> class LDLT : m_matrix(matrix.rows(), matrix.cols()), m_transpositions(matrix.rows()), m_temporary(matrix.rows()), + m_sign(internal::ZeroSign), m_isInitialized(false) { compute(matrix); @@ -139,7 +149,7 @@ template<typename _MatrixType, int _UpLo> class LDLT inline bool isPositive() const { eigen_assert(m_isInitialized && "LDLT is not initialized."); - return m_sign == 1; + return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign; } #ifdef EIGEN2_SUPPORT @@ -153,7 +163,7 @@ template<typename _MatrixType, int _UpLo> class LDLT inline bool isNegative(void) const { eigen_assert(m_isInitialized && "LDLT is not initialized."); - return m_sign == -1; + return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign; } /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A. @@ -235,7 +245,7 @@ template<typename _MatrixType, int _UpLo> class LDLT MatrixType m_matrix; TranspositionType m_transpositions; TmpMatrixType m_temporary; - int m_sign; + internal::SignMatrix m_sign; bool m_isInitialized; }; @@ -246,7 +256,7 @@ template<int UpLo> struct ldlt_inplace; template<> struct ldlt_inplace<Lower> { template<typename MatrixType, typename TranspositionType, typename Workspace> - static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0) + static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) { using std::abs; typedef typename MatrixType::Scalar Scalar; @@ -258,8 +268,9 @@ template<> struct ldlt_inplace<Lower> if (size <= 1) { transpositions.setIdentity(); - if(sign) - *sign = numext::real(mat.coeff(0,0))>0 ? 1:-1; + if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef; + else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef; + else sign = ZeroSign; return true; } @@ -284,7 +295,6 @@ template<> struct ldlt_inplace<Lower> if(biggest_in_corner < cutoff) { for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i; - if(sign) *sign = 0; break; } @@ -325,15 +335,15 @@ template<> struct ldlt_inplace<Lower> } if((rs>0) && (abs(mat.coeffRef(k,k)) > cutoff)) A21 /= mat.coeffRef(k,k); - - if(sign) - { - // LDLT is not guaranteed to work for indefinite matrices, but let's try to get the sign right - int newSign = numext::real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0; - if(k == 0) - *sign = newSign; - else if(*sign != newSign) - *sign = 0; + + RealScalar realAkk = numext::real(mat.coeffRef(k,k)); + if (sign == PositiveSemiDef) { + if (realAkk < 0) sign = Indefinite; + } else if (sign == NegativeSemiDef) { + if (realAkk > 0) sign = Indefinite; + } else if (sign == ZeroSign) { + if (realAkk > 0) sign = PositiveSemiDef; + else if (realAkk < 0) sign = NegativeSemiDef; } } @@ -399,7 +409,7 @@ template<> struct ldlt_inplace<Lower> template<> struct ldlt_inplace<Upper> { template<typename MatrixType, typename TranspositionType, typename Workspace> - static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, int* sign=0) + static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) { Transpose<MatrixType> matt(mat); return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign); @@ -445,7 +455,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a) m_isInitialized = false; m_temporary.resize(size); - internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign); + internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign); m_isInitialized = true; return *this; @@ -473,7 +483,7 @@ LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Deri for (Index i = 0; i < size; i++) m_transpositions.coeffRef(i) = i; m_temporary.resize(size); - m_sign = sigma>=0 ? 1 : -1; + m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef; m_isInitialized = true; } diff --git a/Eigen/src/Core/CwiseNullaryOp.h b/Eigen/src/Core/CwiseNullaryOp.h index 1d4ee50a8..124383114 100644 --- a/Eigen/src/Core/CwiseNullaryOp.h +++ b/Eigen/src/Core/CwiseNullaryOp.h @@ -138,6 +138,9 @@ DenseBase<Derived>::NullaryExpr(Index rows, Index cols, const CustomNullaryOp& f * * The template parameter \a CustomNullaryOp is the type of the functor. * + * Here is an example with C++11 random generators: \include random_cpp11.cpp + * Output: \verbinclude random_cpp11.out + * * \sa class CwiseNullaryOp */ template<typename Derived> diff --git a/Eigen/src/Core/Diagonal.h b/Eigen/src/Core/Diagonal.h index 4436c6a69..b160479ab 100644 --- a/Eigen/src/Core/Diagonal.h +++ b/Eigen/src/Core/Diagonal.h @@ -77,7 +77,12 @@ template<typename MatrixType, int _DiagIndex> class Diagonal EIGEN_DEVICE_FUNC inline Index rows() const - { return m_index.value()<0 ? (std::min<Index>)(m_matrix.cols(),m_matrix.rows()+m_index.value()) : (std::min<Index>)(m_matrix.rows(),m_matrix.cols()-m_index.value()); } + { + EIGEN_USING_STD_MATH(min); + return m_index.value()<0 ? (min)(Index(m_matrix.cols()),Index(m_matrix.rows()+m_index.value())) + : (min)(Index(m_matrix.rows()),Index(m_matrix.cols()-m_index.value())); + + } EIGEN_DEVICE_FUNC inline Index cols() const { return 1; } diff --git a/Eigen/src/Core/EigenBase.h b/Eigen/src/Core/EigenBase.h index a25e823ab..1a577c2dc 100644 --- a/Eigen/src/Core/EigenBase.h +++ b/Eigen/src/Core/EigenBase.h @@ -141,36 +141,6 @@ Derived& DenseBase<Derived>::operator-=(const EigenBase<OtherDerived> &other) return derived(); } -/** replaces \c *this by \c *this * \a other. - * - * \returns a reference to \c *this - */ -template<typename Derived> -template<typename OtherDerived> -inline Derived& -MatrixBase<Derived>::operator*=(const EigenBase<OtherDerived> &other) -{ - other.derived().applyThisOnTheRight(derived()); - return derived(); -} - -/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=(). - */ -template<typename Derived> -template<typename OtherDerived> -inline void MatrixBase<Derived>::applyOnTheRight(const EigenBase<OtherDerived> &other) -{ - other.derived().applyThisOnTheRight(derived()); -} - -/** replaces \c *this by \c *this * \a other. */ -template<typename Derived> -template<typename OtherDerived> -inline void MatrixBase<Derived>::applyOnTheLeft(const EigenBase<OtherDerived> &other) -{ - other.derived().applyThisOnTheLeft(derived()); -} - } // end namespace Eigen #endif // EIGEN_EIGENBASE_H diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 6453145d0..851fa28d5 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -575,6 +575,51 @@ template<typename Derived> class MatrixBase {EIGEN_STATIC_ASSERT(std::ptrdiff_t(sizeof(typename OtherDerived::Scalar))==-1,YOU_CANNOT_MIX_ARRAYS_AND_MATRICES); return *this;} }; + +/*************************************************************************** +* Implementation of matrix base methods +***************************************************************************/ + +/** replaces \c *this by \c *this * \a other. + * + * \returns a reference to \c *this + * + * Example: \include MatrixBase_applyOnTheRight.cpp + * Output: \verbinclude MatrixBase_applyOnTheRight.out + */ +template<typename Derived> +template<typename OtherDerived> +inline Derived& +MatrixBase<Derived>::operator*=(const EigenBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheRight(derived()); + return derived(); +} + +/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=(). + * + * Example: \include MatrixBase_applyOnTheRight.cpp + * Output: \verbinclude MatrixBase_applyOnTheRight.out + */ +template<typename Derived> +template<typename OtherDerived> +inline void MatrixBase<Derived>::applyOnTheRight(const EigenBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheRight(derived()); +} + +/** replaces \c *this by \a other * \c *this. + * + * Example: \include MatrixBase_applyOnTheLeft.cpp + * Output: \verbinclude MatrixBase_applyOnTheLeft.out + */ +template<typename Derived> +template<typename OtherDerived> +inline void MatrixBase<Derived>::applyOnTheLeft(const EigenBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheLeft(derived()); +} + } // end namespace Eigen #endif // EIGEN_MATRIXBASE_H diff --git a/Eigen/src/Core/Random.h b/Eigen/src/Core/Random.h index a0bd3039e..fdd43ed0c 100644 --- a/Eigen/src/Core/Random.h +++ b/Eigen/src/Core/Random.h @@ -34,6 +34,8 @@ struct functor_traits<scalar_random_op<Scalar> > * The parameters \a rows and \a cols are the number of rows and of columns of * the returned matrix. Must be compatible with this MatrixBase type. * + * \not_reentrant + * * This variant is meant to be used for dynamic-size matrix types. For fixed-size types, * it is redundant to pass \a rows and \a cols as arguments, so Random() should be used * instead. @@ -45,8 +47,10 @@ struct functor_traits<scalar_random_op<Scalar> > * This expression has the "evaluate before nesting" flag so that it will be evaluated into * a temporary matrix whenever it is nested in a larger expression. This prevents unexpected * behavior with expressions involving random matrices. + * + * See DenseBase::NullaryExpr(Index, const CustomNullaryOp&) for an example using C++11 random generators. * - * \sa MatrixBase::setRandom(), MatrixBase::Random(Index), MatrixBase::Random() + * \sa DenseBase::setRandom(), DenseBase::Random(Index), DenseBase::Random() */ template<typename Derived> inline const CwiseNullaryOp<internal::scalar_random_op<typename internal::traits<Derived>::Scalar>, Derived> @@ -64,6 +68,7 @@ DenseBase<Derived>::Random(Index rows, Index cols) * Must be compatible with this MatrixBase type. * * \only_for_vectors + * \not_reentrant * * This variant is meant to be used for dynamic-size vector types. For fixed-size types, * it is redundant to pass \a size as argument, so Random() should be used @@ -76,7 +81,7 @@ DenseBase<Derived>::Random(Index rows, Index cols) * a temporary vector whenever it is nested in a larger expression. This prevents unexpected * behavior with expressions involving random matrices. * - * \sa MatrixBase::setRandom(), MatrixBase::Random(Index,Index), MatrixBase::Random() + * \sa DenseBase::setRandom(), DenseBase::Random(Index,Index), DenseBase::Random() */ template<typename Derived> inline const CwiseNullaryOp<internal::scalar_random_op<typename internal::traits<Derived>::Scalar>, Derived> @@ -99,8 +104,10 @@ DenseBase<Derived>::Random(Index size) * This expression has the "evaluate before nesting" flag so that it will be evaluated into * a temporary matrix whenever it is nested in a larger expression. This prevents unexpected * behavior with expressions involving random matrices. + * + * \not_reentrant * - * \sa MatrixBase::setRandom(), MatrixBase::Random(Index,Index), MatrixBase::Random(Index) + * \sa DenseBase::setRandom(), DenseBase::Random(Index,Index), DenseBase::Random(Index) */ template<typename Derived> inline const CwiseNullaryOp<internal::scalar_random_op<typename internal::traits<Derived>::Scalar>, Derived> @@ -114,6 +121,8 @@ DenseBase<Derived>::Random() * Numbers are uniformly spread through their whole definition range for integer types, * and in the [-1:1] range for floating point scalar types. * + * \not_reentrant + * * Example: \include MatrixBase_setRandom.cpp * Output: \verbinclude MatrixBase_setRandom.out * @@ -131,11 +140,12 @@ inline Derived& DenseBase<Derived>::setRandom() * and in the [-1:1] range for floating point scalar types. * * \only_for_vectors + * \not_reentrant * * Example: \include Matrix_setRandom_int.cpp * Output: \verbinclude Matrix_setRandom_int.out * - * \sa MatrixBase::setRandom(), setRandom(Index,Index), class CwiseNullaryOp, MatrixBase::Random() + * \sa DenseBase::setRandom(), setRandom(Index,Index), class CwiseNullaryOp, DenseBase::Random() */ template<typename Derived> EIGEN_STRONG_INLINE Derived& @@ -150,13 +160,15 @@ PlainObjectBase<Derived>::setRandom(Index newSize) * Numbers are uniformly spread through their whole definition range for integer types, * and in the [-1:1] range for floating point scalar types. * + * \not_reentrant + * * \param nbRows the new number of rows * \param nbCols the new number of columns * * Example: \include Matrix_setRandom_int_int.cpp * Output: \verbinclude Matrix_setRandom_int_int.out * - * \sa MatrixBase::setRandom(), setRandom(Index), class CwiseNullaryOp, MatrixBase::Random() + * \sa DenseBase::setRandom(), setRandom(Index), class CwiseNullaryOp, DenseBase::Random() */ template<typename Derived> EIGEN_STRONG_INLINE Derived& diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index c219e2f53..525620bad 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h @@ -17,16 +17,29 @@ namespace internal { template<typename ExpressionType, typename Scalar> inline void stable_norm_kernel(const ExpressionType& bl, Scalar& ssq, Scalar& scale, Scalar& invScale) { - Scalar max = bl.cwiseAbs().maxCoeff(); - if (max>scale) + using std::max; + Scalar maxCoeff = bl.cwiseAbs().maxCoeff(); + + if (maxCoeff>scale) { - ssq = ssq * numext::abs2(scale/max); - scale = max; - invScale = Scalar(1)/scale; + ssq = ssq * numext::abs2(scale/maxCoeff); + Scalar tmp = Scalar(1)/maxCoeff; + if(tmp > NumTraits<Scalar>::highest()) + { + invScale = NumTraits<Scalar>::highest(); + scale = Scalar(1)/invScale; + } + else + { + scale = maxCoeff; + invScale = tmp; + } } - // TODO if the max is much much smaller than the current scale, + + // TODO if the maxCoeff is much much smaller than the current scale, // then we can neglect this sub vector - ssq += (bl*invScale).squaredNorm(); + if(scale>Scalar(0)) // if scale==0, then bl is 0 + ssq += (bl*invScale).squaredNorm(); } template<typename Derived> diff --git a/Eigen/src/Core/arch/SSE/PacketMath.h b/Eigen/src/Core/arch/SSE/PacketMath.h index f85d2e06e..f5a3dab52 100644 --- a/Eigen/src/Core/arch/SSE/PacketMath.h +++ b/Eigen/src/Core/arch/SSE/PacketMath.h @@ -504,13 +504,18 @@ template<> EIGEN_STRONG_INLINE double predux_min<Packet2d>(const Packet2d& a) } template<> EIGEN_STRONG_INLINE int predux_min<Packet4i>(const Packet4i& a) { +#ifdef EIGEN_VECTORIZE_SSE4_1 + Packet4i tmp = _mm_min_epi32(a, _mm_shuffle_epi32(a, _MM_SHUFFLE(0,0,3,2))); + return pfirst(_mm_min_epi32(tmp,_mm_shuffle_epi32(tmp, 1))); +#else // after some experiments, it is seems this is the fastest way to implement it // for GCC (eg., it does not like using std::min after the pstore !!) EIGEN_ALIGN16 int aux[4]; pstore(aux, a); - register int aux0 = aux[0]<aux[1] ? aux[0] : aux[1]; - register int aux2 = aux[2]<aux[3] ? aux[2] : aux[3]; + int aux0 = aux[0]<aux[1] ? aux[0] : aux[1]; + int aux2 = aux[2]<aux[3] ? aux[2] : aux[3]; return aux0<aux2 ? aux0 : aux2; +#endif // EIGEN_VECTORIZE_SSE4_1 } // max @@ -525,13 +530,18 @@ template<> EIGEN_STRONG_INLINE double predux_max<Packet2d>(const Packet2d& a) } template<> EIGEN_STRONG_INLINE int predux_max<Packet4i>(const Packet4i& a) { +#ifdef EIGEN_VECTORIZE_SSE4_1 + Packet4i tmp = _mm_max_epi32(a, _mm_shuffle_epi32(a, _MM_SHUFFLE(0,0,3,2))); + return pfirst(_mm_max_epi32(tmp,_mm_shuffle_epi32(tmp, 1))); +#else // after some experiments, it is seems this is the fastest way to implement it // for GCC (eg., it does not like using std::min after the pstore !!) EIGEN_ALIGN16 int aux[4]; pstore(aux, a); - register int aux0 = aux[0]>aux[1] ? aux[0] : aux[1]; - register int aux2 = aux[2]>aux[3] ? aux[2] : aux[3]; + int aux0 = aux[0]>aux[1] ? aux[0] : aux[1]; + int aux2 = aux[2]>aux[3] ? aux[2] : aux[3]; return aux0>aux2 ? aux0 : aux2; +#endif // EIGEN_VECTORIZE_SSE4_1 } #if (defined __GNUC__) diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h index c40e80f53..f698f67f9 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -79,8 +79,8 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_product<Scalar,Index,StorageOrd for (Index j=FirstTriangular ? bound : 0; j<(FirstTriangular ? size : bound);j+=2) { - register const Scalar* EIGEN_RESTRICT A0 = lhs + j*lhsStride; - register const Scalar* EIGEN_RESTRICT A1 = lhs + (j+1)*lhsStride; + const Scalar* EIGEN_RESTRICT A0 = lhs + j*lhsStride; + const Scalar* EIGEN_RESTRICT A1 = lhs + (j+1)*lhsStride; Scalar t0 = cjAlpha * rhs[j]; Packet ptmp0 = pset1<Packet>(t0); @@ -147,7 +147,7 @@ EIGEN_DONT_INLINE void selfadjoint_matrix_vector_product<Scalar,Index,StorageOrd } for (Index j=FirstTriangular ? 0 : bound;j<(FirstTriangular ? bound : size);j++) { - register const Scalar* EIGEN_RESTRICT A0 = lhs + j*lhsStride; + const Scalar* EIGEN_RESTRICT A0 = lhs + j*lhsStride; Scalar t1 = cjAlpha * rhs[j]; Scalar t2(0); diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index a292495e5..bc472d720 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -608,7 +608,7 @@ template<typename T> class aligned_stack_memory_handler */ #ifdef EIGEN_ALLOCA - #ifdef __arm__ + #if defined(__arm__) || defined(_WIN32) #define EIGEN_ALIGNED_ALLOCA(SIZE) reinterpret_cast<void*>((reinterpret_cast<size_t>(EIGEN_ALLOCA(SIZE+16)) & ~(size_t(15))) + 16) #else #define EIGEN_ALIGNED_ALLOCA EIGEN_ALLOCA @@ -664,7 +664,9 @@ template<typename T> class aligned_stack_memory_handler /* memory allocated we can safely let the default implementation handle */ \ /* this particular case. */ \ static void *operator new(size_t size, void *ptr) { return ::operator new(size,ptr); } \ + static void *operator new[](size_t size, void* ptr) { return ::operator new[](size,ptr); } \ void operator delete(void * memory, void *ptr) throw() { return ::operator delete(memory,ptr); } \ + void operator delete[](void * memory, void *ptr) throw() { return ::operator delete[](memory,ptr); } \ /* nothrow-new (returns zero instead of std::bad_alloc) */ \ EIGEN_MAKE_ALIGNED_OPERATOR_NEW_NOTHROW(NeedsToAlign) \ void operator delete(void *ptr, const std::nothrow_t&) throw() { \ @@ -759,11 +761,27 @@ public: ::new( p ) T( value ); } +#if (__cplusplus >= 201103L) + template <typename U, typename... Args> + void construct( U* u, Args&&... args) + { + ::new( static_cast<void*>(u) ) U( std::forward<Args>( args )... ); + } +#endif + void destroy( pointer p ) { p->~T(); } +#if (__cplusplus >= 201103L) + template <typename U> + void destroy( U* u ) + { + u->~U(); + } +#endif + void deallocate( pointer p, size_type /*num*/ ) { internal::aligned_free( p ); diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 0c3425069..b8146d04d 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -756,7 +756,9 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false> EIGEN_DEVICE_FUNC static inline void run(SolverType& solver, const MatrixType& mat, int options) { - using std::sqrt; + EIGEN_USING_STD_MATH(max) + EIGEN_USING_STD_MATH(sqrt); + eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); eigen_assert((options&~(EigVecMask|GenEigMask))==0 && (options&EigVecMask)!=EigVecMask @@ -768,7 +770,7 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false> // map the matrix coefficients to [-1:1] to avoid over- and underflow. Scalar scale = mat.cwiseAbs().maxCoeff(); - scale = (std::max)(scale,Scalar(1)); + scale = (max)(scale,Scalar(1)); MatrixType scaledMat = mat / scale; // Compute the eigenvalues diff --git a/Eigen/src/Geometry/EulerAngles.h b/Eigen/src/Geometry/EulerAngles.h index 97984d590..82802fb43 100644 --- a/Eigen/src/Geometry/EulerAngles.h +++ b/Eigen/src/Geometry/EulerAngles.h @@ -28,7 +28,7 @@ namespace Eigen { * * AngleAxisf(ea[2], Vector3f::UnitZ()); \endcode * This corresponds to the right-multiply conventions (with right hand side frames). * - * The returned angles are in the ranges [0:pi]x[0:pi]x[-pi:pi]. + * The returned angles are in the ranges [0:pi]x[-pi:pi]x[-pi:pi]. * * \sa class AngleAxis */ diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index d036c018a..803001272 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -150,10 +150,6 @@ public: /** \returns the conjugated quaternion */ Quaternion<Scalar> conjugate() const; - /** \returns an interpolation for a constant motion between \a other and \c *this - * \a t in [0;1] - * see http://en.wikipedia.org/wiki/Slerp - */ template<class OtherDerived> Quaternion<Scalar> slerp(const Scalar& t, const QuaternionBase<OtherDerived>& other) const; /** \returns \c true if \c *this is approximately equal to \a other, within the precision @@ -194,11 +190,11 @@ public: * \brief The quaternion class used to represent 3D orientations and rotations * * \tparam _Scalar the scalar type, i.e., the type of the coefficients - * \tparam _Options controls the memory alignement of the coeffecients. Can be \# AutoAlign or \# DontAlign. Default is AutoAlign. + * \tparam _Options controls the memory alignment of the coefficients. Can be \# AutoAlign or \# DontAlign. Default is AutoAlign. * * This class represents a quaternion \f$ w+xi+yj+zk \f$ that is a convenient representation of * orientations and rotations of objects in three dimensions. Compared to other representations - * like Euler angles or 3x3 matrices, quatertions offer the following advantages: + * like Euler angles or 3x3 matrices, quaternions offer the following advantages: * \li \b compact storage (4 scalars) * \li \b efficient to compose (28 flops), * \li \b stable spherical interpolation @@ -385,7 +381,7 @@ class Map<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. */ @@ -399,16 +395,16 @@ class Map<Quaternion<_Scalar>, _Options > }; /** \ingroup Geometry_Module - * Map an unaligned array of single precision scalar as a quaternion */ + * Map an unaligned array of single precision scalars as a quaternion */ typedef Map<Quaternion<float>, 0> QuaternionMapf; /** \ingroup Geometry_Module - * Map an unaligned array of double precision scalar as a quaternion */ + * Map an unaligned array of double precision scalars as a quaternion */ typedef Map<Quaternion<double>, 0> QuaternionMapd; /** \ingroup Geometry_Module - * Map a 16-bits aligned array of double precision scalars as a quaternion */ + * Map a 16-byte aligned array of single precision scalars as a quaternion */ typedef Map<Quaternion<float>, Aligned> QuaternionMapAlignedf; /** \ingroup Geometry_Module - * Map a 16-bits aligned array of double precision scalars as a quaternion */ + * Map a 16-byte aligned array of double precision scalars as a quaternion */ typedef Map<Quaternion<double>, Aligned> QuaternionMapAlignedd; /*************************************************************************** @@ -579,7 +575,7 @@ inline Derived& QuaternionBase<Derived>::setFromTwoVectors(const MatrixBase<Deri Scalar c = v1.dot(v0); // if dot == -1, vectors are nearly opposites - // => accuraletly compute the rotation axis by computing the + // => accurately compute the rotation axis by computing the // intersection of the two planes. This is done by solving: // x^T v0 = 0 // x^T v1 = 0 @@ -677,8 +673,13 @@ QuaternionBase<Derived>::angularDistance(const QuaternionBase<OtherDerived>& oth return static_cast<Scalar>(2 * acos(d)); } + + /** \returns the spherical linear interpolation between the two quaternions - * \c *this and \a other at the parameter \a t + * \c *this and \a other at the parameter \a t in [0;1]. + * + * This represents an interpolation for a constant motion between \c *this and \a other, + * see also http://en.wikipedia.org/wiki/Slerp. */ template <class Derived> template <class OtherDerived> diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 419f90148..0a5012cfe 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -530,9 +530,9 @@ public: inline Transform& operator=(const UniformScaling<Scalar>& t); inline Transform& operator*=(const UniformScaling<Scalar>& s) { return scale(s.factor()); } - inline Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Isometry)> operator*(const UniformScaling<Scalar>& s) const + inline Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode)> operator*(const UniformScaling<Scalar>& s) const { - Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Isometry),Options> res = *this; + Transform<Scalar,Dim,(int(Mode)==int(Isometry)?Affine:Mode),Options> res = *this; res.scale(s.factor()); return res; } @@ -699,9 +699,13 @@ template<typename Scalar, int Dim, int Mode,int Options> Transform<Scalar,Dim,Mode,Options>& Transform<Scalar,Dim,Mode,Options>::operator=(const QMatrix& other) { EIGEN_STATIC_ASSERT(Dim==2, YOU_MADE_A_PROGRAMMING_MISTAKE) - m_matrix << other.m11(), other.m21(), other.dx(), - other.m12(), other.m22(), other.dy(), - 0, 0, 1; + if (Mode == int(AffineCompact)) + m_matrix << other.m11(), other.m21(), other.dx(), + other.m12(), other.m22(), other.dy(); + else + m_matrix << other.m11(), other.m21(), other.dx(), + other.m12(), other.m22(), other.dy(), + 0, 0, 1; return *this; } diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index e6ab172b3..44eaa1b1a 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -33,13 +33,13 @@ struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> > * * \param MatrixType the type of the matrix of which we are computing the QR decomposition * - * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R + * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b P', \b Q and \b R * such that * \f[ - * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R} + * \mathbf{P} \, \mathbf{A} \, \mathbf{P}' = \mathbf{Q} \, \mathbf{R} * \f] - * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an - * upper triangular matrix. + * by using Householder transformations. Here, \b P and \b P' are permutation matrices, \b Q a unitary matrix + * and \b R an upper triangular matrix. * * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR. diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index ad156396a..352dbf3f0 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -251,56 +251,62 @@ void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename } /** \internal */ -template<typename MatrixQR, typename HCoeffs> -void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs, - typename MatrixQR::Index maxBlockSize=32, - typename MatrixQR::Scalar* tempData = 0) +template<typename MatrixQR, typename HCoeffs, + typename MatrixQRScalar = typename MatrixQR::Scalar, + bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)> +struct householder_qr_inplace_blocked { - typedef typename MatrixQR::Index Index; - typedef typename MatrixQR::Scalar Scalar; - typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; - - Index rows = mat.rows(); - Index cols = mat.cols(); - Index size = (std::min)(rows, cols); - - typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; - TempType tempVector; - if(tempData==0) + // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h + static void run(MatrixQR& mat, HCoeffs& hCoeffs, + typename MatrixQR::Index maxBlockSize=32, + typename MatrixQR::Scalar* tempData = 0) { - tempVector.resize(cols); - tempData = tempVector.data(); - } - - Index blockSize = (std::min)(maxBlockSize,size); - - Index k = 0; - for (k = 0; k < size; k += blockSize) - { - Index bs = (std::min)(size-k,blockSize); // actual size of the block - Index tcols = cols - k - bs; // trailing columns - Index brows = rows-k; // rows of the block + typedef typename MatrixQR::Index Index; + typedef typename MatrixQR::Scalar Scalar; + typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; - // partition the matrix: - // A00 | A01 | A02 - // mat = A10 | A11 | A12 - // A20 | A21 | A22 - // and performs the qr dec of [A11^T A12^T]^T - // and update [A21^T A22^T]^T using level 3 operations. - // Finally, the algorithm continue on A22 + Index rows = mat.rows(); + Index cols = mat.cols(); + Index size = (std::min)(rows, cols); - BlockType A11_21 = mat.block(k,k,brows,bs); - Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); + typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; + TempType tempVector; + if(tempData==0) + { + tempVector.resize(cols); + tempData = tempVector.data(); + } - householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); + Index blockSize = (std::min)(maxBlockSize,size); - if(tcols) + Index k = 0; + for (k = 0; k < size; k += blockSize) { - BlockType A21_22 = mat.block(k,k+bs,brows,tcols); - apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint()); + Index bs = (std::min)(size-k,blockSize); // actual size of the block + Index tcols = cols - k - bs; // trailing columns + Index brows = rows-k; // rows of the block + + // partition the matrix: + // A00 | A01 | A02 + // mat = A10 | A11 | A12 + // A20 | A21 | A22 + // and performs the qr dec of [A11^T A12^T]^T + // and update [A21^T A22^T]^T using level 3 operations. + // Finally, the algorithm continue on A22 + + BlockType A11_21 = mat.block(k,k,brows,bs); + Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); + + householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); + + if(tcols) + { + BlockType A21_22 = mat.block(k,k+bs,brows,tcols); + apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint()); + } } } -} +}; template<typename _MatrixType, typename Rhs> struct solve_retval<HouseholderQR<_MatrixType>, Rhs> @@ -352,7 +358,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& m_temp.resize(cols); - internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data()); + internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data()); m_isInitialized = true; return *this; diff --git a/Eigen/src/QR/HouseholderQR_MKL.h b/Eigen/src/QR/HouseholderQR_MKL.h index 5313de604..8a3a7e406 100644 --- a/Eigen/src/QR/HouseholderQR_MKL.h +++ b/Eigen/src/QR/HouseholderQR_MKL.h @@ -34,7 +34,7 @@ #ifndef EIGEN_QR_MKL_H #define EIGEN_QR_MKL_H -#include "Eigen/src/Core/util/MKL_support.h" +#include "../Core/util/MKL_support.h" namespace Eigen { @@ -44,18 +44,20 @@ namespace internal { #define EIGEN_MKL_QR_NOPIV(EIGTYPE, MKLTYPE, MKLPREFIX) \ template<typename MatrixQR, typename HCoeffs> \ -void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs, \ - typename MatrixQR::Index maxBlockSize=32, \ - EIGTYPE* tempData = 0) \ +struct householder_qr_inplace_blocked<MatrixQR, HCoeffs, EIGTYPE, true> \ { \ - lapack_int m = mat.rows(); \ - lapack_int n = mat.cols(); \ - lapack_int lda = mat.outerStride(); \ - lapack_int matrix_order = (MatrixQR::IsRowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \ - LAPACKE_##MKLPREFIX##geqrf( matrix_order, m, n, (MKLTYPE*)mat.data(), lda, (MKLTYPE*)hCoeffs.data()); \ - hCoeffs.adjointInPlace(); \ -\ -} + static void run(MatrixQR& mat, HCoeffs& hCoeffs, \ + typename MatrixQR::Index = 32, \ + typename MatrixQR::Scalar* = 0) \ + { \ + lapack_int m = (lapack_int) mat.rows(); \ + lapack_int n = (lapack_int) mat.cols(); \ + lapack_int lda = (lapack_int) mat.outerStride(); \ + lapack_int matrix_order = (MatrixQR::IsRowMajor) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; \ + LAPACKE_##MKLPREFIX##geqrf( matrix_order, m, n, (MKLTYPE*)mat.data(), lda, (MKLTYPE*)hCoeffs.data()); \ + hCoeffs.adjointInPlace(); \ + } \ +}; EIGEN_MKL_QR_NOPIV(double, double, d) EIGEN_MKL_QR_NOPIV(float, float, s) diff --git a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h index 4b13f08d4..5c320e2d2 100644 --- a/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h +++ b/Eigen/src/SparseCore/ConservativeSparseSparseProduct.h @@ -66,9 +66,9 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r } // unordered insertion - for(int k=0; k<nnz; ++k) + for(Index k=0; k<nnz; ++k) { - int i = indices[k]; + Index i = indices[k]; res.insertBackByOuterInnerUnordered(j,i) = values[i]; mask[i] = false; } @@ -76,8 +76,8 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r #if 0 // alternative ordered insertion code: - int t200 = rows/(log2(200)*1.39); - int t = (rows*100)/139; + Index t200 = rows/(log2(200)*1.39); + Index t = (rows*100)/139; // FIXME reserve nnz non zeros // FIXME implement fast sort algorithms for very small nnz @@ -90,9 +90,9 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r if(true) { if(nnz>1) std::sort(indices.data(),indices.data()+nnz); - for(int k=0; k<nnz; ++k) + for(Index k=0; k<nnz; ++k) { - int i = indices[k]; + Index i = indices[k]; res.insertBackByOuterInner(j,i) = values[i]; mask[i] = false; } @@ -100,7 +100,7 @@ static void conservative_sparse_sparse_product_impl(const Lhs& lhs, const Rhs& r else { // dense path - for(int i=0; i<rows; ++i) + for(Index i=0; i<rows; ++i) { if(mask[i]) { @@ -134,8 +134,8 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,C static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) { - typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix; - typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix; ColMajorMatrix resCol(lhs.rows(),rhs.cols()); internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol); // sort the non zeros: @@ -149,7 +149,7 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,C { static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) { - typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix; RowMajorMatrix rhsRow = rhs; RowMajorMatrix resRow(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl<RowMajorMatrix,Lhs,RowMajorMatrix>(rhsRow, lhs, resRow); @@ -162,7 +162,7 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,R { static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) { - typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix; RowMajorMatrix lhsRow = lhs; RowMajorMatrix resRow(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl<Rhs,RowMajorMatrix,RowMajorMatrix>(rhs, lhsRow, resRow); @@ -175,7 +175,7 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,R { static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) { - typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix; RowMajorMatrix resRow(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow); res = resRow; @@ -190,7 +190,7 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,C static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) { - typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix; ColMajorMatrix resCol(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl<Lhs,Rhs,ColMajorMatrix>(lhs, rhs, resCol); res = resCol; @@ -202,7 +202,7 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,C { static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) { - typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix; ColMajorMatrix lhsCol = lhs; ColMajorMatrix resCol(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl<ColMajorMatrix,Rhs,ColMajorMatrix>(lhsCol, rhs, resCol); @@ -215,7 +215,7 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,ColMajor,R { static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) { - typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix; ColMajorMatrix rhsCol = rhs; ColMajorMatrix resCol(lhs.rows(), rhs.cols()); internal::conservative_sparse_sparse_product_impl<Lhs,ColMajorMatrix,ColMajorMatrix>(lhs, rhsCol, resCol); @@ -228,8 +228,8 @@ struct conservative_sparse_sparse_product_selector<Lhs,Rhs,ResultType,RowMajor,R { static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res) { - typedef SparseMatrix<typename ResultType::Scalar,RowMajor> RowMajorMatrix; - typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,RowMajor,typename ResultType::Index> RowMajorMatrix; + typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> ColMajorMatrix; RowMajorMatrix resRow(lhs.rows(),rhs.cols()); internal::conservative_sparse_sparse_product_impl<Rhs,Lhs,RowMajorMatrix>(rhs, lhs, resRow); // sort the non zeros: diff --git a/Eigen/src/SparseCore/SparseBlock.h b/Eigen/src/SparseCore/SparseBlock.h index d0436aa7c..6f615e250 100644 --- a/Eigen/src/SparseCore/SparseBlock.h +++ b/Eigen/src/SparseCore/SparseBlock.h @@ -335,6 +335,14 @@ const Block<const Derived,Dynamic,Dynamic,true> SparseMatrixBase<Derived>::inner }
+namespace internal {
+
+template< typename XprType, int BlockRows, int BlockCols, bool InnerPanel,
+ bool OuterVector = (BlockCols==1 && XprType::IsRowMajor) || (BlockRows==1 && !XprType::IsRowMajor)>
+class GenericSparseBlockInnerIteratorImpl;
+
+}
+
/** Generic implementation of sparse Block expression.
* Real-only.
*/
@@ -342,11 +350,12 @@ template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel> class BlockImpl<XprType,BlockRows,BlockCols,InnerPanel,Sparse>
: public SparseMatrixBase<Block<XprType,BlockRows,BlockCols,InnerPanel> >, internal::no_assignment_operator
{
- typedef typename internal::remove_all<typename XprType::Nested>::type _MatrixTypeNested;
typedef Block<XprType, BlockRows, BlockCols, InnerPanel> BlockType;
public:
enum { IsRowMajor = internal::traits<BlockType>::IsRowMajor };
EIGEN_SPARSE_PUBLIC_INTERFACE(BlockType)
+
+ typedef typename internal::remove_all<typename XprType::Nested>::type _MatrixTypeNested;
/** Column or Row constructor
*/
@@ -354,8 +363,8 @@ public: : m_matrix(xpr),
m_startRow( (BlockRows==1) && (BlockCols==XprType::ColsAtCompileTime) ? i : 0),
m_startCol( (BlockRows==XprType::RowsAtCompileTime) && (BlockCols==1) ? i : 0),
- m_blockRows(xpr.rows()),
- m_blockCols(xpr.cols())
+ m_blockRows(BlockRows==1 ? 1 : xpr.rows()),
+ m_blockCols(BlockCols==1 ? 1 : xpr.cols())
{}
/** Dynamic-size constructor
@@ -394,29 +403,8 @@ public: inline const _MatrixTypeNested& nestedExpression() const { return m_matrix; }
- class InnerIterator : public _MatrixTypeNested::InnerIterator
- {
- typedef typename _MatrixTypeNested::InnerIterator Base;
- const BlockType& m_block;
- Index m_end;
- public:
-
- EIGEN_STRONG_INLINE InnerIterator(const BlockType& block, Index outer)
- : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())),
- m_block(block),
- m_end(IsRowMajor ? block.m_startCol.value()+block.m_blockCols.value() : block.m_startRow.value()+block.m_blockRows.value())
- {
- while( (Base::operator bool()) && (Base::index() < (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value())) )
- Base::operator++();
- }
-
- inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); }
- inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); }
- inline Index row() const { return Base::row() - m_block.m_startRow.value(); }
- inline Index col() const { return Base::col() - m_block.m_startCol.value(); }
-
- inline operator bool() const { return Base::operator bool() && Base::index() < m_end; }
- };
+ typedef internal::GenericSparseBlockInnerIteratorImpl<XprType,BlockRows,BlockCols,InnerPanel> InnerIterator;
+
class ReverseInnerIterator : public _MatrixTypeNested::ReverseInnerIterator
{
typedef typename _MatrixTypeNested::ReverseInnerIterator Base;
@@ -441,7 +429,7 @@ public: inline operator bool() const { return Base::operator bool() && Base::index() >= m_begin; }
};
protected:
- friend class InnerIterator;
+ friend class internal::GenericSparseBlockInnerIteratorImpl<XprType,BlockRows,BlockCols,InnerPanel>;
friend class ReverseInnerIterator;
EIGEN_INHERIT_ASSIGNMENT_OPERATORS(BlockImpl)
@@ -454,6 +442,100 @@ public: };
+namespace internal {
+ template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
+ class GenericSparseBlockInnerIteratorImpl<XprType,BlockRows,BlockCols,InnerPanel,false> : public Block<XprType, BlockRows, BlockCols, InnerPanel>::_MatrixTypeNested::InnerIterator
+ {
+ typedef Block<XprType, BlockRows, BlockCols, InnerPanel> BlockType;
+ enum {
+ IsRowMajor = BlockType::IsRowMajor
+ };
+ typedef typename BlockType::_MatrixTypeNested _MatrixTypeNested;
+ typedef typename BlockType::Index Index;
+ typedef typename _MatrixTypeNested::InnerIterator Base;
+ const BlockType& m_block;
+ Index m_end;
+ public:
+
+ EIGEN_STRONG_INLINE GenericSparseBlockInnerIteratorImpl(const BlockType& block, Index outer)
+ : Base(block.derived().nestedExpression(), outer + (IsRowMajor ? block.m_startRow.value() : block.m_startCol.value())),
+ m_block(block),
+ m_end(IsRowMajor ? block.m_startCol.value()+block.m_blockCols.value() : block.m_startRow.value()+block.m_blockRows.value())
+ {
+ while( (Base::operator bool()) && (Base::index() < (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value())) )
+ Base::operator++();
+ }
+
+ inline Index index() const { return Base::index() - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); }
+ inline Index outer() const { return Base::outer() - (IsRowMajor ? m_block.m_startRow.value() : m_block.m_startCol.value()); }
+ inline Index row() const { return Base::row() - m_block.m_startRow.value(); }
+ inline Index col() const { return Base::col() - m_block.m_startCol.value(); }
+
+ inline operator bool() const { return Base::operator bool() && Base::index() < m_end; }
+ };
+
+ // Row vector of a column-major sparse matrix or column of a row-major one.
+ template<typename XprType, int BlockRows, int BlockCols, bool InnerPanel>
+ class GenericSparseBlockInnerIteratorImpl<XprType,BlockRows,BlockCols,InnerPanel,true>
+ {
+ typedef Block<XprType, BlockRows, BlockCols, InnerPanel> BlockType;
+ enum {
+ IsRowMajor = BlockType::IsRowMajor
+ };
+ typedef typename BlockType::_MatrixTypeNested _MatrixTypeNested;
+ typedef typename BlockType::Index Index;
+ typedef typename BlockType::Scalar Scalar;
+ const BlockType& m_block;
+ Index m_outerPos;
+ Index m_innerIndex;
+ Scalar m_value;
+ Index m_end;
+ public:
+
+ EIGEN_STRONG_INLINE GenericSparseBlockInnerIteratorImpl(const BlockType& block, Index outer = 0)
+ :
+ m_block(block),
+ m_outerPos( (IsRowMajor ? block.m_startCol.value() : block.m_startRow.value()) - 1), // -1 so that operator++ finds the first non-zero entry
+ m_innerIndex(IsRowMajor ? block.m_startRow.value() : block.m_startCol.value()),
+ m_end(IsRowMajor ? block.m_startCol.value()+block.m_blockCols.value() : block.m_startRow.value()+block.m_blockRows.value())
+ {
+ EIGEN_UNUSED_VARIABLE(outer);
+ eigen_assert(outer==0);
+
+ ++(*this);
+ }
+
+ inline Index index() const { return m_outerPos - (IsRowMajor ? m_block.m_startCol.value() : m_block.m_startRow.value()); }
+ inline Index outer() const { return 0; }
+ inline Index row() const { return IsRowMajor ? 0 : index(); }
+ inline Index col() const { return IsRowMajor ? index() : 0; }
+
+ inline Scalar value() const { return m_value; }
+
+ inline GenericSparseBlockInnerIteratorImpl& operator++()
+ {
+ // search next non-zero entry
+ while(m_outerPos<m_end)
+ {
+ m_outerPos++;
+ typename XprType::InnerIterator it(m_block.m_matrix, m_outerPos);
+ // search for the key m_innerIndex in the current outer-vector
+ while(it && it.index() < m_innerIndex) ++it;
+ if(it && it.index()==m_innerIndex)
+ {
+ m_value = it.value();
+ break;
+ }
+ }
+ return *this;
+ }
+
+ inline operator bool() const { return m_outerPos < m_end; }
+ };
+
+} // end namespace internal
+
+
} // end namespace Eigen
#endif // EIGEN_SPARSE_BLOCK_H
diff --git a/Eigen/src/SparseCore/SparseDenseProduct.h b/Eigen/src/SparseCore/SparseDenseProduct.h index 6132a4880..610833f3b 100644 --- a/Eigen/src/SparseCore/SparseDenseProduct.h +++ b/Eigen/src/SparseCore/SparseDenseProduct.h @@ -125,7 +125,7 @@ class SparseDenseOuterProduct<Lhs,Rhs,Transpose>::InnerIterator : public _LhsNes inline Scalar value() const { return Base::value() * m_factor; } protected: - int m_outer; + Index m_outer; Scalar m_factor; }; @@ -156,7 +156,7 @@ struct sparse_time_dense_product_impl<SparseLhsType,DenseRhsType,DenseResType, t { for(Index c=0; c<rhs.cols(); ++c) { - int n = lhs.outerSize(); + Index n = lhs.outerSize(); for(Index j=0; j<n; ++j) { typename Res::Scalar tmp(0); diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h index f2e234eee..0e7e760fb 100644 --- a/Eigen/src/SparseCore/SparseMatrix.h +++ b/Eigen/src/SparseCore/SparseMatrix.h @@ -402,7 +402,7 @@ class SparseMatrix * \sa insertBack, insertBackByOuterInner */ inline void startVec(Index outer) { - eigen_assert(m_outerIndex[outer]==int(m_data.size()) && "You must call startVec for each inner vector sequentially"); + eigen_assert(m_outerIndex[outer]==Index(m_data.size()) && "You must call startVec for each inner vector sequentially"); eigen_assert(m_outerIndex[outer+1]==0 && "You must call startVec for each inner vector sequentially"); m_outerIndex[outer+1] = m_outerIndex[outer]; } @@ -480,7 +480,7 @@ class SparseMatrix if(m_innerNonZeros != 0) return; m_innerNonZeros = static_cast<Index*>(std::malloc(m_outerSize * sizeof(Index))); - for (int i = 0; i < m_outerSize; i++) + for (Index i = 0; i < m_outerSize; i++) { m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; } @@ -752,8 +752,8 @@ class SparseMatrix else for (Index i=0; i<m.outerSize(); ++i) { - int p = m.m_outerIndex[i]; - int pe = m.m_outerIndex[i]+m.m_innerNonZeros[i]; + Index p = m.m_outerIndex[i]; + Index pe = m.m_outerIndex[i]+m.m_innerNonZeros[i]; Index k=p; for (; k<pe; ++k) s << "(" << m.m_data.value(k) << "," << m.m_data.index(k) << ") "; @@ -1022,7 +1022,7 @@ void SparseMatrix<Scalar,_Options,_Index>::sumupDuplicates() wi.fill(-1); Index count = 0; // for each inner-vector, wi[inner_index] will hold the position of first element into the index/value buffers - for(int j=0; j<outerSize(); ++j) + for(Index j=0; j<outerSize(); ++j) { Index start = count; Index oldEnd = m_outerIndex[j]+m_innerNonZeros[j]; diff --git a/Eigen/src/SparseCore/SparseMatrixBase.h b/Eigen/src/SparseCore/SparseMatrixBase.h index 706f699b8..bbcf7fb1c 100644 --- a/Eigen/src/SparseCore/SparseMatrixBase.h +++ b/Eigen/src/SparseCore/SparseMatrixBase.h @@ -302,8 +302,8 @@ template<typename Derived> class SparseMatrixBase : public EigenBase<Derived> } else { - SparseMatrix<Scalar, RowMajorBit> trans = m; - s << static_cast<const SparseMatrixBase<SparseMatrix<Scalar, RowMajorBit> >&>(trans); + SparseMatrix<Scalar, RowMajorBit, Index> trans = m; + s << static_cast<const SparseMatrixBase<SparseMatrix<Scalar, RowMajorBit, Index> >&>(trans); } } return s; diff --git a/Eigen/src/SparseCore/SparseProduct.h b/Eigen/src/SparseCore/SparseProduct.h index 70b6480ef..cf7663070 100644 --- a/Eigen/src/SparseCore/SparseProduct.h +++ b/Eigen/src/SparseCore/SparseProduct.h @@ -16,6 +16,7 @@ template<typename Lhs, typename Rhs> struct SparseSparseProductReturnType { typedef typename internal::traits<Lhs>::Scalar Scalar; + typedef typename internal::traits<Lhs>::Index Index; enum { LhsRowMajor = internal::traits<Lhs>::Flags & RowMajorBit, RhsRowMajor = internal::traits<Rhs>::Flags & RowMajorBit, @@ -24,11 +25,11 @@ struct SparseSparseProductReturnType }; typedef typename internal::conditional<TransposeLhs, - SparseMatrix<Scalar,0>, + SparseMatrix<Scalar,0,Index>, typename internal::nested<Lhs,Rhs::RowsAtCompileTime>::type>::type LhsNested; typedef typename internal::conditional<TransposeRhs, - SparseMatrix<Scalar,0>, + SparseMatrix<Scalar,0,Index>, typename internal::nested<Rhs,Lhs::RowsAtCompileTime>::type>::type RhsNested; typedef SparseSparseProduct<LhsNested, RhsNested> Type; diff --git a/Eigen/src/SparseCore/SparseSparseProductWithPruning.h b/Eigen/src/SparseCore/SparseSparseProductWithPruning.h index 70857c7b6..fcc18f5c9 100644 --- a/Eigen/src/SparseCore/SparseSparseProductWithPruning.h +++ b/Eigen/src/SparseCore/SparseSparseProductWithPruning.h @@ -27,7 +27,7 @@ static void sparse_sparse_product_with_pruning_impl(const Lhs& lhs, const Rhs& r // make sure to call innerSize/outerSize since we fake the storage order. Index rows = lhs.innerSize(); Index cols = rhs.outerSize(); - //int size = lhs.outerSize(); + //Index size = lhs.outerSize(); eigen_assert(lhs.outerSize() == rhs.innerSize()); // allocate a temporary buffer @@ -100,7 +100,7 @@ struct sparse_sparse_product_with_pruning_selector<Lhs,Rhs,ResultType,ColMajor,C static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, const RealScalar& tolerance) { // we need a col-major matrix to hold the result - typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType; + typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename ResultType::Index> SparseTemporaryType; SparseTemporaryType _res(res.rows(), res.cols()); internal::sparse_sparse_product_with_pruning_impl<Lhs,Rhs,SparseTemporaryType>(lhs, rhs, _res, tolerance); res = _res; @@ -126,10 +126,11 @@ struct sparse_sparse_product_with_pruning_selector<Lhs,Rhs,ResultType,RowMajor,R typedef typename ResultType::RealScalar RealScalar; static void run(const Lhs& lhs, const Rhs& rhs, ResultType& res, const RealScalar& tolerance) { - typedef SparseMatrix<typename ResultType::Scalar,ColMajor> ColMajorMatrix; - ColMajorMatrix colLhs(lhs); - ColMajorMatrix colRhs(rhs); - internal::sparse_sparse_product_with_pruning_impl<ColMajorMatrix,ColMajorMatrix,ResultType>(colLhs, colRhs, res, tolerance); + typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename Lhs::Index> ColMajorMatrixLhs; + typedef SparseMatrix<typename ResultType::Scalar,ColMajor,typename Lhs::Index> ColMajorMatrixRhs; + ColMajorMatrixLhs colLhs(lhs); + ColMajorMatrixRhs colRhs(rhs); + internal::sparse_sparse_product_with_pruning_impl<ColMajorMatrixLhs,ColMajorMatrixRhs,ResultType>(colLhs, colRhs, res, tolerance); // let's transpose the product to get a column x column product // typedef SparseMatrix<typename ResultType::Scalar> SparseTemporaryType; diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h index d068d62df..1ffa7d54e 100644 --- a/Eigen/src/SparseLU/SparseLU_Memory.h +++ b/Eigen/src/SparseLU/SparseLU_Memory.h @@ -70,7 +70,7 @@ Index SparseLUImpl<Scalar,Index>::expand(VectorType& vec, Index& length, Index if(num_expansions == 0 || keep_prev) new_len = length ; // First time allocate requested else - new_len = Index(alpha * length); + new_len = (std::max)(length+1,Index(alpha * length)); VectorType old_vec; // Temporary vector to hold the previous values if (nbElts > 0 ) @@ -107,7 +107,7 @@ Index SparseLUImpl<Scalar,Index>::expand(VectorType& vec, Index& length, Index do { alpha = (alpha + 1)/2; - new_len = Index(alpha * length); + new_len = (std::max)(length+1,Index(alpha * length)); #ifdef EIGEN_EXCEPTIONS try #endif diff --git a/Eigen/src/misc/SparseSolve.h b/Eigen/src/misc/SparseSolve.h index 244bb8ec7..05caa9266 100644 --- a/Eigen/src/misc/SparseSolve.h +++ b/Eigen/src/misc/SparseSolve.h @@ -54,8 +54,10 @@ template<typename _DecompositionType, typename Rhs> struct sparse_solve_retval_b static const int NbColsAtOnce = 4; int rhsCols = m_rhs.cols(); int size = m_rhs.rows(); - Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols); - Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,rhsCols); + // the temporary matrices do not need more columns than NbColsAtOnce: + int tmpCols = (std::min)(rhsCols, NbColsAtOnce); + Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,tmpCols); + Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmpX(size,tmpCols); for(int k=0; k<rhsCols; k+=NbColsAtOnce) { int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); diff --git a/Eigen/src/plugins/BlockMethods.h b/Eigen/src/plugins/BlockMethods.h index 3bc345211..9b7fdc4aa 100644 --- a/Eigen/src/plugins/BlockMethods.h +++ b/Eigen/src/plugins/BlockMethods.h @@ -119,13 +119,13 @@ inline const Block<const Derived, CRows, CCols> topRightCorner() const /** \returns an expression of a top-right corner of *this. * - * \tparam CRows number of rows in corner as specified at compile time - * \tparam CCols number of columns in corner as specified at compile time - * \param cRows number of rows in corner as specified at run time - * \param cCols number of columns in corner as specified at run time + * \tparam CRows number of rows in corner as specified at compile-time + * \tparam CCols number of columns in corner as specified at compile-time + * \param cRows number of rows in corner as specified at run-time + * \param cCols number of columns in corner as specified at run-time * - * This function is mainly useful for corners where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for corners where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a cRows should equal \a CRows unless * \a CRows is \a Dynamic, and the same for the number of columns. * @@ -198,13 +198,13 @@ inline const Block<const Derived, CRows, CCols> topLeftCorner() const /** \returns an expression of a top-left corner of *this. * - * \tparam CRows number of rows in corner as specified at compile time - * \tparam CCols number of columns in corner as specified at compile time - * \param cRows number of rows in corner as specified at run time - * \param cCols number of columns in corner as specified at run time + * \tparam CRows number of rows in corner as specified at compile-time + * \tparam CCols number of columns in corner as specified at compile-time + * \param cRows number of rows in corner as specified at run-time + * \param cCols number of columns in corner as specified at run-time * - * This function is mainly useful for corners where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for corners where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a cRows should equal \a CRows unless * \a CRows is \a Dynamic, and the same for the number of columns. * @@ -277,13 +277,13 @@ inline const Block<const Derived, CRows, CCols> bottomRightCorner() const /** \returns an expression of a bottom-right corner of *this. * - * \tparam CRows number of rows in corner as specified at compile time - * \tparam CCols number of columns in corner as specified at compile time - * \param cRows number of rows in corner as specified at run time - * \param cCols number of columns in corner as specified at run time + * \tparam CRows number of rows in corner as specified at compile-time + * \tparam CCols number of columns in corner as specified at compile-time + * \param cRows number of rows in corner as specified at run-time + * \param cCols number of columns in corner as specified at run-time * - * This function is mainly useful for corners where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for corners where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a cRows should equal \a CRows unless * \a CRows is \a Dynamic, and the same for the number of columns. * @@ -356,13 +356,13 @@ inline const Block<const Derived, CRows, CCols> bottomLeftCorner() const /** \returns an expression of a bottom-left corner of *this. * - * \tparam CRows number of rows in corner as specified at compile time - * \tparam CCols number of columns in corner as specified at compile time - * \param cRows number of rows in corner as specified at run time - * \param cCols number of columns in corner as specified at run time + * \tparam CRows number of rows in corner as specified at compile-time + * \tparam CCols number of columns in corner as specified at compile-time + * \param cRows number of rows in corner as specified at run-time + * \param cCols number of columns in corner as specified at run-time * - * This function is mainly useful for corners where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for corners where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a cRows should equal \a CRows unless * \a CRows is \a Dynamic, and the same for the number of columns. * @@ -410,7 +410,11 @@ inline ConstRowsBlockXpr topRows(Index n) const /** \returns a block consisting of the top rows of *this. * - * \tparam N the number of rows in the block + * \tparam N the number of rows in the block as specified at compile-time + * \param n the number of rows in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_topRows.cpp * Output: \verbinclude MatrixBase_template_int_topRows.out @@ -419,17 +423,17 @@ inline ConstRowsBlockXpr topRows(Index n) const */ template<int N> EIGEN_DEVICE_FUNC -inline typename NRowsBlockXpr<N>::Type topRows() +inline typename NRowsBlockXpr<N>::Type topRows(Index n = N) { - return typename NRowsBlockXpr<N>::Type(derived(), 0, 0, N, cols()); + return typename NRowsBlockXpr<N>::Type(derived(), 0, 0, n, cols()); } /** This is the const version of topRows<int>().*/ template<int N> EIGEN_DEVICE_FUNC -inline typename ConstNRowsBlockXpr<N>::Type topRows() const +inline typename ConstNRowsBlockXpr<N>::Type topRows(Index n = N) const { - return typename ConstNRowsBlockXpr<N>::Type(derived(), 0, 0, N, cols()); + return typename ConstNRowsBlockXpr<N>::Type(derived(), 0, 0, n, cols()); } @@ -458,7 +462,11 @@ inline ConstRowsBlockXpr bottomRows(Index n) const /** \returns a block consisting of the bottom rows of *this. * - * \tparam N the number of rows in the block + * \tparam N the number of rows in the block as specified at compile-time + * \param n the number of rows in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_bottomRows.cpp * Output: \verbinclude MatrixBase_template_int_bottomRows.out @@ -467,17 +475,17 @@ inline ConstRowsBlockXpr bottomRows(Index n) const */ template<int N> EIGEN_DEVICE_FUNC -inline typename NRowsBlockXpr<N>::Type bottomRows() +inline typename NRowsBlockXpr<N>::Type bottomRows(Index n = N) { - return typename NRowsBlockXpr<N>::Type(derived(), rows() - N, 0, N, cols()); + return typename NRowsBlockXpr<N>::Type(derived(), rows() - n, 0, n, cols()); } /** This is the const version of bottomRows<int>().*/ template<int N> EIGEN_DEVICE_FUNC -inline typename ConstNRowsBlockXpr<N>::Type bottomRows() const +inline typename ConstNRowsBlockXpr<N>::Type bottomRows(Index n = N) const { - return typename ConstNRowsBlockXpr<N>::Type(derived(), rows() - N, 0, N, cols()); + return typename ConstNRowsBlockXpr<N>::Type(derived(), rows() - n, 0, n, cols()); } @@ -485,7 +493,7 @@ inline typename ConstNRowsBlockXpr<N>::Type bottomRows() const /** \returns a block consisting of a range of rows of *this. * * \param startRow the index of the first row in the block - * \param numRows the number of rows in the block + * \param n the number of rows in the block * * Example: \include DenseBase_middleRows_int.cpp * Output: \verbinclude DenseBase_middleRows_int.out @@ -493,22 +501,26 @@ inline typename ConstNRowsBlockXpr<N>::Type bottomRows() const * \sa class Block, block(Index,Index,Index,Index) */ EIGEN_DEVICE_FUNC -inline RowsBlockXpr middleRows(Index startRow, Index numRows) +inline RowsBlockXpr middleRows(Index startRow, Index n) { - return RowsBlockXpr(derived(), startRow, 0, numRows, cols()); + return RowsBlockXpr(derived(), startRow, 0, n, cols()); } /** This is the const version of middleRows(Index,Index).*/ EIGEN_DEVICE_FUNC -inline ConstRowsBlockXpr middleRows(Index startRow, Index numRows) const +inline ConstRowsBlockXpr middleRows(Index startRow, Index n) const { - return ConstRowsBlockXpr(derived(), startRow, 0, numRows, cols()); + return ConstRowsBlockXpr(derived(), startRow, 0, n, cols()); } /** \returns a block consisting of a range of rows of *this. * - * \tparam N the number of rows in the block + * \tparam N the number of rows in the block as specified at compile-time * \param startRow the index of the first row in the block + * \param n the number of rows in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include DenseBase_template_int_middleRows.cpp * Output: \verbinclude DenseBase_template_int_middleRows.out @@ -517,17 +529,17 @@ inline ConstRowsBlockXpr middleRows(Index startRow, Index numRows) const */ template<int N> EIGEN_DEVICE_FUNC -inline typename NRowsBlockXpr<N>::Type middleRows(Index startRow) +inline typename NRowsBlockXpr<N>::Type middleRows(Index startRow, Index n = N) { - return typename NRowsBlockXpr<N>::Type(derived(), startRow, 0, N, cols()); + return typename NRowsBlockXpr<N>::Type(derived(), startRow, 0, n, cols()); } /** This is the const version of middleRows<int>().*/ template<int N> EIGEN_DEVICE_FUNC -inline typename ConstNRowsBlockXpr<N>::Type middleRows(Index startRow) const +inline typename ConstNRowsBlockXpr<N>::Type middleRows(Index startRow, Index n = N) const { - return typename ConstNRowsBlockXpr<N>::Type(derived(), startRow, 0, N, cols()); + return typename ConstNRowsBlockXpr<N>::Type(derived(), startRow, 0, n, cols()); } @@ -556,7 +568,11 @@ inline ConstColsBlockXpr leftCols(Index n) const /** \returns a block consisting of the left columns of *this. * - * \tparam N the number of columns in the block + * \tparam N the number of columns in the block as specified at compile-time + * \param n the number of columns in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_leftCols.cpp * Output: \verbinclude MatrixBase_template_int_leftCols.out @@ -565,17 +581,17 @@ inline ConstColsBlockXpr leftCols(Index n) const */ template<int N> EIGEN_DEVICE_FUNC -inline typename NColsBlockXpr<N>::Type leftCols() +inline typename NColsBlockXpr<N>::Type leftCols(Index n = N) { - return typename NColsBlockXpr<N>::Type(derived(), 0, 0, rows(), N); + return typename NColsBlockXpr<N>::Type(derived(), 0, 0, rows(), n); } /** This is the const version of leftCols<int>().*/ template<int N> EIGEN_DEVICE_FUNC -inline typename ConstNColsBlockXpr<N>::Type leftCols() const +inline typename ConstNColsBlockXpr<N>::Type leftCols(Index n = N) const { - return typename ConstNColsBlockXpr<N>::Type(derived(), 0, 0, rows(), N); + return typename ConstNColsBlockXpr<N>::Type(derived(), 0, 0, rows(), n); } @@ -604,7 +620,11 @@ inline ConstColsBlockXpr rightCols(Index n) const /** \returns a block consisting of the right columns of *this. * - * \tparam N the number of columns in the block + * \tparam N the number of columns in the block as specified at compile-time + * \param n the number of columns in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_rightCols.cpp * Output: \verbinclude MatrixBase_template_int_rightCols.out @@ -613,17 +633,17 @@ inline ConstColsBlockXpr rightCols(Index n) const */ template<int N> EIGEN_DEVICE_FUNC -inline typename NColsBlockXpr<N>::Type rightCols() +inline typename NColsBlockXpr<N>::Type rightCols(Index n = N) { - return typename NColsBlockXpr<N>::Type(derived(), 0, cols() - N, rows(), N); + return typename NColsBlockXpr<N>::Type(derived(), 0, cols() - n, rows(), n); } /** This is the const version of rightCols<int>().*/ template<int N> EIGEN_DEVICE_FUNC -inline typename ConstNColsBlockXpr<N>::Type rightCols() const +inline typename ConstNColsBlockXpr<N>::Type rightCols(Index n = N) const { - return typename ConstNColsBlockXpr<N>::Type(derived(), 0, cols() - N, rows(), N); + return typename ConstNColsBlockXpr<N>::Type(derived(), 0, cols() - n, rows(), n); } @@ -653,8 +673,12 @@ inline ConstColsBlockXpr middleCols(Index startCol, Index numCols) const /** \returns a block consisting of a range of columns of *this. * - * \tparam N the number of columns in the block + * \tparam N the number of columns in the block as specified at compile-time * \param startCol the index of the first column in the block + * \param n the number of columns in the block as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include DenseBase_template_int_middleCols.cpp * Output: \verbinclude DenseBase_template_int_middleCols.out @@ -663,17 +687,17 @@ inline ConstColsBlockXpr middleCols(Index startCol, Index numCols) const */ template<int N> EIGEN_DEVICE_FUNC -inline typename NColsBlockXpr<N>::Type middleCols(Index startCol) +inline typename NColsBlockXpr<N>::Type middleCols(Index startCol, Index n = N) { - return typename NColsBlockXpr<N>::Type(derived(), 0, startCol, rows(), N); + return typename NColsBlockXpr<N>::Type(derived(), 0, startCol, rows(), n); } /** This is the const version of middleCols<int>().*/ template<int N> EIGEN_DEVICE_FUNC -inline typename ConstNColsBlockXpr<N>::Type middleCols(Index startCol) const +inline typename ConstNColsBlockXpr<N>::Type middleCols(Index startCol, Index n = N) const { - return typename ConstNColsBlockXpr<N>::Type(derived(), 0, startCol, rows(), N); + return typename ConstNColsBlockXpr<N>::Type(derived(), 0, startCol, rows(), n); } @@ -711,15 +735,15 @@ inline const Block<const Derived, BlockRows, BlockCols> block(Index startRow, In /** \returns an expression of a block in *this. * - * \tparam BlockRows number of rows in block as specified at compile time - * \tparam BlockCols number of columns in block as specified at compile time + * \tparam BlockRows number of rows in block as specified at compile-time + * \tparam BlockCols number of columns in block as specified at compile-time * \param startRow the first row in the block * \param startCol the first column in the block - * \param blockRows number of rows in block as specified at run time - * \param blockCols number of columns in block as specified at run time + * \param blockRows number of rows in block as specified at run-time + * \param blockCols number of columns in block as specified at run-time * - * This function is mainly useful for blocks where the number of rows is specified at compile time - * and the number of columns is specified at run time, or vice versa. The compile-time and run-time + * This function is mainly useful for blocks where the number of rows is specified at compile-time + * and the number of columns is specified at run-time, or vice versa. The compile-time and run-time * information should not contradict. In other words, \a blockRows should equal \a BlockRows unless * \a BlockRows is \a Dynamic, and the same for the number of columns. * @@ -786,7 +810,7 @@ inline ConstRowXpr row(Index i) const * \only_for_vectors * * \param start the first coefficient in the segment - * \param vecSize the number of coefficients in the segment + * \param n the number of coefficients in the segment * * Example: \include MatrixBase_segment_int_int.cpp * Output: \verbinclude MatrixBase_segment_int_int.out @@ -798,26 +822,26 @@ inline ConstRowXpr row(Index i) const * \sa class Block, segment(Index) */ EIGEN_DEVICE_FUNC -inline SegmentReturnType segment(Index start, Index vecSize) +inline SegmentReturnType segment(Index start, Index n) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return SegmentReturnType(derived(), start, vecSize); + return SegmentReturnType(derived(), start, n); } /** This is the const version of segment(Index,Index).*/ EIGEN_DEVICE_FUNC -inline ConstSegmentReturnType segment(Index start, Index vecSize) const +inline ConstSegmentReturnType segment(Index start, Index n) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return ConstSegmentReturnType(derived(), start, vecSize); + return ConstSegmentReturnType(derived(), start, n); } /** \returns a dynamic-size expression of the first coefficients of *this. * * \only_for_vectors * - * \param vecSize the number of coefficients in the block + * \param n the number of coefficients in the segment * * Example: \include MatrixBase_start_int.cpp * Output: \verbinclude MatrixBase_start_int.out @@ -829,26 +853,25 @@ inline ConstSegmentReturnType segment(Index start, Index vecSize) const * \sa class Block, block(Index,Index) */ EIGEN_DEVICE_FUNC -inline SegmentReturnType head(Index vecSize) +inline SegmentReturnType head(Index n) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return SegmentReturnType(derived(), 0, vecSize); + return SegmentReturnType(derived(), 0, n); } /** This is the const version of head(Index).*/ EIGEN_DEVICE_FUNC -inline ConstSegmentReturnType - head(Index vecSize) const +inline ConstSegmentReturnType head(Index n) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return ConstSegmentReturnType(derived(), 0, vecSize); + return ConstSegmentReturnType(derived(), 0, n); } /** \returns a dynamic-size expression of the last coefficients of *this. * * \only_for_vectors * - * \param vecSize the number of coefficients in the block + * \param n the number of coefficients in the segment * * Example: \include MatrixBase_end_int.cpp * Output: \verbinclude MatrixBase_end_int.out @@ -860,102 +883,113 @@ inline ConstSegmentReturnType * \sa class Block, block(Index,Index) */ EIGEN_DEVICE_FUNC -inline SegmentReturnType tail(Index vecSize) +inline SegmentReturnType tail(Index n) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return SegmentReturnType(derived(), this->size() - vecSize, vecSize); + return SegmentReturnType(derived(), this->size() - n, n); } /** This is the const version of tail(Index).*/ EIGEN_DEVICE_FUNC -inline ConstSegmentReturnType tail(Index vecSize) const +inline ConstSegmentReturnType tail(Index n) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return ConstSegmentReturnType(derived(), this->size() - vecSize, vecSize); + return ConstSegmentReturnType(derived(), this->size() - n, n); } /** \returns a fixed-size expression of a segment (i.e. a vector block) in \c *this * * \only_for_vectors * - * The template parameter \a Size is the number of coefficients in the block + * \tparam N the number of coefficients in the segment as specified at compile-time + * \param start the index of the first element in the segment + * \param n the number of coefficients in the segment as specified at compile-time * - * \param start the index of the first element of the sub-vector + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_segment.cpp * Output: \verbinclude MatrixBase_template_int_segment.out * * \sa class Block */ -template<int Size> +template<int N> EIGEN_DEVICE_FUNC -inline typename FixedSegmentReturnType<Size>::Type segment(Index start) +inline typename FixedSegmentReturnType<N>::Type segment(Index start, Index n = N) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename FixedSegmentReturnType<Size>::Type(derived(), start); + return typename FixedSegmentReturnType<N>::Type(derived(), start, n); } /** This is the const version of segment<int>(Index).*/ -template<int Size> +template<int N> EIGEN_DEVICE_FUNC -inline typename ConstFixedSegmentReturnType<Size>::Type segment(Index start) const +inline typename ConstFixedSegmentReturnType<N>::Type segment(Index start, Index n = N) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename ConstFixedSegmentReturnType<Size>::Type(derived(), start); + return typename ConstFixedSegmentReturnType<N>::Type(derived(), start, n); } /** \returns a fixed-size expression of the first coefficients of *this. * * \only_for_vectors * - * The template parameter \a Size is the number of coefficients in the block + * \tparam N the number of coefficients in the segment as specified at compile-time + * \param n the number of coefficients in the segment as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_start.cpp * Output: \verbinclude MatrixBase_template_int_start.out * * \sa class Block */ -template<int Size> +template<int N> EIGEN_DEVICE_FUNC -inline typename FixedSegmentReturnType<Size>::Type head() +inline typename FixedSegmentReturnType<N>::Type head(Index n = N) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename FixedSegmentReturnType<Size>::Type(derived(), 0); + return typename FixedSegmentReturnType<N>::Type(derived(), 0, n); } /** This is the const version of head<int>().*/ -template<int Size> +template<int N> EIGEN_DEVICE_FUNC -inline typename ConstFixedSegmentReturnType<Size>::Type head() const +inline typename ConstFixedSegmentReturnType<N>::Type head(Index n = N) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename ConstFixedSegmentReturnType<Size>::Type(derived(), 0); + return typename ConstFixedSegmentReturnType<N>::Type(derived(), 0, n); } /** \returns a fixed-size expression of the last coefficients of *this. * * \only_for_vectors * - * The template parameter \a Size is the number of coefficients in the block + * \tparam N the number of coefficients in the segment as specified at compile-time + * \param n the number of coefficients in the segment as specified at run-time + * + * The compile-time and run-time information should not contradict. In other words, + * \a n should equal \a N unless \a N is \a Dynamic. * * Example: \include MatrixBase_template_int_end.cpp * Output: \verbinclude MatrixBase_template_int_end.out * * \sa class Block */ -template<int Size> +template<int N> EIGEN_DEVICE_FUNC -inline typename FixedSegmentReturnType<Size>::Type tail() +inline typename FixedSegmentReturnType<N>::Type tail(Index n = N) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename FixedSegmentReturnType<Size>::Type(derived(), size() - Size); + return typename FixedSegmentReturnType<N>::Type(derived(), size() - n); } /** This is the const version of tail<int>.*/ -template<int Size> +template<int N> EIGEN_DEVICE_FUNC -inline typename ConstFixedSegmentReturnType<Size>::Type tail() const +inline typename ConstFixedSegmentReturnType<N>::Type tail(Index n = N) const { EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename ConstFixedSegmentReturnType<Size>::Type(derived(), size() - Size); + return typename ConstFixedSegmentReturnType<N>::Type(derived(), size() - n); } |