diff options
Diffstat (limited to 'Eigen/src/Core')
-rw-r--r-- | Eigen/src/Core/Block.h | 249 | ||||
-rw-r--r-- | Eigen/src/Core/CwiseUnaryView.h | 15 | ||||
-rw-r--r-- | Eigen/src/Core/Dot.h | 138 | ||||
-rw-r--r-- | Eigen/src/Core/MathFunctions.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/Matrix.h | 4 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 40 | ||||
-rw-r--r-- | Eigen/src/Core/NumTraits.h | 8 | ||||
-rw-r--r-- | Eigen/src/Core/Product.h | 17 | ||||
-rw-r--r-- | Eigen/src/Core/SolveTriangular.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/VectorBlock.h | 311 | ||||
-rw-r--r-- | Eigen/src/Core/util/Constants.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 3 | ||||
-rw-r--r-- | Eigen/src/Core/util/Macros.h | 16 | ||||
-rw-r--r-- | Eigen/src/Core/util/Memory.h | 12 | ||||
-rw-r--r-- | Eigen/src/Core/util/XprHelper.h | 6 |
15 files changed, 530 insertions, 295 deletions
diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h index de4268344..382518696 100644 --- a/Eigen/src/Core/Block.h +++ b/Eigen/src/Core/Block.h @@ -271,13 +271,19 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess> inline int stride(void) const { return m_matrix.stride(); } + #ifndef __SUNPRO_CC + // FIXME sunstudio is not friendly with the above friend... protected: + #endif + #ifndef EIGEN_PARSED_BY_DOXYGEN /** \internal used by allowAligned() */ inline Block(const MatrixType& matrix, const Scalar* data, int blockRows, int blockCols) : Base(data, blockRows, blockCols), m_matrix(matrix) {} + #endif + protected: const typename MatrixType::Nested m_matrix; }; @@ -314,249 +320,6 @@ inline const typename BlockReturnType<Derived>::Type MatrixBase<Derived> return typename BlockReturnType<Derived>::Type(derived(), startRow, startCol, blockRows, blockCols); } -/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this. - * - * \only_for_vectors - * - * \addexample SegmentIntInt \label How to reference a sub-vector (dynamic size) - * - * \param start the first coefficient in the segment - * \param size the number of coefficients in the segment - * - * Example: \include MatrixBase_segment_int_int.cpp - * Output: \verbinclude MatrixBase_segment_int_int.out - * - * \note Even though the returned expression has dynamic size, in the case - * when it is applied to a fixed-size vector, it inherits a fixed maximal size, - * which means that evaluating it does not cause a dynamic memory allocation. - * - * \sa class Block, segment(int) - */ -template<typename Derived> -inline typename BlockReturnType<Derived>::SubVectorType MatrixBase<Derived> - ::segment(int start, int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename BlockReturnType<Derived>::SubVectorType(derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** This is the const version of segment(int,int).*/ -template<typename Derived> -inline const typename BlockReturnType<Derived>::SubVectorType -MatrixBase<Derived>::segment(int start, int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename BlockReturnType<Derived>::SubVectorType(derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** \returns a dynamic-size expression of the first coefficients of *this. - * - * \only_for_vectors - * - * \param size the number of coefficients in the block - * - * \addexample BlockInt \label How to reference a sub-vector (fixed-size) - * - * Example: \include MatrixBase_start_int.cpp - * Output: \verbinclude MatrixBase_start_int.out - * - * \note Even though the returned expression has dynamic size, in the case - * when it is applied to a fixed-size vector, it inherits a fixed maximal size, - * which means that evaluating it does not cause a dynamic memory allocation. - * - * \sa class Block, block(int,int) - */ -template<typename Derived> -inline typename BlockReturnType<Derived,Dynamic>::SubVectorType -MatrixBase<Derived>::start(int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, - RowsAtCompileTime == 1 ? 1 : Dynamic, - ColsAtCompileTime == 1 ? 1 : Dynamic> - (derived(), 0, 0, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** This is the const version of start(int).*/ -template<typename Derived> -inline const typename BlockReturnType<Derived,Dynamic>::SubVectorType -MatrixBase<Derived>::start(int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, - RowsAtCompileTime == 1 ? 1 : Dynamic, - ColsAtCompileTime == 1 ? 1 : Dynamic> - (derived(), 0, 0, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** \returns a dynamic-size expression of the last coefficients of *this. - * - * \only_for_vectors - * - * \param size the number of coefficients in the block - * - * \addexample BlockEnd \label How to reference the end of a vector (fixed-size) - * - * Example: \include MatrixBase_end_int.cpp - * Output: \verbinclude MatrixBase_end_int.out - * - * \note Even though the returned expression has dynamic size, in the case - * when it is applied to a fixed-size vector, it inherits a fixed maximal size, - * which means that evaluating it does not cause a dynamic memory allocation. - * - * \sa class Block, block(int,int) - */ -template<typename Derived> -inline typename BlockReturnType<Derived,Dynamic>::SubVectorType -MatrixBase<Derived>::end(int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, - RowsAtCompileTime == 1 ? 1 : Dynamic, - ColsAtCompileTime == 1 ? 1 : Dynamic> - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - size, - ColsAtCompileTime == 1 ? 0 : cols() - size, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** This is the const version of end(int).*/ -template<typename Derived> -inline const typename BlockReturnType<Derived,Dynamic>::SubVectorType -MatrixBase<Derived>::end(int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, - RowsAtCompileTime == 1 ? 1 : Dynamic, - ColsAtCompileTime == 1 ? 1 : Dynamic> - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - size, - ColsAtCompileTime == 1 ? 0 : cols() - size, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** \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 - * - * \param start the index of the first element of the sub-vector - * - * Example: \include MatrixBase_template_int_segment.cpp - * Output: \verbinclude MatrixBase_template_int_segment.out - * - * \sa class Block - */ -template<typename Derived> -template<int Size> -inline typename BlockReturnType<Derived,Size>::SubVectorType -MatrixBase<Derived>::segment(int start) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size), - (ColsAtCompileTime == 1 ? 1 : Size)> - (derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start); -} - -/** This is the const version of segment<int>(int).*/ -template<typename Derived> -template<int Size> -inline const typename BlockReturnType<Derived,Size>::SubVectorType -MatrixBase<Derived>::segment(int start) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size), - (ColsAtCompileTime == 1 ? 1 : Size)> - (derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start); -} - -/** \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 - * - * \addexample BlockStart \label How to reference the start of a vector (fixed-size) - * - * Example: \include MatrixBase_template_int_start.cpp - * Output: \verbinclude MatrixBase_template_int_start.out - * - * \sa class Block - */ -template<typename Derived> -template<int Size> -inline typename BlockReturnType<Derived,Size>::SubVectorType -MatrixBase<Derived>::start() -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size), - (ColsAtCompileTime == 1 ? 1 : Size)>(derived(), 0, 0); -} - -/** This is the const version of start<int>().*/ -template<typename Derived> -template<int Size> -inline const typename BlockReturnType<Derived,Size>::SubVectorType -MatrixBase<Derived>::start() const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, (RowsAtCompileTime == 1 ? 1 : Size), - (ColsAtCompileTime == 1 ? 1 : Size)>(derived(), 0, 0); -} - -/** \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 - * - * Example: \include MatrixBase_template_int_end.cpp - * Output: \verbinclude MatrixBase_template_int_end.out - * - * \sa class Block - */ -template<typename Derived> -template<int Size> -inline typename BlockReturnType<Derived,Size>::SubVectorType -MatrixBase<Derived>::end() -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size, - ColsAtCompileTime == 1 ? 1 : Size> - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - Size, - ColsAtCompileTime == 1 ? 0 : cols() - Size); -} - -/** This is the const version of end<int>.*/ -template<typename Derived> -template<int Size> -inline const typename BlockReturnType<Derived,Size>::SubVectorType -MatrixBase<Derived>::end() const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block<Derived, RowsAtCompileTime == 1 ? 1 : Size, - ColsAtCompileTime == 1 ? 1 : Size> - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - Size, - ColsAtCompileTime == 1 ? 0 : cols() - Size); -} - /** \returns a dynamic-size expression of a corner of *this. * * \param type the type of corner. Can be \a Eigen::TopLeft, \a Eigen::TopRight, diff --git a/Eigen/src/Core/CwiseUnaryView.h b/Eigen/src/Core/CwiseUnaryView.h index c8fd98ea1..4785f01f4 100644 --- a/Eigen/src/Core/CwiseUnaryView.h +++ b/Eigen/src/Core/CwiseUnaryView.h @@ -42,13 +42,13 @@ struct ei_traits<CwiseUnaryView<ViewOp, MatrixType> > : ei_traits<MatrixType> { typedef typename ei_result_of< - ViewOp(typename MatrixType::Scalar) + ViewOp(typename ei_traits<MatrixType>::Scalar) >::type Scalar; typedef typename MatrixType::Nested MatrixTypeNested; - typedef typename ei_unref<MatrixTypeNested>::type _MatrixTypeNested; + typedef typename ei_cleantype<MatrixTypeNested>::type _MatrixTypeNested; enum { - Flags = (_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)), - CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits<ViewOp>::Cost + Flags = (ei_traits<_MatrixTypeNested>::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)), + CoeffReadCost = ei_traits<_MatrixTypeNested>::CoeffReadCost + ei_functor_traits<ViewOp>::Cost }; }; @@ -62,7 +62,7 @@ class CwiseUnaryView : ei_no_assignment_operator, inline CwiseUnaryView(const MatrixType& mat, const ViewOp& func = ViewOp()) : m_matrix(mat), m_functor(func) {} - + EIGEN_INHERIT_ASSIGNMENT_OPERATORS(CwiseUnaryView) EIGEN_STRONG_INLINE int rows() const { return m_matrix.rows(); } @@ -77,7 +77,7 @@ class CwiseUnaryView : ei_no_assignment_operator, { return m_functor(m_matrix.coeff(index)); } - + EIGEN_STRONG_INLINE Scalar& coeffRef(int row, int col) { return m_functor(m_matrix.const_cast_derived().coeffRef(row, col)); @@ -89,7 +89,8 @@ class CwiseUnaryView : ei_no_assignment_operator, } protected: - const typename MatrixType::Nested m_matrix; + // FIXME changed from MatrixType::Nested because of a weird compilation error with sun CC + const typename ei_nested<MatrixType>::type m_matrix; const ViewOp m_functor; }; diff --git a/Eigen/src/Core/Dot.h b/Eigen/src/Core/Dot.h index 3861984da..4f185ea5b 100644 --- a/Eigen/src/Core/Dot.h +++ b/Eigen/src/Core/Dot.h @@ -295,7 +295,7 @@ inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real MatrixBase< /** \returns the \em l2 norm of \c *this using a numerically more stable * algorithm. * - * \sa norm(), dot(), squaredNorm() + * \sa norm(), dot(), squaredNorm(), blueNorm() */ template<typename Derived> inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real @@ -304,6 +304,142 @@ MatrixBase<Derived>::stableNorm() const return this->cwise().abs().redux(ei_scalar_hypot_op<RealScalar>()); } +/** \internal Computes ibeta^iexp by binary expansion of iexp, + * exact if ibeta is the machine base */ +template<typename T> inline T bexp(int ibeta, int iexp) +{ + T tbeta = T(ibeta); + T res = 1.0; + int n = iexp; + if (n<0) + { + n = - n; + tbeta = 1.0/tbeta; + } + for(;;) + { + if ((n % 2)==0) + res = res * tbeta; + n = n/2; + if (n==0) return res; + tbeta = tbeta*tbeta; + } + return res; +} + +/** \returns the \em l2 norm of \c *this using the Blue's algorithm. + * A Portable Fortran Program to Find the Euclidean Norm of a Vector, + * ACM TOMS, Vol 4, Issue 1, 1978. + * + * \sa norm(), dot(), squaredNorm(), stableNorm() + */ +template<typename Derived> +inline typename NumTraits<typename ei_traits<Derived>::Scalar>::Real +MatrixBase<Derived>::blueNorm() const +{ + static int nmax; + static Scalar b1, b2, s1m, s2m, overfl, rbig, relerr; + int n; + Scalar ax, abig, amed, asml; + + if(nmax <= 0) + { + int nbig, ibeta, it, iemin, iemax, iexp; + Scalar abig, eps; + // This program calculates the machine-dependent constants + // bl, b2, slm, s2m, relerr overfl, nmax + // from the "basic" machine-dependent numbers + // nbig, ibeta, it, iemin, iemax, rbig. + // The following define the basic machine-dependent constants. + // For portability, the PORT subprograms "ilmaeh" and "rlmach" + // are used. For any specific computer, each of the assignment + // statements can be replaced + nbig = std::numeric_limits<int>::max(); // largest integer + ibeta = NumTraits<Scalar>::Base; // base for floating-point numbers + it = NumTraits<Scalar>::Mantissa; // number of base-beta digits in mantissa + iemin = std::numeric_limits<Scalar>::min_exponent; // minimum exponent + iemax = std::numeric_limits<Scalar>::max_exponent; // maximum exponent + rbig = std::numeric_limits<Scalar>::max(); // largest floating-point number + + // Check the basic machine-dependent constants. + if(iemin > 1 - 2*it || 1+it>iemax || (it==2 && ibeta<5) + || (it<=4 && ibeta <= 3 ) || it<2) + { + ei_assert(false && "the algorithm cannot be guaranteed on this computer"); + } + iexp = -((1-iemin)/2); + b1 = bexp<Scalar>(ibeta, iexp); // lower boundary of midrange + iexp = (iemax + 1 - it)/2; + b2 = bexp<Scalar>(ibeta,iexp); // upper boundary of midrange + + iexp = (2-iemin)/2; + s1m = bexp<Scalar>(ibeta,iexp); // scaling factor for lower range + iexp = - ((iemax+it)/2); + s2m = bexp<Scalar>(ibeta,iexp); // scaling factor for upper range + + overfl = rbig*s2m; // overfow boundary for abig + eps = bexp<Scalar>(ibeta, 1-it); + relerr = ei_sqrt(eps); // tolerance for neglecting asml + abig = 1.0/eps - 1.0; + if (Scalar(nbig)>abig) nmax = abig; // largest safe n + else nmax = nbig; + } + n = size(); + if(n==0) + return 0; + asml = Scalar(0); + amed = Scalar(0); + abig = Scalar(0); + for(int j=0; j<n; ++j) + { + ax = ei_abs(coeff(j)); + if(ax > b2) abig += ei_abs2(ax*s2m); + else if(ax < b2) asml += ei_abs2(ax*s1m); + else amed += ei_abs2(ax); + } + if(abig > Scalar(0)) + { + abig = ei_sqrt(abig); + if(abig > overfl) + { + ei_assert(false && "overflow"); + return rbig; + } + if(amed > Scalar(0)) + { + abig = abig/s2m; + amed = ei_sqrt(amed); + } + else + { + return abig/s2m; + } + + } + else if(asml > Scalar(0)) + { + if (amed > Scalar(0)) + { + abig = ei_sqrt(amed); + amed = ei_sqrt(asml) / s1m; + } + else + { + return ei_sqrt(asml)/s1m; + } + } + else + { + return ei_sqrt(amed); + } + asml = std::min(abig, amed); + abig = std::max(abig, amed); + if(asml <= abig*relerr) + return abig; + else + return abig * ei_sqrt(Scalar(1) + ei_abs2(asml/abig)); +} + /** \returns an expression of the quotient of *this by its own norm. * * \only_for_vectors diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 8cb1f7b40..3527042f9 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -151,7 +151,7 @@ inline bool ei_isApproxOrLessThan(float a, float b, float prec = precision<float *** double *** **************/ -template<> inline double precision<double>() { return 1e-11; } +template<> inline double precision<double>() { return 1e-12; } template<> inline double machine_epsilon<double>() { return 2.220e-16; } inline double ei_real(double x) { return x; } diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 18c8f969a..2b4c4634a 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -124,7 +124,7 @@ class Matrix { public: EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix) - + enum { Options = _Options }; friend class Eigen::Map<Matrix, Unaligned>; typedef class Eigen::Map<Matrix, Unaligned> UnalignedMapType; @@ -218,7 +218,7 @@ class Matrix * * This method is intended for dynamic-size matrices, although it is legal to call it on any * matrix as long as fixed dimensions are left unchanged. If you only want to change the number - * of rows and/or of columns, you can use resize(NoChange_t, int), resize(int, NoChange_t). + * of rows and/or of columns, you can use resize(NoChange_t, int), resize(int, NoChange_t). * * If the current number of coefficients of \c *this exactly matches the * product \a rows * \a cols, then no memory allocation is performed and diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index cfa41978e..0a3342758 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -163,10 +163,14 @@ template<typename Derived> class MatrixBase * constructed from this one. See the \ref flags "list of flags". */ - CoeffReadCost = ei_traits<Derived>::CoeffReadCost + CoeffReadCost = ei_traits<Derived>::CoeffReadCost, /**< This is a rough measure of how expensive it is to read one coefficient from * this expression. */ + +#ifndef EIGEN_PARSED_BY_DOXYGEN + _HasDirectAccess = (int(Flags)&DirectAccessBit) ? 1 : 0 // workaround sunCC +#endif }; /** Default constructor. Just checks at compile-time for self-consistency of the flags. */ @@ -230,7 +234,7 @@ template<typename Derived> class MatrixBase /** \internal the return type of coeff() */ - typedef typename ei_meta_if<bool(int(Flags)&DirectAccessBit), const Scalar&, Scalar>::ret CoeffReturnType; + typedef typename ei_meta_if<_HasDirectAccess, const Scalar&, Scalar>::ret CoeffReturnType; /** \internal Represents a matrix with all coefficients equal to one another*/ typedef CwiseNullaryOp<ei_scalar_constant_op<Scalar>,Derived> ConstantReturnType; @@ -426,8 +430,9 @@ template<typename Derived> class MatrixBase template<typename OtherDerived> Scalar dot(const MatrixBase<OtherDerived>& other) const; RealScalar squaredNorm() const; - RealScalar norm() const; - RealScalar stableNorm() const; + RealScalar norm() const; + RealScalar stableNorm() const; + RealScalar blueNorm() const; const PlainMatrixType normalized() const; void normalize(); @@ -450,14 +455,14 @@ template<typename Derived> class MatrixBase const typename BlockReturnType<Derived>::Type block(int startRow, int startCol, int blockRows, int blockCols) const; - typename BlockReturnType<Derived>::SubVectorType segment(int start, int size); - const typename BlockReturnType<Derived>::SubVectorType segment(int start, int size) const; + VectorBlock<Derived> segment(int start, int size); + const VectorBlock<Derived> segment(int start, int size) const; - typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size); - const typename BlockReturnType<Derived,Dynamic>::SubVectorType start(int size) const; + VectorBlock<Derived> start(int size); + const VectorBlock<Derived> start(int size) const; - typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size); - const typename BlockReturnType<Derived,Dynamic>::SubVectorType end(int size) const; + VectorBlock<Derived> end(int size); + const VectorBlock<Derived> end(int size) const; typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols); const typename BlockReturnType<Derived>::Type corner(CornerType type, int cRows, int cCols) const; @@ -472,14 +477,14 @@ template<typename Derived> class MatrixBase template<int CRows, int CCols> const typename BlockReturnType<Derived, CRows, CCols>::Type corner(CornerType type) const; - template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType start(void); - template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType start() const; + template<int Size> VectorBlock<Derived,Size> start(void); + template<int Size> const VectorBlock<Derived,Size> start() const; - template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType end(); - template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType end() const; + template<int Size> VectorBlock<Derived,Size> end(); + template<int Size> const VectorBlock<Derived,Size> end() const; - template<int Size> typename BlockReturnType<Derived,Size>::SubVectorType segment(int start); - template<int Size> const typename BlockReturnType<Derived,Size>::SubVectorType segment(int start) const; + template<int Size> VectorBlock<Derived,Size> segment(int start); + template<int Size> const VectorBlock<Derived,Size> segment(int start) const; Diagonal<Derived,0> diagonal(); const Diagonal<Derived,0> diagonal() const; @@ -696,6 +701,7 @@ template<typename Derived> class MatrixBase const PartialLU<PlainMatrixType> partialLu() const; const PlainMatrixType inverse() const; void computeInverse(PlainMatrixType *result) const; + bool computeInverseWithCheck( PlainMatrixType *result ) const; Scalar determinant() const; /////////// Cholesky module /////////// @@ -705,7 +711,7 @@ template<typename Derived> class MatrixBase /////////// QR module /////////// - const QR<PlainMatrixType> qr() const; + const HouseholderQR<PlainMatrixType> householderQr() const; EigenvaluesReturnType eigenvalues() const; RealScalar operatorNorm() const; diff --git a/Eigen/src/Core/NumTraits.h b/Eigen/src/Core/NumTraits.h index 24afe54c5..dec07a036 100644 --- a/Eigen/src/Core/NumTraits.h +++ b/Eigen/src/Core/NumTraits.h @@ -70,7 +70,9 @@ template<> struct NumTraits<float> HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1 + MulCost = 1, + Base = 2, + Mantissa = 23 }; }; @@ -83,7 +85,9 @@ template<> struct NumTraits<double> HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1 + MulCost = 1, + Base = 2, + Mantissa = 52 }; }; diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 44e3f606e..f0c4480de 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -155,13 +155,14 @@ template<typename Lhs, typename Rhs> struct ei_product_mode typedef typename ei_blas_traits<Lhs>::ActualXprType ActualLhs; typedef typename ei_blas_traits<Rhs>::ActualXprType ActualRhs; enum{ - - value = Lhs::MaxColsAtCompileTime == Dynamic - && ( Lhs::MaxRowsAtCompileTime == Dynamic - || Rhs::MaxColsAtCompileTime == Dynamic ) - && (!(Rhs::IsVectorAtCompileTime && (Lhs::Flags&RowMajorBit) && (!(ActualLhs::Flags&DirectAccessBit)))) - && (!(Lhs::IsVectorAtCompileTime && (!(Rhs::Flags&RowMajorBit)) && (!(ActualRhs::Flags&DirectAccessBit)))) - && (ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret) + // workaround sun studio: + LhsIsVectorAtCompileTime = ei_traits<Lhs>::ColsAtCompileTime==1 || ei_traits<Rhs>::ColsAtCompileTime==1, + value = ei_traits<Lhs>::MaxColsAtCompileTime == Dynamic + && ( ei_traits<Lhs>::MaxRowsAtCompileTime == Dynamic + || ei_traits<Rhs>::MaxColsAtCompileTime == Dynamic ) + && (!(Rhs::IsVectorAtCompileTime && (ei_traits<Lhs>::Flags&RowMajorBit) && (!(ei_traits<Lhs>::Flags&DirectAccessBit)))) + && (!(LhsIsVectorAtCompileTime && (!(ei_traits<Rhs>::Flags&RowMajorBit)) && (!(ei_traits<Rhs>::Flags&DirectAccessBit)))) + && (ei_is_same_type<typename ei_traits<Lhs>::Scalar, typename ei_traits<Rhs>::Scalar>::ret) ? CacheFriendlyProduct : NormalProduct }; }; @@ -577,7 +578,7 @@ struct ei_product_packet_impl<ColMajor, Dynamic, Lhs, Rhs, PacketScalar, LoadMod // Forward declarations template<typename Scalar, bool ConjugateLhs, bool ConjugateRhs> -void ei_cache_friendly_product( +static void ei_cache_friendly_product( int _rows, int _cols, int depth, bool _lhsRowMajor, const Scalar* _lhs, int _lhsStride, bool _rhsRowMajor, const Scalar* _rhs, int _rhsStride, diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index 76e4289de..cb162ca91 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -180,7 +180,7 @@ struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,false> { template<typename Lhs, typename Rhs, int Mode, int Index, int Size> struct ei_triangular_solver_unroller<Lhs,Rhs,Mode,Index,Size,true> { - static void run(const Lhs& lhs, Rhs& rhs) {} + static void run(const Lhs&, Rhs&) {} }; template<typename Lhs, typename Rhs, int Mode, int StorageOrder> diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h new file mode 100644 index 000000000..7ce5977f6 --- /dev/null +++ b/Eigen/src/Core/VectorBlock.h @@ -0,0 +1,311 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_VECTORBLOCK_H +#define EIGEN_VECTORBLOCK_H + +/** \class VectorBlock + * + * \brief Expression of a fixed-size or dynamic-size sub-vector + * + * \param VectorType the type of the object in which we are taking a sub-vector + * \param Size size of the sub-vector we are taking at compile time (optional) + * \param _PacketAccess allows to enforce aligned loads and stores if set to ForceAligned. + * The default is AsRequested. This parameter is internaly used by Eigen + * in expressions such as \code mat.segment() += other; \endcode and most of + * the time this is the only way it is used. + * + * This class represents an expression of either a fixed-size or dynamic-size sub-vector. + * It is the return type of MatrixBase::segment(int,int) and MatrixBase::segment<int>(int) and + * most of the time this is the only way it is used. + * + * However, if you want to directly maniputate sub-vector expressions, + * for instance if you want to write a function returning such an expression, you + * will need to use this class. + * + * Here is an example illustrating the dynamic case: + * \include class_VectorBlock.cpp + * Output: \verbinclude class_VectorBlock.out + * + * \note Even though this expression has dynamic size, in the case where \a VectorType + * has fixed size, this expression inherits a fixed maximal size which means that evaluating + * it does not cause a dynamic memory allocation. + * + * Here is an example illustrating the fixed-size case: + * \include class_FixedVectorBlock.cpp + * Output: \verbinclude class_FixedVectorBlock.out + * + * \sa class Block, MatrixBase::segment(int,int,int,int), MatrixBase::segment(int,int) + */ +template<typename VectorType, int Size, int _PacketAccess> +struct ei_traits<VectorBlock<VectorType, Size, _PacketAccess> > + : public ei_traits<Block<VectorType, + ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, + ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size, + _PacketAccess> > +{ +}; + +template<typename VectorType, int Size, int PacketAccess> class VectorBlock + : public Block<VectorType, + ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, + ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size, + PacketAccess> +{ + typedef Block<VectorType, + ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, + ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size, + PacketAccess> Base; + enum { + IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1 + }; + public: + + using Base::operator=; + using Base::operator+=; + using Base::operator-=; + using Base::operator*=; + using Base::operator/=; + + /** Dynamic-size constructor + */ + inline VectorBlock(const VectorType& vector, int start, int size) + : Base(vector, + IsColVector ? start : 0, IsColVector ? 0 : start, + IsColVector ? size : 1, IsColVector ? 1 : size) + { + + EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock); + } + + /** Fixed-size constructor + */ + inline VectorBlock(const VectorType& vector, int start) + : Base(vector, IsColVector ? start : 0, IsColVector ? 0 : start) + { + EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock); + } +}; + + +/** \returns a dynamic-size expression of a segment (i.e. a vector block) in *this. + * + * \only_for_vectors + * + * \addexample VectorBlockIntInt \label How to reference a sub-vector (dynamic size) + * + * \param start the first coefficient in the segment + * \param size the number of coefficients in the segment + * + * Example: \include MatrixBase_segment_int_int.cpp + * Output: \verbinclude MatrixBase_segment_int_int.out + * + * \note Even though the returned expression has dynamic size, in the case + * when it is applied to a fixed-size vector, it inherits a fixed maximal size, + * which means that evaluating it does not cause a dynamic memory allocation. + * + * \sa class Block, segment(int) + */ +template<typename Derived> +inline VectorBlock<Derived> MatrixBase<Derived> + ::segment(int start, int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived>(derived(), start, size); +} + +/** This is the const version of segment(int,int).*/ +template<typename Derived> +inline const VectorBlock<Derived> +MatrixBase<Derived>::segment(int start, int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived>(derived(), start, size); +} + +/** \returns a dynamic-size expression of the first coefficients of *this. + * + * \only_for_vectors + * + * \param size the number of coefficients in the block + * + * \addexample BlockInt \label How to reference a sub-vector (fixed-size) + * + * Example: \include MatrixBase_start_int.cpp + * Output: \verbinclude MatrixBase_start_int.out + * + * \note Even though the returned expression has dynamic size, in the case + * when it is applied to a fixed-size vector, it inherits a fixed maximal size, + * which means that evaluating it does not cause a dynamic memory allocation. + * + * \sa class Block, block(int,int) + */ +template<typename Derived> +inline VectorBlock<Derived> +MatrixBase<Derived>::start(int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived>(derived(), 0, size); +} + +/** This is the const version of start(int).*/ +template<typename Derived> +inline const VectorBlock<Derived> +MatrixBase<Derived>::start(int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived>(derived(), 0, size); +} + +/** \returns a dynamic-size expression of the last coefficients of *this. + * + * \only_for_vectors + * + * \param size the number of coefficients in the block + * + * \addexample BlockEnd \label How to reference the end of a vector (fixed-size) + * + * Example: \include MatrixBase_end_int.cpp + * Output: \verbinclude MatrixBase_end_int.out + * + * \note Even though the returned expression has dynamic size, in the case + * when it is applied to a fixed-size vector, it inherits a fixed maximal size, + * which means that evaluating it does not cause a dynamic memory allocation. + * + * \sa class Block, block(int,int) + */ +template<typename Derived> +inline VectorBlock<Derived> +MatrixBase<Derived>::end(int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived>(derived(), this->size() - size, size); +} + +/** This is the const version of end(int).*/ +template<typename Derived> +inline const VectorBlock<Derived> +MatrixBase<Derived>::end(int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived>(derived(), this->size() - size, size); +} + +/** \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 + * + * \param start the index of the first element of the sub-vector + * + * Example: \include MatrixBase_template_int_segment.cpp + * Output: \verbinclude MatrixBase_template_int_segment.out + * + * \sa class Block + */ +template<typename Derived> +template<int Size> +inline VectorBlock<Derived,Size> +MatrixBase<Derived>::segment(int start) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived,Size>(derived(), start); +} + +/** This is the const version of segment<int>(int).*/ +template<typename Derived> +template<int Size> +inline const VectorBlock<Derived,Size> +MatrixBase<Derived>::segment(int start) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived,Size>(derived(), start); +} + +/** \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 + * + * \addexample BlockStart \label How to reference the start of a vector (fixed-size) + * + * Example: \include MatrixBase_template_int_start.cpp + * Output: \verbinclude MatrixBase_template_int_start.out + * + * \sa class Block + */ +template<typename Derived> +template<int Size> +inline VectorBlock<Derived,Size> +MatrixBase<Derived>::start() +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived,Size>(derived(), 0); +} + +/** This is the const version of start<int>().*/ +template<typename Derived> +template<int Size> +inline const VectorBlock<Derived,Size> +MatrixBase<Derived>::start() const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived,Size>(derived(), 0); +} + +/** \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 + * + * Example: \include MatrixBase_template_int_end.cpp + * Output: \verbinclude MatrixBase_template_int_end.out + * + * \sa class Block + */ +template<typename Derived> +template<int Size> +inline VectorBlock<Derived,Size> +MatrixBase<Derived>::end() +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived, Size>(derived(), size() - Size); +} + +/** This is the const version of end<int>.*/ +template<typename Derived> +template<int Size> +inline const VectorBlock<Derived,Size> +MatrixBase<Derived>::end() const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock<Derived, Size>(derived(), size() - Size); +} + + +#endif // EIGEN_VECTORBLOCK_H diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 454b44e52..0c251022b 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -185,7 +185,7 @@ const unsigned int HereditaryBits = RowMajorBit | EvalBeforeNestingBit | EvalBeforeAssigningBit | SparseBit; - + // Possible values for the Mode parameter of part() const unsigned int UpperTriangular = UpperTriangularBit; const unsigned int StrictlyUpperTriangular = UpperTriangularBit | ZeroDiagBit; diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 3e752c84b..363972b60 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -42,6 +42,7 @@ template<typename MatrixType> class Minor; template<typename MatrixType, int BlockRows=Dynamic, int BlockCols=Dynamic, int PacketAccess=AsRequested, int _DirectAccessStatus = ei_traits<MatrixType>::Flags&DirectAccessBit ? DirectAccessBit : ei_traits<MatrixType>::Flags&SparseBit> class Block; +template<typename MatrixType, int Size=Dynamic, int PacketAccess=AsRequested> class VectorBlock; template<typename MatrixType> class Transpose; template<typename MatrixType> class Conjugate; template<typename NullaryOp, typename MatrixType> class CwiseNullaryOp; @@ -111,7 +112,7 @@ template<typename MatrixType, int Direction = BothDirections> class Reverse; template<typename MatrixType> class LU; template<typename MatrixType> class PartialLU; -template<typename MatrixType> class QR; +template<typename MatrixType> class HouseholderQR; template<typename MatrixType> class SVD; template<typename MatrixType, int UpLo = LowerTriangular> class LLT; template<typename MatrixType> class LDLT; diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index 0cc5ae9aa..c9f2f4d40 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -51,7 +51,8 @@ #define EIGEN_GCC3_OR_OLDER 0 #endif -#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER +// FIXME vectorization + alignment is completely disabled with sun studio +#if !EIGEN_GCC_AND_ARCH_DOESNT_WANT_ALIGNMENT && !EIGEN_GCC3_OR_OLDER && !defined(__SUNPRO_CC) #define EIGEN_ARCH_WANTS_ALIGNMENT 1 #else #define EIGEN_ARCH_WANTS_ALIGNMENT 0 @@ -104,7 +105,7 @@ /** 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: - * - single precision Cwise::sin() and Cwise::cos() when SSE vectorization is enabled. + * - single precision Cwise::sin() and Cwise::cos() when SSE vectorization is enabled. */ #ifndef EIGEN_FAST_MATH #define EIGEN_FAST_MATH 1 @@ -206,13 +207,16 @@ using Eigen::ei_cos; * vectorized and non-vectorized code. */ #if !EIGEN_ALIGN -#define EIGEN_ALIGN_128 + #define EIGEN_ALIGN_128 #elif (defined __GNUC__) -#define EIGEN_ALIGN_128 __attribute__((aligned(16))) + #define EIGEN_ALIGN_128 __attribute__((aligned(16))) #elif (defined _MSC_VER) -#define EIGEN_ALIGN_128 __declspec(align(16)) + #define EIGEN_ALIGN_128 __declspec(align(16)) +#elif (defined __SUNPRO_CC) + // FIXME not sure about this one: + #define EIGEN_ALIGN_128 __attribute__((aligned(16))) #else -#error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler + #error Please tell me what is the equivalent of __attribute__((aligned(16))) for your compiler #endif #define EIGEN_RESTRICT __restrict diff --git a/Eigen/src/Core/util/Memory.h b/Eigen/src/Core/util/Memory.h index cc3aa4fac..f8581eebc 100644 --- a/Eigen/src/Core/util/Memory.h +++ b/Eigen/src/Core/util/Memory.h @@ -27,7 +27,17 @@ #ifndef EIGEN_MEMORY_H #define EIGEN_MEMORY_H -#if defined(__APPLE__) || defined(_WIN64) || defined (__FreeBSD__) +// FreeBSD 6 seems to have 16-byte aligned malloc +// See http://svn.freebsd.org/viewvc/base/stable/6/lib/libc/stdlib/malloc.c?view=markup +// FreeBSD 7 seems to have 16-byte aligned malloc except on ARM and MIPS architectures +// See http://svn.freebsd.org/viewvc/base/stable/7/lib/libc/stdlib/malloc.c?view=markup +#if defined(__FreeBSD__) && !defined(__arm__) && !defined(__mips__) +#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 1 +#else +#define EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED 0 +#endif + +#if defined(__APPLE__) || defined(_WIN64) || EIGEN_FREEBSD_MALLOC_ALREADY_ALIGNED #define EIGEN_MALLOC_ALREADY_ALIGNED 1 #else #define EIGEN_MALLOC_ALREADY_ALIGNED 0 diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 7b139b0c1..a4ffb8a20 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -217,7 +217,7 @@ template<typename Derived,typename Scalar,typename OtherScalar, bool EnableIt = !ei_is_same_type<Scalar,OtherScalar>::ret > struct ei_special_scalar_op_base { - // dummy operator* so that the + // dummy operator* so that the // "using ei_special_scalar_op_base::operator*" compiles void operator*() const; }; @@ -237,8 +237,6 @@ struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true> * TODO: could be a good idea to define a big ReturnType struct ?? */ template<typename ExpressionType, int RowsOrSize=Dynamic, int Cols=Dynamic> struct BlockReturnType { - typedef Block<ExpressionType, (ei_traits<ExpressionType>::RowsAtCompileTime == 1 ? 1 : RowsOrSize), - (ei_traits<ExpressionType>::ColsAtCompileTime == 1 ? 1 : RowsOrSize)> SubVectorType; typedef Block<ExpressionType, RowsOrSize, Cols> Type; }; @@ -251,7 +249,7 @@ template<typename ExpressionType> struct HNormalizedReturnType { typedef Block<ExpressionType, ei_traits<ExpressionType>::ColsAtCompileTime==1 ? SizeMinusOne : 1, ei_traits<ExpressionType>::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne; - typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<ExpressionType>::Scalar>, + typedef CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<ExpressionType>::Scalar>, NestByValue<StartMinusOne> > Type; }; |