diff --git a/Eigen/Core b/Eigen/Core index 87e7e9a01..2bec68f9c 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -88,6 +88,8 @@ #include #include #include +// for min/max: +#include #if (defined(_CPPUNWIND) || defined(__EXCEPTIONS)) && !defined(EIGEN_NO_EXCEPTIONS) #define EIGEN_EXCEPTIONS @@ -162,6 +164,7 @@ namespace Eigen { #include "src/Core/MapBase.h" #include "src/Core/Map.h" #include "src/Core/Block.h" +#include "src/Core/VectorBlock.h" #include "src/Core/Minor.h" #include "src/Core/Transpose.h" #include "src/Core/DiagonalMatrix.h" @@ -182,7 +185,6 @@ namespace Eigen { #include "src/Core/SolveTriangular.h" #include "src/Core/products/SelfadjointRank2Update.h" #include "src/Core/products/TriangularMatrixVector.h" -#include "src/Core/BandMatrix.h" } // namespace Eigen diff --git a/Eigen/Geometry b/Eigen/Geometry index 7860d8fca..9cae3459c 100644 --- a/Eigen/Geometry +++ b/Eigen/Geometry @@ -6,6 +6,7 @@ #include "src/Core/util/DisableMSVCWarnings.h" #include "Array" +#include "SVD" #include #ifndef M_PI diff --git a/Eigen/Sparse b/Eigen/Sparse index 364fd50c9..a8888daa3 100644 --- a/Eigen/Sparse +++ b/Eigen/Sparse @@ -84,6 +84,7 @@ namespace Eigen { #include "src/Sparse/SparseUtil.h" #include "src/Sparse/SparseMatrixBase.h" +#include "src/Sparse/SparseNestByValue.h" #include "src/Sparse/CompressedStorage.h" #include "src/Sparse/AmbiVector.h" #include "src/Sparse/RandomSetter.h" 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::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 -template -const Replicate -VectorwiseOp::replicate(int factor) const -{ - return Replicate - (_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 class VectorwiseOp > Type; }; + enum { + IsVertical = (Direction==Vertical) ? 1 : 0, + IsHorizontal = (Direction==Horizontal) ? 1 : 0 + }; + protected: /** \internal @@ -222,9 +227,17 @@ template 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 const typename ReduxReturnType::Type - redux(const BinaryOp& func = BinaryOp()) const; + redux(const BinaryOp& func = BinaryOp()) const + { return typename ReduxReturnType::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 class VectorwiseOp * * \sa MatrixBase::reverse() */ const Reverse reverse() const - { - return Reverse( _expression() ); - } + { return Reverse( _expression() ); } const Replicate replicate(int factor) const; - template - const Replicate - 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 const Replicate + replicate(int factor = Factor) const + { + return Replicate + (_expression(),Direction==Vertical?factor:1,Direction==Horizontal?factor:1); + } /////////// Artithmetic operators /////////// @@ -466,19 +489,4 @@ MatrixBase::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 -template -const typename VectorwiseOp::template ReduxReturnType::Type -VectorwiseOp::redux(const BinaryOp& func) const -{ - return typename ReduxReturnType::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 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::Type MatrixBase return typename BlockReturnType::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 -inline typename BlockReturnType::SubVectorType MatrixBase - ::segment(int start, int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename BlockReturnType::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 -inline const typename BlockReturnType::SubVectorType -MatrixBase::segment(int start, int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return typename BlockReturnType::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 -inline typename BlockReturnType::SubVectorType -MatrixBase::start(int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), 0, 0, - RowsAtCompileTime == 1 ? 1 : size, - ColsAtCompileTime == 1 ? 1 : size); -} - -/** This is the const version of start(int).*/ -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::start(int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (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 -inline typename BlockReturnType::SubVectorType -MatrixBase::end(int size) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (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 -inline const typename BlockReturnType::SubVectorType -MatrixBase::end(int size) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (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 -template -inline typename BlockReturnType::SubVectorType -MatrixBase::segment(int start) -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), RowsAtCompileTime == 1 ? 0 : start, - ColsAtCompileTime == 1 ? 0 : start); -} - -/** This is the const version of segment(int).*/ -template -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::segment(int start) const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (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 -template -inline typename BlockReturnType::SubVectorType -MatrixBase::start() -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block(derived(), 0, 0); -} - -/** This is the const version of start().*/ -template -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::start() const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block(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 -template -inline typename BlockReturnType::SubVectorType -MatrixBase::end() -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (derived(), - RowsAtCompileTime == 1 ? 0 : rows() - Size, - ColsAtCompileTime == 1 ? 0 : cols() - Size); -} - -/** This is the const version of end.*/ -template -template -inline const typename BlockReturnType::SubVectorType -MatrixBase::end() const -{ - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) - return Block - (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 > : ei_traits { typedef typename ei_result_of< - ViewOp(typename MatrixType::Scalar) + ViewOp(typename ei_traits::Scalar) >::type Scalar; typedef typename MatrixType::Nested MatrixTypeNested; - typedef typename ei_unref::type _MatrixTypeNested; + typedef typename ei_cleantype::type _MatrixTypeNested; enum { - Flags = (_MatrixTypeNested::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)), - CoeffReadCost = _MatrixTypeNested::CoeffReadCost + ei_functor_traits::Cost + Flags = (ei_traits<_MatrixTypeNested>::Flags & (HereditaryBits | LinearAccessBit | AlignedBit)), + CoeffReadCost = ei_traits<_MatrixTypeNested>::CoeffReadCost + ei_functor_traits::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::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::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 inline typename NumTraits::Scalar>::Real @@ -304,6 +304,142 @@ MatrixBase::stableNorm() const return this->cwise().abs().redux(ei_scalar_hypot_op()); } +/** \internal Computes ibeta^iexp by binary expansion of iexp, + * exact if ibeta is the machine base */ +template 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 +inline typename NumTraits::Scalar>::Real +MatrixBase::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::max(); // largest integer + ibeta = NumTraits::Base; // base for floating-point numbers + it = NumTraits::Mantissa; // number of base-beta digits in mantissa + iemin = std::numeric_limits::min_exponent; // minimum exponent + iemax = std::numeric_limits::max_exponent; // maximum exponent + rbig = std::numeric_limits::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(ibeta, iexp); // lower boundary of midrange + iexp = (iemax + 1 - it)/2; + b2 = bexp(ibeta,iexp); // upper boundary of midrange + + iexp = (2-iemin)/2; + s1m = bexp(ibeta,iexp); // scaling factor for lower range + iexp = - ((iemax+it)/2); + s2m = bexp(ibeta,iexp); // scaling factor for upper range + + overfl = rbig*s2m; // overfow boundary for abig + eps = bexp(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 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 inline double precision() { return 1e-11; } +template<> inline double precision() { return 1e-12; } template<> inline double machine_epsilon() { 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; typedef class Eigen::Map 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 class MatrixBase * constructed from this one. See the \ref flags "list of flags". */ - CoeffReadCost = ei_traits::CoeffReadCost + CoeffReadCost = ei_traits::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 class MatrixBase /** \internal the return type of coeff() */ - typedef typename ei_meta_if::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,Derived> ConstantReturnType; @@ -426,8 +430,9 @@ template class MatrixBase template Scalar dot(const MatrixBase& 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 class MatrixBase const typename BlockReturnType::Type block(int startRow, int startCol, int blockRows, int blockCols) const; - typename BlockReturnType::SubVectorType segment(int start, int size); - const typename BlockReturnType::SubVectorType segment(int start, int size) const; + VectorBlock segment(int start, int size); + const VectorBlock segment(int start, int size) const; - typename BlockReturnType::SubVectorType start(int size); - const typename BlockReturnType::SubVectorType start(int size) const; + VectorBlock start(int size); + const VectorBlock start(int size) const; - typename BlockReturnType::SubVectorType end(int size); - const typename BlockReturnType::SubVectorType end(int size) const; + VectorBlock end(int size); + const VectorBlock end(int size) const; typename BlockReturnType::Type corner(CornerType type, int cRows, int cCols); const typename BlockReturnType::Type corner(CornerType type, int cRows, int cCols) const; @@ -472,14 +477,14 @@ template class MatrixBase template const typename BlockReturnType::Type corner(CornerType type) const; - template typename BlockReturnType::SubVectorType start(void); - template const typename BlockReturnType::SubVectorType start() const; + template VectorBlock start(void); + template const VectorBlock start() const; - template typename BlockReturnType::SubVectorType end(); - template const typename BlockReturnType::SubVectorType end() const; + template VectorBlock end(); + template const VectorBlock end() const; - template typename BlockReturnType::SubVectorType segment(int start); - template const typename BlockReturnType::SubVectorType segment(int start) const; + template VectorBlock segment(int start); + template const VectorBlock segment(int start) const; Diagonal diagonal(); const Diagonal diagonal() const; @@ -696,6 +701,7 @@ template class MatrixBase const PartialLU 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 class MatrixBase /////////// QR module /////////// - const QR qr() const; + const HouseholderQR 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 HasFloatingPoint = 1, ReadCost = 1, AddCost = 1, - MulCost = 1 + MulCost = 1, + Base = 2, + Mantissa = 23 }; }; @@ -83,7 +85,9 @@ template<> struct NumTraits 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 struct ei_product_mode typedef typename ei_blas_traits::ActualXprType ActualLhs; typedef typename ei_blas_traits::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::ret) + // workaround sun studio: + LhsIsVectorAtCompileTime = ei_traits::ColsAtCompileTime==1 || ei_traits::ColsAtCompileTime==1, + value = ei_traits::MaxColsAtCompileTime == Dynamic + && ( ei_traits::MaxRowsAtCompileTime == Dynamic + || ei_traits::MaxColsAtCompileTime == Dynamic ) + && (!(Rhs::IsVectorAtCompileTime && (ei_traits::Flags&RowMajorBit) && (!(ei_traits::Flags&DirectAccessBit)))) + && (!(LhsIsVectorAtCompileTime && (!(ei_traits::Flags&RowMajorBit)) && (!(ei_traits::Flags&DirectAccessBit)))) + && (ei_is_same_type::Scalar, typename ei_traits::Scalar>::ret) ? CacheFriendlyProduct : NormalProduct }; }; @@ -577,7 +578,7 @@ struct ei_product_packet_impl -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 { template struct ei_triangular_solver_unroller { - static void run(const Lhs& lhs, Rhs& rhs) {} + static void run(const Lhs&, Rhs&) {} }; template 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 +// Copyright (C) 2006-2008 Benoit Jacob +// +// 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 . + +#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) 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 +struct ei_traits > + : public ei_traits::RowsAtCompileTime==1 ? 1 : Size, + ei_traits::ColsAtCompileTime==1 ? 1 : Size, + _PacketAccess> > +{ +}; + +template class VectorBlock + : public Block::RowsAtCompileTime==1 ? 1 : Size, + ei_traits::ColsAtCompileTime==1 ? 1 : Size, + PacketAccess> +{ + typedef Block::RowsAtCompileTime==1 ? 1 : Size, + ei_traits::ColsAtCompileTime==1 ? 1 : Size, + PacketAccess> Base; + enum { + IsColVector = ei_traits::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 +inline VectorBlock MatrixBase + ::segment(int start, int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), start, size); +} + +/** This is the const version of segment(int,int).*/ +template +inline const VectorBlock +MatrixBase::segment(int start, int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(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 +inline VectorBlock +MatrixBase::start(int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), 0, size); +} + +/** This is the const version of start(int).*/ +template +inline const VectorBlock +MatrixBase::start(int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(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 +inline VectorBlock +MatrixBase::end(int size) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), this->size() - size, size); +} + +/** This is the const version of end(int).*/ +template +inline const VectorBlock +MatrixBase::end(int size) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(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 +template +inline VectorBlock +MatrixBase::segment(int start) +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), start); +} + +/** This is the const version of segment(int).*/ +template +template +inline const VectorBlock +MatrixBase::segment(int start) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(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 +template +inline VectorBlock +MatrixBase::start() +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), 0); +} + +/** This is the const version of start().*/ +template +template +inline const VectorBlock +MatrixBase::start() const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(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 +template +inline VectorBlock +MatrixBase::end() +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(derived(), size() - Size); +} + +/** This is the const version of end.*/ +template +template +inline const VectorBlock +MatrixBase::end() const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + return VectorBlock(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 class Minor; template::Flags&DirectAccessBit ? DirectAccessBit : ei_traits::Flags&SparseBit> class Block; +template class VectorBlock; template class Transpose; template class Conjugate; template class CwiseNullaryOp; @@ -111,7 +112,7 @@ template class Reverse; template class LU; template class PartialLU; -template class QR; +template class HouseholderQR; template class SVD; template class LLT; template 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::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 * TODO: could be a good idea to define a big ReturnType struct ?? */ template struct BlockReturnType { - typedef Block::RowsAtCompileTime == 1 ? 1 : RowsOrSize), - (ei_traits::ColsAtCompileTime == 1 ? 1 : RowsOrSize)> SubVectorType; typedef Block Type; }; @@ -251,7 +249,7 @@ template struct HNormalizedReturnType { typedef Block::ColsAtCompileTime==1 ? SizeMinusOne : 1, ei_traits::ColsAtCompileTime==1 ? 1 : SizeMinusOne> StartMinusOne; - typedef CwiseUnaryOp::Scalar>, + typedef CwiseUnaryOp::Scalar>, NestByValue > 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& Quaternion::operator=(const MatrixBase inline typename Quaternion::Matrix3 Quaternion::toRotationMatrix(void) const @@ -340,11 +342,15 @@ Quaternion::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 template @@ -354,21 +360,26 @@ inline Quaternion& Quaternion::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()) { - // set to identity - this->w() = 1; this->vec().setZero(); - return *this; - } - // if dot == -1, vectors are opposites - if (ei_isApprox(c,Scalar(-1))) - { - this->vec() = v0.unitOrthogonal(); - this->w() = 0; - return *this; - } + c = std::max(c,-1); + Matrix m; m << v0.transpose(), v1.transpose(); + SVD > 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; + } 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 -void ei_compute_inverse_in_size2_case(const MatrixType& matrix, MatrixType* result) +template +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 +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 -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 -void ei_compute_inverse_in_size3_case(const MatrixType& matrix, MatrixType* result) +template +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 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 -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 XprBlock22; typedef typename MatrixBase::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 -void ei_compute_inverse_in_size4_case(const MatrixType& matrix, MatrixType* result) +template +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 { 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 { static inline void run(const MatrixType& matrix, MatrixType* result) { - ei_compute_inverse_in_size3_case(matrix, result); + ei_compute_inverse_size3(matrix, result); } }; @@ -207,7 +237,7 @@ struct ei_compute_inverse { 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 * * 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 inline void MatrixBase::computeInverse(PlainMatrixType *result) const @@ -236,7 +267,8 @@ inline void MatrixBase::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::computeInverse(PlainMatrixType *result) const * Example: \include MatrixBase_inverse.cpp * Output: \verbinclude MatrixBase_inverse.out * - * \sa computeInverse() + * \sa computeInverse(), computeInverseWithCheck() */ template inline const typename MatrixBase::PlainMatrixType MatrixBase::inverse() const @@ -254,4 +286,81 @@ inline const typename MatrixBase::PlainMatrixType MatrixBase:: return result; } + +/******************************************** + * Compute inverse with invertibility check * + *******************************************/ + +template +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + typedef typename MatrixType::Scalar Scalar; + LU lu( matrix ); + if( !lu.isInvertible() ) return false; + lu.computeInverse(result); + return true; + } +}; + +template +struct ei_compute_inverse_with_check +{ + 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 +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + return ei_compute_inverse_size2_with_check(matrix, result); + } +}; + +template +struct ei_compute_inverse_with_check +{ + static inline bool run(const MatrixType& matrix, MatrixType* result) + { + return ei_compute_inverse_size3(matrix, result); + } +}; + +template +struct ei_compute_inverse_with_check +{ + 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 +inline bool MatrixBase::computeInverseWithCheck(PlainMatrixType *result) const +{ + ei_assert(rows() == cols()); + EIGEN_STATIC_ASSERT(NumTraits::HasFloatingPoint,NUMERIC_TYPE_MUST_BE_FLOATING_POINT) + return ei_compute_inverse_with_check::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 class QR +template class HouseholderQR { public: @@ -53,88 +55,23 @@ template 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, 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(); } @@ -148,58 +85,35 @@ template 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 - bool solve(const MatrixBase& b, ResultType *result) const; + void solve(const MatrixBase& 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 -int QR::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 -void QR::compute(const MatrixType& matrix) +void HouseholderQR::compute(const MatrixType& matrix) { - m_rankIsUptodate = false; m_qr = matrix; m_hCoeffs.resize(matrix.cols()); @@ -262,12 +176,12 @@ void QR::compute(const MatrixType& matrix) template template -bool QR::solve( +void HouseholderQR::solve( const MatrixBase& 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::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() - .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols())); - - return true; + .solveTriangularInPlace(result->corner(TopLeft, rank, result->cols())); } /** \returns the matrix Q */ template -MatrixType QR::matrixQ() const +MatrixType HouseholderQR::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::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 -const QR::PlainMatrixType> -MatrixBase::qr() const +const HouseholderQR::PlainMatrixType> +MatrixBase::householderQr() const { - return QR(eval()); + return HouseholderQR(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 +// Copyright (C) 2008-2009 Gael Guennebaud // // 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 class SVD typedef Matrix ColVector; typedef Matrix RowVector; - typedef Matrix MatrixUType; + typedef Matrix MatrixUType; typedef Matrix MatrixVType; - typedef Matrix SingularValuesType; + typedef Matrix SingularValuesType; public: - /** + /** * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to @@ -70,9 +68,9 @@ template 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 class SVD template bool solve(const MatrixBase &b, ResultType* result) const; - const MatrixUType& matrixU() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matU; - } - - const SingularValuesType& singularValues() const + const MatrixUType& matrixU() const { ei_assert(m_isInitialized && "SVD is not initialized."); - return m_sigma; + return m_matU; } - const MatrixVType& matrixV() const + const SingularValuesType& singularValues() const { ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matV; + return m_sigma; + } + + const MatrixVType& matrixV() const + { + ei_assert(m_isInitialized && "SVD is not initialized."); + return m_matV; } void compute(const MatrixType& matrix); @@ -111,6 +109,23 @@ template class SVD template 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; @@ -123,380 +138,271 @@ template 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 void SVD::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; + int max_iters = 30; - // 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) + MatrixVType& V = m_matV; + MatrixType A = matrix; + SingularValuesType& W = m_sigma; + + 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(); + + Matrix rv1(n); + g = scale = anorm = 0; + // Householder reduction to bidiagonal form. + for (i=0; i= 0; k--) - { - if (m_sigma[k] != 0.0) - { - for (j = k+1; j < nu; ++j) + for (k=i; k0) - m_matU.col(k).start(k-1).setZero(); - } - else - { - m_matU.col(k).setZero(); - m_matU(k,k) = 1.0; + f = A(i, i); + g = -sign( ei_sqrt(s), f ); + h = f*g - s; + A(i, i)=f-g; + for (j=l-1; j= 0; k--) + W[i] = scale * g; + g = s = scale = 0.0; + if (i+1 <= m && i+1 != n) { - if ((k < nrt) & (e[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=0; i--) + { + //Accumulation of right-hand transformations. + if (i < n-1) + { + if (g != Scalar(0.0)) + { + for (j=l; j::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

= -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=0; k--) + { + for (its=0; its= 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) + if (flag) { - kase = 3; + c = 0.0; //Cancellation of rv1[l], if l > 0. + s = 1.0; + for (i=l ;i= k; --j) - { - 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) - { - 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; - } - } - } - } - break; - - // Split at negligible s(k). - case 2: - { - 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; - } - } - } - } - break; - - // Perform one qr step. - case 3: - { - // 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)) - { - shift = ei_sqrt(b*b + c); - if (b < 0.0) - shift = -shift; - shift = c/(b + shift); - } - Scalar f = (sk + sp)*(sk - sp) + shift; - Scalar g = sk*ek; - - // Chase zeros. - - for (j = k; j < p-1; ++j) - { - 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; - } - } - } - e[p-2] = f; - iter = iter + 1; - } - break; - - // Convergence. - case 4: - { - // Make the singular values positive. - if (m_sigma[k] <= 0.0) - { - 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); - } - - // 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--; - } - 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::solve(const MatrixBase &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::type OtherCopy; - typedef typename ei_eval::type OtherCopy; + // Since each coeff of the rhs has to be evaluated twice, let's evaluate it if needed + typedef typename ei_nested::type OtherCopy; typedef typename ei_cleantype::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 class SparseMatrixBase /** \internal the return type of MatrixBase::imag() */ typedef SparseCwiseUnaryOp, Derived> ImagReturnType; /** \internal the return type of MatrixBase::adjoint() */ - typedef SparseTranspose::type> /*>*/ - AdjointReturnType; + typedef typename ei_meta_if::IsComplex, + SparseCwiseUnaryOp, SparseNestByValue > >, + SparseTranspose + >::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 class SparseMatrixBase SparseTranspose transpose() { return derived(); } const SparseTranspose transpose() const { return derived(); } // void transposeInPlace(); - const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; } + const AdjointReturnType adjoint() const { return transpose().nestByValue(); } // sub-vector SparseInnerVectorSet row(int i); @@ -529,7 +531,7 @@ template class SparseMatrixBase */ // inline int stride(void) const { return derived().stride(); } -// inline const NestByValue nestByValue() const; + inline const SparseNestByValue 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 +// Copyright (C) 2006-2008 Benoit Jacob +// +// 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 . + +#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 +struct ei_traits > : public ei_traits +{}; + +template class SparseNestByValue + : public SparseMatrixBase > +{ + 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 +inline const SparseNestByValue +SparseMatrixBase::nestByValue() const +{ + return SparseNestByValue(derived()); +} + +// template +// class SparseNestByValue::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 class DynamicSparseMatrix; template class SparseVector; template class MappedSparseMatrix; +template class SparseNestByValue; template class SparseTranspose; template class SparseInnerVectorSet; template class SparseCwise; @@ -146,4 +147,6 @@ template class ei_eval typedef SparseMatrix<_Scalar, _Flags> type; }; +template struct ei_must_nest_by_value > { enum { ret = true }; }; + #endif // EIGEN_SPARSEUTIL_H diff --git a/cmake/FindEigen2.cmake b/cmake/FindEigen2.cmake index 0e86aaf29..ee054b0da 100644 --- a/cmake/FindEigen2.cmake +++ b/cmake/FindEigen2.cmake @@ -46,7 +46,7 @@ macro(_eigen2_check_version) endif(${EIGEN2_VERSION} VERSION_LESS ${Eigen2_FIND_VERSION}) if(NOT EIGEN2_VERSION_OK) - + message(STATUS "Eigen2 version ${EIGEN2_VERSION} found in ${EIGEN2_INCLUDE_DIR}, " "but at least version ${Eigen2_FIND_VERSION} is required") endif(NOT EIGEN2_VERSION_OK) diff --git a/doc/D09_StructHavingEigenMembers.dox b/doc/D09_StructHavingEigenMembers.dox index 4ff2b23aa..07b6570dc 100644 --- a/doc/D09_StructHavingEigenMembers.dox +++ b/doc/D09_StructHavingEigenMembers.dox @@ -88,7 +88,7 @@ The solution is to let class Foo have an aligned "operator new", as we showed in \section movetotop Should I then put all the members of Eigen types at the beginning of my class? -No, that's not needed. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So when you have code like +That's not required. Since Eigen takes care of declaring 128-bit alignment, all members that need it are automatically 128-bit aligned relatively to the class. So code like this works fine: \code class Foo @@ -100,25 +100,13 @@ public: }; \endcode -it will work just fine. You do \b not need to rewrite it as - -\code -class Foo -{ - Eigen::Vector2d v; - double x; -public: - EIGEN_MAKE_ALIGNED_OPERATOR_NEW -}; -\endcode - \section dynamicsize What about dynamic-size matrices and vectors? -Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with fixed-size matrices and vectors. +Dynamic-size matrices and vectors, such as Eigen::VectorXd, allocate dynamically their own array of coefficients, so they take care of requiring absolute alignment automatically. So they don't cause this issue. The issue discussed here is only with \ref FixedSizeVectorizable "fixed-size vectorizable matrices and vectors". \section bugineigen So is this a bug in Eigen? -No, it's not our bug. It's more like an inherent problem of the C++ language -- though it must be said that any other existing language probably has the same problem. The problem is that there is no way that you can specify an aligned "operator new" that would propagate to classes having you as member data. +No, it's not our bug. It's more like an inherent problem of the C++98 language specification, and seems to be taken care of in the upcoming language revision: see this document. \section conditional What if I want to do this conditionnally (depending on template parameters) ? diff --git a/doc/D11_UnalignedArrayAssert.dox b/doc/D11_UnalignedArrayAssert.dox index da1b49a08..e95c767b4 100644 --- a/doc/D11_UnalignedArrayAssert.dox +++ b/doc/D11_UnalignedArrayAssert.dox @@ -15,12 +15,24 @@ is explained here: http://eigen.tuxfamily.org/dox/UnalignedArrayAssert.html There are 4 known causes for this issue. Please read on to understand them and learn how to fix them. \b Table \b of \b contents + - \ref where - \ref c1 - \ref c2 - \ref c3 - \ref c4 - \ref explanation +\section where Where in my own code is the cause of the problem? + +First of all, you need to find out where in your own code this assertion was triggered from. At first glance, the error message doesn't look helpful, as it refers to a file inside Eigen! However, since your program crashed, if you can reproduce the crash, you can get a backtrace using any debugger. For example, if you're using GCC, you can use the GDB debugger as follows: +\code +$ gdb ./my_program # Start GDB on your program +> run # Start running your program +... # Now reproduce the crash! +> bt # Obtain the backtrace +\endcode +Now that you know precisely where in your own code the problem is happening, read on to understand what you need to change. + \section c1 Cause 1: Structures having Eigen objects as members If you have code like this, diff --git a/doc/eigendoxy_header.html.in b/doc/eigendoxy_header.html.in index 97cafd41d..572c47158 100644 --- a/doc/eigendoxy_header.html.in +++ b/doc/eigendoxy_header.html.in @@ -6,5 +6,5 @@ \ No newline at end of file diff --git a/doc/examples/class_FixedVectorBlock.cpp b/doc/examples/class_FixedVectorBlock.cpp new file mode 100644 index 000000000..c176be495 --- /dev/null +++ b/doc/examples/class_FixedVectorBlock.cpp @@ -0,0 +1,26 @@ +#include +USING_PART_OF_NAMESPACE_EIGEN +using namespace std; + +template +Eigen::VectorBlock +firstTwo(MatrixBase& v) +{ + return Eigen::VectorBlock(v.derived(), 0); +} + +template +const Eigen::VectorBlock +firstTwo(const MatrixBase& v) +{ + return Eigen::VectorBlock(v.derived(), 0); +} + +int main(int, char**) +{ + Matrix v; v << 1,2,3,4,5,6; + cout << firstTwo(4*v) << endl; // calls the const version + firstTwo(v) *= 2; // calls the non-const version + cout << "Now the vector v is:" << endl << v << endl; + return 0; +} diff --git a/doc/examples/class_VectorBlock.cpp b/doc/examples/class_VectorBlock.cpp new file mode 100644 index 000000000..d979f973a --- /dev/null +++ b/doc/examples/class_VectorBlock.cpp @@ -0,0 +1,26 @@ +#include +USING_PART_OF_NAMESPACE_EIGEN +using namespace std; + +template +Eigen::VectorBlock +segmentFromRange(MatrixBase& v, int start, int end) +{ + return Eigen::VectorBlock(v.derived(), start, end-start); +} + +template +const Eigen::VectorBlock +segmentFromRange(const MatrixBase& v, int start, int end) +{ + return Eigen::VectorBlock(v.derived(), start, end-start); +} + +int main(int, char**) +{ + Matrix v; v << 1,2,3,4,5,6; + cout << segmentFromRange(2*v, 2, 4) << endl; // calls the const version + segmentFromRange(v, 1, 3) *= 5; // calls the non-const version + cout << "Now the vector v is:" << endl << v << endl; + return 0; +} diff --git a/doc/snippets/QR_solve.cpp b/doc/snippets/HouseholderQR_solve.cpp similarity index 52% rename from doc/snippets/QR_solve.cpp rename to doc/snippets/HouseholderQR_solve.cpp index 7e629f851..aa9404951 100644 --- a/doc/snippets/QR_solve.cpp +++ b/doc/snippets/HouseholderQR_solve.cpp @@ -4,11 +4,6 @@ Matrix3f y = Matrix3f::Random(); cout << "Here is the matrix m:" << endl << m << endl; cout << "Here is the matrix y:" << endl << y << endl; Matrix3f x; -if(m.qr().solve(y, &x)) -{ - assert(y.isApprox(m*x)); - cout << "Here is a solution x to the equation mx=y:" << endl << x << endl; -} -else - cout << "The equation mx=y does not have any solution." << endl; - +m.householderQr().solve(y, &x)); +assert(y.isApprox(m*x)); +cout << "Here is a solution x to the equation mx=y:" << endl << x << endl; diff --git a/doc/snippets/MatrixBase_computeInverseWithCheck.cpp b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp new file mode 100644 index 000000000..19e24c90b --- /dev/null +++ b/doc/snippets/MatrixBase_computeInverseWithCheck.cpp @@ -0,0 +1,9 @@ +Matrix3d m = Matrix3d::Random(); +cout << "Here is the matrix m:" << endl << m << endl; +Matrix3d inv; +if(m.computeInverseWithCheck(&inv)) { + cout << "It is invertible, and its inverse is:" << endl << inv << endl; +} +else { + cout << "It is not invertible." << endl; +} diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 965e4d256..403f16a45 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -76,6 +76,7 @@ template void adjoint(const MatrixType& m) { VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast(1)); VERIFY_IS_APPROX(v1.norm(), v1.stableNorm()); + VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm()); } // check compatibility of dot and adjoint diff --git a/test/array.cpp b/test/array.cpp index 37deeaa4f..a308f7366 100644 --- a/test/array.cpp +++ b/test/array.cpp @@ -1,7 +1,7 @@ // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // -// Copyright (C) 2008 Gael Guennebaud +// Copyright (C) 2008-2009 Gael Guennebaud // // Eigen is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index 2f088988f..52aff93ee 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -107,7 +107,7 @@ template void basicStuff(const MatrixType& m) { VERIFY_IS_NOT_APPROX(m3, m1); } - + m3.real() = m1.real(); VERIFY_IS_APPROX(static_cast(m3).real(), static_cast(m1).real()); VERIFY_IS_APPROX(static_cast(m3).real(), m1.real()); @@ -121,16 +121,16 @@ template void basicStuffComplex(const MatrixType& m) int rows = m.rows(); int cols = m.cols(); - + Scalar s1 = ei_random(), s2 = ei_random(); - + VERIFY(ei_real(s1)==ei_real_ref(s1)); VERIFY(ei_imag(s1)==ei_imag_ref(s1)); ei_real_ref(s1) = ei_real(s2); ei_imag_ref(s1) = ei_imag(s2); VERIFY(s1==s2); - + RealMatrixType rm1 = RealMatrixType::Random(rows,cols), rm2 = RealMatrixType::Random(rows,cols); MatrixType cm(rows,cols); @@ -162,7 +162,7 @@ void test_basicstuff() CALL_SUBTEST( basicStuff(MatrixXcd(20, 20)) ); CALL_SUBTEST( basicStuff(Matrix()) ); CALL_SUBTEST( basicStuff(Matrix(10,10)) ); - + CALL_SUBTEST( basicStuffComplex(MatrixXcf(21, 17)) ); CALL_SUBTEST( basicStuffComplex(MatrixXcd(2, 3)) ); } diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp index c5d2715db..f8281a16b 100644 --- a/test/geo_hyperplane.cpp +++ b/test/geo_hyperplane.cpp @@ -64,7 +64,7 @@ template void hyperplane(const HyperplaneType& _plane) // transform if (!NumTraits::IsComplex) { - MatrixType rot = MatrixType::Random(dim,dim).qr().matrixQ(); + MatrixType rot = MatrixType::Random(dim,dim).householderQr().matrixQ(); DiagonalMatrix scaling(VectorType::Random()); Translation translation(VectorType::Random()); diff --git a/test/inverse.cpp b/test/inverse.cpp index 90b897197..5ac39e35a 100644 --- a/test/inverse.cpp +++ b/test/inverse.cpp @@ -65,6 +65,21 @@ template void inverse(const MatrixType& m) // since for the general case we implement separately row-major and col-major, test that VERIFY_IS_APPROX(m1.transpose().inverse(), m1.inverse().transpose()); + + //computeInverseWithCheck tests + //First: an invertible matrix + bool invertible = m1.computeInverseWithCheck(&m2); + VERIFY(invertible); + VERIFY_IS_APPROX(identity, m1*m2); + + //Second: a rank one matrix (not invertible, except for 1x1 matrices) + VectorType v3 = VectorType::Random(rows); + MatrixType m3 = v3*v3.transpose(), m4(rows,cols); + invertible = m3.computeInverseWithCheck( &m4 ); + if( 1 == rows ){ + VERIFY( invertible ); } + else{ + VERIFY( !invertible ); } } void test_inverse() @@ -77,7 +92,7 @@ void test_inverse() CALL_SUBTEST( inverse(MatrixXf(8,8)) ); CALL_SUBTEST( inverse(MatrixXcd(7,7)) ); } - + // test some tricky cases for 4x4 matrices VERIFY_IS_APPROX((Matrix4f() << 0,0,1,0, 1,0,0,0, 0,1,0,0, 0,0,0,1).finished().inverse(), (Matrix4f() << 0,1,0,0, 0,0,1,0, 1,0,0,0, 0,0,0,1).finished()); diff --git a/test/lu.cpp b/test/lu.cpp index 60e45ff2b..4ad92bb11 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -57,6 +57,11 @@ template void lu_non_invertible() VERIFY_IS_APPROX(m3, m1*m2); m3 = MatrixType::Random(rows,cols2); VERIFY(!lu.solve(m3, &m2)); + + typedef Matrix SquareMatrixType; + SquareMatrixType m4(rows, rows), m5(rows, rows); + createRandomMatrixOfRank(rows/2, rows, rows, m4); + VERIFY(!m4.computeInverseWithCheck(&m5)); } template void lu_invertible() diff --git a/test/main.h b/test/main.h index 76572d9b3..3947451cc 100644 --- a/test/main.h +++ b/test/main.h @@ -242,8 +242,8 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, Eigen::Matri const int diag_size = std::min(d.rows(),d.cols()); d.diagonal().segment(desired_rank, diag_size-desired_rank) = VectorType::Zero(diag_size-desired_rank); - QR qra(a); - QR qrb(b); + HouseholderQR qra(a); + HouseholderQR qrb(b); m = (qra.matrixQ() * d * qrb.matrixQ()).lazy(); } diff --git a/test/qr.cpp b/test/qr.cpp index 6e96f1e97..88a447c4b 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -36,7 +36,7 @@ template void qr(const MatrixType& m) typedef Matrix VectorType; MatrixType a = MatrixType::Random(rows,cols); - QR qrOfA(a); + HouseholderQR qrOfA(a); VERIFY_IS_APPROX(a, qrOfA.matrixQ() * qrOfA.matrixR()); VERIFY_IS_NOT_APPROX(a+MatrixType::Identity(rows, cols), qrOfA.matrixQ() * qrOfA.matrixR()); @@ -55,40 +55,6 @@ template void qr(const MatrixType& m) VERIFY_IS_APPROX(b, hess.matrixQ() * hess.matrixH() * hess.matrixQ().adjoint()); } -template void qr_non_invertible() -{ - /* this test covers the following files: QR.h */ - int rows = ei_random(20,200), cols = ei_random(20,rows), cols2 = ei_random(20,rows); - int rank = ei_random(1, std::min(rows, cols)-1); - - MatrixType m1(rows, cols), m2(cols, cols2), m3(rows, cols2), k(1,1); - createRandomMatrixOfRank(rank, rows, cols, m1); - - QR lu(m1); -// typename LU::KernelResultType m1kernel = lu.kernel(); -// typename LU::ImageResultType m1image = lu.image(); - std::cerr << rows << "x" << cols << " " << rank << " " << lu.rank() << "\n"; - if (rank != lu.rank()) - std::cerr << lu.matrixR().diagonal().transpose() << "\n"; - VERIFY(rank == lu.rank()); - VERIFY(cols - lu.rank() == lu.dimensionOfKernel()); - VERIFY(!lu.isInjective()); - VERIFY(!lu.isInvertible()); - VERIFY(lu.isSurjective() == (lu.rank() == rows)); -// VERIFY((m1 * m1kernel).isMuchSmallerThan(m1)); -// VERIFY(m1image.lu().rank() == rank); -// MatrixType sidebyside(m1.rows(), m1.cols() + m1image.cols()); -// sidebyside << m1, m1image; -// VERIFY(sidebyside.lu().rank() == rank); - m2 = MatrixType::Random(cols,cols2); - m3 = m1*m2; - m2 = MatrixType::Random(cols,cols2); - lu.solve(m3, &m2); - VERIFY_IS_APPROX(m3, m1*m2); - m3 = MatrixType::Random(rows,cols2); - VERIFY(!lu.solve(m3, &m2)); -} - template void qr_invertible() { /* this test covers the following files: QR.h */ @@ -105,33 +71,18 @@ template void qr_invertible() m1 += a * a.adjoint(); } - QR lu(m1); - VERIFY(0 == lu.dimensionOfKernel()); - VERIFY(size == lu.rank()); - VERIFY(lu.isInjective()); - VERIFY(lu.isSurjective()); - VERIFY(lu.isInvertible()); -// VERIFY(lu.image().lu().isInvertible()); + HouseholderQR qr(m1); m3 = MatrixType::Random(size,size); - lu.solve(m3, &m2); + qr.solve(m3, &m2); //std::cerr << m3 - m1*m2 << "\n\n"; VERIFY_IS_APPROX(m3, m1*m2); -// VERIFY_IS_APPROX(m2, lu.inverse()*m3); - m3 = MatrixType::Random(size,size); - VERIFY(lu.solve(m3, &m2)); } template void qr_verify_assert() { MatrixType tmp; - QR qr; - VERIFY_RAISES_ASSERT(qr.isFullRank()) - VERIFY_RAISES_ASSERT(qr.rank()) - VERIFY_RAISES_ASSERT(qr.dimensionOfKernel()) - VERIFY_RAISES_ASSERT(qr.isInjective()) - VERIFY_RAISES_ASSERT(qr.isSurjective()) - VERIFY_RAISES_ASSERT(qr.isInvertible()) + HouseholderQR qr; VERIFY_RAISES_ASSERT(qr.matrixR()) VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) VERIFY_RAISES_ASSERT(qr.matrixQ()) @@ -149,11 +100,6 @@ void test_qr() } for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( qr_non_invertible() ); - CALL_SUBTEST( qr_non_invertible() ); - // TODO fix issue with complex -// CALL_SUBTEST( qr_non_invertible() ); -// CALL_SUBTEST( qr_non_invertible() ); CALL_SUBTEST( qr_invertible() ); CALL_SUBTEST( qr_invertible() ); // TODO fix issue with complex diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 9b7d7c4e6..666eea872 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp @@ -299,6 +299,8 @@ template void sparse_basic(const SparseMatrixType& re initSparse(density, refMat2, m2); VERIFY_IS_APPROX(m2.transpose().eval(), refMat2.transpose().eval()); VERIFY_IS_APPROX(m2.transpose(), refMat2.transpose()); + + VERIFY_IS_APPROX(SparseMatrixType(m2.adjoint()), refMat2.adjoint()); } // test prune diff --git a/test/svd.cpp b/test/svd.cpp index 2eb7881a0..93e5ea017 100644 --- a/test/svd.cpp +++ b/test/svd.cpp @@ -50,7 +50,7 @@ template void svd(const MatrixType& m) MatrixType sigma = MatrixType::Zero(rows,cols); MatrixType matU = MatrixType::Zero(rows,rows); sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal(); - matU.block(0,0,rows,cols) = svd.matrixU(); + matU = svd.matrixU(); VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); }