diff options
Diffstat (limited to 'Eigen/src')
25 files changed, 1095 insertions, 849 deletions
diff --git a/Eigen/src/Array/Replicate.h b/Eigen/src/Array/Replicate.h index df3afbbdb..02f9c0601 100644 --- a/Eigen/src/Array/Replicate.h +++ b/Eigen/src/Array/Replicate.h @@ -140,21 +140,4 @@ VectorwiseOp<ExpressionType,Direction>::replicate(int factor) const (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1); } -/** \nonstableyet - * \return an expression of the replication of each column (or row) of \c *this - * - * Example: \include DirectionWise_replicate.cpp - * Output: \verbinclude DirectionWise_replicate.out - * - * \sa VectorwiseOp::replicate(int), MatrixBase::replicate(), class Replicate - */ -template<typename ExpressionType, int Direction> -template<int Factor> -const Replicate<ExpressionType,(Direction==Vertical?Factor:1),(Direction==Horizontal?Factor:1)> -VectorwiseOp<ExpressionType,Direction>::replicate(int factor) const -{ - return Replicate<ExpressionType,Direction==Vertical?Factor:1,Direction==Horizontal?Factor:1> - (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1); -} - #endif // EIGEN_REPLICATE_H diff --git a/Eigen/src/Array/VectorwiseOp.h b/Eigen/src/Array/VectorwiseOp.h index 50302bba4..3684d7cd1 100644 --- a/Eigen/src/Array/VectorwiseOp.h +++ b/Eigen/src/Array/VectorwiseOp.h @@ -179,6 +179,11 @@ template<typename ExpressionType, int Direction> class VectorwiseOp > Type; }; + enum { + IsVertical = (Direction==Vertical) ? 1 : 0, + IsHorizontal = (Direction==Horizontal) ? 1 : 0 + }; + protected: /** \internal @@ -222,9 +227,17 @@ template<typename ExpressionType, int Direction> class VectorwiseOp /** \internal */ inline const ExpressionType& _expression() const { return m_matrix; } + /** \returns a row or column vector expression of \c *this reduxed by \a func + * + * The template parameter \a BinaryOp is the type of the functor + * of the custom redux operator. Note that func must be an associative operator. + * + * \sa class VectorwiseOp, MatrixBase::colwise(), MatrixBase::rowwise() + */ template<typename BinaryOp> const typename ReduxReturnType<BinaryOp>::Type - redux(const BinaryOp& func = BinaryOp()) const; + redux(const BinaryOp& func = BinaryOp()) const + { return typename ReduxReturnType<BinaryOp>::Type(_expression(), func); } /** \returns a row (or column) vector expression of the smallest coefficient * of each column (or row) of the referenced expression. @@ -319,16 +332,26 @@ template<typename ExpressionType, int Direction> class VectorwiseOp * * \sa MatrixBase::reverse() */ const Reverse<ExpressionType, Direction> reverse() const - { - return Reverse<ExpressionType, Direction>( _expression() ); - } + { return Reverse<ExpressionType, Direction>( _expression() ); } const Replicate<ExpressionType,Direction==Vertical?Dynamic:1,Direction==Horizontal?Dynamic:1> replicate(int factor) const; - template<int Factor> - const Replicate<ExpressionType,(Direction==Vertical?Factor:1),(Direction==Horizontal?Factor:1)> - replicate(int factor = Factor) const; + /** \nonstableyet + * \return an expression of the replication of each column (or row) of \c *this + * + * Example: \include DirectionWise_replicate.cpp + * Output: \verbinclude DirectionWise_replicate.out + * + * \sa VectorwiseOp::replicate(int), MatrixBase::replicate(), class Replicate + */ + // NOTE implemented here because of sunstudio's compilation errors + template<int Factor> const Replicate<ExpressionType,(IsVertical?Factor:1),(IsHorizontal?Factor:1)> + replicate(int factor = Factor) const + { + return Replicate<ExpressionType,Direction==Vertical?Factor:1,Direction==Horizontal?Factor:1> + (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1); + } /////////// Artithmetic operators /////////// @@ -466,19 +489,4 @@ MatrixBase<Derived>::rowwise() return derived(); } -/** \returns a row or column vector expression of \c *this reduxed by \a func - * - * The template parameter \a BinaryOp is the type of the functor - * of the custom redux operator. Note that func must be an associative operator. - * - * \sa class VectorwiseOp, MatrixBase::colwise(), MatrixBase::rowwise() - */ -template<typename ExpressionType, int Direction> -template<typename BinaryOp> -const typename VectorwiseOp<ExpressionType,Direction>::template ReduxReturnType<BinaryOp>::Type -VectorwiseOp<ExpressionType,Direction>::redux(const BinaryOp& func) const -{ - return typename ReduxReturnType<BinaryOp>::Type(_expression(), func); -} - #endif // EIGEN_PARTIAL_REDUX_H 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; }; diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index 66580cf1b..9385f259d 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -303,7 +303,9 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::operator=(const MatrixBase<Derive return *this; } -/** Convert the quaternion to a 3x3 rotation matrix */ +/** Convert the quaternion to a 3x3 rotation matrix. The quaternion is required to + * be normalized, otherwise the result is undefined. + */ template<typename Scalar> inline typename Quaternion<Scalar>::Matrix3 Quaternion<Scalar>::toRotationMatrix(void) const @@ -340,11 +342,15 @@ Quaternion<Scalar>::toRotationMatrix(void) const return res; } -/** Sets *this to be a quaternion representing a rotation sending the vector \a a to the vector \a b. +/** Sets \c *this to be a quaternion representing a rotation between + * the two arbitrary vectors \a a and \a b. In other words, the built + * rotation represent a rotation sending the line of direction \a a + * to the line of direction \a b, both lines passing through the origin. * - * \returns a reference to *this. + * \returns a reference to \c *this. * - * Note that the two input vectors do \b not have to be normalized. + * Note that the two input vectors do \b not have to be normalized, and + * do not need to have the same norm. */ template<typename Scalar> template<typename Derived1, typename Derived2> @@ -354,21 +360,26 @@ inline Quaternion<Scalar>& Quaternion<Scalar>::setFromTwoVectors(const MatrixBas Vector3 v1 = b.normalized(); Scalar c = v0.dot(v1); - // if dot == 1, vectors are the same - if (ei_isApprox(c,Scalar(1))) + // if dot == -1, vectors are nearly opposites + // => accuraletly 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 + // under the constraint: + // ||x|| = 1 + // which yields a singular value problem + if (c < Scalar(-1)+precision<Scalar>()) { - // set to identity - this->w() = 1; this->vec().setZero(); + c = std::max<Scalar>(c,-1); + Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose(); + SVD<Matrix<Scalar,2,3> > svd(m); + Vector3 axis = svd.matrixV().col(2); + + Scalar w2 = (Scalar(1)+c)*Scalar(0.5); + this->w() = ei_sqrt(w2); + this->vec() = axis * ei_sqrt(Scalar(1) - w2); return *this; } - // if dot == -1, vectors are opposites - if (ei_isApprox(c,Scalar(-1))) - { - this->vec() = v0.unitOrthogonal(); - this->w() = 0; - return *this; - } - Vector3 axis = v0.cross(v1); Scalar s = ei_sqrt((Scalar(1)+c)*Scalar(2)); Scalar invs = Scalar(1)/s; diff --git a/Eigen/src/LU/Inverse.h b/Eigen/src/LU/Inverse.h index b29efb60c..5af14813d 100644 --- a/Eigen/src/LU/Inverse.h +++ b/Eigen/src/LU/Inverse.h @@ -29,41 +29,45 @@ *** Part 1 : optimized implementations for fixed-size 2,3,4 cases *** ********************************************************************/ -template<typename MatrixType> -void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result) +template<typename XprType, typename MatrixType> +inline void ei_compute_inverse_size2_helper( + const XprType& matrix, const typename MatrixType::Scalar& invdet, + MatrixType* result) { - typedef typename MatrixType::Scalar Scalar; - const Scalar invdet = Scalar(1) / matrix.determinant(); result->coeffRef(0,0) = matrix.coeff(1,1) * invdet; result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet; result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet; result->coeffRef(1,1) = matrix.coeff(0,0) * invdet; } +template<typename MatrixType> +inline void ei_compute_inverse_size2(const MatrixType& matrix, MatrixType* result) +{ + typedef typename MatrixType::Scalar Scalar; + const Scalar invdet = Scalar(1) / matrix.determinant(); + ei_compute_inverse_size2_helper( matrix, invdet, result ); +} + template<typename XprType, typename MatrixType> -bool ei_compute_inverse_in_size2_case_with_check(const XprType& matrix, MatrixType* result) +bool ei_compute_inverse_size2_with_check(const XprType& matrix, MatrixType* result) { typedef typename MatrixType::Scalar Scalar; const Scalar det = matrix.determinant(); if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false; const Scalar invdet = Scalar(1) / det; - result->coeffRef(0,0) = matrix.coeff(1,1) * invdet; - result->coeffRef(1,0) = -matrix.coeff(1,0) * invdet; - result->coeffRef(0,1) = -matrix.coeff(0,1) * invdet; - result->coeffRef(1,1) = matrix.coeff(0,0) * invdet; + ei_compute_inverse_size2_helper( matrix, invdet, result ); return true; } -template<typename MatrixType> -void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result) +template<typename XprType, typename MatrixType> +void ei_compute_inverse_size3_helper( + const XprType& matrix, + const typename MatrixType::Scalar& invdet, + const typename MatrixType::Scalar& det_minor00, + const typename MatrixType::Scalar& det_minor10, + const typename MatrixType::Scalar& det_minor20, + MatrixType* result) { - typedef typename MatrixType::Scalar Scalar; - const Scalar det_minor00 = matrix.minor(0,0).determinant(); - const Scalar det_minor10 = matrix.minor(1,0).determinant(); - const Scalar det_minor20 = matrix.minor(2,0).determinant(); - const Scalar invdet = Scalar(1) / ( det_minor00 * matrix.coeff(0,0) - - det_minor10 * matrix.coeff(1,0) - + det_minor20 * matrix.coeff(2,0) ); result->coeffRef(0, 0) = det_minor00 * invdet; result->coeffRef(0, 1) = -det_minor10 * invdet; result->coeffRef(0, 2) = det_minor20 * invdet; @@ -75,8 +79,24 @@ void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* resu result->coeffRef(2, 2) = matrix.minor(2,2).determinant() * invdet; } +template<bool Check, typename XprType, typename MatrixType> +bool ei_compute_inverse_size3(const XprType& matrix, MatrixType* result) +{ + typedef typename MatrixType::Scalar Scalar; + const Scalar det_minor00 = matrix.minor(0,0).determinant(); + const Scalar det_minor10 = matrix.minor(1,0).determinant(); + const Scalar det_minor20 = matrix.minor(2,0).determinant(); + const Scalar det = ( det_minor00 * matrix.coeff(0,0) + - det_minor10 * matrix.coeff(1,0) + + det_minor20 * matrix.coeff(2,0) ); + if(Check) if(ei_isMuchSmallerThan(det, matrix.cwise().abs().maxCoeff())) return false; + const Scalar invdet = Scalar(1) / det; + ei_compute_inverse_size3_helper( matrix, invdet, det_minor00, det_minor10, det_minor20, result ); + return true; +} + template<typename MatrixType> -bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixType* result) +bool ei_compute_inverse_size4_helper(const MatrixType& matrix, MatrixType* result) { /* Let's split M into four 2x2 blocks: * (P Q) @@ -94,7 +114,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp typedef Block<MatrixType,2,2> XprBlock22; typedef typename MatrixBase<XprBlock22>::PlainMatrixType Block22; Block22 P_inverse; - if(ei_compute_inverse_in_size2_case_with_check(matrix.template block<2,2>(0,0), &P_inverse)) + if(ei_compute_inverse_size2_with_check(matrix.template block<2,2>(0,0), &P_inverse)) { const Block22 Q = matrix.template block<2,2>(0,2); const Block22 P_inverse_times_Q = P_inverse * Q; @@ -104,7 +124,7 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp const XprBlock22 S = matrix.template block<2,2>(2,2); const Block22 X = S - R_times_P_inverse_times_Q; Block22 Y; - ei_compute_inverse_in_size2_case(X, &Y); + ei_compute_inverse_size2(X, &Y); result->template block<2,2>(2,2) = Y; result->template block<2,2>(2,0) = - Y * R_times_P_inverse; const Block22 Z = P_inverse_times_Q * Y; @@ -118,13 +138,13 @@ bool ei_compute_inverse_in_size4_case_helper(const MatrixType& matrix, MatrixTyp } } -template<typename MatrixType> -void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result) +template<typename XprType, typename MatrixType> +bool ei_compute_inverse_size4_with_check(const XprType& matrix, MatrixType* result) { - if(ei_compute_inverse_in_size4_case_helper(matrix, result)) + if(ei_compute_inverse_size4_helper(matrix, result)) { // good ! The topleft 2x2 block was invertible, so the 2x2 blocks approach is successful. - return; + return true; } else { @@ -134,17 +154,17 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu MatrixType m(matrix); m.row(0).swap(m.row(2)); m.row(1).swap(m.row(3)); - if(ei_compute_inverse_in_size4_case_helper(m, result)) + if(ei_compute_inverse_size4_helper(m, result)) { // good, the topleft 2x2 block of m is invertible. Since m is different from matrix in that some // rows were permuted, the actual inverse of matrix is derived from the inverse of m by permuting // the corresponding columns. result->col(0).swap(result->col(2)); result->col(1).swap(result->col(3)); + return true; } else { - // last possible case. Since matrix is assumed to be invertible, this last case has to work. // first, undo the swaps previously made m.row(0).swap(m.row(2)); m.row(1).swap(m.row(3)); @@ -154,13 +174,23 @@ void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* resu // swap row 1 with the the row among 2 and 3 that has the biggest 2 first coeffs int swap1with = ei_abs(m.coeff(2,0))+ei_abs(m.coeff(2,1))>ei_abs(m.coeff(3,0))+ei_abs(m.coeff(3,1)) ? 2 : 3; m.row(1).swap(m.row(swap1with)); - ei_compute_inverse_in_size4_case_helper(m, result); - result->col(1).swap(result->col(swap1with)); - result->col(0).swap(result->col(swap0with)); + if( ei_compute_inverse_size4_helper(m, result) ) + { + result->col(1).swap(result->col(swap1with)); + result->col(0).swap(result->col(swap0with)); + return true; + } + else + { + // non-invertible matrix + return false; + } } } } + + /*********************************************** *** Part 2 : selector and MatrixBase methods *** ***********************************************/ @@ -189,7 +219,7 @@ struct ei_compute_inverse<MatrixType, 2> { static inline void run(const MatrixType& matrix, MatrixType* result) { - ei_compute_inverse_in_size2_case(matrix, result); + ei_compute_inverse_size2(matrix, result); } }; @@ -198,7 +228,7 @@ struct ei_compute_inverse<MatrixType, 3> { static inline void run(const MatrixType& matrix, MatrixType* result) { - ei_compute_inverse_in_size3_case(matrix, result); + ei_compute_inverse_size3<false, MatrixType, MatrixType>(matrix, result); } }; @@ -207,7 +237,7 @@ struct ei_compute_inverse<MatrixType, 4> { static inline void run(const MatrixType& matrix, MatrixType* result) { - ei_compute_inverse_in_size4_case(matrix, result); + ei_compute_inverse_size4_with_check(matrix, result); } }; @@ -215,14 +245,15 @@ struct ei_compute_inverse<MatrixType, 4> * * Computes the matrix inverse of this matrix. * - * \note This matrix must be invertible, otherwise the result is undefined. + * \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use + * computeInverseWithCheck(). * * \param result Pointer to the matrix in which to store the result. * * Example: \include MatrixBase_computeInverse.cpp * Output: \verbinclude MatrixBase_computeInverse.out * - * \sa inverse() + * \sa inverse(), computeInverseWithCheck() */ template<typename Derived> inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const @@ -236,7 +267,8 @@ inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const * * \returns the matrix inverse of this matrix. * - * \note This matrix must be invertible, otherwise the result is undefined. + * \note This matrix must be invertible, otherwise the result is undefined. If you need an invertibility check, use + * computeInverseWithCheck(). * * \note This method returns a matrix by value, which can be inefficient. To avoid that overhead, * use computeInverse() instead. @@ -244,7 +276,7 @@ inline void MatrixBase<Derived>::computeInverse(PlainMatrixType *result) const * Example: \include MatrixBase_inverse.cpp * Output: \verbinclude MatrixBase_inverse.out * - * \sa computeInverse() + * \sa computeInverse(), computeInverseWithCheck() */ template<typename Derived> inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>::inverse() const @@ -254,4 +286,81 @@ inline const typename MatrixBase<Derived>::PlainMatrixType MatrixBase<Derived>:: return result; } + +/******************************************** + * Compute inverse with invertibility check * + *******************************************/ + +template<typename MatrixType, int Size = MatrixType::RowsAtCompileTime> +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + typedef typename MatrixType::Scalar Scalar; + LU<MatrixType> lu( matrix ); + if( !lu.isInvertible() ) return false; + lu.computeInverse(result); + return true; + } +}; + +template<typename MatrixType> +struct ei_compute_inverse_with_check<MatrixType, 1> +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + if( 0 == result->coeffRef(0,0) ) return false; + + typedef typename MatrixType::Scalar Scalar; + result->coeffRef(0,0) = Scalar(1) / matrix.coeff(0,0); + return true; + } +}; + +template<typename MatrixType> +struct ei_compute_inverse_with_check<MatrixType, 2> +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + return ei_compute_inverse_size2_with_check(matrix, result); + } +}; + +template<typename MatrixType> +struct ei_compute_inverse_with_check<MatrixType, 3> +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + return ei_compute_inverse_size3<true, MatrixType, MatrixType>(matrix, result); + } +}; + +template<typename MatrixType> +struct ei_compute_inverse_with_check<MatrixType, 4> +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + return ei_compute_inverse_size4_with_check(matrix, result); + } +}; + +/** \lu_module + * + * Computation of matrix inverse, with invertibility check. + * + * \returns true if the matrix is invertible, false otherwise. + * + * \param result Pointer to the matrix in which to store the result. + * + * \sa inverse(), computeInverse() + */ +template<typename Derived> +inline bool MatrixBase<Derived>::computeInverseWithCheck(PlainMatrixType *result) const +{ + ei_assert(rows() == cols()); + EIGEN_STATIC_ASSERT(NumTraits<Scalar>::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT) + return ei_compute_inverse_with_check<PlainMatrixType>::run(eval(), result); +} + + #endif // EIGEN_INVERSE_H diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 5883eed65..a23a3a035 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -28,18 +28,20 @@ /** \ingroup QR_Module * \nonstableyet * - * \class QR + * \class HouseholderQR * - * \brief QR decomposition of a matrix + * \brief Householder QR decomposition of a matrix * * \param MatrixType the type of the matrix of which we are computing the QR decomposition * * This class performs a QR decomposition using Householder transformations. The result is * stored in a compact way compatible with LAPACK. * + * Note that no pivoting is performed. This is \b not a rank-revealing decomposition. + * * \sa MatrixBase::qr() */ -template<typename MatrixType> class QR +template<typename MatrixType> class HouseholderQR { public: @@ -53,88 +55,23 @@ template<typename MatrixType> class QR * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to - * perform decompositions via QR::compute(const MatrixType&). + * perform decompositions via HouseholderQR::compute(const MatrixType&). */ - QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + HouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} - QR(const MatrixType& matrix) + HouseholderQR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), m_hCoeffs(matrix.cols()), m_isInitialized(false) { compute(matrix); } - - /** \deprecated use isInjective() - * \returns whether or not the matrix is of full rank - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - EIGEN_DEPRECATED bool isFullRank() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.cols(); - } - - /** \returns the rank of the matrix of which *this is the QR decomposition. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - int rank() const; - - /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline int dimensionOfKernel() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return m_qr.cols() - rank(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents an injective - * linear map, i.e. has trivial kernel; false otherwise. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isInjective() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.cols(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents a surjective - * linear map; false otherwise. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isSurjective() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return rank() == m_qr.rows(); - } - - /** \returns true if the matrix of which *this is the QR decomposition is invertible. - * - * \note Since the rank is computed only once, i.e. the first time it is needed, this - * method almost does not perform any further computation. - */ - inline bool isInvertible() const - { - ei_assert(m_isInitialized && "QR is not initialized."); - return isInjective() && isSurjective(); - } - + /** \returns a read-only expression of the matrix R of the actual the QR decomposition */ const TriangularView<NestByValue<MatrixRBlockType>, UpperTriangular> matrixR(void) const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); int cols = m_qr.cols(); return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<UpperTriangular>(); } @@ -148,58 +85,35 @@ template<typename MatrixType> class QR * Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols(). * If no solution exists, *result is left with undefined coefficients. * - * \returns true if any solution exists, false if no solution exists. - * - * \note If there exist more than one solution, this method will arbitrarily choose one. - * If you need a complete analysis of the space of solutions, take the one solution obtained - * by this method and add to it elements of the kernel, as determined by kernel(). - * * \note The case where b is a matrix is not yet implemented. Also, this * code is space inefficient. * - * Example: \include QR_solve.cpp - * Output: \verbinclude QR_solve.out - * - * \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse() + * Example: \include HouseholderQR_solve.cpp + * Output: \verbinclude HouseholderQR_solve.out */ template<typename OtherDerived, typename ResultType> - bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; + void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; MatrixType matrixQ(void) const; + + /** \returns a reference to the matrix where the Householder QR decomposition is stored + * in a LAPACK-compatible way. + */ + const MatrixType& matrixQR() const { return m_qr; } void compute(const MatrixType& matrix); protected: MatrixType m_qr; VectorType m_hCoeffs; - mutable int m_rank; - mutable bool m_rankIsUptodate; bool m_isInitialized; }; -/** \returns the rank of the matrix of which *this is the QR decomposition. */ -template<typename MatrixType> -int QR<MatrixType>::rank() const -{ - ei_assert(m_isInitialized && "QR is not initialized."); - if (!m_rankIsUptodate) - { - RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff(); - int n = m_qr.cols(); - m_rank = 0; - while(m_rank<n && !ei_isMuchSmallerThan(m_qr.diagonal().coeff(m_rank), maxCoeff)) - ++m_rank; - m_rankIsUptodate = true; - } - return m_rank; -} - #ifndef EIGEN_HIDE_HEAVY_CODE template<typename MatrixType> -void QR<MatrixType>::compute(const MatrixType& matrix) +void HouseholderQR<MatrixType>::compute(const MatrixType& matrix) { - m_rankIsUptodate = false; m_qr = matrix; m_hCoeffs.resize(matrix.cols()); @@ -262,12 +176,12 @@ void QR<MatrixType>::compute(const MatrixType& matrix) template<typename MatrixType> template<typename OtherDerived, typename ResultType> -bool QR<MatrixType>::solve( +void HouseholderQR<MatrixType>::solve( const MatrixBase<OtherDerived>& b, ResultType *result ) const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); const int rows = m_qr.rows(); ei_assert(b.rows() == rows); result->resize(rows, b.cols()); @@ -276,27 +190,17 @@ bool QR<MatrixType>::solve( // Q^T without explicitly forming matrixQ(). Investigate. *result = matrixQ().transpose()*b; - if(!isSurjective()) - { - // is result is in the image of R ? - RealScalar biggest_in_res = result->corner(TopLeft, m_rank, result->cols()).cwise().abs().maxCoeff(); - for(int col = 0; col < result->cols(); ++col) - for(int row = m_rank; row < result->rows(); ++row) - if(!ei_isMuchSmallerThan(result->coeff(row,col), biggest_in_res)) - return false; - } - m_qr.corner(TopLeft, m_rank, m_rank) + const int rank = std::min(result->rows(), result->cols()); + m_qr.corner(TopLeft, rank, rank) .template marked<UpperTriangular>() - .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols())); - - return true; + .solveTriangularInPlace(result->corner(TopLeft, rank, result->cols())); } /** \returns the matrix Q */ template<typename MatrixType> -MatrixType QR<MatrixType>::matrixQ() const +MatrixType HouseholderQR<MatrixType>::matrixQ() const { - ei_assert(m_isInitialized && "QR is not initialized."); + ei_assert(m_isInitialized && "HouseholderQR is not initialized."); // compute the product Q_0 Q_1 ... Q_n-1, // where Q_k is the k-th Householder transformation I - h_k v_k v_k' // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] @@ -319,15 +223,15 @@ MatrixType QR<MatrixType>::matrixQ() const #endif // EIGEN_HIDE_HEAVY_CODE -/** \return the QR decomposition of \c *this. +/** \return the Householder QR decomposition of \c *this. * - * \sa class QR + * \sa class HouseholderQR */ template<typename Derived> -const QR<typename MatrixBase<Derived>::PlainMatrixType> -MatrixBase<Derived>::qr() const +const HouseholderQR<typename MatrixBase<Derived>::PlainMatrixType> +MatrixBase<Derived>::householderQr() const { - return QR<PlainMatrixType>(eval()); + return HouseholderQR<PlainMatrixType>(eval()); } diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 9733fbd21..f9f9feb89 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr> // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -34,9 +34,7 @@ * * \param MatrixType the type of the matrix of which we are computing the SVD decomposition * - * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N - * with \c M \>= \c N. - * + * This class performs a standard SVD decomposition of a real matrix A of size \c M x \c N. * * \sa MatrixBase::SVD() */ @@ -55,13 +53,13 @@ template<typename MatrixType> class SVD typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> ColVector; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> RowVector; - typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MinSize> MatrixUType; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> MatrixUType; typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixVType; - typedef Matrix<Scalar, MinSize, 1> SingularValuesType; + typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> SingularValuesType; public: - /** + /** * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to @@ -70,9 +68,9 @@ template<typename MatrixType> class SVD SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {} SVD(const MatrixType& matrix) - : m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())), + : m_matU(matrix.rows(), matrix.rows()), m_matV(matrix.cols(),matrix.cols()), - m_sigma(std::min(matrix.rows(),matrix.cols())), + m_sigma(matrix.cols()), m_isInitialized(false) { compute(matrix); @@ -81,22 +79,22 @@ template<typename MatrixType> class SVD template<typename OtherDerived, typename ResultType> bool solve(const MatrixBase<OtherDerived> &b, ResultType* result) const; - const MatrixUType& matrixU() const - { + const MatrixUType& matrixU() const + { ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matU; + return m_matU; } - const SingularValuesType& singularValues() const + const SingularValuesType& singularValues() const { ei_assert(m_isInitialized && "SVD is not initialized."); - return m_sigma; + return m_sigma; } - const MatrixVType& matrixV() const + const MatrixVType& matrixV() const { ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matV; + return m_matV; } void compute(const MatrixType& matrix); @@ -112,6 +110,23 @@ template<typename MatrixType> class SVD void computeScalingRotation(ScalingType *positive, RotationType *unitary) const; protected: + // Computes (a^2 + b^2)^(1/2) without destructive underflow or overflow. + inline static Scalar pythag(Scalar a, Scalar b) + { + Scalar abs_a = ei_abs(a); + Scalar abs_b = ei_abs(b); + if (abs_a > abs_b) + return abs_a*ei_sqrt(Scalar(1.0)+ei_abs2(abs_b/abs_a)); + else + return (abs_b == Scalar(0.0) ? Scalar(0.0) : abs_b*ei_sqrt(Scalar(1.0)+ei_abs2(abs_a/abs_b))); + } + + inline static Scalar sign(Scalar a, Scalar b) + { + return (b >= Scalar(0.0) ? ei_abs(a) : -ei_abs(a)); + } + + protected: /** \internal */ MatrixUType m_matU; /** \internal */ @@ -123,380 +138,271 @@ template<typename MatrixType> class SVD /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix * - * \note this code has been adapted from JAMA (public domain) + * \note this code has been adapted from Numerical Recipes, third edition. */ template<typename MatrixType> void SVD<MatrixType>::compute(const MatrixType& matrix) { const int m = matrix.rows(); const int n = matrix.cols(); - const int nu = std::min(m,n); - m_matU.resize(m, nu); + m_matU.resize(m, m); m_matU.setZero(); - m_sigma.resize(std::min(m,n)); + m_sigma.resize(n); m_matV.resize(n,n); - RowVector e(n); - ColVector work(m); - MatrixType matA(matrix); - const bool wantu = true; - const bool wantv = true; - int i=0, j=0, k=0; - - // Reduce A to bidiagonal form, storing the diagonal elements - // in s and the super-diagonal elements in e. - int nct = std::min(m-1,n); - int nrt = std::max(0,std::min(n-2,m)); - for (k = 0; k < std::max(nct,nrt); ++k) - { - if (k < nct) - { - // Compute the transformation for the k-th column and - // place the k-th diagonal in m_sigma[k]. - m_sigma[k] = matA.col(k).end(m-k).norm(); - if (m_sigma[k] != 0.0) // FIXME - { - if (matA(k,k) < 0.0) - m_sigma[k] = -m_sigma[k]; - matA.col(k).end(m-k) /= m_sigma[k]; - matA(k,k) += 1.0; - } - m_sigma[k] = -m_sigma[k]; - } - - for (j = k+1; j < n; ++j) - { - if ((k < nct) && (m_sigma[k] != 0.0)) - { - // Apply the transformation. - Scalar t = matA.col(k).end(m-k).dot(matA.col(j).end(m-k)); // FIXME dot product or cwise prod + .sum() ?? - t = -t/matA(k,k); - matA.col(j).end(m-k) += t * matA.col(k).end(m-k); - } + int max_iters = 30; - // Place the k-th row of A into e for the - // subsequent calculation of the row transformation. - e[j] = matA(k,j); - } + MatrixVType& V = m_matV; + MatrixType A = matrix; + SingularValuesType& W = m_sigma; - // Place the transformation in U for subsequent back multiplication. - if (wantu & (k < nct)) - m_matU.col(k).end(m-k) = matA.col(k).end(m-k); + bool flag; + int i,its,j,jj,k,l,nm; + Scalar anorm, c, f, g, h, s, scale, x, y, z; + bool convergence = true; + Scalar eps = precision<Scalar>(); - if (k < nrt) + Matrix<Scalar,Dynamic,1> rv1(n); + g = scale = anorm = 0; + // Householder reduction to bidiagonal form. + for (i=0; i<n; i++) + { + l = i+2; + rv1[i] = scale*g; + g = s = scale = 0.0; + if (i < m) { - // Compute the k-th row transformation and place the - // k-th super-diagonal in e[k]. - e[k] = e.end(n-k-1).norm(); - if (e[k] != 0.0) - { - if (e[k+1] < 0.0) - e[k] = -e[k]; - e.end(n-k-1) /= e[k]; - e[k+1] += 1.0; - } - e[k] = -e[k]; - if ((k+1 < m) & (e[k] != 0.0)) + scale = A.col(i).end(m-i).cwise().abs().sum(); + if (scale != Scalar(0)) { - // Apply the transformation. - work.end(m-k-1) = matA.corner(BottomRight,m-k-1,n-k-1) * e.end(n-k-1); - for (j = k+1; j < n; ++j) - matA.col(j).end(m-k-1) += (-e[j]/e[k+1]) * work.end(m-k-1); + for (k=i; k<m; k++) + { + A(k, i) /= scale; + s += A(k, i)*A(k, i); + } + f = A(i, i); + g = -sign( ei_sqrt(s), f ); + h = f*g - s; + A(i, i)=f-g; + for (j=l-1; j<n; j++) + { + s = A.col(i).end(m-i).dot(A.col(j).end(m-i)); + f = s/h; + A.col(j).end(m-i) += f*A.col(i).end(m-i); + } + A.col(i).end(m-i) *= scale; } - - // Place the transformation in V for subsequent back multiplication. - if (wantv) - m_matV.col(k).end(n-k-1) = e.end(n-k-1); - } - } - - - // Set up the final bidiagonal matrix or order p. - int p = std::min(n,m+1); - if (nct < n) - m_sigma[nct] = matA(nct,nct); - if (m < p) - m_sigma[p-1] = 0.0; - if (nrt+1 < p) - e[nrt] = matA(nrt,p-1); - e[p-1] = 0.0; - - // If required, generate U. - if (wantu) - { - for (j = nct; j < nu; ++j) - { - m_matU.col(j).setZero(); - m_matU(j,j) = 1.0; } - for (k = nct-1; k >= 0; k--) + W[i] = scale * g; + g = s = scale = 0.0; + if (i+1 <= m && i+1 != n) { - if (m_sigma[k] != 0.0) + scale = A.row(i).end(n-l+1).cwise().abs().sum(); + if (scale != Scalar(0)) { - for (j = k+1; j < nu; ++j) + for (k=l-1; k<n; k++) { - Scalar t = m_matU.col(k).end(m-k).dot(m_matU.col(j).end(m-k)); // FIXME is it really a dot product we want ? - t = -t/m_matU(k,k); - m_matU.col(j).end(m-k) += t * m_matU.col(k).end(m-k); + A(i, k) /= scale; + s += A(i, k)*A(i, k); } - m_matU.col(k).end(m-k) = - m_matU.col(k).end(m-k); - m_matU(k,k) = Scalar(1) + m_matU(k,k); - if (k-1>0) - m_matU.col(k).start(k-1).setZero(); - } - else - { - m_matU.col(k).setZero(); - m_matU(k,k) = 1.0; + f = A(i,l-1); + g = -sign(ei_sqrt(s),f); + h = f*g - s; + A(i,l-1) = f-g; + rv1.end(n-l+1) = A.row(i).end(n-l+1)/h; + for (j=l-1; j<m; j++) + { + s = A.row(j).end(n-l+1).dot(A.row(i).end(n-l+1)); + A.row(j).end(n-l+1) += s*rv1.end(n-l+1).transpose(); + } + A.row(i).end(n-l+1) *= scale; } } + anorm = std::max( anorm, (ei_abs(W[i])+ei_abs(rv1[i])) ); } - - // If required, generate V. - if (wantv) + // Accumulation of right-hand transformations. + for (i=n-1; i>=0; i--) { - for (k = n-1; k >= 0; k--) + //Accumulation of right-hand transformations. + if (i < n-1) { - if ((k < nrt) & (e[k] != 0.0)) + if (g != Scalar(0.0)) { - for (j = k+1; j < nu; ++j) + for (j=l; j<n; j++) //Double division to avoid possible underflow. + V(j, i) = (A(i, j)/A(i, l))/g; + for (j=l; j<n; j++) { - Scalar t = m_matV.col(k).end(n-k-1).dot(m_matV.col(j).end(n-k-1)); // FIXME is it really a dot product we want ? - t = -t/m_matV(k+1,k); - m_matV.col(j).end(n-k-1) += t * m_matV.col(k).end(n-k-1); + s = A.row(i).end(n-l).dot(V.col(j).end(n-l)); + V.col(j).end(n-l) += s * V.col(i).end(n-l); } } - m_matV.col(k).setZero(); - m_matV(k,k) = 1.0; + V.row(i).end(n-l).setZero(); + V.col(i).end(n-l).setZero(); } + V(i, i) = 1.0; + g = rv1[i]; + l = i; } - - // Main iteration loop for the singular values. - int pp = p-1; - int iter = 0; - Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52)); - while (p > 0) + // Accumulation of left-hand transformations. + for (i=std::min(m,n)-1; i>=0; i--) { - int k=0; - int kase=0; - - // Here is where a test for too many iterations would go. - - // This section of the program inspects for - // negligible elements in the s and e arrays. On - // completion the variables kase and k are set as follows. - - // kase = 1 if s(p) and e[k-1] are negligible and k<p - // kase = 2 if s(k) is negligible and k<p - // kase = 3 if e[k-1] is negligible, k<p, and - // s(k), ..., s(p) are not negligible (qr step). - // kase = 4 if e(p-1) is negligible (convergence). - - for (k = p-2; k >= -1; --k) + l = i+1; + g = W[i]; + if (n-l>0) + A.row(i).end(n-l).setZero(); + if (g != Scalar(0.0)) { - if (k == -1) - break; - if (ei_abs(e[k]) <= eps*(ei_abs(m_sigma[k]) + ei_abs(m_sigma[k+1]))) + g = Scalar(1.0)/g; + for (j=l; j<n; j++) { - e[k] = 0.0; - break; + s = A.col(i).end(m-l).dot(A.col(j).end(m-l)); + f = (s/A(i,i))*g; + A.col(j).end(m-i) += f * A.col(i).end(m-i); } - } - if (k == p-2) - { - kase = 4; + A.col(i).end(m-i) *= g; } else + A.col(i).end(m-i).setZero(); + ++A(i,i); + } + // Diagonalization of the bidiagonal form: Loop over + // singular values, and over allowed iterations. + for (k=n-1; k>=0; k--) + { + for (its=0; its<max_iters; its++) { - int ks; - for (ks = p-1; ks >= k; --ks) + flag = true; + for (l=k; l>=0; l--) { - if (ks == k) - break; - Scalar t = (ks != p ? ei_abs(e[ks]) : Scalar(0)) + (ks != k+1 ? ei_abs(e[ks-1]) : Scalar(0)); - if (ei_abs(m_sigma[ks]) <= eps*t) + // Test for splitting. + nm = l-1; + // Note that rv1[1] is always zero. + //if ((double)(ei_abs(rv1[l])+anorm) == anorm) + if (l==0 || ei_abs(rv1[l]) <= eps*anorm) { - m_sigma[ks] = 0.0; + flag = false; break; } + //if ((double)(ei_abs(W[nm])+anorm) == anorm) + if (ei_abs(W[nm]) <= eps*anorm) + break; } - if (ks == k) - { - kase = 3; - } - else if (ks == p-1) - { - kase = 1; - } - else - { - kase = 2; - k = ks; - } - } - ++k; - - // Perform the task indicated by kase. - switch (kase) - { - - // Deflate negligible s(p). - case 1: + if (flag) { - Scalar f(e[p-2]); - e[p-2] = 0.0; - for (j = p-2; j >= k; --j) + c = 0.0; //Cancellation of rv1[l], if l > 0. + s = 1.0; + for (i=l ;i<k+1; i++) { - Scalar t(ei_hypot(m_sigma[j],f)); - Scalar cs(m_sigma[j]/t); - Scalar sn(f/t); - m_sigma[j] = t; - if (j != k) - { - f = -sn*e[j-1]; - e[j-1] = cs*e[j-1]; - } - if (wantv) + f = s*rv1[i]; + rv1[i] = c*rv1[i]; + //if ((double)(ei_abs(f)+anorm) == anorm) + if (ei_abs(f) <= eps*anorm) + break; + g = W[i]; + h = pythag(f,g); + W[i] = h; + h = Scalar(1.0)/h; + c = g*h; + s = -f*h; + for (j=0; j<m; j++) { - for (i = 0; i < n; ++i) - { - t = cs*m_matV(i,j) + sn*m_matV(i,p-1); - m_matV(i,p-1) = -sn*m_matV(i,j) + cs*m_matV(i,p-1); - m_matV(i,j) = t; - } + y = A(j,nm); + z = A(j,i); + A(j,nm) = y*c + z*s; + A(j,i) = z*c - y*s; } } } - break; - - // Split at negligible s(k). - case 2: + z = W[k]; + if (l == k) //Convergence. { - Scalar f(e[k-1]); - e[k-1] = 0.0; - for (j = k; j < p; ++j) - { - Scalar t(ei_hypot(m_sigma[j],f)); - Scalar cs( m_sigma[j]/t); - Scalar sn(f/t); - m_sigma[j] = t; - f = -sn*e[j]; - e[j] = cs*e[j]; - if (wantu) - { - for (i = 0; i < m; ++i) - { - t = cs*m_matU(i,j) + sn*m_matU(i,k-1); - m_matU(i,k-1) = -sn*m_matU(i,j) + cs*m_matU(i,k-1); - m_matU(i,j) = t; - } - } + if (z < 0.0) { // Singular value is made nonnegative. + W[k] = -z; + V.col(k) = -V.col(k); } + break; } - break; - - // Perform one qr step. - case 3: + if (its+1 == max_iters) { - // Calculate the shift. - Scalar scale = std::max(std::max(std::max(std::max( - ei_abs(m_sigma[p-1]),ei_abs(m_sigma[p-2])),ei_abs(e[p-2])), - ei_abs(m_sigma[k])),ei_abs(e[k])); - Scalar sp = m_sigma[p-1]/scale; - Scalar spm1 = m_sigma[p-2]/scale; - Scalar epm1 = e[p-2]/scale; - Scalar sk = m_sigma[k]/scale; - Scalar ek = e[k]/scale; - Scalar b = ((spm1 + sp)*(spm1 - sp) + epm1*epm1)/Scalar(2); - Scalar c = (sp*epm1)*(sp*epm1); - Scalar shift = 0.0; - if ((b != 0.0) || (c != 0.0)) + convergence = false; + } + x = W[l]; // Shift from bottom 2-by-2 minor. + nm = k-1; + y = W[nm]; + g = rv1[nm]; + h = rv1[k]; + f = ((y-z)*(y+z) + (g-h)*(g+h))/(Scalar(2.0)*h*y); + g = pythag(f,1.0); + f = ((x-z)*(x+z) + h*((y/(f+sign(g,f)))-h))/x; + c = s = 1.0; + //Next QR transformation: + for (j=l; j<=nm; j++) + { + i = j+1; + g = rv1[i]; + y = W[i]; + h = s*g; + g = c*g; + z = pythag(f,h); + rv1[j] = z; + c = f/z; + s = h/z; + f = x*c + g*s; + g = g*c - x*s; + h = y*s; + y *= c; + for (jj=0; jj<n; jj++) { - shift = ei_sqrt(b*b + c); - if (b < 0.0) - shift = -shift; - shift = c/(b + shift); + x = V(jj,j); + z = V(jj,i); + V(jj,j) = x*c + z*s; + V(jj,i) = z*c - x*s; } - Scalar f = (sk + sp)*(sk - sp) + shift; - Scalar g = sk*ek; - - // Chase zeros. - - for (j = k; j < p-1; ++j) + z = pythag(f,h); + W[j] = z; + // Rotation can be arbitrary if z = 0. + if (z!=Scalar(0)) { - Scalar t = ei_hypot(f,g); - Scalar cs = f/t; - Scalar sn = g/t; - if (j != k) - e[j-1] = t; - f = cs*m_sigma[j] + sn*e[j]; - e[j] = cs*e[j] - sn*m_sigma[j]; - g = sn*m_sigma[j+1]; - m_sigma[j+1] = cs*m_sigma[j+1]; - if (wantv) - { - for (i = 0; i < n; ++i) - { - t = cs*m_matV(i,j) + sn*m_matV(i,j+1); - m_matV(i,j+1) = -sn*m_matV(i,j) + cs*m_matV(i,j+1); - m_matV(i,j) = t; - } - } - t = ei_hypot(f,g); - cs = f/t; - sn = g/t; - m_sigma[j] = t; - f = cs*e[j] + sn*m_sigma[j+1]; - m_sigma[j+1] = -sn*e[j] + cs*m_sigma[j+1]; - g = sn*e[j+1]; - e[j+1] = cs*e[j+1]; - if (wantu && (j < m-1)) - { - for (i = 0; i < m; ++i) - { - t = cs*m_matU(i,j) + sn*m_matU(i,j+1); - m_matU(i,j+1) = -sn*m_matU(i,j) + cs*m_matU(i,j+1); - m_matU(i,j) = t; - } - } + z = Scalar(1.0)/z; + c = f*z; + s = h*z; } - e[p-2] = f; - iter = iter + 1; - } - break; - - // Convergence. - case 4: - { - // Make the singular values positive. - if (m_sigma[k] <= 0.0) + f = c*g + s*y; + x = c*y - s*g; + for (jj=0; jj<m; jj++) { - m_sigma[k] = m_sigma[k] < Scalar(0) ? -m_sigma[k] : Scalar(0); - if (wantv) - m_matV.col(k).start(pp+1) = -m_matV.col(k).start(pp+1); + y = A(jj,j); + z = A(jj,i); + A(jj,j) = y*c + z*s; + A(jj,i) = z*c - y*s; } + } + rv1[l] = 0.0; + rv1[k] = f; + W[k] = x; + } + } - // Order the singular values. - while (k < pp) - { - if (m_sigma[k] >= m_sigma[k+1]) - break; - Scalar t = m_sigma[k]; - m_sigma[k] = m_sigma[k+1]; - m_sigma[k+1] = t; - if (wantv && (k < n-1)) - m_matV.col(k).swap(m_matV.col(k+1)); - if (wantu && (k < m-1)) - m_matU.col(k).swap(m_matU.col(k+1)); - ++k; - } - iter = 0; - p--; + // sort the singular values: + { + for (int i=0; i<n; i++) + { + int k; + W.end(n-i).minCoeff(&k); + if (k != i) + { + std::swap(W[k],W[i]); + A.col(i).swap(A.col(k)); + V.col(i).swap(V.col(k)); } - break; - } // end big switch - } // end iterations + } + } + m_matU.setZero(); + if (m>=n) + m_matU.block(0,0,m,n) = A; + else + m_matU = A.block(0,0,m,m); m_isInitialized = true; } @@ -554,6 +460,8 @@ bool SVD<MatrixType>::solve(const MatrixBase<OtherDerived> &b, ResultType* resul const int rows = m_matU.rows(); ei_assert(b.rows() == rows); + result->resize(m_matV.rows(), b.cols()); + Scalar maxVal = m_sigma.cwise().abs().maxCoeff(); for (int j=0; j<b.cols(); ++j) { diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h index b751b5cb9..af3b5e5eb 100644 --- a/Eigen/src/Sparse/SparseMatrix.h +++ b/Eigen/src/Sparse/SparseMatrix.h @@ -443,9 +443,8 @@ class SparseMatrix // two passes algorithm: // 1 - compute the number of coeffs per dest inner vector // 2 - do the actual copy/eval - // Since each coeff of the rhs has to be evaluated twice, let's evauluate it if needed - //typedef typename ei_nested<OtherDerived,2>::type OtherCopy; - typedef typename ei_eval<OtherDerived>::type OtherCopy; + // Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed + typedef typename ei_nested<OtherDerived,2>::type OtherCopy; typedef typename ei_cleantype<OtherCopy>::type _OtherCopy; OtherCopy otherCopy(other.derived()); diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h index d62d23a78..1af452251 100644 --- a/Eigen/src/Sparse/SparseMatrixBase.h +++ b/Eigen/src/Sparse/SparseMatrixBase.h @@ -102,8 +102,10 @@ template<typename Derived> class SparseMatrixBase /** \internal the return type of MatrixBase::imag() */ typedef SparseCwiseUnaryOp<ei_scalar_imag_op<Scalar>, Derived> ImagReturnType; /** \internal the return type of MatrixBase::adjoint() */ - typedef SparseTranspose</*NestByValue<*/typename ei_cleantype<ConjugateReturnType>::type> /*>*/ - AdjointReturnType; + typedef typename ei_meta_if<NumTraits<Scalar>::IsComplex, + SparseCwiseUnaryOp<ei_scalar_conjugate_op<Scalar>, SparseNestByValue<Eigen::SparseTranspose<Derived> > >, + SparseTranspose<Derived> + >::ret AdjointReturnType; #ifndef EIGEN_PARSED_BY_DOXYGEN /** This is the "real scalar" type; if the \a Scalar type is already real numbers @@ -357,7 +359,7 @@ template<typename Derived> class SparseMatrixBase SparseTranspose<Derived> transpose() { return derived(); } const SparseTranspose<Derived> transpose() const { return derived(); } // void transposeInPlace(); - const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; } + const AdjointReturnType adjoint() const { return transpose().nestByValue(); } // sub-vector SparseInnerVectorSet<Derived,1> row(int i); @@ -529,7 +531,7 @@ template<typename Derived> class SparseMatrixBase */ // inline int stride(void) const { return derived().stride(); } -// inline const NestByValue<Derived> nestByValue() const; + inline const SparseNestByValue<Derived> nestByValue() const; ConjugateReturnType conjugate() const; diff --git a/Eigen/src/Sparse/SparseNestByValue.h b/Eigen/src/Sparse/SparseNestByValue.h new file mode 100644 index 000000000..b48277232 --- /dev/null +++ b/Eigen/src/Sparse/SparseNestByValue.h @@ -0,0 +1,84 @@ +// 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_SPARSENESTBYVALUE_H +#define EIGEN_SPARSENESTBYVALUE_H + +/** \class SparseNestByValue + * + * \brief Expression which must be nested by value + * + * \param ExpressionType the type of the object of which we are requiring nesting-by-value + * + * This class is the return type of MatrixBase::nestByValue() + * and most of the time this is the only way it is used. + * + * \sa SparseMatrixBase::nestByValue(), class NestByValue + */ +template<typename ExpressionType> +struct ei_traits<SparseNestByValue<ExpressionType> > : public ei_traits<ExpressionType> +{}; + +template<typename ExpressionType> class SparseNestByValue + : public SparseMatrixBase<SparseNestByValue<ExpressionType> > +{ + public: + + typedef typename ExpressionType::InnerIterator InnerIterator; + + EIGEN_SPARSE_GENERIC_PUBLIC_INTERFACE(SparseNestByValue) + + inline SparseNestByValue(const ExpressionType& matrix) : m_expression(matrix) {} + + EIGEN_STRONG_INLINE int rows() const { return m_expression.rows(); } + EIGEN_STRONG_INLINE int cols() const { return m_expression.cols(); } + + operator const ExpressionType&() const { return m_expression; } + + protected: + const ExpressionType m_expression; +}; + +/** \returns an expression of the temporary version of *this. + */ +template<typename Derived> +inline const SparseNestByValue<Derived> +SparseMatrixBase<Derived>::nestByValue() const +{ + return SparseNestByValue<Derived>(derived()); +} + +// template<typename MatrixType> +// class SparseNestByValue<MatrixType>::InnerIterator : public MatrixType::InnerIterator +// { +// typedef typename MatrixType::InnerIterator Base; +// public: +// +// EIGEN_STRONG_INLINE InnerIterator(const SparseNestByValue& expr, int outer) +// : Base(expr.m_expression, outer) +// {} +// }; + +#endif // EIGEN_SPARSENESTBYVALUE_H diff --git a/Eigen/src/Sparse/SparseUtil.h b/Eigen/src/Sparse/SparseUtil.h index 52781aa46..b5fc7c7b7 100644 --- a/Eigen/src/Sparse/SparseUtil.h +++ b/Eigen/src/Sparse/SparseUtil.h @@ -106,6 +106,7 @@ template<typename _Scalar, int _Flags = 0> class DynamicSparseMatrix; template<typename _Scalar, int _Flags = 0> class SparseVector; template<typename _Scalar, int _Flags = 0> class MappedSparseMatrix; +template<typename MatrixType> class SparseNestByValue; template<typename MatrixType> class SparseTranspose; template<typename MatrixType, int Size> class SparseInnerVectorSet; template<typename Derived> class SparseCwise; @@ -146,4 +147,6 @@ template<typename T> class ei_eval<T,IsSparse> typedef SparseMatrix<_Scalar, _Flags> type; }; +template<typename T> struct ei_must_nest_by_value<SparseNestByValue<T> > { enum { ret = true }; }; + #endif // EIGEN_SPARSEUTIL_H |