diff options
| author |  Thomas Capricelli <orzel@freehackers.org> | 2009-09-10 02:37:08 +0200 |
|---|---|---|
| committer | Thomas Capricelli <orzel@freehackers.org> | 2009-09-10 02:37:08 +0200 |
| commit | 72746838ad381ea974143d92b30296c2597d26c3 (patch) | |
| tree | 9ed2642a8fa75a3a5e022503dc7d7a182d35a122 | |
| parent | e3db42611b79761cf7d45b94bcb503a66efd7e77 (diff) | |
| parent | 2a6db40f1064c71400bd0be35da440aa82591331 (diff) | |
merge with tip
91 files changed, 2590 insertions, 893 deletions
diff --git a/Eigen/CMakeLists.txt b/Eigen/CMakeLists.txt index b692cdc51..931cc6e20 100644 --- a/Eigen/CMakeLists.txt +++ b/Eigen/CMakeLists.txt @@ -1,4 +1,4 @@ -set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi) +set(Eigen_HEADERS Core LU Cholesky QR Geometry Sparse Array SVD LeastSquares QtAlignedMalloc StdVector Householder Jacobi Eigenvalues) if(EIGEN_BUILD_LIB) set(Eigen_SRCS diff --git a/Eigen/Eigenvalues b/Eigen/Eigenvalues new file mode 100644 index 000000000..9a6443f39 --- /dev/null +++ b/Eigen/Eigenvalues @@ -0,0 +1,74 @@ +#ifndef EIGEN_EIGENVALUES_MODULE_H +#define EIGEN_EIGENVALUES_MODULE_H + +#include "Core" + +#include "src/Core/util/DisableMSVCWarnings.h" + +#include "Cholesky" +#include "Jacobi" +#include "Householder" + +// Note that EIGEN_HIDE_HEAVY_CODE has to be defined per module +#if (defined EIGEN_EXTERN_INSTANTIATIONS) && (EIGEN_EXTERN_INSTANTIATIONS>=2) + #ifndef EIGEN_HIDE_HEAVY_CODE + #define EIGEN_HIDE_HEAVY_CODE + #endif +#elif defined EIGEN_HIDE_HEAVY_CODE + #undef EIGEN_HIDE_HEAVY_CODE +#endif + +namespace Eigen { + +/** \defgroup Eigenvalues_Module Eigenvalues module + * + * \nonstableyet + * + * This module mainly provides various eigenvalue solvers. + * This module also provides some MatrixBase methods, including: + * - MatrixBase::eigenvalues(), + * - MatrixBase::operatorNorm() + * + * \code + * #include <Eigen/Eigenvalues> + * \endcode + */ + +#include "src/Eigenvalues/Tridiagonalization.h" +#include "src/Eigenvalues/EigenSolver.h" +#include "src/Eigenvalues/SelfAdjointEigenSolver.h" +#include "src/Eigenvalues/HessenbergDecomposition.h" +#include "src/Eigenvalues/ComplexSchur.h" +#include "src/Eigenvalues/ComplexEigenSolver.h" + +// declare all classes for a given matrix type +#define EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \ + PREFIX template class Tridiagonalization<MATRIXTYPE>; \ + PREFIX template class HessenbergDecomposition<MATRIXTYPE>; \ + PREFIX template class SelfAdjointEigenSolver<MATRIXTYPE> + +// removed because it does not support complex yet +// PREFIX template class EigenSolver<MATRIXTYPE> + +// declare all class for all types +#define EIGEN_EIGENVALUES_MODULE_INSTANTIATE(PREFIX) \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix2f,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix2d,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix3f,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix3d,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix4f,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(Matrix4d,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXf,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXd,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXcf,PREFIX); \ + EIGEN_EIGENVALUES_MODULE_INSTANTIATE_TYPE(MatrixXcd,PREFIX) + +#ifdef EIGEN_EXTERN_INSTANTIATIONS + EIGEN_EIGENVALUES_MODULE_INSTANTIATE(extern); +#endif // EIGEN_EXTERN_INSTANTIATIONS + +} // namespace Eigen + +#include "src/Core/util/EnableMSVCWarnings.h" + +#endif // EIGEN_EIGENVALUES_MODULE_H diff --git a/Eigen/Jacobi b/Eigen/Jacobi index 33a6b757e..da67d2a6a 100644 --- a/Eigen/Jacobi +++ b/Eigen/Jacobi @@ -8,11 +8,15 @@ namespace Eigen { /** \defgroup Jacobi_Module Jacobi module - * This module provides Jacobi rotations. + * This module provides Jacobi and Givens rotations. * * \code * #include <Eigen/Jacobi> * \endcode + * + * In addition to listed classes, it defines the two following MatrixBase methods to apply a Jacobi or Givens rotation: + * - MatrixBase::applyOnTheLeft() + * - MatrixBase::applyOnTheRight(). */ #include "src/Jacobi/Jacobi.h" diff --git a/Eigen/LeastSquares b/Eigen/LeastSquares index 573a13cb4..75620a349 100644 --- a/Eigen/LeastSquares +++ b/Eigen/LeastSquares @@ -5,7 +5,7 @@ #include "src/Core/util/DisableMSVCWarnings.h" -#include "QR" +#include "Eigenvalues" #include "Geometry" namespace Eigen { @@ -24,11 +24,9 @@ namespace Eigen { * * \nonstableyet * - * This module mainly provides QR decomposition and an eigen value solver. + * This module provides various QR decompositions * This module also provides some MatrixBase methods, including: * - MatrixBase::qr(), - * - MatrixBase::eigenvalues(), - * - MatrixBase::operatorNorm() * * \code * #include <Eigen/QR> @@ -38,20 +36,10 @@ namespace Eigen { #include "src/QR/HouseholderQR.h" #include "src/QR/FullPivotingHouseholderQR.h" #include "src/QR/ColPivotingHouseholderQR.h" -#include "src/QR/Tridiagonalization.h" -#include "src/QR/EigenSolver.h" -#include "src/QR/SelfAdjointEigenSolver.h" -#include "src/QR/HessenbergDecomposition.h" // declare all classes for a given matrix type #define EIGEN_QR_MODULE_INSTANTIATE_TYPE(MATRIXTYPE,PREFIX) \ PREFIX template class HouseholderQR<MATRIXTYPE>; \ - PREFIX template class Tridiagonalization<MATRIXTYPE>; \ - PREFIX template class HessenbergDecomposition<MATRIXTYPE>; \ - PREFIX template class SelfAdjointEigenSolver<MATRIXTYPE> - -// removed because it does not support complex yet -// PREFIX template class EigenSolver<MATRIXTYPE> // declare all class for all types #define EIGEN_QR_MODULE_INSTANTIATE(PREFIX) \ @@ -74,4 +62,7 @@ namespace Eigen { #include "src/Core/util/EnableMSVCWarnings.h" +// FIXME for compatibility we include Eigenvalues here: +#include "Eigenvalues" + #endif // EIGEN_QR_MODULE_H diff --git a/Eigen/QtAlignedMalloc b/Eigen/QtAlignedMalloc index de4da752c..fce7edf9b 100644 --- a/Eigen/QtAlignedMalloc +++ b/Eigen/QtAlignedMalloc @@ -6,8 +6,10 @@ #if (!EIGEN_MALLOC_ALREADY_ALIGNED) +#include "src/Core/util/DisableMSVCWarnings.h" + void *qMalloc(size_t size) -{std::cerr << "ok\n"; +{ return Eigen::ei_aligned_malloc(size); } @@ -24,6 +26,8 @@ void *qRealloc(void *ptr, size_t size) return newPtr; } +#include "src/Core/util/EnableMSVCWarnings.h" + #endif #endif // EIGEN_QTMALLOC_MODULE_H @@ -1,7 +1,7 @@ #ifndef EIGEN_SVD_MODULE_H #define EIGEN_SVD_MODULE_H -#include "Core" +#include "QR" #include "Householder" #include "Jacobi" @@ -23,7 +23,7 @@ namespace Eigen { */ #include "src/SVD/SVD.h" -#include "src/SVD/JacobiSquareSVD.h" +#include "src/SVD/JacobiSVD.h" } // namespace Eigen diff --git a/Eigen/src/Array/BooleanRedux.h b/Eigen/src/Array/BooleanRedux.h index 6a5e208de..ab6e06d56 100644 --- a/Eigen/src/Array/BooleanRedux.h +++ b/Eigen/src/Array/BooleanRedux.h @@ -78,10 +78,8 @@ struct ei_any_unroller<Derived, Dynamic> }; /** \array_module - * - * \returns true if all coefficients are true * - * \addexample CwiseAll \label How to check whether a point is inside a box (using operator< and all()) + * \returns true if all coefficients are true * * Example: \include MatrixBase_all.cpp * Output: \verbinclude MatrixBase_all.out @@ -107,7 +105,7 @@ inline bool MatrixBase<Derived>::all() const } /** \array_module - * + * * \returns true if at least one coefficient is true * * \sa MatrixBase::all() @@ -131,7 +129,7 @@ inline bool MatrixBase<Derived>::any() const } /** \array_module - * + * * \returns the number of coefficients which evaluate to true * * \sa MatrixBase::all(), MatrixBase::any() diff --git a/Eigen/src/Array/Random.h b/Eigen/src/Array/Random.h index 15cc6ae7c..adadf99e3 100644 --- a/Eigen/src/Array/Random.h +++ b/Eigen/src/Array/Random.h @@ -34,7 +34,7 @@ struct ei_functor_traits<ei_scalar_random_op<Scalar> > { enum { Cost = 5 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = false }; }; /** \array_module - * + * * \returns a random matrix expression * * The parameters \a rows and \a cols are the number of rows and of columns of @@ -44,8 +44,6 @@ struct ei_functor_traits<ei_scalar_random_op<Scalar> > * it is redundant to pass \a rows and \a cols as arguments, so Random() should be used * instead. * - * \addexample RandomExample \label How to create a matrix with random coefficients - * * Example: \include MatrixBase_random_int_int.cpp * Output: \verbinclude MatrixBase_random_int_int.out * @@ -63,7 +61,7 @@ MatrixBase<Derived>::Random(int rows, int cols) } /** \array_module - * + * * \returns a random vector expression * * The parameter \a size is the size of the returned vector. @@ -92,7 +90,7 @@ MatrixBase<Derived>::Random(int size) } /** \array_module - * + * * \returns a fixed-size random matrix or vector expression * * This variant is only for fixed-size MatrixBase types. For dynamic-size types, you @@ -115,7 +113,7 @@ MatrixBase<Derived>::Random() } /** \array_module - * + * * Sets all coefficients in this expression to random values. * * Example: \include MatrixBase_setRandom.cpp diff --git a/Eigen/src/CMakeLists.txt b/Eigen/src/CMakeLists.txt index 2f6a83a1b..0df8273d1 100644 --- a/Eigen/src/CMakeLists.txt +++ b/Eigen/src/CMakeLists.txt @@ -9,3 +9,4 @@ ADD_SUBDIRECTORY(LeastSquares) ADD_SUBDIRECTORY(Sparse) ADD_SUBDIRECTORY(Jacobi) ADD_SUBDIRECTORY(Householder) +ADD_SUBDIRECTORY(Eigenvalues) diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index adca9fe2e..c8d92f3c0 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -218,8 +218,8 @@ LDLT<MatrixType>& LDLT<MatrixType>::compute(const MatrixType& a) int endSize = size - j - 1; if (endSize > 0) { - _temporary.end(endSize) = ( m_matrix.block(j+1,0, endSize, j) - * m_matrix.col(j).start(j).conjugate() ).lazy(); + _temporary.end(endSize).noalias() = m_matrix.block(j+1,0, endSize, j) + * m_matrix.col(j).start(j).conjugate(); m_matrix.row(j).end(endSize) = m_matrix.row(j).end(endSize).conjugate() - _temporary.end(endSize).transpose(); diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 3ce85b2d9..ec7b8123c 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -141,7 +141,7 @@ template<> struct ei_llt_inplace<LowerTriangular> if (x<=RealScalar(0)) return false; mat.coeffRef(k,k) = x = ei_sqrt(x); - if (k>0 && rs>0) A21 -= (A20 * A10.adjoint()).lazy(); + if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); if (rs>0) A21 *= RealScalar(1)/x; } return true; diff --git a/Eigen/src/Core/Block.h b/Eigen/src/Core/Block.h index cf7730170..cebfeaf75 100644 --- a/Eigen/src/Core/Block.h +++ b/Eigen/src/Core/Block.h @@ -201,6 +201,13 @@ template<typename MatrixType, int BlockRows, int BlockCols, int PacketAccess, in m_startCol.value() + (RowsAtCompileTime == 1 ? index : 0), x); } + #ifdef EIGEN_PARSED_BY_DOXYGEN + /** \sa MapBase::data() */ + inline const Scalar* data() const; + /** \sa MapBase::stride() */ + inline int stride() const; + #endif + protected: const typename MatrixType::Nested m_matrix; @@ -269,7 +276,14 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess> && startCol >= 0 && blockCols >= 0 && startCol + blockCols <= matrix.cols()); } - inline int stride(void) const { return m_matrix.stride(); } + /** \sa MapBase::stride() */ + inline int stride() const + { + return ((!Base::IsVectorAtCompileTime) + || (BlockRows==1 && ((Flags&RowMajorBit)==0)) + || (BlockCols==1 && ((Flags&RowMajorBit)==RowMajorBit))) + ? m_matrix.stride() : 1; + } #ifndef __SUNPRO_CC // FIXME sunstudio is not friendly with the above friend... @@ -294,8 +308,6 @@ class Block<MatrixType,BlockRows,BlockCols,PacketAccess,HasDirectAccess> * \param blockRows the number of rows in the block * \param blockCols the number of columns in the block * - * \addexample BlockIntIntIntInt \label How to reference a sub-matrix (dynamic-size) - * * Example: \include MatrixBase_block_int_int_int_int.cpp * Output: \verbinclude MatrixBase_block_int_int_int_int.out * @@ -327,8 +339,6 @@ inline const typename BlockReturnType<Derived>::Type MatrixBase<Derived> * \param cRows the number of rows in the corner * \param cCols the number of columns in the corner * - * \addexample BlockCornerDynamicSize \label How to reference a sub-corner of a matrix - * * Example: \include MatrixBase_corner_enum_int_int.cpp * Output: \verbinclude MatrixBase_corner_enum_int_int.out * @@ -438,8 +448,6 @@ MatrixBase<Derived>::corner(CornerType type) const * \param startRow the first row in the block * \param startCol the first column in the block * - * \addexample BlockSubMatrixFixedSize \label How to reference a sub-matrix (fixed-size) - * * Example: \include MatrixBase_block_int_int.cpp * Output: \verbinclude MatrixBase_block_int_int.out * @@ -467,8 +475,6 @@ MatrixBase<Derived>::block(int startRow, int startCol) const /** \returns an expression of the \a i-th column of *this. Note that the numbering starts at 0. * - * \addexample BlockColumn \label How to reference a single column of a matrix - * * Example: \include MatrixBase_col.cpp * Output: \verbinclude MatrixBase_col.out * @@ -490,8 +496,6 @@ MatrixBase<Derived>::col(int i) const /** \returns an expression of the \a i-th row of *this. Note that the numbering starts at 0. * - * \addexample BlockRow \label How to reference a single row of a matrix - * * Example: \include MatrixBase_row.cpp * Output: \verbinclude MatrixBase_row.out * diff --git a/Eigen/src/Core/CommaInitializer.h b/Eigen/src/Core/CommaInitializer.h index c0dd701f2..e86f47ad0 100644 --- a/Eigen/src/Core/CommaInitializer.h +++ b/Eigen/src/Core/CommaInitializer.h @@ -124,8 +124,6 @@ struct CommaInitializer * The coefficients must be provided in a row major order and exactly match * the size of the matrix. Otherwise an assertion is raised. * - * \addexample CommaInit \label How to easily set all the coefficients of a matrix - * * Example: \include MatrixBase_set.cpp * Output: \verbinclude MatrixBase_set.out * diff --git a/Eigen/src/Core/CwiseBinaryOp.h b/Eigen/src/Core/CwiseBinaryOp.h index 78de9fbce..875bc9aa5 100644 --- a/Eigen/src/Core/CwiseBinaryOp.h +++ b/Eigen/src/Core/CwiseBinaryOp.h @@ -291,8 +291,6 @@ Cwise<ExpressionType>::max(const MatrixBase<OtherDerived> &other) const * The template parameter \a CustomBinaryOp is the type of the functor * of the custom operator (see class CwiseBinaryOp for an example) * - * \addexample CustomCwiseBinaryFunctors \label How to use custom coeff wise binary functors - * * Here is an example illustrating the use of custom functors: * \include class_CwiseBinaryOp.cpp * Output: \verbinclude class_CwiseBinaryOp.out diff --git a/Eigen/src/Core/CwiseNullaryOp.h b/Eigen/src/Core/CwiseNullaryOp.h index 804bc2e74..61ce51885 100644 --- a/Eigen/src/Core/CwiseNullaryOp.h +++ b/Eigen/src/Core/CwiseNullaryOp.h @@ -317,8 +317,6 @@ Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setConstant(int row * it is redundant to pass \a rows and \a cols as arguments, so Zero() should be used * instead. * - * \addexample Zero \label How to take get a zero matrix - * * Example: \include MatrixBase_zero_int_int.cpp * Output: \verbinclude MatrixBase_zero_int_int.out * @@ -448,8 +446,6 @@ Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setZero(int rows, i * it is redundant to pass \a rows and \a cols as arguments, so Ones() should be used * instead. * - * \addexample One \label How to get a matrix with all coefficients equal one - * * Example: \include MatrixBase_ones_int_int.cpp * Output: \verbinclude MatrixBase_ones_int_int.out * @@ -576,8 +572,6 @@ Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::setOnes(int rows, i * it is redundant to pass \a rows and \a cols as arguments, so Identity() should be used * instead. * - * \addexample Identity \label How to get an identity matrix - * * Example: \include MatrixBase_identity_int_int.cpp * Output: \verbinclude MatrixBase_identity_int_int.out * diff --git a/Eigen/src/Core/CwiseUnaryOp.h b/Eigen/src/Core/CwiseUnaryOp.h index d701c25d9..03011800c 100644 --- a/Eigen/src/Core/CwiseUnaryOp.h +++ b/Eigen/src/Core/CwiseUnaryOp.h @@ -109,8 +109,6 @@ class CwiseUnaryOp : ei_no_assignment_operator, * The template parameter \a CustomUnaryOp is the type of the functor * of the custom unary operator. * - * \addexample CustomCwiseUnaryFunctors \label How to use custom coeff wise unary functors - * * Example: * \include class_CwiseUnaryOp.cpp * Output: \verbinclude class_CwiseUnaryOp.out @@ -234,7 +232,7 @@ Cwise<ExpressionType>::log() const } -/** \relates MatrixBase */ +/** \returns an expression of \c *this scaled by the scalar factor \a scalar */ template<typename Derived> EIGEN_STRONG_INLINE const typename MatrixBase<Derived>::ScalarMultipleReturnType MatrixBase<Derived>::operator*(const Scalar& scalar) const @@ -243,7 +241,17 @@ MatrixBase<Derived>::operator*(const Scalar& scalar) const (derived(), ei_scalar_multiple_op<Scalar>(scalar)); } -/** \relates MatrixBase */ +/** Overloaded for efficient real matrix times complex scalar value */ +template<typename Derived> +EIGEN_STRONG_INLINE const CwiseUnaryOp<ei_scalar_multiple2_op<typename ei_traits<Derived>::Scalar, + std::complex<typename ei_traits<Derived>::Scalar> >, Derived> +MatrixBase<Derived>::operator*(const std::complex<Scalar>& scalar) const +{ + return CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,std::complex<Scalar> >, Derived> + (*static_cast<const Derived*>(this), ei_scalar_multiple2_op<Scalar,std::complex<Scalar> >(scalar)); +} + +/** \returns an expression of \c *this divided by the scalar value \a scalar */ template<typename Derived> EIGEN_STRONG_INLINE const CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived> MatrixBase<Derived>::operator/(const Scalar& scalar) const diff --git a/Eigen/src/Core/CwiseUnaryView.h b/Eigen/src/Core/CwiseUnaryView.h index 4785f01f4..580344379 100644 --- a/Eigen/src/Core/CwiseUnaryView.h +++ b/Eigen/src/Core/CwiseUnaryView.h @@ -53,8 +53,7 @@ struct ei_traits<CwiseUnaryView<ViewOp, MatrixType> > }; template<typename ViewOp, typename MatrixType> -class CwiseUnaryView : ei_no_assignment_operator, - public MatrixBase<CwiseUnaryView<ViewOp, MatrixType> > +class CwiseUnaryView : public MatrixBase<CwiseUnaryView<ViewOp, MatrixType> > { public: @@ -99,8 +98,6 @@ class CwiseUnaryView : ei_no_assignment_operator, * The template parameter \a CustomUnaryOp is the type of the functor * of the custom unary operator. * - * \addexample CustomCwiseUnaryFunctors \label How to use custom coeff wise unary functors - * * Example: * \include class_CwiseUnaryOp.cpp * Output: \verbinclude class_CwiseUnaryOp.out diff --git a/Eigen/src/Core/DiagonalMatrix.h b/Eigen/src/Core/DiagonalMatrix.h index ebbed15d4..a74695921 100644 --- a/Eigen/src/Core/DiagonalMatrix.h +++ b/Eigen/src/Core/DiagonalMatrix.h @@ -32,7 +32,7 @@ class DiagonalBase : public AnyMatrixBase<Derived> public: typedef typename ei_traits<Derived>::DiagonalVectorType DiagonalVectorType; typedef typename DiagonalVectorType::Scalar Scalar; - + enum { RowsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, ColsAtCompileTime = DiagonalVectorType::SizeAtCompileTime, @@ -41,14 +41,14 @@ class DiagonalBase : public AnyMatrixBase<Derived> IsVectorAtCompileTime = 0, Flags = 0 }; - + typedef Matrix<Scalar, RowsAtCompileTime, ColsAtCompileTime, 0, MaxRowsAtCompileTime, MaxColsAtCompileTime> DenseMatrixType; - + #ifndef EIGEN_PARSED_BY_DOXYGEN inline const Derived& derived() const { return *static_cast<const Derived*>(this); } inline Derived& derived() { return *static_cast<Derived*>(this); } #endif // not EIGEN_PARSED_BY_DOXYGEN - + DenseMatrixType toDenseMatrix() const { return derived(); } template<typename DenseDerived> void evalToDense(MatrixBase<DenseDerived> &other) const; @@ -64,7 +64,7 @@ class DiagonalBase : public AnyMatrixBase<Derived> inline int rows() const { return diagonal().size(); } inline int cols() const { return diagonal().size(); } - + template<typename MatrixDerived> const DiagonalProduct<MatrixDerived, Derived, DiagonalOnTheLeft> operator*(const MatrixBase<MatrixDerived> &matrix) const; @@ -100,20 +100,20 @@ class DiagonalMatrix : public DiagonalBase<DiagonalMatrix<_Scalar,_Size> > { public: - + typedef typename ei_traits<DiagonalMatrix>::DiagonalVectorType DiagonalVectorType; typedef const DiagonalMatrix& Nested; typedef _Scalar Scalar; - + protected: DiagonalVectorType m_diagonal; - + public: inline const DiagonalVectorType& diagonal() const { return m_diagonal; } inline DiagonalVectorType& diagonal() { return m_diagonal; } - + /** Default constructor without initialization */ inline DiagonalMatrix() {} @@ -152,7 +152,7 @@ class DiagonalMatrix m_diagonal = other.m_diagonal(); return *this; } - + inline void resize(int size) { m_diagonal.resize(size); } inline void setZero() { m_diagonal.setZero(); } inline void setZero(int size) { m_diagonal.setZero(size); } @@ -194,10 +194,10 @@ class DiagonalWrapper public: typedef _DiagonalVectorType DiagonalVectorType; typedef DiagonalWrapper Nested; - + inline DiagonalWrapper(const DiagonalVectorType& diagonal) : m_diagonal(diagonal) {} const DiagonalVectorType& diagonal() const { return m_diagonal; } - + protected: const typename DiagonalVectorType::Nested m_diagonal; }; @@ -207,8 +207,6 @@ class DiagonalWrapper * * \only_for_vectors * - * \addexample AsDiagonalExample \label How to build a diagonal matrix from a vector - * * Example: \include MatrixBase_asDiagonal.cpp * Output: \verbinclude MatrixBase_asDiagonal.out * diff --git a/Eigen/src/Core/DiagonalProduct.h b/Eigen/src/Core/DiagonalProduct.h index b838d1b31..7eaa380eb 100644 --- a/Eigen/src/Core/DiagonalProduct.h +++ b/Eigen/src/Core/DiagonalProduct.h @@ -29,7 +29,7 @@ template<typename MatrixType, typename DiagonalType, int ProductOrder> struct ei_traits<DiagonalProduct<MatrixType, DiagonalType, ProductOrder> > { - typedef typename MatrixType::Scalar Scalar; + typedef typename ei_scalar_product_traits<typename MatrixType::Scalar, typename DiagonalType::Scalar>::ReturnType Scalar; enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, @@ -62,7 +62,7 @@ class DiagonalProduct : ei_no_assignment_operator, { return m_diagonal.diagonal().coeff(ProductOrder == DiagonalOnTheLeft ? row : col) * m_matrix.coeff(row, col); } - + template<int LoadMode> EIGEN_STRONG_INLINE PacketScalar packet(int row, int col) const { @@ -72,7 +72,7 @@ class DiagonalProduct : ei_no_assignment_operator, DiagonalVectorPacketLoadMode = (LoadMode == Aligned && ((InnerSize%16) == 0)) ? Aligned : Unaligned }; const int indexInDiagonalVector = ProductOrder == DiagonalOnTheLeft ? row : col; - + if((int(StorageOrder) == RowMajor && int(ProductOrder) == DiagonalOnTheLeft) ||(int(StorageOrder) == ColMajor && int(ProductOrder) == DiagonalOnTheRight)) { diff --git a/Eigen/src/Core/MapBase.h b/Eigen/src/Core/MapBase.h index 9cb085b61..88a3fac1e 100644 --- a/Eigen/src/Core/MapBase.h +++ b/Eigen/src/Core/MapBase.h @@ -62,7 +62,21 @@ template<typename Derived> class MapBase inline int rows() const { return m_rows.value(); } inline int cols() const { return m_cols.value(); } + /** Returns the leading dimension (for matrices) or the increment (for vectors) to be used with data(). + * + * More precisely: + * - for a column major matrix it returns the number of elements between two successive columns + * - for a row major matrix it returns the number of elements between two successive rows + * - for a vector it returns the number of elements between two successive coefficients + * This function has to be used together with the MapBase::data() function. + * + * \sa MapBase::data() */ inline int stride() const { return derived().stride(); } + + /** Returns a pointer to the first coefficient of the matrix or vector. + * This function has to be used together with the stride() function. + * + * \sa MapBase::stride() */ inline const Scalar* data() const { return m_data; } template<bool IsForceAligned,typename Dummy> struct force_aligned_impl { diff --git a/Eigen/src/Core/MathFunctions.h b/Eigen/src/Core/MathFunctions.h index 1570f01e0..40edf4a3c 100644 --- a/Eigen/src/Core/MathFunctions.h +++ b/Eigen/src/Core/MathFunctions.h @@ -116,6 +116,7 @@ inline float ei_imag(float) { return 0.f; } inline float ei_conj(float x) { return x; } inline float ei_abs(float x) { return std::abs(x); } inline float ei_abs2(float x) { return x*x; } +inline float ei_norm1(float x) { return ei_abs(x); } inline float ei_sqrt(float x) { return std::sqrt(x); } inline float ei_exp(float x) { return std::exp(x); } inline float ei_log(float x) { return std::log(x); } @@ -164,6 +165,7 @@ inline double ei_imag(double) { return 0.; } inline double ei_conj(double x) { return x; } inline double ei_abs(double x) { return std::abs(x); } inline double ei_abs2(double x) { return x*x; } +inline double ei_norm1(double x) { return ei_abs(x); } inline double ei_sqrt(double x) { return std::sqrt(x); } inline double ei_exp(double x) { return std::exp(x); } inline double ei_log(double x) { return std::log(x); } @@ -212,6 +214,7 @@ inline float& ei_imag_ref(std::complex<float>& x) { return reinterpret_cast<floa inline std::complex<float> ei_conj(const std::complex<float>& x) { return std::conj(x); } inline float ei_abs(const std::complex<float>& x) { return std::abs(x); } inline float ei_abs2(const std::complex<float>& x) { return std::norm(x); } +inline float ei_norm1(const std::complex<float> &x) { return(ei_abs(x.real()) + ei_abs(x.imag())); } inline std::complex<float> ei_exp(std::complex<float> x) { return std::exp(x); } inline std::complex<float> ei_sin(std::complex<float> x) { return std::sin(x); } inline std::complex<float> ei_cos(std::complex<float> x) { return std::cos(x); } @@ -248,6 +251,7 @@ inline double& ei_imag_ref(std::complex<double>& x) { return reinterpret_cast<do inline std::complex<double> ei_conj(const std::complex<double>& x) { return std::conj(x); } inline double ei_abs(const std::complex<double>& x) { return std::abs(x); } inline double ei_abs2(const std::complex<double>& x) { return std::norm(x); } +inline double ei_norm1(const std::complex<double> &x) { return(ei_abs(x.real()) + ei_abs(x.imag())); } inline std::complex<double> ei_exp(std::complex<double> x) { return std::exp(x); } inline std::complex<double> ei_sin(std::complex<double> x) { return std::sin(x); } inline std::complex<double> ei_cos(std::complex<double> x) { return std::cos(x); } diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index f58424ba2..c08f12491 100644 --- a/Eigen/src/Core/Matrix.h +++ b/Eigen/src/Core/Matrix.h @@ -25,6 +25,7 @@ #ifndef EIGEN_MATRIX_H #define EIGEN_MATRIX_H +template <typename Derived, typename OtherDerived, bool IsVector = static_cast<bool>(Derived::IsVectorAtCompileTime)> struct ei_conservative_resize_like_impl; /** \class Matrix * @@ -126,6 +127,7 @@ class Matrix EIGEN_GENERIC_PUBLIC_INTERFACE(Matrix) enum { Options = _Options }; + typedef typename Base::PlainMatrixType PlainMatrixType; friend class Eigen::Map<Matrix, Unaligned>; typedef class Eigen::Map<Matrix, Unaligned> UnalignedMapType; friend class Eigen::Map<Matrix, Aligned>; @@ -139,15 +141,27 @@ class Matrix && SizeAtCompileTime!=Dynamic && ((sizeof(Scalar)*SizeAtCompileTime)%16)==0 }; EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF(NeedsToAlign) + Base& base() { return *static_cast<Base*>(this); } + const Base& base() const { return *static_cast<const Base*>(this); } + EIGEN_STRONG_INLINE int rows() const { return m_storage.rows(); } EIGEN_STRONG_INLINE int cols() const { return m_storage.cols(); } - EIGEN_STRONG_INLINE int stride(void) const + /** Returns the leading dimension (for matrices) or the increment (for vectors) to be used with data(). + * + * More precisely: + * - for a column major matrix it returns the number of elements between two successive columns + * - for a row major matrix it returns the number of elements between two successive rows + * - for a vector it returns the number of elements between two successive coefficients + * This function has to be used together with the MapBase::data() function. + * + * \sa Matrix::data() */ + EIGEN_STRONG_INLINE int stride() const { - if(Flags & RowMajorBit) - return m_storage.cols(); + if(IsVectorAtCompileTime) + return 1; else - return m_storage.rows(); + return (Flags & RowMajorBit) ? m_storage.cols() : m_storage.rows(); } EIGEN_STRONG_INLINE const Scalar& coeff(int row, int col) const @@ -295,7 +309,7 @@ class Matrix */ template<typename OtherDerived> EIGEN_STRONG_INLINE void resizeLike(const MatrixBase<OtherDerived>& other) - { + { if(RowsAtCompileTime == 1) { ei_assert(other.isVector()); @@ -309,6 +323,51 @@ class Matrix else resize(other.rows(), other.cols()); } + /** Resizes \c *this to a \a rows x \a cols matrix while leaving old values of *this untouched. + * + * This method is intended for dynamic-size matrices. If you only want to change the number + * of rows and/or of columns, you can use conservativeResize(NoChange_t, int), + * conservativeResize(int, NoChange_t). + * + * The top-left part of the resized matrix will be the same as the overlapping top-left corner + * of *this. In case values need to be appended to the matrix they will be uninitialized. + */ + EIGEN_STRONG_INLINE void conservativeResize(int rows, int cols) + { + conservativeResizeLike(PlainMatrixType(rows, cols)); + } + + EIGEN_STRONG_INLINE void conservativeResize(int rows, NoChange_t) + { + // Note: see the comment in conservativeResize(int,int) + conservativeResize(rows, cols()); + } + + EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, int cols) + { + // Note: see the comment in conservativeResize(int,int) + conservativeResize(rows(), cols); + } + + /** Resizes \c *this to a vector of length \a size while retaining old values of *this. + * + * \only_for_vectors. This method does not work for + * partially dynamic matrices when the static dimension is anything other + * than 1. For example it will not work with Matrix<double, 2, Dynamic>. + * + * When values are appended, they will be uninitialized. + */ + EIGEN_STRONG_INLINE void conservativeResize(int size) + { + conservativeResizeLike(PlainMatrixType(size)); + } + + template<typename OtherDerived> + EIGEN_STRONG_INLINE void conservativeResizeLike(const MatrixBase<OtherDerived>& other) + { + ei_conservative_resize_like_impl<Matrix, OtherDerived>::run(*this, other); + } + /** Copies the value of the expression \a other into \c *this with automatic resizing. * * *this might be resized to match the dimensions of \a other. If *this was a null matrix (not already initialized), @@ -329,7 +388,15 @@ class Matrix */ EIGEN_STRONG_INLINE Matrix& operator=(const Matrix& other) { - return _set(other); + return _set(other); + } + + /** \sa MatrixBase::lazyAssign() */ + template<typename OtherDerived> + EIGEN_STRONG_INLINE Matrix& lazyAssign(const MatrixBase<OtherDerived>& other) + { + _resize_to_match(other); + return Base::lazyAssign(other.derived()); } template<typename OtherDerived,typename OtherEvalType> @@ -470,13 +537,8 @@ class Matrix /** Override MatrixBase::swap() since for dynamic-sized matrices of same type it is enough to swap the * data pointers. */ - inline void swap(Matrix& other) - { - if (Base::SizeAtCompileTime==Dynamic) - m_storage.swap(other.m_storage); - else - this->Base::swap(other); - } + template<typename OtherDerived> + void swap(const MatrixBase<OtherDerived>& other); /** \name Map * These are convenience functions returning Map objects. The Map() static functions return unaligned Map objects, @@ -635,8 +697,71 @@ class Matrix m_storage.data()[0] = x; m_storage.data()[1] = y; } + + template<typename MatrixType, typename OtherDerived, bool SwapPointers> + friend struct ei_matrix_swap_impl; +}; + +template <typename Derived, typename OtherDerived, bool IsVector> +struct ei_conservative_resize_like_impl +{ + static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other) + { + // Note: Here is space for improvement. Basically, for conservativeResize(int,int), + // neither RowsAtCompileTime or ColsAtCompileTime must be Dynamic. If only one of the + // dimensions is dynamic, one could use either conservativeResize(int rows, NoChange_t) or + // conservativeResize(NoChange_t, int cols). For these methods new static asserts like + // EIGEN_STATIC_ASSERT_DYNAMIC_ROWS and EIGEN_STATIC_ASSERT_DYNAMIC_COLS would be good. + EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Derived) + EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(OtherDerived) + + typename MatrixBase<Derived>::PlainMatrixType tmp(other); + const int common_rows = std::min(tmp.rows(), _this.rows()); + const int common_cols = std::min(tmp.cols(), _this.cols()); + tmp.block(0,0,common_rows,common_cols) = _this.block(0,0,common_rows,common_cols); + _this.derived().swap(tmp); + } }; +template <typename Derived, typename OtherDerived> +struct ei_conservative_resize_like_impl<Derived,OtherDerived,true> +{ + static void run(MatrixBase<Derived>& _this, const MatrixBase<OtherDerived>& other) + { + // segment(...) will check whether Derived/OtherDerived are vectors! + typename MatrixBase<Derived>::PlainMatrixType tmp(other); + const int common_size = std::min<int>(_this.size(),tmp.size()); + tmp.segment(0,common_size) = _this.segment(0,common_size); + _this.derived().swap(tmp); + } +}; + +template<typename MatrixType, typename OtherDerived, bool SwapPointers> +struct ei_matrix_swap_impl +{ + static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other) + { + matrix.base().swap(other); + } +}; + +template<typename MatrixType, typename OtherDerived> +struct ei_matrix_swap_impl<MatrixType, OtherDerived, true> +{ + static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other) + { + matrix.m_storage.swap(other.derived().m_storage); + } +}; + +template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols> +template<typename OtherDerived> +inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other) +{ + enum { SwapPointers = ei_is_same_type<Matrix, OtherDerived>::ret && Base::SizeAtCompileTime==Dynamic }; + ei_matrix_swap_impl<Matrix, OtherDerived, bool(SwapPointers)>::run(*this, *const_cast<MatrixBase<OtherDerived>*>(&other)); +} + /** \defgroup matrixtypedefs Global matrix typedefs * * \ingroup Core_Module diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 7cb8853e6..25a0545c6 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -91,9 +91,8 @@ template<typename Derived> struct AnyMatrixBase */ template<typename Derived> class MatrixBase #ifndef EIGEN_PARSED_BY_DOXYGEN - : public AnyMatrixBase<Derived>, - public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar, - typename NumTraits<typename ei_traits<Derived>::Scalar>::Real> + : public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar, + typename NumTraits<typename ei_traits<Derived>::Scalar>::Real> #endif // not EIGEN_PARSED_BY_DOXYGEN { public: @@ -419,10 +418,17 @@ template<typename Derived> class MatrixBase const CwiseUnaryOp<ei_scalar_quotient1_op<typename ei_traits<Derived>::Scalar>, Derived> operator/(const Scalar& scalar) const; - inline friend const CwiseUnaryOp<ei_scalar_multiple_op<typename ei_traits<Derived>::Scalar>, Derived> + const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,std::complex<Scalar> >, Derived> + operator*(const std::complex<Scalar>& scalar) const; + + inline friend const ScalarMultipleReturnType operator*(const Scalar& scalar, const MatrixBase& matrix) { return matrix*scalar; } + inline friend const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,std::complex<Scalar> >, Derived> + operator*(const std::complex<Scalar>& scalar, const MatrixBase& matrix) + { return matrix*scalar; } + template<typename OtherDerived> const typename ProductReturnType<Derived,OtherDerived>::Type operator*(const MatrixBase<OtherDerived> &other) const; @@ -803,11 +809,10 @@ template<typename Derived> class MatrixBase ///////// Jacobi module ///////// - void applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s); - void applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s); - bool makeJacobi(int p, int q, Scalar *c, Scalar *s) const; - bool makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const; - bool makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const; + template<typename OtherScalar> + void applyOnTheLeft(int p, int q, const PlanarRotation<OtherScalar>& j); + template<typename OtherScalar> + void applyOnTheRight(int p, int q, const PlanarRotation<OtherScalar>& j); #ifdef EIGEN_MATRIXBASE_PLUGIN #include EIGEN_MATRIXBASE_PLUGIN diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index 71203a362..e7227d4f6 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -63,7 +63,7 @@ template<typename Lhs, typename Rhs> struct ei_product_type Cols = Rhs::ColsAtCompileTime, Depth = EIGEN_ENUM_MIN(Lhs::ColsAtCompileTime,Rhs::RowsAtCompileTime) }; - + // the splitting into different lines of code here, introducing the _select enums and the typedef below, // is to work around an internal compiler error with gcc 4.1 and 4.2. private: @@ -73,7 +73,7 @@ private: depth_select = Depth>=EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ? Large : (Depth==1 ? 1 : Small) }; typedef ei_product_type_selector<rows_select, cols_select, depth_select> product_type_selector; - + public: enum { value = product_type_selector::ret @@ -84,18 +84,18 @@ public: * based on the three dimensions of the product. * This is a compile time mapping from {1,Small,Large}^3 -> {product types} */ // FIXME I'm not sure the current mapping is the ideal one. -template<int Rows, int Cols> struct ei_product_type_selector<Rows,Cols,1> { enum { ret = OuterProduct }; }; -template<int Depth> struct ei_product_type_selector<1,1,Depth> { enum { ret = InnerProduct }; }; -template<> struct ei_product_type_selector<1,1,1> { enum { ret = InnerProduct }; }; -template<> struct ei_product_type_selector<Small,1,Small> { enum { ret = UnrolledProduct }; }; -template<> struct ei_product_type_selector<1,Small,Small> { enum { ret = UnrolledProduct }; }; +template<int Rows, int Cols> struct ei_product_type_selector<Rows, Cols, 1> { enum { ret = OuterProduct }; }; +template<int Depth> struct ei_product_type_selector<1, 1, Depth> { enum { ret = InnerProduct }; }; +template<> struct ei_product_type_selector<1, 1, 1> { enum { ret = InnerProduct }; }; +template<> struct ei_product_type_selector<Small,1, Small> { enum { ret = UnrolledProduct }; }; +template<> struct ei_product_type_selector<1, Small,Small> { enum { ret = UnrolledProduct }; }; template<> struct ei_product_type_selector<Small,Small,Small> { enum { ret = UnrolledProduct }; }; -template<> struct ei_product_type_selector<1,Large,Small> { enum { ret = GemvProduct }; }; -template<> struct ei_product_type_selector<1,Large,Large> { enum { ret = GemvProduct }; }; -template<> struct ei_product_type_selector<1,Small,Large> { enum { ret = GemvProduct }; }; -template<> struct ei_product_type_selector<Large,1,Small> { enum { ret = GemvProduct }; }; -template<> struct ei_product_type_selector<Large,1,Large> { enum { ret = GemvProduct }; }; -template<> struct ei_product_type_selector<Small,1,Large> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<1, Large,Small> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<1, Large,Large> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<1, Small,Large> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<Large,1, Small> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<Large,1, Large> { enum { ret = GemvProduct }; }; +template<> struct ei_product_type_selector<Small,1, Large> { enum { ret = GemvProduct }; }; template<> struct ei_product_type_selector<Small,Small,Large> { enum { ret = GemmProduct }; }; template<> struct ei_product_type_selector<Large,Small,Large> { enum { ret = GemmProduct }; }; template<> struct ei_product_type_selector<Small,Large,Large> { enum { ret = GemmProduct }; }; @@ -144,7 +144,7 @@ struct ProductReturnType<Lhs,Rhs,UnrolledProduct> ***********************************************************************/ // FIXME : maybe the "inner product" could return a Scalar -// instead of a 1x1 matrix ?? +// instead of a 1x1 matrix ?? // Pro: more natural for the user // Cons: this could be a problem if in a meta unrolled algorithm a matrix-matrix // product ends up to a row-vector times col-vector product... To tackle this use @@ -162,7 +162,11 @@ class GeneralProduct<Lhs, Rhs, InnerProduct> public: EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) - GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} + GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + { + EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + } EIGEN_STRONG_INLINE Scalar value() const { @@ -197,11 +201,15 @@ class GeneralProduct<Lhs, Rhs, OuterProduct> public: EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) - GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} + GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + { + EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::RealScalar, typename Rhs::RealScalar>::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + } template<typename Dest> void scaleAndAddTo(Dest& dest, Scalar alpha) const { - ei_outer_product_selector<Dest::Flags&RowMajorBit>::run(*this, dest, alpha); + ei_outer_product_selector<(int(Dest::Flags)&RowMajorBit) ? RowMajor : ColMajor>::run(*this, dest, alpha); } }; @@ -209,6 +217,7 @@ template<> struct ei_outer_product_selector<ColMajor> { template<typename ProductType, typename Dest> EIGEN_DONT_INLINE static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) { // FIXME make sure lhs is sequentially stored + // FIXME not very good if rhs is real and lhs complex while alpha is real too const int cols = dest.cols(); for (int j=0; j<cols; ++j) dest.col(j) += (alpha * prod.rhs().coeff(j)) * prod.lhs(); @@ -219,6 +228,7 @@ template<> struct ei_outer_product_selector<RowMajor> { template<typename ProductType, typename Dest> EIGEN_DONT_INLINE static void run(const ProductType& prod, Dest& dest, typename ProductType::Scalar alpha) { // FIXME make sure rhs is sequentially stored + // FIXME not very good if lhs is real and rhs complex while alpha is real too const int rows = dest.rows(); for (int i=0; i<rows; ++i) dest.row(i) += (alpha * prod.lhs().coeff(i)) * prod.rhs(); @@ -251,7 +261,11 @@ class GeneralProduct<Lhs, Rhs, GemvProduct> public: EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) - GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} + GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + { + EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + } enum { Side = Lhs::IsVectorAtCompileTime ? OnTheLeft : OnTheRight }; typedef typename ei_meta_if<int(Side)==OnTheRight,_LhsNested,_RhsNested>::ret MatrixType; @@ -259,8 +273,8 @@ class GeneralProduct<Lhs, Rhs, GemvProduct> template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const { ei_assert(m_lhs.rows() == dst.rows() && m_rhs.cols() == dst.cols()); - ei_gemv_selector<Side,int(MatrixType::Flags)&RowMajorBit, - ei_blas_traits<MatrixType>::ActualAccess>::run(*this, dst, alpha); + ei_gemv_selector<Side,(int(MatrixType::Flags)&RowMajorBit) ? RowMajor : ColMajor, + bool(ei_blas_traits<MatrixType>::ActualAccess)>::run(*this, dst, alpha); } }; @@ -339,7 +353,7 @@ template<> struct ei_gemv_selector<OnTheRight,RowMajor,true> Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(prod.lhs()) * RhsBlasTraits::extractScalarFactor(prod.rhs()); - + enum { DirectlyUseRhs = ((ei_packet_traits<Scalar>::size==1) || (_ActualRhsType::Flags&ActualPacketAccessBit)) && (!(_ActualRhsType::Flags & RowMajorBit)) diff --git a/Eigen/src/Core/ProductBase.h b/Eigen/src/Core/ProductBase.h index b2c4cd989..764dc4d8e 100644 --- a/Eigen/src/Core/ProductBase.h +++ b/Eigen/src/Core/ProductBase.h @@ -33,7 +33,7 @@ struct ei_traits<ProductBase<Derived,_Lhs,_Rhs> > { typedef typename ei_cleantype<_Lhs>::type Lhs; typedef typename ei_cleantype<_Rhs>::type Rhs; - typedef typename ei_traits<Lhs>::Scalar Scalar; + typedef typename ei_scalar_product_traits<typename Lhs::Scalar, typename Rhs::Scalar>::ReturnType Scalar; enum { RowsAtCompileTime = ei_traits<Lhs>::RowsAtCompileTime, ColsAtCompileTime = ei_traits<Rhs>::ColsAtCompileTime, @@ -146,7 +146,7 @@ class ScaledProduct; // functions of ProductBase, because, otherwise we would have to // define all overloads defined in MatrixBase. Furthermore, Using // "using Base::operator*" would not work with MSVC. -// +// // Also note that here we accept any compatible scalar types template<typename Derived,typename Lhs,typename Rhs> const ScaledProduct<Derived> diff --git a/Eigen/src/Core/SolveTriangular.h b/Eigen/src/Core/SolveTriangular.h index 810b08240..c7f0cd227 100644 --- a/Eigen/src/Core/SolveTriangular.h +++ b/Eigen/src/Core/SolveTriangular.h @@ -30,7 +30,7 @@ template<typename Lhs, typename Rhs, int Side, // can be OnTheLeft/OnTheRight int Unrolling = Rhs::IsVectorAtCompileTime && Rhs::SizeAtCompileTime <= 8 // FIXME ? CompleteUnrolling : NoUnrolling, - int StorageOrder = int(Lhs::Flags) & RowMajorBit, + int StorageOrder = (int(Lhs::Flags) & RowMajorBit) ? RowMajor : ColMajor, int RhsCols = Rhs::ColsAtCompileTime > struct ei_triangular_solver_selector; @@ -157,7 +157,7 @@ struct ei_triangular_solver_selector<Lhs,Rhs,Side,Mode,NoUnrolling,StorageOrder, { const ActualLhsType actualLhs = LhsProductTraits::extract(lhs); ei_triangular_solve_matrix<Scalar,Side,Mode,LhsProductTraits::NeedToConjugate,StorageOrder, - Rhs::Flags&RowMajorBit> + (Rhs::Flags&RowMajorBit) ? RowMajor : ColMajor> ::run(lhs.rows(), Side==OnTheLeft? rhs.cols() : rhs.rows(), &actualLhs.coeff(0,0), actualLhs.stride(), &rhs.coeffRef(0,0), rhs.stride()); } }; diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h index 7ce5977f6..b291f7b1a 100644 --- a/Eigen/src/Core/VectorBlock.h +++ b/Eigen/src/Core/VectorBlock.h @@ -88,7 +88,7 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock using Base::operator-=; using Base::operator*=; using Base::operator/=; - + /** Dynamic-size constructor */ inline VectorBlock(const VectorType& vector, int start, int size) @@ -96,7 +96,7 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock IsColVector ? start : 0, IsColVector ? 0 : start, IsColVector ? size : 1, IsColVector ? 1 : size) { - + EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorBlock); } @@ -114,8 +114,6 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock * * \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 * @@ -151,8 +149,6 @@ MatrixBase<Derived>::segment(int start, int size) const * * \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 * @@ -185,8 +181,6 @@ MatrixBase<Derived>::start(int size) const * * \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 * @@ -251,8 +245,6 @@ MatrixBase<Derived>::segment(int start) const * * 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 * diff --git a/Eigen/src/Core/products/GeneralMatrixMatrix.h b/Eigen/src/Core/products/GeneralMatrixMatrix.h index 8b3b13266..beec17ee4 100644 --- a/Eigen/src/Core/products/GeneralMatrixMatrix.h +++ b/Eigen/src/Core/products/GeneralMatrixMatrix.h @@ -135,7 +135,11 @@ class GeneralProduct<Lhs, Rhs, GemmProduct> public: EIGEN_PRODUCT_PUBLIC_INTERFACE(GeneralProduct) - GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) {} + GeneralProduct(const Lhs& lhs, const Rhs& rhs) : Base(lhs,rhs) + { + EIGEN_STATIC_ASSERT((ei_is_same_type<typename Lhs::Scalar, typename Rhs::Scalar>::ret), + YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) + } template<typename Dest> void scaleAndAddTo(Dest& dst, Scalar alpha) const { @@ -149,9 +153,9 @@ class GeneralProduct<Lhs, Rhs, GemmProduct> ei_general_matrix_matrix_product< Scalar, - (_ActualLhsType::Flags&RowMajorBit)?RowMajor:ColMajor, bool(LhsBlasTraits::NeedToConjugate), - (_ActualRhsType::Flags&RowMajorBit)?RowMajor:ColMajor, bool(RhsBlasTraits::NeedToConjugate), - (Dest::Flags&RowMajorBit)?RowMajor:ColMajor> + (_ActualLhsType::Flags&RowMajorBit) ? RowMajor : ColMajor, bool(LhsBlasTraits::NeedToConjugate), + (_ActualRhsType::Flags&RowMajorBit) ? RowMajor : ColMajor, bool(RhsBlasTraits::NeedToConjugate), + (Dest::Flags&RowMajorBit) ? RowMajor : ColMajor> ::run( this->rows(), this->cols(), lhs.cols(), (const Scalar*)&(lhs.const_cast_derived().coeffRef(0,0)), lhs.stride(), diff --git a/Eigen/src/Core/products/SelfadjointMatrixVector.h b/Eigen/src/Core/products/SelfadjointMatrixVector.h index c2c33d5b8..c27454bee 100644 --- a/Eigen/src/Core/products/SelfadjointMatrixVector.h +++ b/Eigen/src/Core/products/SelfadjointMatrixVector.h @@ -41,7 +41,7 @@ static EIGEN_DONT_INLINE void ei_product_selfadjoint_vector( const int PacketSize = sizeof(Packet)/sizeof(Scalar); enum { - IsRowMajor = StorageOrder==RowMajorBit ? 1 : 0, + IsRowMajor = StorageOrder==RowMajor ? 1 : 0, IsLower = UpLo == LowerTriangularBit ? 1 : 0, FirstTriangular = IsRowMajor == IsLower }; @@ -185,14 +185,14 @@ struct SelfadjointProductMatrix<Lhs,LhsMode,false,Rhs,0,true> Scalar actualAlpha = alpha * LhsBlasTraits::extractScalarFactor(m_lhs) * RhsBlasTraits::extractScalarFactor(m_rhs); - ei_assert((&dst.coeff(1))-(&dst.coeff(0))==1 && "not implemented yet"); - ei_product_selfadjoint_vector<Scalar, ei_traits<_ActualLhsType>::Flags&RowMajorBit, int(LhsUpLo), bool(LhsBlasTraits::NeedToConjugate), bool(RhsBlasTraits::NeedToConjugate)> + ei_assert(dst.stride()==1 && "not implemented yet"); + ei_product_selfadjoint_vector<Scalar, (ei_traits<_ActualLhsType>::Flags&RowMajorBit) ? RowMajor : ColMajor, int(LhsUpLo), bool(LhsBlasTraits::NeedToConjugate), bool(RhsBlasTraits::NeedToConjugate)> ( - lhs.rows(), // size - &lhs.coeff(0,0), lhs.stride(), // lhs info - &rhs.coeff(0), (&rhs.coeff(1))-(&rhs.coeff(0)), // rhs info - &dst.coeffRef(0), // result info - actualAlpha // scale factor + lhs.rows(), // size + &lhs.coeff(0,0), lhs.stride(), // lhs info + &rhs.coeff(0), rhs.stride(), // rhs info + &dst.coeffRef(0), // result info + actualAlpha // scale factor ); } }; diff --git a/Eigen/src/Core/products/SelfadjointProduct.h b/Eigen/src/Core/products/SelfadjointProduct.h index d3fc962d9..f9cdda637 100644 --- a/Eigen/src/Core/products/SelfadjointProduct.h +++ b/Eigen/src/Core/products/SelfadjointProduct.h @@ -132,7 +132,7 @@ SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived()); - enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit)?1:0 }; + enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 }; ei_selfadjoint_product<Scalar, _ActualUType::Flags&RowMajorBit ? RowMajor : ColMajor, diff --git a/Eigen/src/Core/products/SelfadjointRank2Update.h b/Eigen/src/Core/products/SelfadjointRank2Update.h index f7dcc003f..64fcf8660 100644 --- a/Eigen/src/Core/products/SelfadjointRank2Update.h +++ b/Eigen/src/Core/products/SelfadjointRank2Update.h @@ -85,7 +85,7 @@ SelfAdjointView<MatrixType,UpLo>& SelfAdjointView<MatrixType,UpLo> Scalar actualAlpha = alpha * UBlasTraits::extractScalarFactor(u.derived()) * VBlasTraits::extractScalarFactor(v.derived()); - enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit)?1:0 }; + enum { IsRowMajor = (ei_traits<MatrixType>::Flags&RowMajorBit) ? 1 : 0 }; ei_selfadjoint_rank2_update_selector<Scalar, typename ei_cleantype<typename ei_conj_expr_if<IsRowMajor ^ UBlasTraits::NeedToConjugate,_ActualUType>::ret>::type, typename ei_cleantype<typename ei_conj_expr_if<IsRowMajor ^ VBlasTraits::NeedToConjugate,_ActualVType>::ret>::type, diff --git a/Eigen/src/Core/products/TriangularMatrixVector.h b/Eigen/src/Core/products/TriangularMatrixVector.h index 620b090b9..291177445 100644 --- a/Eigen/src/Core/products/TriangularMatrixVector.h +++ b/Eigen/src/Core/products/TriangularMatrixVector.h @@ -145,7 +145,7 @@ struct TriangularProduct<Mode,true,Lhs,false,Rhs,true> Mode, LhsBlasTraits::NeedToConjugate, RhsBlasTraits::NeedToConjugate, - ei_traits<Lhs>::Flags&RowMajorBit> + (int(ei_traits<Lhs>::Flags)&RowMajorBit) ? RowMajor : ColMajor> ::run(lhs,rhs,dst,actualAlpha); } }; diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 21e7f62c3..affc1d478 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -238,6 +238,15 @@ enum { OnTheRight = 2 }; +// options for SVD decomposition +enum { + SkipU = 0x1, + SkipV = 0x2, + AtLeastAsManyRowsAsCols = 0x4, + AtLeastAsManyColsAsRows = 0x8, + Square = AtLeastAsManyRowsAsCols | AtLeastAsManyColsAsRows +}; + /* the following could as well be written: * enum NoChange_t { NoChange }; * but it feels dangerous to disambiguate overloaded functions on enum/integer types. diff --git a/Eigen/src/Core/util/DisableMSVCWarnings.h b/Eigen/src/Core/util/DisableMSVCWarnings.h index 765ddecc5..c08d0389f 100644 --- a/Eigen/src/Core/util/DisableMSVCWarnings.h +++ b/Eigen/src/Core/util/DisableMSVCWarnings.h @@ -1,5 +1,6 @@ #ifdef _MSC_VER + // 4273 - QtAlignedMalloc, inconsistent dll linkage #pragma warning( push ) - #pragma warning( disable : 4181 4244 4127 4211 4717 ) + #pragma warning( disable : 4181 4244 4127 4211 4273 4522 4717 ) #endif diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index 414aaceca..c5f27d80b 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -120,8 +120,10 @@ template<typename MatrixType> class HouseholderQR; template<typename MatrixType> class ColPivotingHouseholderQR; template<typename MatrixType> class FullPivotingHouseholderQR; template<typename MatrixType> class SVD; +template<typename MatrixType, unsigned int Options = 0> class JacobiSVD; template<typename MatrixType, int UpLo = LowerTriangular> class LLT; template<typename MatrixType> class LDLT; +template<typename Scalar> class PlanarRotation; // Geometry module: template<typename Derived, int _Dim> class RotationBase; diff --git a/Eigen/src/Core/util/Macros.h b/Eigen/src/Core/util/Macros.h index dc3b3ee0a..ec8337e33 100644 --- a/Eigen/src/Core/util/Macros.h +++ b/Eigen/src/Core/util/Macros.h @@ -246,6 +246,14 @@ using Eigen::ei_cos; // needed to define it here as escaping characters in CMake add_definition's argument seems very problematic. #define EIGEN_DOCS_IO_FORMAT IOFormat(3, 0, " ", "\n", "", "") +#ifdef _MSC_VER +#define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \ +using Base::operator =; \ +using Base::operator +=; \ +using Base::operator -=; \ +using Base::operator *=; \ +using Base::operator /=; +#else #define EIGEN_INHERIT_ASSIGNMENT_OPERATORS(Derived) \ using Base::operator =; \ using Base::operator +=; \ @@ -256,6 +264,7 @@ EIGEN_STRONG_INLINE Derived& operator=(const Derived& other) \ { \ return Base::operator=(other); \ } +#endif #define _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, BaseClass) \ typedef BaseClass Base; \ @@ -278,6 +287,9 @@ enum { RowsAtCompileTime = Eigen::ei_traits<Derived>::RowsAtCompileTime, \ _EIGEN_GENERIC_PUBLIC_INTERFACE(Derived, Eigen::MatrixBase<Derived>) #define EIGEN_ENUM_MIN(a,b) (((int)a <= (int)b) ? (int)a : (int)b) +#define EIGEN_SIZE_MIN(a,b) (((int)a == 1 || (int)b == 1) ? 1 \ + : ((int)a == Dynamic || (int)b == Dynamic) ? Dynamic \ + : ((int)a <= (int)b) ? (int)a : (int)b) #define EIGEN_ENUM_MAX(a,b) (((int)a >= (int)b) ? (int)a : (int)b) #define EIGEN_LOGICAL_XOR(a,b) (((a) || (b)) && !((a) && (b))) diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h index e1ccd58e7..73d302fda 100644 --- a/Eigen/src/Core/util/StaticAssert.h +++ b/Eigen/src/Core/util/StaticAssert.h @@ -62,6 +62,7 @@ THIS_METHOD_IS_ONLY_FOR_MATRICES_OF_A_SPECIFIC_SIZE, YOU_MADE_A_PROGRAMMING_MISTAKE, YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR, + YOU_CALLED_A_DYNAMIC_SIZE_METHOD_ON_A_FIXED_SIZE_MATRIX_OR_VECTOR, UNALIGNED_LOAD_AND_STORE_OPERATIONS_UNIMPLEMENTED_ON_ALTIVEC, NUMERIC_TYPE_MUST_BE_FLOATING_POINT, NUMERIC_TYPE_MUST_BE_REAL, @@ -114,6 +115,11 @@ EIGEN_STATIC_ASSERT(TYPE::SizeAtCompileTime!=Eigen::Dynamic, \ YOU_CALLED_A_FIXED_SIZE_METHOD_ON_A_DYNAMIC_SIZE_MATRIX_OR_VECTOR) +// static assertion failing if the type \a TYPE is not dynamic-size +#define EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(TYPE) \ + EIGEN_STATIC_ASSERT(TYPE::SizeAtCompileTime==Eigen::Dynamic, \ + YOU_CALLED_A_DYNAMIC_SIZE_METHOD_ON_A_FIXED_SIZE_MATRIX_OR_VECTOR) + // static assertion failing if the type \a TYPE is not a vector type of the given size #define EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(TYPE, SIZE) \ EIGEN_STATIC_ASSERT(TYPE::IsVectorAtCompileTime && TYPE::SizeAtCompileTime==SIZE, \ diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 871259b08..cea2faaa8 100644 --- a/Eigen/src/Core/util/XprHelper.h +++ b/Eigen/src/Core/util/XprHelper.h @@ -217,7 +217,7 @@ template<unsigned int Flags> struct ei_are_flags_consistent * overloads for complex types */ template<typename Derived,typename Scalar,typename OtherScalar, bool EnableIt = !ei_is_same_type<Scalar,OtherScalar>::ret > -struct ei_special_scalar_op_base +struct ei_special_scalar_op_base : public AnyMatrixBase<Derived> { // dummy operator* so that the // "using ei_special_scalar_op_base::operator*" compiles @@ -225,7 +225,7 @@ struct ei_special_scalar_op_base }; template<typename Derived,typename Scalar,typename OtherScalar> -struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true> +struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true> : public AnyMatrixBase<Derived> { const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,OtherScalar>, Derived> operator*(const OtherScalar& scalar) const @@ -233,6 +233,10 @@ struct ei_special_scalar_op_base<Derived,Scalar,OtherScalar,true> return CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,OtherScalar>, Derived> (*static_cast<const Derived*>(this), ei_scalar_multiple2_op<Scalar,OtherScalar>(scalar)); } + + inline friend const CwiseUnaryOp<ei_scalar_multiple2_op<Scalar,OtherScalar>, Derived> + operator*(const OtherScalar& scalar, const Derived& matrix) + { return matrix*scalar; } }; /** \internal Gives the type of a sub-matrix or sub-vector of a matrix of type \a ExpressionType and size \a Size diff --git a/Eigen/src/Eigenvalues/CMakeLists.txt b/Eigen/src/Eigenvalues/CMakeLists.txt new file mode 100644 index 000000000..193e02685 --- /dev/null +++ b/Eigen/src/Eigenvalues/CMakeLists.txt @@ -0,0 +1,6 @@ +FILE(GLOB Eigen_EIGENVALUES_SRCS "*.h") + +INSTALL(FILES + ${Eigen_EIGENVALUES_SRCS} + DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Eigenvalues COMPONENT Devel + ) diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h new file mode 100644 index 000000000..666381949 --- /dev/null +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -0,0 +1,148 @@ +// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Claire Maurice
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
+#define EIGEN_COMPLEX_EIGEN_SOLVER_H
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+ * \nonstableyet
+ *
+ * \class ComplexEigenSolver
+ *
+ * \brief Eigen values/vectors solver for general complex matrices
+ *
+ * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
+ *
+ * \sa class EigenSolver, class SelfAdjointEigenSolver
+ */
+template<typename _MatrixType> class ComplexEigenSolver
+{
+ public:
+ typedef _MatrixType MatrixType;
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ typedef std::complex<RealScalar> Complex;
+ typedef Matrix<Complex, MatrixType::ColsAtCompileTime,1> EigenvalueType;
+ typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> EigenvectorType;
+
+ /**
+ * \brief Default Constructor.
+ *
+ * The default constructor is useful in cases in which the user intends to
+ * perform decompositions via ComplexEigenSolver::compute(const MatrixType&).
+ */
+ ComplexEigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false)
+ {}
+
+ ComplexEigenSolver(const MatrixType& matrix)
+ : m_eivec(matrix.rows(),matrix.cols()),
+ m_eivalues(matrix.cols()),
+ m_isInitialized(false)
+ {
+ compute(matrix);
+ }
+
+ EigenvectorType eigenvectors(void) const
+ {
+ ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+ return m_eivec;
+ }
+
+ EigenvalueType eigenvalues() const
+ {
+ ei_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
+ return m_eivalues;
+ }
+
+ void compute(const MatrixType& matrix);
+
+ protected:
+ MatrixType m_eivec;
+ EigenvalueType m_eivalues;
+ bool m_isInitialized;
+};
+
+
+template<typename MatrixType>
+void ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix)
+{
+ // this code is inspired from Jampack
+ assert(matrix.cols() == matrix.rows());
+ int n = matrix.cols();
+ m_eivalues.resize(n,1);
+
+ RealScalar eps = epsilon<RealScalar>();
+
+ // Reduce to complex Schur form
+ ComplexSchur<MatrixType> schur(matrix);
+
+ m_eivalues = schur.matrixT().diagonal();
+
+ m_eivec.setZero();
+
+ Scalar d2, z;
+ RealScalar norm = matrix.norm();
+
+ // compute the (normalized) eigenvectors
+ for(int k=n-1 ; k>=0 ; k--)
+ {
+ d2 = schur.matrixT().coeff(k,k);
+ m_eivec.coeffRef(k,k) = Scalar(1.0,0.0);
+ for(int i=k-1 ; i>=0 ; i--)
+ {
+ m_eivec.coeffRef(i,k) = -schur.matrixT().coeff(i,k);
+ if(k-i-1>0)
+ m_eivec.coeffRef(i,k) -= (schur.matrixT().row(i).segment(i+1,k-i-1) * m_eivec.col(k).segment(i+1,k-i-1)).value();
+ z = schur.matrixT().coeff(i,i) - d2;
+ if(z==Scalar(0))
+ ei_real_ref(z) = eps * norm;
+ m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) / z;
+
+ }
+ m_eivec.col(k).normalize();
+ }
+
+ m_eivec = schur.matrixU() * m_eivec;
+ m_isInitialized = true;
+
+ // sort the eigenvalues
+ {
+ for (int i=0; i<n; i++)
+ {
+ int k;
+ m_eivalues.cwise().abs().end(n-i).minCoeff(&k);
+ if (k != 0)
+ {
+ k += i;
+ std::swap(m_eivalues[k],m_eivalues[i]);
+ m_eivec.col(i).swap(m_eivec.col(k));
+ }
+ }
+ }
+}
+
+
+
+#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h new file mode 100644 index 000000000..58e2ea440 --- /dev/null +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -0,0 +1,237 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Claire Maurice +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_COMPLEX_SCHUR_H +#define EIGEN_COMPLEX_SCHUR_H + +/** \eigenvalues_module \ingroup Eigenvalues_Module + * \nonstableyet + * + * \class ComplexShur + * + * \brief Performs a complex Shur decomposition of a real or complex square matrix + * + */ +template<typename _MatrixType> class ComplexSchur +{ + public: + typedef _MatrixType MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef std::complex<RealScalar> Complex; + typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType; + + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via ComplexSchur::compute(const MatrixType&). + */ + ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false) + {} + + ComplexSchur(const MatrixType& matrix) + : m_matT(matrix.rows(),matrix.cols()), + m_matU(matrix.rows(),matrix.cols()), + m_isInitialized(false) + { + compute(matrix); + } + + ComplexMatrixType matrixU() const + { + ei_assert(m_isInitialized && "ComplexSchur is not initialized."); + return m_matU; + } + + ComplexMatrixType matrixT() const + { + ei_assert(m_isInitialized && "ComplexShur is not initialized."); + return m_matT; + } + + void compute(const MatrixType& matrix); + + protected: + ComplexMatrixType m_matT, m_matU; + bool m_isInitialized; +}; + +/** Computes the principal value of the square root of the complex \a z. */ +template<typename RealScalar> +std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z) +{ + RealScalar t, tre, tim; + + t = ei_abs(z); + + if (ei_abs(ei_real(z)) <= ei_abs(ei_imag(z))) + { + // No cancellation in these formulas + tre = ei_sqrt(0.5*(t + ei_real(z))); + tim = ei_sqrt(0.5*(t - ei_real(z))); + } + else + { + // Stable computation of the above formulas + if (z.real() > 0) + { + tre = t + z.real(); + tim = ei_abs(ei_imag(z))*ei_sqrt(0.5/tre); + tre = ei_sqrt(0.5*tre); + } + else + { + tim = t - z.real(); + tre = ei_abs(ei_imag(z))*ei_sqrt(0.5/tim); + tim = ei_sqrt(0.5*tim); + } + } + if(z.imag() < 0) + tim = -tim; + + return (std::complex<RealScalar>(tre,tim)); + +} + +template<typename MatrixType> +void ComplexSchur<MatrixType>::compute(const MatrixType& matrix) +{ + // this code is inspired from Jampack + assert(matrix.cols() == matrix.rows()); + int n = matrix.cols(); + + // Reduce to Hessenberg form + HessenbergDecomposition<MatrixType> hess(matrix); + + m_matT = hess.matrixH(); + m_matU = hess.matrixQ(); + + int iu = m_matT.cols() - 1; + int il; + RealScalar d,sd,sf; + Complex c,b,disc,r1,r2,kappa; + + RealScalar eps = epsilon<RealScalar>(); + + int iter = 0; + while(true) + { + //locate the range in which to iterate + while(iu > 0) + { + d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1)); + sd = ei_norm1(m_matT.coeffRef(iu,iu-1)); + + if(sd >= eps * d) break; // FIXME : precision criterion ?? + + m_matT.coeffRef(iu,iu-1) = Complex(0); + iter = 0; + --iu; + } + if(iu==0) break; + iter++; + + if(iter >= 30) + { + // FIXME : what to do when iter==MAXITER ?? + std::cerr << "MAXITER" << std::endl; + return; + } + + il = iu-1; + while( il > 0 ) + { + // check if the current 2x2 block on the diagonal is upper triangular + d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1)); + sd = ei_norm1(m_matT.coeffRef(il,il-1)); + + if(sd < eps * d) break; // FIXME : precision criterion ?? + + --il; + } + + if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0); + + // compute the shift (the normalization by sf is to avoid under/overflow) + Matrix<Scalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); + sf = t.cwise().abs().sum(); + t /= sf; + + c = t.determinant(); + b = t.diagonal().sum(); + + disc = ei_sqrt(b*b - RealScalar(4)*c); + + r1 = (b+disc)/RealScalar(2); + r2 = (b-disc)/RealScalar(2); + + if(ei_norm1(r1) > ei_norm1(r2)) + r2 = c/r1; + else + r1 = c/r2; + + if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1))) + kappa = sf * r1; + else + kappa = sf * r2; + + // perform the QR step using Givens rotations + PlanarRotation<Complex> rot; + rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il)); + + for(int i=il ; i<iu ; i++) + { + m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint()); + m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot); + m_matU.applyOnTheRight(i, i+1, rot); + + if(i != iu-1) + { + int i1 = i+1; + int i2 = i+2; + + rot.makeGivens(m_matT.coeffRef(i1,i), m_matT.coeffRef(i2,i), &m_matT.coeffRef(i1,i)); + m_matT.coeffRef(i2,i) = Complex(0); + } + } + } + + // FIXME : is it necessary ? + /* + for(int i=0 ; i<n ; i++) + for(int j=0 ; j<n ; j++) + { + if(ei_abs(ei_real(m_matT.coeff(i,j))) < eps) + ei_real_ref(m_matT.coeffRef(i,j)) = 0; + if(ei_imag(ei_abs(m_matT.coeff(i,j))) < eps) + ei_imag_ref(m_matT.coeffRef(i,j)) = 0; + } + */ + + m_isInitialized = true; +} + +#endif // EIGEN_COMPLEX_SCHUR_H diff --git a/Eigen/src/QR/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index bd7a76c49..3fc36c080 100644 --- a/Eigen/src/QR/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -25,7 +25,7 @@ #ifndef EIGEN_EIGENSOLVER_H #define EIGEN_EIGENSOLVER_H -/** \ingroup QR_Module +/** \eigenvalues_module \ingroup Eigenvalues_Module * \nonstableyet * * \class EigenSolver @@ -53,7 +53,7 @@ template<typename _MatrixType> class EigenSolver typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType; typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX; - /** + /** * \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to @@ -103,19 +103,19 @@ template<typename _MatrixType> class EigenSolver * * \sa pseudoEigenvalueMatrix() */ - const MatrixType& pseudoEigenvectors() const - { + const MatrixType& pseudoEigenvectors() const + { ei_assert(m_isInitialized && "EigenSolver is not initialized."); - return m_eivec; + return m_eivec; } MatrixType pseudoEigenvalueMatrix() const; /** \returns the eigenvalues as a column vector */ - EigenvalueType eigenvalues() const - { + EigenvalueType eigenvalues() const + { ei_assert(m_isInitialized && "EigenSolver is not initialized."); - return m_eivalues; + return m_eivalues; } EigenSolver& compute(const MatrixType& matrix); @@ -244,11 +244,11 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort) // Apply Householder similarity transformation // H = (I-u*u'/h)*H*(I-u*u')/h) int bSize = high-m+1; - matH.block(m, m, bSize, n-m) -= ((ort.segment(m, bSize)/h) - * (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m))).lazy(); + matH.block(m, m, bSize, n-m).noalias() -= ((ort.segment(m, bSize)/h) + * (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m))); - matH.block(0, m, high+1, bSize) -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize)) - * (ort.segment(m, bSize)/h).transpose()).lazy(); + matH.block(0, m, high+1, bSize).noalias() -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize)) + * (ort.segment(m, bSize)/h).transpose()); ort.coeffRef(m) = scale*ort.coeff(m); matH.coeffRef(m,m-1) = scale*g; @@ -265,8 +265,8 @@ void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort) ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m); int bSize = high-m+1; - m_eivec.block(m, m, bSize, bSize) += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m))) - * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) ).lazy(); + m_eivec.block(m, m, bSize, bSize).noalias() += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m))) + * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)) ); } } } diff --git a/Eigen/src/QR/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index 6b261a8a7..b1e21d4ee 100644 --- a/Eigen/src/QR/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -25,7 +25,7 @@ #ifndef EIGEN_HESSENBERGDECOMPOSITION_H #define EIGEN_HESSENBERGDECOMPOSITION_H -/** \ingroup QR_Module +/** \eigenvalues_module \ingroup Eigenvalues_Module * \nonstableyet * * \class HessenbergDecomposition @@ -156,7 +156,7 @@ void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVector // Apply similarity transformation to remaining columns, // i.e., compute A = H A H' - + // A = H A matA.corner(BottomRight, remainingSize, remainingSize) .applyHouseholderOnTheLeft(matA.col(i).end(remainingSize-1), h, &temp.coeffRef(0)); diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 809a717b9..84856aa66 100644 --- a/Eigen/src/QR/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -25,7 +25,7 @@ #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H #define EIGEN_SELFADJOINTEIGENSOLVER_H -/** \qr_module \ingroup QR_Module +/** \eigenvalues_module \ingroup Eigenvalues_Module * \nonstableyet * * \class SelfAdjointEigenSolver @@ -135,31 +135,9 @@ template<typename _MatrixType> class SelfAdjointEigenSolver #ifndef EIGEN_HIDE_HEAVY_CODE -// from Golub's "Matrix Computations", algorithm 5.1.3 -template<typename Scalar> -static void ei_givens_rotation(Scalar a, Scalar b, Scalar& c, Scalar& s) -{ - if (b==0) - { - c = 1; s = 0; - } - else if (ei_abs(b)>ei_abs(a)) - { - Scalar t = -a/b; - s = Scalar(1)/ei_sqrt(1+t*t); - c = s * t; - } - else - { - Scalar t = -b/a; - c = Scalar(1)/ei_sqrt(1+t*t); - s = c * t; - } -} - /** \internal * - * \qr_module + * \eigenvalues_module \ingroup Eigenvalues_Module * * Performs a QR step on a tridiagonal symmetric matrix represented as a * pair of two vectors \a diag and \a subdiag. @@ -288,7 +266,7 @@ compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors #endif // EIGEN_HIDE_HEAVY_CODE -/** \qr_module +/** \eigenvalues_module * * \returns a vector listing the eigenvalues of this matrix. */ @@ -329,7 +307,7 @@ template<typename Derived> struct ei_operatorNorm_selector<Derived, false> } }; -/** \qr_module +/** \eigenvalues_module * * \returns the matrix norm of this matrix. */ @@ -353,34 +331,33 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, int st for (int k = start; k < end; ++k) { - RealScalar c, s; - ei_givens_rotation(x, z, c, s); + PlanarRotation<RealScalar> rot; + rot.makeGivens(x, z); // do T = G' T G - RealScalar sdk = s * diag[k] + c * subdiag[k]; - RealScalar dkp1 = s * subdiag[k] + c * diag[k+1]; + RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k]; + RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1]; - diag[k] = c * (c * diag[k] - s * subdiag[k]) - s * (c * subdiag[k] - s * diag[k+1]); - diag[k+1] = s * sdk + c * dkp1; - subdiag[k] = c * sdk - s * dkp1; + diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]); + diag[k+1] = rot.s() * sdk + rot.c() * dkp1; + subdiag[k] = rot.c() * sdk - rot.s() * dkp1; if (k > start) - subdiag[k - 1] = c * subdiag[k-1] - s * z; + subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; x = subdiag[k]; if (k < end - 1) { - z = -s * subdiag[k+1]; - subdiag[k + 1] = c * subdiag[k+1]; + z = -rot.s() * subdiag[k+1]; + subdiag[k + 1] = rot.c() * subdiag[k+1]; } // apply the givens rotation to the unit matrix Q = Q * G - // G only modifies the two columns k and k+1 if (matrixQ) { Map<Matrix<Scalar,Dynamic,Dynamic> > q(matrixQ,n,n); - q.applyJacobiOnTheRight(k,k+1,c,s); + q.applyOnTheRight(k,k+1,rot); } } } diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index 2053505f6..5f891bfa6 100644 --- a/Eigen/src/QR/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -25,7 +25,7 @@ #ifndef EIGEN_TRIDIAGONALIZATION_H #define EIGEN_TRIDIAGONALIZATION_H -/** \ingroup QR_Module +/** \eigenvalues_module \ingroup Eigenvalues_Module * \nonstableyet * * \class Tridiagonalization diff --git a/Eigen/src/Geometry/AngleAxis.h b/Eigen/src/Geometry/AngleAxis.h index 9fd239354..b9dfa6972 100644 --- a/Eigen/src/Geometry/AngleAxis.h +++ b/Eigen/src/Geometry/AngleAxis.h @@ -39,8 +39,6 @@ * \li \c AngleAxisf for \c float * \li \c AngleAxisd for \c double * - * \addexample AngleAxisForEuler \label How to define a rotation from Euler-angles - * * Combined with MatrixBase::Unit{X,Y,Z}, AngleAxis can be used to easily * mimic Euler-angles. Here is an example: * \include AngleAxis_mimic_euler.cpp @@ -85,7 +83,7 @@ public: AngleAxis() {} /** Constructs and initialize the angle-axis rotation from an \a angle in radian * and an \a axis which \b must \b be \b normalized. - * + * * \warning If the \a axis vector is not normalized, then the angle-axis object * represents an invalid rotation. */ template<typename Derived> diff --git a/Eigen/src/Geometry/Umeyama.h b/Eigen/src/Geometry/Umeyama.h index 64c06fe66..7652aa92e 100644 --- a/Eigen/src/Geometry/Umeyama.h +++ b/Eigen/src/Geometry/Umeyama.h @@ -144,7 +144,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo // const Scalar dst_var = dst_demean.rowwise().squaredNorm().sum() * one_over_n; // Eq. (38) - const MatrixType sigma = (one_over_n * dst_demean * src_demean.transpose()).lazy(); + const MatrixType sigma = one_over_n * dst_demean * src_demean.transpose(); SVD<MatrixType> svd(sigma); @@ -160,14 +160,14 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo int rank = 0; for (int i=0; i<m; ++i) if (!ei_isMuchSmallerThan(d.coeff(i),d.coeff(0))) ++rank; if (rank == m-1) { if ( svd.matrixU().determinant() * svd.matrixV().determinant() > 0 ) { - Rt.block(0,0,m,m) = (svd.matrixU()*svd.matrixV().transpose()).lazy(); + Rt.block(0,0,m,m).noalias() = svd.matrixU()*svd.matrixV().transpose(); } else { const Scalar s = S(m-1); S(m-1) = -1; - Rt.block(0,0,m,m) = (svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose()).lazy(); + Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose(); S(m-1) = s; } } else { - Rt.block(0,0,m,m) = (svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose()).lazy(); + Rt.block(0,0,m,m).noalias() = svd.matrixU() * S.asDiagonal() * svd.matrixV().transpose(); } // Eq. (42) @@ -177,7 +177,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo // Note that we first assign dst_mean to the destination so that there no need // for a temporary. Rt.col(m).start(m) = dst_mean; - Rt.col(m).start(m) -= (c*Rt.corner(TopLeft,m,m)*src_mean).lazy(); + Rt.col(m).start(m).noalias() -= c*Rt.corner(TopLeft,m,m)*src_mean; if (with_scaling) Rt.block(0,0,m,m) *= c; diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 76f0800fe..eeb81c178 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -26,108 +26,282 @@ #ifndef EIGEN_JACOBI_H #define EIGEN_JACOBI_H -/** Applies the counter clock wise 2D rotation of angle \c theta given by its - * cosine \a c and sine \a s to the set of 2D vectors of cordinates \a x and \a y: - * \f$ x = c x - s' y \f$ - * \f$ y = s x + c y \f$ +/** \ingroup Jacobi_Module + * \jacobi_module + * \class PlanarRotation + * \brief Represents a rotation in the plane from a cosine-sine pair. * - * \sa MatrixBase::applyJacobiOnTheLeft(), MatrixBase::applyJacobiOnTheRight() + * This class represents a Jacobi or Givens rotation. + * This is a 2D rotation in the plane \c J of angle \f$ \theta \f$ defined by + * its cosine \c c and sine \c s as follow: + * \f$ J = \left ( \begin{array}{cc} c & \overline s \\ -s & \overline c \end{array} \right ) \f$ + * + * You can apply the respective counter-clockwise rotation to a column vector \c v by + * applying its adjoint on the left: \f$ v = J^* v \f$ that translates to the following Eigen code: + * \code + * v.applyOnTheLeft(J.adjoint()); + * \endcode + * + * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ -template<typename VectorX, typename VectorY> -void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s); +template<typename Scalar> class PlanarRotation +{ + public: + typedef typename NumTraits<Scalar>::Real RealScalar; + + /** Default constructor without any initialization. */ + PlanarRotation() {} + + /** Construct a planar rotation from a cosine-sine pair (\a c, \c s). */ + PlanarRotation(const Scalar& c, const Scalar& s) : m_c(c), m_s(s) {} + + Scalar& c() { return m_c; } + Scalar c() const { return m_c; } + Scalar& s() { return m_s; } + Scalar s() const { return m_s; } + + /** Concatenates two planar rotation */ + PlanarRotation operator*(const PlanarRotation& other) + { + return PlanarRotation(m_c * other.m_c - ei_conj(m_s) * other.m_s, + ei_conj(m_c * ei_conj(other.m_s) + ei_conj(m_s) * ei_conj(other.m_c))); + } + + /** Returns the transposed transformation */ + PlanarRotation transpose() const { return PlanarRotation(m_c, -ei_conj(m_s)); } + + /** Returns the adjoint transformation */ + PlanarRotation adjoint() const { return PlanarRotation(ei_conj(m_c), -m_s); } + + template<typename Derived> + bool makeJacobi(const MatrixBase<Derived>&, int p, int q); + bool makeJacobi(RealScalar x, Scalar y, RealScalar z); -/** Applies a rotation in the plane defined by \a c, \a s to the rows \a p and \a q of \c *this. - * More precisely, it computes B = J' * B, with J = [c s ; -s' c] and B = [ *this.row(p) ; *this.row(q) ] - * \sa MatrixBase::applyJacobiOnTheRight(), ei_apply_rotation_in_the_plane() + void makeGivens(const Scalar& p, const Scalar& q, Scalar* z=0); + + protected: + void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, ei_meta_true); + void makeGivens(const Scalar& p, const Scalar& q, Scalar* z, ei_meta_false); + + Scalar m_c, m_s; +}; + +/** Makes \c *this as a Jacobi rotation \a J such that applying \a J on both the right and left sides of the selfadjoint 2x2 matrix + * \f$ B = \left ( \begin{array}{cc} x & y \\ \overline y & z \end{array} \right )\f$ yields a diagonal matrix \f$ A = J^* B J \f$ + * + * \sa MatrixBase::makeJacobi(const MatrixBase<Derived>&, int, int), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ -template<typename Derived> -inline void MatrixBase<Derived>::applyJacobiOnTheLeft(int p, int q, Scalar c, Scalar s) +template<typename Scalar> +bool PlanarRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z) { - RowXpr x(row(p)); - RowXpr y(row(q)); - ei_apply_rotation_in_the_plane(x, y, ei_conj(c), ei_conj(s)); + typedef typename NumTraits<Scalar>::Real RealScalar; + if(y == Scalar(0)) + { + m_c = Scalar(1); + m_s = Scalar(0); + return false; + } + else + { + RealScalar tau = (x-z)/(RealScalar(2)*ei_abs(y)); + RealScalar w = ei_sqrt(ei_abs2(tau) + 1); + RealScalar t; + if(tau>0) + { + t = RealScalar(1) / (tau + w); + } + else + { + t = RealScalar(1) / (tau - w); + } + RealScalar sign_t = t > 0 ? 1 : -1; + RealScalar n = RealScalar(1) / ei_sqrt(ei_abs2(t)+1); + m_s = - sign_t * (ei_conj(y) / ei_abs(y)) * ei_abs(t) * n; + m_c = n; + return true; + } } -/** Applies a rotation in the plane defined by \a c, \a s to the columns \a p and \a q of \c *this. - * More precisely, it computes B = B * J, with J = [c s ; -s' c] and B = [ *this.col(p) ; *this.col(q) ] - * \sa MatrixBase::applyJacobiOnTheLeft(), ei_apply_rotation_in_the_plane() +/** Makes \c *this as a Jacobi rotation \c J such that applying \a J on both the right and left sides of the 2x2 selfadjoint matrix + * \f$ B = \left ( \begin{array}{cc} \text{this}_{pp} & \text{this}_{pq} \\ (\text{this}_{pq})^* & \text{this}_{qq} \end{array} \right )\f$ yields + * a diagonal matrix \f$ A = J^* B J \f$ + * + * Example: \include Jacobi_makeJacobi.cpp + * Output: \verbinclude Jacobi_makeJacobi.out + * + * \sa PlanarRotation::makeJacobi(RealScalar, Scalar, RealScalar), MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ +template<typename Scalar> template<typename Derived> -inline void MatrixBase<Derived>::applyJacobiOnTheRight(int p, int q, Scalar c, Scalar s) +inline bool PlanarRotation<Scalar>::makeJacobi(const MatrixBase<Derived>& m, int p, int q) { - ColXpr x(col(p)); - ColXpr y(col(q)); - ei_apply_rotation_in_the_plane(x, y, c, s); + return makeJacobi(ei_real(m.coeff(p,p)), m.coeff(p,q), ei_real(m.coeff(q,q))); } -/** Computes the cosine-sine pair (\a c, \a s) such that its associated - * rotation \f$ J = ( \begin{array}{cc} c & s \\ -s' c \end{array} )\f$ - * applied to both the right and left of the 2x2 matrix - * \f$ B = ( \begin{array}{cc} x & y \\ * & z \end{array} )\f$ yields - * a diagonal matrix A: \f$ A = J' B J \f$ +/** Makes \c *this as a Givens rotation \c G such that applying \f$ G^* \f$ to the left of the vector + * \f$ V = \left ( \begin{array}{c} p \\ q \end{array} \right )\f$ yields: + * \f$ G^* V = \left ( \begin{array}{c} r \\ 0 \end{array} \right )\f$. + * + * The value of \a z is returned if \a z is not null (the default is null). + * Also note that G is built such that the cosine is always real. + * + * Example: \include Jacobi_makeGivens.cpp + * Output: \verbinclude Jacobi_makeGivens.out + * + * This function implements the continuous Givens rotation generation algorithm + * found in Anderson (2000), Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem. + * LAPACK Working Note 150, University of Tennessee, UT-CS-00-454, December 4, 2000. + * + * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() */ template<typename Scalar> -bool ei_makeJacobi(Scalar x, Scalar y, Scalar z, Scalar *c, Scalar *s) +void PlanarRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* z) { - if(y == 0) + makeGivens(p, q, z, typename ei_meta_if<NumTraits<Scalar>::IsComplex, ei_meta_true, ei_meta_false>::ret()); +} + + +// specialization for complexes +template<typename Scalar> +void PlanarRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, ei_meta_true) +{ + if(q==Scalar(0)) { - *c = Scalar(1); - *s = Scalar(0); - return false; + m_c = ei_real(p)<0 ? Scalar(-1) : Scalar(1); + m_s = 0; + if(r) *r = m_c * p; + } + else if(p==Scalar(0)) + { + m_c = 0; + m_s = -q/ei_abs(q); + if(r) *r = ei_abs(q); } else { - Scalar tau = (z - x) / (2 * y); - Scalar w = ei_sqrt(1 + ei_abs2(tau)); - Scalar t; - if(tau>0) - t = Scalar(1) / (tau + w); + RealScalar p1 = ei_norm1(p); + RealScalar q1 = ei_norm1(q); + if(p1>=q1) + { + Scalar ps = p / p1; + RealScalar p2 = ei_abs2(ps); + Scalar qs = q / p1; + RealScalar q2 = ei_abs2(qs); + + RealScalar u = ei_sqrt(RealScalar(1) + q2/p2); + if(ei_real(p)<RealScalar(0)) + u = -u; + + m_c = Scalar(1)/u; + m_s = -qs*ei_conj(ps)*(m_c/p2); + if(r) *r = p * u; + } else - t = Scalar(1) / (tau - w); - *c = Scalar(1) / ei_sqrt(1 + ei_abs2(t)); - *s = *c * t; - return true; + { + Scalar ps = p / q1; + RealScalar p2 = ei_abs2(ps); + Scalar qs = q / q1; + RealScalar q2 = ei_abs2(qs); + + RealScalar u = q1 * ei_sqrt(p2 + q2); + if(ei_real(p)<RealScalar(0)) + u = -u; + + p1 = ei_abs(p); + ps = p/p1; + m_c = p1/u; + m_s = -ei_conj(ps) * (q/u); + if(r) *r = ps * u; + } } } -template<typename Derived> -inline bool MatrixBase<Derived>::makeJacobi(int p, int q, Scalar *c, Scalar *s) const +// specialization for reals +template<typename Scalar> +void PlanarRotation<Scalar>::makeGivens(const Scalar& p, const Scalar& q, Scalar* r, ei_meta_false) { - return ei_makeJacobi(coeff(p,p), coeff(p,q), coeff(q,q), c, s); + + if(q==0) + { + m_c = p<Scalar(0) ? Scalar(-1) : Scalar(1); + m_s = 0; + if(r) *r = ei_abs(p); + } + else if(p==0) + { + m_c = 0; + m_s = q<Scalar(0) ? Scalar(1) : Scalar(-1); + if(r) *r = ei_abs(q); + } + else if(ei_abs(p) > ei_abs(q)) + { + Scalar t = q/p; + Scalar u = ei_sqrt(Scalar(1) + ei_abs2(t)); + if(p<Scalar(0)) + u = -u; + m_c = Scalar(1)/u; + m_s = -t * m_c; + if(r) *r = p * u; + } + else + { + Scalar t = p/q; + Scalar u = ei_sqrt(Scalar(1) + ei_abs2(t)); + if(q<Scalar(0)) + u = -u; + m_s = -Scalar(1)/u; + m_c = -t * m_s; + if(r) *r = q * u; + } + } +/**************************************************************************************** +* Implementation of MatrixBase methods +****************************************************************************************/ + +/** \jacobi_module + * Applies the clock wise 2D rotation \a j to the set of 2D vectors of cordinates \a x and \a y: + * \f$ \left ( \begin{array}{cc} x \\ y \end{array} \right ) = J \left ( \begin{array}{cc} x \\ y \end{array} \right ) \f$ + * + * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() + */ +template<typename VectorX, typename VectorY, typename OtherScalar> +void ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const PlanarRotation<OtherScalar>& j); + +/** \jacobi_module + * Applies the rotation in the plane \a j to the rows \a p and \a q of \c *this, i.e., it computes B = J * B, + * with \f$ B = \left ( \begin{array}{cc} \text{*this.row}(p) \\ \text{*this.row}(q) \end{array} \right ) \f$. + * + * \sa class PlanarRotation, MatrixBase::applyOnTheRight(), ei_apply_rotation_in_the_plane() + */ template<typename Derived> -inline bool MatrixBase<Derived>::makeJacobiForAtA(int p, int q, Scalar *c, Scalar *s) const +template<typename OtherScalar> +inline void MatrixBase<Derived>::applyOnTheLeft(int p, int q, const PlanarRotation<OtherScalar>& j) { - return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(q,p)), - ei_conj(coeff(p,p))*coeff(p,q) + ei_conj(coeff(q,p))*coeff(q,q), - ei_abs2(coeff(p,q)) + ei_abs2(coeff(q,q)), - c,s); + RowXpr x(row(p)); + RowXpr y(row(q)); + ei_apply_rotation_in_the_plane(x, y, j); } +/** \ingroup Jacobi_Module + * Applies the rotation in the plane \a j to the columns \a p and \a q of \c *this, i.e., it computes B = B * J + * with \f$ B = \left ( \begin{array}{cc} \text{*this.col}(p) & \text{*this.col}(q) \end{array} \right ) \f$. + * + * \sa class PlanarRotation, MatrixBase::applyOnTheLeft(), ei_apply_rotation_in_the_plane() + */ template<typename Derived> -inline bool MatrixBase<Derived>::makeJacobiForAAt(int p, int q, Scalar *c, Scalar *s) const +template<typename OtherScalar> +inline void MatrixBase<Derived>::applyOnTheRight(int p, int q, const PlanarRotation<OtherScalar>& j) { - return ei_makeJacobi(ei_abs2(coeff(p,p)) + ei_abs2(coeff(p,q)), - ei_conj(coeff(q,p))*coeff(p,p) + ei_conj(coeff(q,q))*coeff(p,q), - ei_abs2(coeff(q,p)) + ei_abs2(coeff(q,q)), - c,s); + ColXpr x(col(p)); + ColXpr y(col(q)); + ei_apply_rotation_in_the_plane(x, y, j.transpose()); } -template<typename Scalar> -inline void ei_normalizeJacobi(Scalar *c, Scalar *s, const Scalar& x, const Scalar& y) -{ - Scalar a = x * *c - y * *s; - Scalar b = x * *s + y * *c; - if(ei_abs(b)>ei_abs(a)) { - Scalar x = *c; - *c = -*s; - *s = x; - } -} -template<typename VectorX, typename VectorY> -void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, typename VectorX::Scalar c, typename VectorY::Scalar s) +template<typename VectorX, typename VectorY, typename OtherScalar> +void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& _y, const PlanarRotation<OtherScalar>& j) { typedef typename VectorX::Scalar Scalar; ei_assert(_x.size() == _y.size()); @@ -138,7 +312,7 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& Scalar* EIGEN_RESTRICT x = &_x.coeffRef(0); Scalar* EIGEN_RESTRICT y = &_y.coeffRef(0); - if (incrx==1 && incry==1) + if((VectorX::Flags & VectorY::Flags & PacketAccessBit) && incrx==1 && incry==1) { // both vectors are sequentially stored in memory => vectorization typedef typename ei_packet_traits<Scalar>::type Packet; @@ -147,16 +321,16 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& int alignedStart = ei_alignmentOffset(y, size); int alignedEnd = alignedStart + ((size-alignedStart)/PacketSize)*PacketSize; - const Packet pc = ei_pset1(c); - const Packet ps = ei_pset1(s); + const Packet pc = ei_pset1(Scalar(j.c())); + const Packet ps = ei_pset1(Scalar(j.s())); ei_conj_helper<NumTraits<Scalar>::IsComplex,false> cj; for(int i=0; i<alignedStart; ++i) { Scalar xi = x[i]; Scalar yi = y[i]; - x[i] = c * xi - ei_conj(s) * yi; - y[i] = s * xi + c * yi; + x[i] = j.c() * xi + ei_conj(j.s()) * yi; + y[i] = -j.s() * xi + ei_conj(j.c()) * yi; } Scalar* px = x + alignedStart; @@ -168,8 +342,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Packet xi = ei_pload(px); Packet yi = ei_pload(py); - ei_pstore(px, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi))); - ei_pstore(py, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi))); + ei_pstore(px, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi))); + ei_pstore(py, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi))); px += PacketSize; py += PacketSize; } @@ -183,10 +357,10 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& Packet xi1 = ei_ploadu(px+PacketSize); Packet yi = ei_pload (py); Packet yi1 = ei_pload (py+PacketSize); - ei_pstoreu(px, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi))); - ei_pstoreu(px+PacketSize, ei_psub(ei_pmul(pc,xi1),cj.pmul(ps,yi1))); - ei_pstore (py, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi))); - ei_pstore (py+PacketSize, ei_padd(ei_pmul(ps,xi1),ei_pmul(pc,yi1))); + ei_pstoreu(px, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi))); + ei_pstoreu(px+PacketSize, ei_padd(ei_pmul(pc,xi1),cj.pmul(ps,yi1))); + ei_pstore (py, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi))); + ei_pstore (py+PacketSize, ei_psub(ei_pmul(pc,yi1),ei_pmul(ps,xi1))); px += Peeling*PacketSize; py += Peeling*PacketSize; } @@ -194,8 +368,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Packet xi = ei_ploadu(x+peelingEnd); Packet yi = ei_pload (y+peelingEnd); - ei_pstoreu(x+peelingEnd, ei_psub(ei_pmul(pc,xi),cj.pmul(ps,yi))); - ei_pstore (y+peelingEnd, ei_padd(ei_pmul(ps,xi),ei_pmul(pc,yi))); + ei_pstoreu(x+peelingEnd, ei_padd(ei_pmul(pc,xi),cj.pmul(ps,yi))); + ei_pstore (y+peelingEnd, ei_psub(ei_pmul(pc,yi),ei_pmul(ps,xi))); } } @@ -203,8 +377,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Scalar xi = x[i]; Scalar yi = y[i]; - x[i] = c * xi - ei_conj(s) * yi; - y[i] = s * xi + c * yi; + x[i] = j.c() * xi + ei_conj(j.s()) * yi; + y[i] = -j.s() * xi + ei_conj(j.c()) * yi; } } else @@ -213,8 +387,8 @@ void /*EIGEN_DONT_INLINE*/ ei_apply_rotation_in_the_plane(VectorX& _x, VectorY& { Scalar xi = *x; Scalar yi = *y; - *x = c * xi - ei_conj(s) * yi; - *y = s * xi + c * yi; + *x = j.c() * xi + ei_conj(j.s()) * yi; + *y = -j.s() * xi + ei_conj(j.c()) * yi; x += incrx; y += incry; } diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h index cee41b152..0d542cf7a 100644 --- a/Eigen/src/QR/FullPivotingHouseholderQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -67,12 +67,10 @@ template<typename MatrixType> class FullPivotingHouseholderQR * The default constructor is useful in cases in which the user intends to * perform decompositions via FullPivotingHouseholderQR::compute(const MatrixType&). */ - FullPivotingHouseholderQR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + FullPivotingHouseholderQR() : m_isInitialized(false) {} FullPivotingHouseholderQR(const MatrixType& matrix) - : m_qr(matrix.rows(), matrix.cols()), - m_hCoeffs(std::min(matrix.rows(),matrix.cols())), - m_isInitialized(false) + : m_isInitialized(false) { compute(matrix); } @@ -286,7 +284,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; - RealScalar biggest; + RealScalar biggest(0); for (int k = 0; k < size; ++k) { diff --git a/Eigen/src/QR/QrInstantiations.cpp b/Eigen/src/QR/QrInstantiations.cpp index 38b0a57da..695377d69 100644 --- a/Eigen/src/QR/QrInstantiations.cpp +++ b/Eigen/src/QR/QrInstantiations.cpp @@ -33,11 +33,6 @@ namespace Eigen { -template static void ei_tridiagonal_qr_step(float* , float* , int, int, float* , int); -template static void ei_tridiagonal_qr_step(double* , double* , int, int, double* , int); -template static void ei_tridiagonal_qr_step(float* , float* , int, int, std::complex<float>* , int); -template static void ei_tridiagonal_qr_step(double* , double* , int, int, std::complex<double>* , int); - EIGEN_QR_MODULE_INSTANTIATE(); } diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h new file mode 100644 index 000000000..2801ee077 --- /dev/null +++ b/Eigen/src/SVD/JacobiSVD.h @@ -0,0 +1,352 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#ifndef EIGEN_JACOBISVD_H +#define EIGEN_JACOBISVD_H + +/** \ingroup SVD_Module + * \nonstableyet + * + * \class JacobiSVD + * + * \brief Jacobi SVD decomposition of a square matrix + * + * \param MatrixType the type of the matrix of which we are computing the SVD decomposition + * \param Options a bit field of flags offering the following options: \c SkipU and \c SkipV allow to skip the computation of + * the unitaries \a U and \a V respectively; \c AtLeastAsManyRowsAsCols and \c AtLeastAsManyRowsAsCols allows + * to hint the compiler to only generate the corresponding code paths; \c Square is equivantent to the combination of + * the latter two bits and is useful when you know that the matrix is square. Note that when this information can + * be automatically deduced from the matrix type (e.g. a Matrix3f is always square), Eigen does it for you. + * + * \sa MatrixBase::jacobiSvd() + */ +template<typename MatrixType, unsigned int Options> class JacobiSVD +{ + private: + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; + enum { + ComputeU = (Options & SkipU) == 0, + ComputeV = (Options & SkipV) == 0, + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + DiagSizeAtCompileTime = EIGEN_SIZE_MIN(RowsAtCompileTime,ColsAtCompileTime), + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + MaxDiagSizeAtCompileTime = EIGEN_SIZE_MIN(MaxRowsAtCompileTime,MaxColsAtCompileTime), + MatrixOptions = MatrixType::Options + }; + + typedef Matrix<Scalar, Dynamic, Dynamic, MatrixOptions> DummyMatrixType; + typedef typename ei_meta_if<ComputeU, + Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, + MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime>, + DummyMatrixType>::ret MatrixUType; + typedef typename ei_meta_if<ComputeV, + Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, + MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime>, + DummyMatrixType>::ret MatrixVType; + typedef Matrix<RealScalar, DiagSizeAtCompileTime, 1, + MatrixOptions, MaxDiagSizeAtCompileTime, 1> SingularValuesType; + typedef Matrix<Scalar, 1, RowsAtCompileTime, MatrixOptions, 1, MaxRowsAtCompileTime> RowType; + typedef Matrix<Scalar, RowsAtCompileTime, 1, MatrixOptions, MaxRowsAtCompileTime, 1> ColType; + typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, + MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime> + WorkMatrixType; + + public: + + JacobiSVD() : m_isInitialized(false) {} + + JacobiSVD(const MatrixType& matrix) : m_isInitialized(false) + { + compute(matrix); + } + + JacobiSVD& compute(const MatrixType& matrix); + + const MatrixUType& matrixU() const + { + ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + return m_matrixU; + } + + const SingularValuesType& singularValues() const + { + ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + return m_singularValues; + } + + const MatrixVType& matrixV() const + { + ei_assert(m_isInitialized && "JacobiSVD is not initialized."); + return m_matrixV; + } + + protected: + MatrixUType m_matrixU; + MatrixVType m_matrixV; + SingularValuesType m_singularValues; + bool m_isInitialized; + + template<typename _MatrixType, unsigned int _Options, bool _IsComplex> + friend struct ei_svd_precondition_2x2_block_to_be_real; + template<typename _MatrixType, unsigned int _Options, bool _PossiblyMoreRowsThanCols> + friend struct ei_svd_precondition_if_more_rows_than_cols; + template<typename _MatrixType, unsigned int _Options, bool _PossiblyMoreRowsThanCols> + friend struct ei_svd_precondition_if_more_cols_than_rows; +}; + +template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> +struct ei_svd_precondition_2x2_block_to_be_real +{ + typedef JacobiSVD<MatrixType, Options> SVD; + static void run(typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&, int, int) {} +}; + +template<typename MatrixType, unsigned int Options> +struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options, true> +{ + typedef JacobiSVD<MatrixType, Options> SVD; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV }; + static void run(typename SVD::WorkMatrixType& work_matrix, JacobiSVD<MatrixType, Options>& svd, int p, int q) + { + Scalar z; + PlanarRotation<Scalar> rot; + RealScalar n = ei_sqrt(ei_abs2(work_matrix.coeff(p,p)) + ei_abs2(work_matrix.coeff(q,p))); + if(n==0) + { + z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); + work_matrix.row(p) *= z; + if(ComputeU) svd.m_matrixU.col(p) *= ei_conj(z); + z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + work_matrix.row(q) *= z; + if(ComputeU) svd.m_matrixU.col(q) *= ei_conj(z); + } + else + { + rot.c() = ei_conj(work_matrix.coeff(p,p)) / n; + rot.s() = work_matrix.coeff(q,p) / n; + work_matrix.applyOnTheLeft(p,q,rot); + if(ComputeU) svd.m_matrixU.applyOnTheRight(p,q,rot.adjoint()); + if(work_matrix.coeff(p,q) != Scalar(0)) + { + Scalar z = ei_abs(work_matrix.coeff(p,q)) / work_matrix.coeff(p,q); + work_matrix.col(q) *= z; + if(ComputeV) svd.m_matrixV.col(q) *= z; + } + if(work_matrix.coeff(q,q) != Scalar(0)) + { + z = ei_abs(work_matrix.coeff(q,q)) / work_matrix.coeff(q,q); + work_matrix.row(q) *= z; + if(ComputeU) svd.m_matrixU.col(q) *= ei_conj(z); + } + } + } +}; + +template<typename MatrixType, typename RealScalar> +void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q, + PlanarRotation<RealScalar> *j_left, + PlanarRotation<RealScalar> *j_right) +{ + Matrix<RealScalar,2,2> m; + m << ei_real(matrix.coeff(p,p)), ei_real(matrix.coeff(p,q)), + ei_real(matrix.coeff(q,p)), ei_real(matrix.coeff(q,q)); + PlanarRotation<RealScalar> rot1; + RealScalar t = m.coeff(0,0) + m.coeff(1,1); + RealScalar d = m.coeff(1,0) - m.coeff(0,1); + if(t == RealScalar(0)) + { + rot1.c() = 0; + rot1.s() = d > 0 ? 1 : -1; + } + else + { + RealScalar u = d / t; + rot1.c() = RealScalar(1) / ei_sqrt(1 + ei_abs2(u)); + rot1.s() = rot1.c() * u; + } + m.applyOnTheLeft(0,1,rot1); + j_right->makeJacobi(m,0,1); + *j_left = rot1 * j_right->transpose(); +} + +template<typename MatrixType, unsigned int Options, + bool PossiblyMoreRowsThanCols = (Options & AtLeastAsManyColsAsRows) == 0 + && (MatrixType::RowsAtCompileTime==Dynamic + || MatrixType::RowsAtCompileTime>MatrixType::ColsAtCompileTime)> +struct ei_svd_precondition_if_more_rows_than_cols +{ + typedef JacobiSVD<MatrixType, Options> SVD; + static bool run(const MatrixType&, typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&) { return false; } +}; + +template<typename MatrixType, unsigned int Options> +struct ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options, true> +{ + typedef JacobiSVD<MatrixType, Options> SVD; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + enum { ComputeU = SVD::ComputeU, ComputeV = SVD::ComputeV }; + static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd) + { + int rows = matrix.rows(); + int cols = matrix.cols(); + int diagSize = cols; + if(rows > cols) + { + FullPivotingHouseholderQR<MatrixType> qr(matrix); + work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>(); + if(ComputeU) svd.m_matrixU = qr.matrixQ(); + if(ComputeV) + for(int i = 0; i < cols; i++) + svd.m_matrixV.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); + return true; + } + else return false; + } +}; + +template<typename MatrixType, unsigned int Options, + bool PossiblyMoreColsThanRows = (Options & AtLeastAsManyRowsAsCols) == 0 + && (MatrixType::ColsAtCompileTime==Dynamic + || MatrixType::ColsAtCompileTime>MatrixType::RowsAtCompileTime)> +struct ei_svd_precondition_if_more_cols_than_rows +{ + typedef JacobiSVD<MatrixType, Options> SVD; + static bool run(const MatrixType&, typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&) { return false; } +}; + +template<typename MatrixType, unsigned int Options> +struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true> +{ + typedef JacobiSVD<MatrixType, Options> SVD; + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + enum { + ComputeU = SVD::ComputeU, + ComputeV = SVD::ComputeV, + RowsAtCompileTime = SVD::RowsAtCompileTime, + ColsAtCompileTime = SVD::ColsAtCompileTime, + MaxRowsAtCompileTime = SVD::MaxRowsAtCompileTime, + MaxColsAtCompileTime = SVD::MaxColsAtCompileTime, + MatrixOptions = SVD::MatrixOptions + }; + + static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd) + { + int rows = matrix.rows(); + int cols = matrix.cols(); + int diagSize = rows; + if(cols > rows) + { + typedef Matrix<Scalar,ColsAtCompileTime,RowsAtCompileTime, + MatrixOptions,MaxColsAtCompileTime,MaxRowsAtCompileTime> + TransposeTypeWithSameStorageOrder; + FullPivotingHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint()); + work_matrix = qr.matrixQR().block(0,0,diagSize,diagSize).template triangularView<UpperTriangular>().adjoint(); + if(ComputeV) svd.m_matrixV = qr.matrixQ(); + if(ComputeU) + for(int i = 0; i < rows; i++) + svd.m_matrixU.coeffRef(qr.colsPermutation().coeff(i),i) = Scalar(1); + return true; + } + else return false; + } +}; + +template<typename MatrixType, unsigned int Options> +JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute(const MatrixType& matrix) +{ + WorkMatrixType work_matrix; + int rows = matrix.rows(); + int cols = matrix.cols(); + int diagSize = std::min(rows, cols); + if(ComputeU) m_matrixU = MatrixUType::Zero(rows,rows); + if(ComputeV) m_matrixV = MatrixVType::Zero(cols,cols); + m_singularValues.resize(diagSize); + const RealScalar precision = 2 * epsilon<Scalar>(); + + if(!ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options>::run(matrix, work_matrix, *this) + && !ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options>::run(matrix, work_matrix, *this)) + { + work_matrix = matrix.block(0,0,diagSize,diagSize); + if(ComputeU) m_matrixU.diagonal().setOnes(); + if(ComputeV) m_matrixV.diagonal().setOnes(); + } + + bool finished = false; + while(!finished) + { + finished = true; + for(int p = 1; p < diagSize; ++p) + { + for(int q = 0; q < p; ++q) + { + if(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) + > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision) + { + finished = false; + ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(work_matrix, *this, p, q); + + PlanarRotation<RealScalar> j_left, j_right; + ei_real_2x2_jacobi_svd(work_matrix, p, q, &j_left, &j_right); + + work_matrix.applyOnTheLeft(p,q,j_left); + if(ComputeU) m_matrixU.applyOnTheRight(p,q,j_left.transpose()); + + work_matrix.applyOnTheRight(p,q,j_right); + if(ComputeV) m_matrixV.applyOnTheRight(p,q,j_right); + } + } + } + } + + for(int i = 0; i < diagSize; ++i) + { + RealScalar a = ei_abs(work_matrix.coeff(i,i)); + m_singularValues.coeffRef(i) = a; + if(ComputeU && (a!=RealScalar(0))) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a; + } + + for(int i = 0; i < diagSize; i++) + { + int pos; + m_singularValues.end(diagSize-i).maxCoeff(&pos); + if(pos) + { + pos += i; + std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); + if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i)); + if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); + } + } + + m_isInitialized = true; + return *this; +} +#endif // EIGEN_JACOBISVD_H diff --git a/Eigen/src/SVD/JacobiSquareSVD.h b/Eigen/src/SVD/JacobiSquareSVD.h deleted file mode 100644 index f21c9edca..000000000 --- a/Eigen/src/SVD/JacobiSquareSVD.h +++ /dev/null @@ -1,169 +0,0 @@ -// This file is part of Eigen, a lightweight C++ template library -// for linear algebra. -// -// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> -// -// Eigen is free software; you can redistribute it and/or -// modify it under the terms of the GNU Lesser General Public -// License as published by the Free Software Foundation; either -// version 3 of the License, or (at your option) any later version. -// -// Alternatively, you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation; either version 2 of -// the License, or (at your option) any later version. -// -// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY -// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS -// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU Lesser General Public -// License and a copy of the GNU General Public License along with -// Eigen. If not, see <http://www.gnu.org/licenses/>. - -#ifndef EIGEN_JACOBISQUARESVD_H -#define EIGEN_JACOBISQUARESVD_H - -/** \ingroup SVD_Module - * \nonstableyet - * - * \class JacobiSquareSVD - * - * \brief Jacobi SVD decomposition of a square matrix - * - * \param MatrixType the type of the matrix of which we are computing the SVD decomposition - * \param ComputeU whether the U matrix should be computed - * \param ComputeV whether the V matrix should be computed - * - * \sa MatrixBase::jacobiSvd() - */ -template<typename MatrixType, bool ComputeU, bool ComputeV> class JacobiSquareSVD -{ - private: - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - Options = MatrixType::Options - }; - - typedef Matrix<Scalar, Dynamic, Dynamic, Options> DummyMatrixType; - typedef typename ei_meta_if<ComputeU, - Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, - Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime>, - DummyMatrixType>::ret MatrixUType; - typedef typename Diagonal<MatrixType,0>::PlainMatrixType SingularValuesType; - typedef Matrix<Scalar, 1, RowsAtCompileTime, Options, 1, MaxRowsAtCompileTime> RowType; - typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> ColType; - - public: - - JacobiSquareSVD() : m_isInitialized(false) {} - - JacobiSquareSVD(const MatrixType& matrix) : m_isInitialized(false) - { - compute(matrix); - } - - JacobiSquareSVD& compute(const MatrixType& matrix); - - const MatrixUType& matrixU() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matrixU; - } - - const SingularValuesType& singularValues() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_singularValues; - } - - const MatrixUType& matrixV() const - { - ei_assert(m_isInitialized && "SVD is not initialized."); - return m_matrixV; - } - - protected: - MatrixUType m_matrixU; - MatrixUType m_matrixV; - SingularValuesType m_singularValues; - bool m_isInitialized; -}; - -template<typename MatrixType, bool ComputeU, bool ComputeV> -JacobiSquareSVD<MatrixType, ComputeU, ComputeV>& JacobiSquareSVD<MatrixType, ComputeU, ComputeV>::compute(const MatrixType& matrix) -{ - MatrixType work_matrix(matrix); - int size = matrix.rows(); - if(ComputeU) m_matrixU = MatrixUType::Identity(size,size); - if(ComputeV) m_matrixV = MatrixUType::Identity(size,size); - m_singularValues.resize(size); - const RealScalar precision = 2 * epsilon<Scalar>(); - -sweep_again: - for(int p = 1; p < size; ++p) - { - for(int q = 0; q < p; ++q) - { - Scalar c, s; - while(std::max(ei_abs(work_matrix.coeff(p,q)),ei_abs(work_matrix.coeff(q,p))) - > std::max(ei_abs(work_matrix.coeff(p,p)),ei_abs(work_matrix.coeff(q,q)))*precision) - { - if(work_matrix.makeJacobiForAtA(p,q,&c,&s)) - { - work_matrix.applyJacobiOnTheRight(p,q,c,s); - if(ComputeV) m_matrixV.applyJacobiOnTheRight(p,q,c,s); - } - if(work_matrix.makeJacobiForAAt(p,q,&c,&s)) - { - ei_normalizeJacobi(&c, &s, work_matrix.coeff(p,p), work_matrix.coeff(q,p)), - work_matrix.applyJacobiOnTheLeft(p,q,c,s); - if(ComputeU) m_matrixU.applyJacobiOnTheRight(p,q,c,s); - } - } - } - } - - RealScalar biggestOnDiag = work_matrix.diagonal().cwise().abs().maxCoeff(); - RealScalar maxAllowedOffDiag = biggestOnDiag * precision; - for(int p = 0; p < size; ++p) - { - for(int q = 0; q < p; ++q) - if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) - goto sweep_again; - for(int q = p+1; q < size; ++q) - if(ei_abs(work_matrix.coeff(p,q)) > maxAllowedOffDiag) - goto sweep_again; - } - - m_singularValues = work_matrix.diagonal().cwise().abs(); - RealScalar biggestSingularValue = m_singularValues.maxCoeff(); - - for(int i = 0; i < size; ++i) - { - RealScalar a = ei_abs(work_matrix.coeff(i,i)); - m_singularValues.coeffRef(i) = a; - if(ComputeU && !ei_isMuchSmallerThan(a, biggestSingularValue)) m_matrixU.col(i) *= work_matrix.coeff(i,i)/a; - } - - for(int i = 0; i < size; i++) - { - int pos; - m_singularValues.end(size-i).maxCoeff(&pos); - if(pos) - { - pos += i; - std::swap(m_singularValues.coeffRef(i), m_singularValues.coeffRef(pos)); - if(ComputeU) m_matrixU.col(pos).swap(m_matrixU.col(i)); - if(ComputeV) m_matrixV.col(pos).swap(m_matrixV.col(i)); - } - } - - m_isInitialized = true; - return *this; -} -#endif // EIGEN_JACOBISQUARESVD_H diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 1a7e6c49a..da01cf396 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -309,7 +309,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) h = Scalar(1.0)/h; c = g*h; s = -f*h; - V.applyJacobiOnTheRight(i,nm,c,s); + V.applyOnTheRight(i,nm,PlanarRotation<Scalar>(c,s)); } } z = W[k]; @@ -342,7 +342,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) y = W[i]; h = s*g; g = c*g; - + z = pythag(f,h); rv1[j] = z; c = f/z; @@ -351,8 +351,8 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) g = g*c - x*s; h = y*s; y *= c; - V.applyJacobiOnTheRight(i,j,c,s); - + V.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s)); + z = pythag(f,h); W[j] = z; // Rotation can be arbitrary if z = 0. @@ -364,7 +364,7 @@ SVD<MatrixType>& SVD<MatrixType>::compute(const MatrixType& matrix) } f = c*g + s*y; x = c*y - s*g; - A.applyJacobiOnTheRight(i,j,c,s); + A.applyOnTheRight(i,j,PlanarRotation<Scalar>(c,s)); } rv1[l] = 0.0; rv1[k] = f; diff --git a/Eigen/src/Sparse/AmbiVector.h b/Eigen/src/Sparse/AmbiVector.h index 974e5c6c4..474626848 100644 --- a/Eigen/src/Sparse/AmbiVector.h +++ b/Eigen/src/Sparse/AmbiVector.h @@ -41,7 +41,7 @@ template<typename _Scalar> class AmbiVector resize(size); } - void init(RealScalar estimatedDensity); + void init(double estimatedDensity); void init(int mode); int nonZeros() const; @@ -143,7 +143,7 @@ int AmbiVector<Scalar>::nonZeros() const } template<typename Scalar> -void AmbiVector<Scalar>::init(RealScalar estimatedDensity) +void AmbiVector<Scalar>::init(double estimatedDensity) { if (estimatedDensity>0.1) init(IsDense); diff --git a/Eigen/src/Sparse/CholmodSupport.h b/Eigen/src/Sparse/CholmodSupport.h index ad59c89af..30a33c3dc 100644 --- a/Eigen/src/Sparse/CholmodSupport.h +++ b/Eigen/src/Sparse/CholmodSupport.h @@ -99,7 +99,7 @@ cholmod_dense ei_cholmod_map_eigen_to_dense(MatrixBase<Derived>& mat) res.nrow = mat.rows(); res.ncol = mat.cols(); res.nzmax = res.nrow * res.ncol; - res.d = mat.derived().stride(); + res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().stride(); res.x = mat.derived().data(); res.z = 0; @@ -157,7 +157,7 @@ class SparseLLT<MatrixType,Cholmod> : public SparseLLT<MatrixType> inline const typename Base::CholMatrixType& matrixL(void) const; template<typename Derived> - void solveInPlace(MatrixBase<Derived> &b) const; + bool solveInPlace(MatrixBase<Derived> &b) const; void compute(const MatrixType& matrix); @@ -216,7 +216,7 @@ SparseLLT<MatrixType,Cholmod>::matrixL() const template<typename MatrixType> template<typename Derived> -void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const +bool SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const { const int size = m_cholmodFactor->n; ei_assert(size==b.rows()); @@ -228,9 +228,16 @@ void SparseLLT<MatrixType,Cholmod>::solveInPlace(MatrixBase<Derived> &b) const // as long as our own triangular sparse solver is not fully optimal, // let's use CHOLMOD's one: cholmod_dense cdb = ei_cholmod_map_eigen_to_dense(b); - cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod); + //cholmod_dense* x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod); + cholmod_dense* x = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &cdb, &m_cholmod); + if(!x) + { + std::cerr << "Eigen: cholmod_solve failed\n"; + return false; + } b = Matrix<typename Base::Scalar,Dynamic,1>::Map(reinterpret_cast<typename Base::Scalar*>(x->x),b.rows()); cholmod_free_dense(&x, &m_cholmod); + return true; } #endif // EIGEN_CHOLMODSUPPORT_H diff --git a/Eigen/src/Sparse/SuperLUSupport.h b/Eigen/src/Sparse/SuperLUSupport.h index b4f34f094..708f177e8 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/Eigen/src/Sparse/SuperLUSupport.h @@ -48,6 +48,37 @@ DECL_GSSVX(SuperLU_C,cgssvx,float,std::complex<float>) DECL_GSSVX(SuperLU_D,dgssvx,double,double) DECL_GSSVX(SuperLU_Z,zgssvx,double,std::complex<double>) +#ifdef MILU_ALPHA +#define EIGEN_SUPERLU_HAS_ILU +#endif + +#ifdef EIGEN_SUPERLU_HAS_ILU + +// similarly for the incomplete factorization using gsisx +#define DECL_GSISX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \ + inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \ + int *perm_c, int *perm_r, int *etree, char *equed, \ + FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \ + SuperMatrix *U, void *work, int lwork, \ + SuperMatrix *B, SuperMatrix *X, \ + FLOATTYPE *recip_pivot_growth, \ + FLOATTYPE *rcond, \ + SuperLUStat_t *stats, int *info, KEYTYPE) { \ + using namespace NAMESPACE; \ + mem_usage_t mem_usage; \ + NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \ + U, work, lwork, B, X, recip_pivot_growth, rcond, \ + &mem_usage, stats, info); \ + return mem_usage.for_lu; /* bytes used by the factor storage */ \ + } + +DECL_GSISX(SuperLU_S,sgsisx,float,float) +DECL_GSISX(SuperLU_C,cgsisx,float,std::complex<float>) +DECL_GSISX(SuperLU_D,dgsisx,double,double) +DECL_GSISX(SuperLU_Z,zgsisx,double,std::complex<double>) + +#endif + template<typename MatrixType> struct SluMatrixMapHelper; @@ -71,7 +102,7 @@ struct SluMatrix : SuperMatrix Store = &storage; storage = other.storage; } - + SluMatrix& operator=(const SluMatrix& other) { SuperMatrix::operator=(static_cast<const SuperMatrix&>(other)); @@ -130,7 +161,7 @@ struct SluMatrix : SuperMatrix res.nrow = mat.rows(); res.ncol = mat.cols(); - res.storage.lda = mat.stride(); + res.storage.lda = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.stride(); res.storage.values = mat.data(); return res; } @@ -373,7 +404,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a) m_sluA = m_matrix.asSluMatrix(); memset(&m_sluL,0,sizeof m_sluL); memset(&m_sluU,0,sizeof m_sluU); - m_sluEqued = 'B'; + //m_sluEqued = 'B'; int info = 0; m_p.resize(size); @@ -395,14 +426,44 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a) m_sluX = m_sluB; StatInit(&m_sluStat); - SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], - &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], - &m_sluL, &m_sluU, - NULL, 0, - &m_sluB, &m_sluX, - &recip_pivot_gross, &rcond, - &ferr, &berr, - &m_sluStat, &info, Scalar()); + if (m_flags&IncompleteFactorization) + { + #ifdef EIGEN_SUPERLU_HAS_ILU + ilu_set_default_options(&m_sluOptions); + + // no attempt to preserve column sum + m_sluOptions.ILU_MILU = SILU; + + // only basic ILU(k) support -- no direct control over memory consumption + // better to use ILU_DropRule = DROP_BASIC | DROP_AREA + // and set ILU_FillFactor to max memory growth + m_sluOptions.ILU_DropRule = DROP_BASIC; + m_sluOptions.ILU_DropTol = Base::m_precision; + + SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], + &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_gross, &rcond, + &m_sluStat, &info, Scalar()); + #else + std::cerr << "Incomplete factorization is only available in SuperLU v4\n"; + Base::m_succeeded = false; + return; + #endif + } + else + { + SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], + &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_gross, &rcond, + &ferr, &berr, + &m_sluStat, &info, Scalar()); + } StatFree(&m_sluStat); m_extractedDataAreDirty = true; @@ -440,17 +501,36 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b, StatInit(&m_sluStat); int info = 0; RealScalar recip_pivot_gross, rcond; - SuperLU_gssvx( - &m_sluOptions, &m_sluA, - m_q.data(), m_p.data(), - &m_sluEtree[0], &m_sluEqued, - &m_sluRscale[0], &m_sluCscale[0], - &m_sluL, &m_sluU, - NULL, 0, - &m_sluB, &m_sluX, - &recip_pivot_gross, &rcond, - &m_sluFerr[0], &m_sluBerr[0], - &m_sluStat, &info, Scalar()); + + if (m_flags&IncompleteFactorization) + { + #ifdef EIGEN_SUPERLU_HAS_ILU + SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], + &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_gross, &rcond, + &m_sluStat, &info, Scalar()); + #else + std::cerr << "Incomplete factorization is only available in SuperLU v4\n"; + return false; + #endif + } + else + { + SuperLU_gssvx( + &m_sluOptions, &m_sluA, + m_q.data(), m_p.data(), + &m_sluEtree[0], &m_sluEqued, + &m_sluRscale[0], &m_sluCscale[0], + &m_sluL, &m_sluU, + NULL, 0, + &m_sluB, &m_sluX, + &recip_pivot_gross, &rcond, + &m_sluFerr[0], &m_sluBerr[0], + &m_sluStat, &info, Scalar()); + } StatFree(&m_sluStat); // reset to previous state diff --git a/cmake/EigenTesting.cmake b/cmake/EigenTesting.cmake index e3a87f645..faa75c6f4 100644 --- a/cmake/EigenTesting.cmake +++ b/cmake/EigenTesting.cmake @@ -153,11 +153,17 @@ macro(ei_init_testing) endmacro(ei_init_testing) if(CMAKE_COMPILER_IS_GNUCXX) + option(EIGEN_COVERAGE_TESTING "Enable/disable gcov" OFF) + if(EIGEN_COVERAGE_TESTING) + set(COVERAGE_FLAGS "-fprofile-arcs -ftest-coverage") + else(EIGEN_COVERAGE_TESTING) + set(COVERAGE_FLAGS "") + endif(EIGEN_COVERAGE_TESTING) if(CMAKE_SYSTEM_NAME MATCHES Linux) - set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -g2") - set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -O2 -g2") - set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fno-inline-functions") - set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -g2") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${COVERAGE_FLAGS} -g2") + set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} ${COVERAGE_FLAGS} -O2 -g2") + set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} ${COVERAGE_FLAGS} -fno-inline-functions") + set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${COVERAGE_FLAGS} -O0 -g2") endif(CMAKE_SYSTEM_NAME MATCHES Linux) set(EI_OFLAG "-O2") elseif(MSVC) diff --git a/doc/C01_QuickStartGuide.dox b/doc/C01_QuickStartGuide.dox index d43dbd72d..2943aed80 100644 --- a/doc/C01_QuickStartGuide.dox +++ b/doc/C01_QuickStartGuide.dox @@ -129,7 +129,7 @@ The default constructor leaves coefficients uninitialized. Any dynamic size is s Matrix3f A; // construct 3x3 matrix with uninitialized coefficients A(0,0) = 5; // OK MatrixXf B; // construct 0x0 matrix without allocating anything -A(0,0) = 5; // Error, B is uninitialized, doesn't have any coefficients to address +B(0,0) = 5; // Error, B is uninitialized, doesn't have any coefficients to address \endcode In the above example, B is an uninitialized matrix. What to do with such a matrix? You can call resize() on it, or you can assign another matrix to it. Like this: @@ -261,7 +261,7 @@ v = 6 6 6 \subsection TutorialCasting Casting -In Eigen, any matrices of same size and same scalar type are all naturally compatible. The scalar type can be explicitely casted to another one using the template MatrixBase::cast() function: +In Eigen, any matrices of same size and same scalar type are all naturally compatible. The scalar type can be explicitly casted to another one using the template MatrixBase::cast() function: \code Matrix3d md(1,2,3); Matrix3f mf = md.cast<float>(); @@ -328,7 +328,7 @@ In short, all arithmetic operators can be used right away as in the following ex mat4 -= mat1*1.5 + mat2 * (mat3/4); \endcode which includes two matrix scalar products ("mat1*1.5" and "mat3/4"), a matrix-matrix product ("mat2 * (mat3/4)"), -a matrix addition ("+") and substraction with assignment ("-="). +a matrix addition ("+") and subtraction with assignment ("-="). <table class="tutorial_code"> <tr><td> @@ -464,7 +464,7 @@ mat = 2 7 8 Also note that maxCoeff and minCoeff can takes optional arguments returning the coordinates of the respective min/max coeff: \link MatrixBase::maxCoeff(int*,int*) const maxCoeff(int* i, int* j) \endlink, \link MatrixBase::minCoeff(int*,int*) const minCoeff(int* i, int* j) \endlink. -<span class="note">\b Side \b note: The all() and any() functions are especially useful in combinaison with coeff-wise comparison operators (\ref CwiseAll "example").</span> +<span class="note">\b Side \b note: The all() and any() functions are especially useful in combination with coeff-wise comparison operators (\ref CwiseAll "example").</span> @@ -578,7 +578,7 @@ vec1.normalize();\endcode <a href="#" class="top">top</a>\section TutorialCoreTriangularMatrix Dealing with triangular matrices -Currently, Eigen does not provide any explcit triangular matrix, with storage class. Instead, we +Currently, Eigen does not provide any explicit triangular matrix, with storage class. Instead, we can reference a triangular part of a square matrix or expression to perform special treatment on it. This is achieved by the class TriangularView and the MatrixBase::triangularView template function. Note that the opposite triangular part of the matrix is never referenced, and so it can, e.g., store @@ -595,12 +595,12 @@ m.triangularView<Eigen::LowerTriangular>() m.triangularView<Eigen::UnitLowerTriangular>()\endcode </td></tr> <tr><td> -Writting to a specific triangular part:\n (only the referenced triangular part is evaluated) +Writing to a specific triangular part:\n (only the referenced triangular part is evaluated) </td><td>\code m1.triangularView<Eigen::LowerTriangular>() = m2 + m3 \endcode </td></tr> <tr><td> -Convertion to a dense matrix setting the opposite triangular part to zero: +Conversion to a dense matrix setting the opposite triangular part to zero: </td><td>\code m2 = m1.triangularView<Eigen::UnitUpperTriangular>()\endcode </td></tr> diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index cd67bb07d..5b055ed11 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -208,10 +208,11 @@ ALIASES = "only_for_vectors=This is only for vectors (either row- "lu_module=This is defined in the %LU module. \code #include <Eigen/LU> \endcode" \ "cholesky_module=This is defined in the %Cholesky module. \code #include <Eigen/Cholesky> \endcode" \ "qr_module=This is defined in the %QR module. \code #include <Eigen/QR> \endcode" \ + "jacobi_module=This is defined in the %Jacobi module. \code #include <Eigen/Jacobi> \endcode" \ "svd_module=This is defined in the %SVD module. \code #include <Eigen/SVD> \endcode" \ "geometry_module=This is defined in the %Geometry module. \code #include <Eigen/Geometry> \endcode" \ "leastsquares_module=This is defined in the %LeastSquares module. \code #include <Eigen/LeastSquares> \endcode" \ - "addexample=\anchor" \ + "eigenvalues_module=This is defined in the %Eigenvalues module. \code #include <Eigen/Eigenvalues> \endcode" \ "label=\bug" \ "redstar=<a href='#warningarraymodule' style='color:red;text-decoration: none;'>*</a>" \ "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" diff --git a/doc/snippets/Jacobi_makeGivens.cpp b/doc/snippets/Jacobi_makeGivens.cpp new file mode 100644 index 000000000..3a4defe24 --- /dev/null +++ b/doc/snippets/Jacobi_makeGivens.cpp @@ -0,0 +1,6 @@ +Vector2f v = Vector2f::Random(); +PlanarRotation<float> G; +G.makeGivens(v.x(), v.y()); +cout << "Here is the vector v:" << endl << v << endl; +v.applyOnTheLeft(0, 1, G.adjoint()); +cout << "Here is the vector J' * v:" << endl << v << endl;
\ No newline at end of file diff --git a/doc/snippets/Jacobi_makeJacobi.cpp b/doc/snippets/Jacobi_makeJacobi.cpp new file mode 100644 index 000000000..5c0ab7374 --- /dev/null +++ b/doc/snippets/Jacobi_makeJacobi.cpp @@ -0,0 +1,8 @@ +Matrix2f m = Matrix2f::Random(); +m = (m + m.adjoint()).eval(); +PlanarRotation<float> J; +J.makeJacobi(m, 0, 1); +cout << "Here is the matrix m:" << endl << m << endl; +m.applyOnTheLeft(0, 1, J.adjoint()); +m.applyOnTheRight(0, 1, J); +cout << "Here is the matrix J' * m * J:" << endl << m << endl;
\ No newline at end of file diff --git a/doc/snippets/compile_snippet.cpp.in b/doc/snippets/compile_snippet.cpp.in index d074cac50..3bea1ac8d 100644 --- a/doc/snippets/compile_snippet.cpp.in +++ b/doc/snippets/compile_snippet.cpp.in @@ -4,6 +4,7 @@ #include <Eigen/QR> #include <Eigen/Cholesky> #include <Eigen/Geometry> +#include <Eigen/Jacobi> using namespace Eigen; using namespace std; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index fc43bbb1d..4e279ea47 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -115,6 +115,7 @@ ei_add_test(product_trmv ${EI_OFLAG}) ei_add_test(product_trmm ${EI_OFLAG}) ei_add_test(product_trsm ${EI_OFLAG}) ei_add_test(product_notemporary ${EI_OFLAG}) +ei_add_test(stable_norm) ei_add_test(bandmatrix) ei_add_test(cholesky " " "${GSL_LIBRARIES}") ei_add_test(lu ${EI_OFLAG}) @@ -125,7 +126,9 @@ ei_add_test(qr_colpivoting) ei_add_test(qr_fullpivoting) ei_add_test(eigensolver_selfadjoint " " "${GSL_LIBRARIES}") ei_add_test(eigensolver_generic " " "${GSL_LIBRARIES}") +ei_add_test(eigensolver_complex) ei_add_test(svd) +ei_add_test(jacobisvd ${EI_OFLAG}) ei_add_test(geo_orthomethods) ei_add_test(geo_homogeneous) ei_add_test(geo_quaternion) @@ -146,6 +149,8 @@ ei_add_test(sparse_product) ei_add_test(sparse_solvers " " "${SPARSE_LIBS}") ei_add_test(umeyama) ei_add_test(householder) +ei_add_test(swap) +ei_add_test(conservative_resize) ei_add_property(EIGEN_TESTING_SUMMARY "CXX: ${CMAKE_CXX_COMPILER}\n") if(CMAKE_COMPILER_IS_GNUCXX) diff --git a/test/adjoint.cpp b/test/adjoint.cpp index 964658c65..bebf47ac3 100644 --- a/test/adjoint.cpp +++ b/test/adjoint.cpp @@ -72,13 +72,6 @@ template<typename MatrixType> void adjoint(const MatrixType& m) if(NumTraits<Scalar>::HasFloatingPoint) VERIFY_IS_APPROX(v1.squaredNorm(), v1.norm() * v1.norm()); VERIFY_IS_MUCH_SMALLER_THAN(ei_abs(vzero.dot(v1)), static_cast<RealScalar>(1)); - if(NumTraits<Scalar>::HasFloatingPoint) - { - VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1)); - VERIFY_IS_APPROX(v1.norm(), v1.stableNorm()); - VERIFY_IS_APPROX(v1.blueNorm(), v1.stableNorm()); - VERIFY_IS_APPROX(v1.hypotNorm(), v1.stableNorm()); - } // check compatibility of dot and adjoint VERIFY(ei_isApprox(v1.dot(square * v2), (square.adjoint() * v1).dot(v2), largerEps)); @@ -124,7 +117,7 @@ void test_adjoint() } // test a large matrix only once CALL_SUBTEST( adjoint(Matrix<float, 100, 100>()) ); - + { MatrixXcf a(10,10), b(10,10); VERIFY_RAISES_ASSERT(a = a.transpose()); diff --git a/test/bandmatrix.cpp b/test/bandmatrix.cpp index f677e7b85..2bdc67e28 100644 --- a/test/bandmatrix.cpp +++ b/test/bandmatrix.cpp @@ -44,21 +44,21 @@ template<typename MatrixType> void bandmatrix(const MatrixType& _m) dm1.diagonal().setConstant(123); for (int i=1; i<=m.supers();++i) { - m.diagonal(i).setConstant(i); - dm1.diagonal(i).setConstant(i); + m.diagonal(i).setConstant(static_cast<RealScalar>(i)); + dm1.diagonal(i).setConstant(static_cast<RealScalar>(i)); } for (int i=1; i<=m.subs();++i) { - m.diagonal(-i).setConstant(-i); - dm1.diagonal(-i).setConstant(-i); + m.diagonal(-i).setConstant(-static_cast<RealScalar>(i)); + dm1.diagonal(-i).setConstant(-static_cast<RealScalar>(i)); } //std::cerr << m.m_data << "\n\n" << m.toDense() << "\n\n" << dm1 << "\n\n\n\n"; VERIFY_IS_APPROX(dm1,m.toDense()); for (int i=0; i<cols; ++i) { - m.col(i).setConstant(i+1); - dm1.col(i).setConstant(i+1); + m.col(i).setConstant(static_cast<RealScalar>(i+1)); + dm1.col(i).setConstant(static_cast<RealScalar>(i+1)); } int d = std::min(rows,cols); int a = std::max(0,cols-d-supers); diff --git a/test/basicstuff.cpp b/test/basicstuff.cpp index 52aff93ee..29df99f9e 100644 --- a/test/basicstuff.cpp +++ b/test/basicstuff.cpp @@ -99,15 +99,6 @@ template<typename MatrixType> void basicStuff(const MatrixType& m) MatrixType m4; VERIFY_IS_APPROX(m4 = m1,m1); - // test swap - m3 = m1; - m1.swap(m2); - VERIFY_IS_APPROX(m3, m2); - if(rows*cols>=3) - { - VERIFY_IS_NOT_APPROX(m3, m1); - } - m3.real() = m1.real(); VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), static_cast<const MatrixType&>(m1).real()); VERIFY_IS_APPROX(static_cast<const MatrixType&>(m3).real(), m1.real()); diff --git a/test/cholesky.cpp b/test/cholesky.cpp index df937fd0f..526a9f9d0 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -82,7 +82,7 @@ template<typename MatrixType> void cholesky(const MatrixType& m) // // test gsl itself ! // VERIFY_IS_APPROX(vecB, _vecB); // VERIFY_IS_APPROX(vecX, _vecX); -// +// // Gsl::free(gMatA); // Gsl::free(gSymm); // Gsl::free(gVecB); @@ -149,16 +149,16 @@ void test_cholesky() { for(int i = 0; i < g_repeat; i++) { CALL_SUBTEST( cholesky(Matrix<double,1,1>()) ); -// CALL_SUBTEST( cholesky(MatrixXd(1,1)) ); -// CALL_SUBTEST( cholesky(Matrix2d()) ); -// CALL_SUBTEST( cholesky(Matrix3f()) ); -// CALL_SUBTEST( cholesky(Matrix4d()) ); + CALL_SUBTEST( cholesky(MatrixXd(1,1)) ); + CALL_SUBTEST( cholesky(Matrix2d()) ); + CALL_SUBTEST( cholesky(Matrix3f()) ); + CALL_SUBTEST( cholesky(Matrix4d()) ); CALL_SUBTEST( cholesky(MatrixXd(200,200)) ); CALL_SUBTEST( cholesky(MatrixXcd(100,100)) ); } -// CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() ); -// CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() ); -// CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() ); -// CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() ); + CALL_SUBTEST( cholesky_verify_assert<Matrix3f>() ); + CALL_SUBTEST( cholesky_verify_assert<Matrix3d>() ); + CALL_SUBTEST( cholesky_verify_assert<MatrixXf>() ); + CALL_SUBTEST( cholesky_verify_assert<MatrixXd>() ); } diff --git a/test/conservative_resize.cpp b/test/conservative_resize.cpp new file mode 100644 index 000000000..b92dd5449 --- /dev/null +++ b/test/conservative_resize.cpp @@ -0,0 +1,129 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2009 Hauke Heibel <hauke.heibel@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or1 FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" + +#include <Eigen/Core> +#include <Eigen/Array> + +using namespace Eigen; + +template <typename Scalar, int Storage> +void run_matrix_tests() +{ + typedef Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic, Storage> MatrixType; + + MatrixType m, n; + + // boundary cases ... + m = n = MatrixType::Random(50,50); + m.conservativeResize(1,50); + VERIFY_IS_APPROX(m, n.block(0,0,1,50)); + + m = n = MatrixType::Random(50,50); + m.conservativeResize(50,1); + VERIFY_IS_APPROX(m, n.block(0,0,50,1)); + + m = n = MatrixType::Random(50,50); + m.conservativeResize(50,50); + VERIFY_IS_APPROX(m, n.block(0,0,50,50)); + + // random shrinking ... + for (int i=0; i<25; ++i) + { + const int rows = ei_random<int>(1,50); + const int cols = ei_random<int>(1,50); + m = n = MatrixType::Random(50,50); + m.conservativeResize(rows,cols); + VERIFY_IS_APPROX(m, n.block(0,0,rows,cols)); + } + + // random growing with zeroing ... + for (int i=0; i<25; ++i) + { + const int rows = ei_random<int>(50,75); + const int cols = ei_random<int>(50,75); + m = n = MatrixType::Random(50,50); + m.conservativeResizeLike(MatrixType::Zero(rows,cols)); + VERIFY_IS_APPROX(m.block(0,0,n.rows(),n.cols()), n); + VERIFY( rows<=50 || m.block(50,0,rows-50,cols).sum() == Scalar(0) ); + VERIFY( cols<=50 || m.block(0,50,rows,cols-50).sum() == Scalar(0) ); + } +} + +template <typename Scalar> +void run_vector_tests() +{ + typedef Matrix<Scalar, 1, Eigen::Dynamic> MatrixType; + + MatrixType m, n; + + // boundary cases ... + m = n = MatrixType::Random(50); + m.conservativeResize(1); + VERIFY_IS_APPROX(m, n.segment(0,1)); + + m = n = MatrixType::Random(50); + m.conservativeResize(50); + VERIFY_IS_APPROX(m, n.segment(0,50)); + + // random shrinking ... + for (int i=0; i<50; ++i) + { + const int size = ei_random<int>(1,50); + m = n = MatrixType::Random(50); + m.conservativeResize(size); + VERIFY_IS_APPROX(m, n.segment(0,size)); + } + + // random growing with zeroing ... + for (int i=0; i<50; ++i) + { + const int size = ei_random<int>(50,100); + m = n = MatrixType::Random(50); + m.conservativeResizeLike(MatrixType::Zero(size)); + VERIFY_IS_APPROX(m.segment(0,50), n); + VERIFY( size<=50 || m.segment(50,size-50).sum() == Scalar(0) ); + } +} + +void test_conservative_resize() +{ + run_matrix_tests<int, Eigen::RowMajor>(); + run_matrix_tests<int, Eigen::ColMajor>(); + run_matrix_tests<float, Eigen::RowMajor>(); + run_matrix_tests<float, Eigen::ColMajor>(); + run_matrix_tests<double, Eigen::RowMajor>(); + run_matrix_tests<double, Eigen::ColMajor>(); + run_matrix_tests<std::complex<float>, Eigen::RowMajor>(); + run_matrix_tests<std::complex<float>, Eigen::ColMajor>(); + run_matrix_tests<std::complex<double>, Eigen::RowMajor>(); + run_matrix_tests<std::complex<double>, Eigen::ColMajor>(); + + run_vector_tests<int>(); + run_vector_tests<float>(); + run_vector_tests<double>(); + run_vector_tests<std::complex<float> >(); + run_vector_tests<std::complex<double> >(); +} diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp new file mode 100644 index 000000000..38ede7c4a --- /dev/null +++ b/test/eigensolver_complex.cpp @@ -0,0 +1,61 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. Eigen itself is part of the KDE project. +// +// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/Eigenvalues> +#include <Eigen/LU> + +template<typename MatrixType> void eigensolver(const MatrixType& m) +{ + /* this test covers the following files: + ComplexEigenSolver.h, and indirectly ComplexSchur.h + */ + int rows = m.rows(); + int cols = m.cols(); + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType; + typedef Matrix<RealScalar, MatrixType::RowsAtCompileTime, 1> RealVectorType; + typedef typename std::complex<typename NumTraits<typename MatrixType::Scalar>::Real> Complex; + + MatrixType a = MatrixType::Random(rows,cols); + MatrixType symmA = a.adjoint() * a; + + ComplexEigenSolver<MatrixType> ei0(symmA); + VERIFY_IS_APPROX(symmA * ei0.eigenvectors(), ei0.eigenvectors() * ei0.eigenvalues().asDiagonal()); + + ComplexEigenSolver<MatrixType> ei1(a); + VERIFY_IS_APPROX(a * ei1.eigenvectors(), ei1.eigenvectors() * ei1.eigenvalues().asDiagonal()); + +} + +void test_eigensolver_complex() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( eigensolver(Matrix4cf()) ); + CALL_SUBTEST( eigensolver(MatrixXcd(14,14)) ); + } +} + diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index 2be49faf4..e2b2055b4 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -23,7 +23,7 @@ // Eigen. If not, see <http://www.gnu.org/licenses/>. #include "main.h" -#include <Eigen/QR> +#include <Eigen/Eigenvalues> #ifdef HAS_GSL #include "gsl_helper.h" diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 2571dbf43..3836c074b 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -23,7 +23,7 @@ // Eigen. If not, see <http://www.gnu.org/licenses/>. #include "main.h" -#include <Eigen/QR> +#include <Eigen/Eigenvalues> #ifdef HAS_GSL #include "gsl_helper.h" diff --git a/test/geo_orthomethods.cpp b/test/geo_orthomethods.cpp index 5e1d5bdb4..540a63b82 100644 --- a/test/geo_orthomethods.cpp +++ b/test/geo_orthomethods.cpp @@ -86,10 +86,10 @@ template<typename Scalar, int Size> void orthomethods(int size=Size) VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1)); VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), RealScalar(1)); - if (size>3) + if (size>=3) { - v0.template start<3>().setZero(); - v0.end(size-3).setRandom(); + v0.template start<2>().setZero(); + v0.end(size-2).setRandom(); VERIFY_IS_MUCH_SMALLER_THAN(v0.unitOrthogonal().dot(v0), Scalar(1)); VERIFY_IS_APPROX(v0.unitOrthogonal().norm(), RealScalar(1)); diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp new file mode 100644 index 000000000..5940b8497 --- /dev/null +++ b/test/jacobisvd.cpp @@ -0,0 +1,106 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr> +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" +#include <Eigen/SVD> +#include <Eigen/LU> + +template<typename MatrixType, unsigned int Options> void svd(const MatrixType& m = MatrixType(), bool pickrandom = true) +{ + int rows = m.rows(); + int cols = m.cols(); + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime + }; + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixUType; + typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> MatrixVType; + typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; + typedef Matrix<Scalar, ColsAtCompileTime, 1> InputVectorType; + + MatrixType a; + if(pickrandom) a = MatrixType::Random(rows,cols); + else a = m; + + JacobiSVD<MatrixType,Options> svd(a); + MatrixType sigma = MatrixType::Zero(rows,cols); + sigma.diagonal() = svd.singularValues().template cast<Scalar>(); + MatrixUType u = svd.matrixU(); + MatrixVType v = svd.matrixV(); + + VERIFY_IS_APPROX(a, u * sigma * v.adjoint()); + VERIFY_IS_UNITARY(u); + VERIFY_IS_UNITARY(v); +} + +template<typename MatrixType> void svd_verify_assert() +{ + MatrixType tmp; + + SVD<MatrixType> svd; + //VERIFY_RAISES_ASSERT(svd.solve(tmp, &tmp)) + VERIFY_RAISES_ASSERT(svd.matrixU()) + VERIFY_RAISES_ASSERT(svd.singularValues()) + VERIFY_RAISES_ASSERT(svd.matrixV()) + /*VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp))*/ +} + +void test_jacobisvd() +{ + for(int i = 0; i < g_repeat; i++) { + Matrix2cd m; + m << 0, 1, + 0, 1; + CALL_SUBTEST(( svd<Matrix2cd,0>(m, false) )); + m << 1, 0, + 1, 0; + CALL_SUBTEST(( svd<Matrix2cd,0>(m, false) )); + Matrix2d n; + n << 1, 1, + 1, -1; + CALL_SUBTEST(( svd<Matrix2d,0>(n, false) )); + CALL_SUBTEST(( svd<Matrix3f,0>() )); + CALL_SUBTEST(( svd<Matrix4d,Square>() )); + CALL_SUBTEST(( svd<Matrix<float,3,5> , AtLeastAsManyColsAsRows>() )); + CALL_SUBTEST(( svd<Matrix<double,Dynamic,2> , AtLeastAsManyRowsAsCols>(Matrix<double,Dynamic,2>(10,2)) )); + + CALL_SUBTEST(( svd<MatrixXf,Square>(MatrixXf(50,50)) )); + CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyRowsAsCols>(MatrixXcd(14,7)) )); + } + CALL_SUBTEST(( svd<MatrixXf,0>(MatrixXf(300,200)) )); + CALL_SUBTEST(( svd<MatrixXcd,AtLeastAsManyColsAsRows>(MatrixXcd(100,150)) )); + + CALL_SUBTEST(( svd_verify_assert<Matrix3f>() )); + CALL_SUBTEST(( svd_verify_assert<Matrix3d>() )); + CALL_SUBTEST(( svd_verify_assert<MatrixXf>() )); + CALL_SUBTEST(( svd_verify_assert<MatrixXd>() )); +} diff --git a/test/main.h b/test/main.h index e3866c68b..619fc9e06 100644 --- a/test/main.h +++ b/test/main.h @@ -30,6 +30,10 @@ #include <vector> #include <typeinfo> +#ifdef NDEBUG +#undef NDEBUG +#endif + #ifndef EIGEN_TEST_FUNC #error EIGEN_TEST_FUNC must be defined #endif @@ -153,6 +157,8 @@ namespace Eigen #define VERIFY_IS_APPROX_OR_LESS_THAN(a, b) VERIFY(test_ei_isApproxOrLessThan(a, b)) #define VERIFY_IS_NOT_APPROX_OR_LESS_THAN(a, b) VERIFY(!test_ei_isApproxOrLessThan(a, b)) +#define VERIFY_IS_UNITARY(a) VERIFY(test_isUnitary(a)) + #define CALL_SUBTEST(FUNC) do { \ g_test_stack.push_back(EI_PP_MAKE_STRING(FUNC)); \ FUNC; \ @@ -227,6 +233,12 @@ inline bool test_ei_isMuchSmallerThan(const MatrixBase<Derived>& m, return m.isMuchSmallerThan(s, test_precision<typename ei_traits<Derived>::Scalar>()); } +template<typename Derived> +inline bool test_isUnitary(const MatrixBase<Derived>& m) +{ + return m.isUnitary(test_precision<typename ei_traits<Derived>::Scalar>()); +} + template<typename MatrixType> void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType& m) { diff --git a/test/mixingtypes.cpp b/test/mixingtypes.cpp index d14232bd4..3e322c7fe 100644 --- a/test/mixingtypes.cpp +++ b/test/mixingtypes.cpp @@ -33,6 +33,7 @@ #include "main.h" +using namespace std; template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) { @@ -45,6 +46,104 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf; typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd; + Mat_f mf = Mat_f::Random(size,size); + Mat_d md = mf.template cast<double>(); + Mat_cf mcf = Mat_cf::Random(size,size); + Mat_cd mcd = mcf.template cast<complex<double> >(); + Vec_f vf = Vec_f::Random(size,1); + Vec_d vd = vf.template cast<double>(); + Vec_cf vcf = Vec_cf::Random(size,1); + Vec_cd vcd = vcf.template cast<complex<double> >(); + float sf = ei_random<float>(); + double sd = ei_random<double>(); + complex<float> scf = ei_random<complex<float> >(); + complex<double> scd = ei_random<complex<double> >(); + + + mf+mf; + VERIFY_RAISES_ASSERT(mf+md); + VERIFY_RAISES_ASSERT(mf+mcf); + VERIFY_RAISES_ASSERT(vf=vd); + VERIFY_RAISES_ASSERT(vf+=vd); + VERIFY_RAISES_ASSERT(mcd=md); + + // check scalar products + VERIFY_IS_APPROX(vcf * sf , vcf * complex<float>(sf)); + VERIFY_IS_APPROX(sd * vcd, complex<double>(sd) * vcd); + VERIFY_IS_APPROX(vf * scf , vf.template cast<complex<float> >() * scf); + VERIFY_IS_APPROX(scd * vd, scd * vd.template cast<complex<double> >()); + + // check dot product + vf.dot(vf); + VERIFY_RAISES_ASSERT(vd.dot(vf)); + VERIFY_RAISES_ASSERT(vcf.dot(vf)); // yeah eventually we should allow this but i'm too lazy to make that change now in Dot.h + // especially as that might be rewritten as cwise product .sum() which would make that automatic. + + // check diagonal product + VERIFY_IS_APPROX(vf.asDiagonal() * mcf, vf.template cast<complex<float> >().asDiagonal() * mcf); + VERIFY_IS_APPROX(vcd.asDiagonal() * md, vcd.asDiagonal() * md.template cast<complex<double> >()); + VERIFY_IS_APPROX(mcf * vf.asDiagonal(), mcf * vf.template cast<complex<float> >().asDiagonal()); + VERIFY_IS_APPROX(md * vcd.asDiagonal(), md.template cast<complex<double> >() * vcd.asDiagonal()); +// vd.asDiagonal() * mf; // does not even compile +// vcd.asDiagonal() * mf; // does not even compile + + // check inner product + VERIFY_IS_APPROX((vf.transpose() * vcf).value(), (vf.template cast<complex<float> >().transpose() * vcf).value()); + + // check outer product + VERIFY_IS_APPROX((vf * vcf.transpose()).eval(), (vf.template cast<complex<float> >() * vcf.transpose()).eval()); +} + + +void mixingtypes_large(int size) +{ + static const int SizeAtCompileType = Dynamic; + typedef Matrix<float, SizeAtCompileType, SizeAtCompileType> Mat_f; + typedef Matrix<double, SizeAtCompileType, SizeAtCompileType> Mat_d; + typedef Matrix<std::complex<float>, SizeAtCompileType, SizeAtCompileType> Mat_cf; + typedef Matrix<std::complex<double>, SizeAtCompileType, SizeAtCompileType> Mat_cd; + typedef Matrix<float, SizeAtCompileType, 1> Vec_f; + typedef Matrix<double, SizeAtCompileType, 1> Vec_d; + typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf; + typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd; + + Mat_f mf(size,size); + Mat_d md(size,size); + Mat_cf mcf(size,size); + Mat_cd mcd(size,size); + Vec_f vf(size,1); + Vec_d vd(size,1); + Vec_cf vcf(size,1); + Vec_cd vcd(size,1); + + mf*mf; + // FIXME large products does not allow mixing types + VERIFY_RAISES_ASSERT(md*mcd); + VERIFY_RAISES_ASSERT(mcd*md); + VERIFY_RAISES_ASSERT(mf*vcf); + VERIFY_RAISES_ASSERT(mcf*vf); + VERIFY_RAISES_ASSERT(mcf *= mf); + // VERIFY_RAISES_ASSERT(vcd = md*vcd); // does not even compile (cannot convert complex to double) + VERIFY_RAISES_ASSERT(vcf = mcf*vf); + +// VERIFY_RAISES_ASSERT(mf*md); // does not even compile +// VERIFY_RAISES_ASSERT(mcf*mcd); // does not even compile +// VERIFY_RAISES_ASSERT(mcf*vcd); // does not even compile + VERIFY_RAISES_ASSERT(vcf = mf*vf); +} + +template<int SizeAtCompileType> void mixingtypes_small() +{ + int size = SizeAtCompileType; + typedef Matrix<float, SizeAtCompileType, SizeAtCompileType> Mat_f; + typedef Matrix<double, SizeAtCompileType, SizeAtCompileType> Mat_d; + typedef Matrix<std::complex<float>, SizeAtCompileType, SizeAtCompileType> Mat_cf; + typedef Matrix<std::complex<double>, SizeAtCompileType, SizeAtCompileType> Mat_cd; + typedef Matrix<float, SizeAtCompileType, 1> Vec_f; + typedef Matrix<double, SizeAtCompileType, 1> Vec_d; + typedef Matrix<std::complex<float>, SizeAtCompileType, 1> Vec_cf; + typedef Matrix<std::complex<double>, SizeAtCompileType, 1> Vec_cd; + Mat_f mf(size,size); Mat_d md(size,size); Mat_cf mcf(size,size); @@ -54,14 +153,11 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) Vec_cf vcf(size,1); Vec_cd vcd(size,1); - mf+mf; - VERIFY_RAISES_ASSERT(mf+md); - VERIFY_RAISES_ASSERT(mf+mcf); - VERIFY_RAISES_ASSERT(vf=vd); - VERIFY_RAISES_ASSERT(vf+=vd); - VERIFY_RAISES_ASSERT(mcd=md); mf*mf; + // FIXME shall we discard those products ? + // 1) currently they work only if SizeAtCompileType is small enough + // 2) in case we vectorize complexes this might be difficult to still allow that md*mcd; mcd*md; mf*vcf; @@ -69,20 +165,19 @@ template<int SizeAtCompileType> void mixingtypes(int size = SizeAtCompileType) mcf *= mf; vcd = md*vcd; vcf = mcf*vf; - VERIFY_RAISES_ASSERT(mf*md); - VERIFY_RAISES_ASSERT(mcf*mcd); - VERIFY_RAISES_ASSERT(mcf*vcd); +// VERIFY_RAISES_ASSERT(mf*md); // does not even compile +// VERIFY_RAISES_ASSERT(mcf*mcd); // does not even compile +// VERIFY_RAISES_ASSERT(mcf*vcd); // does not even compile VERIFY_RAISES_ASSERT(vcf = mf*vf); - - vf.dot(vf); - VERIFY_RAISES_ASSERT(vd.dot(vf)); - VERIFY_RAISES_ASSERT(vcf.dot(vf)); // yeah eventually we should allow this but i'm too lazy to make that change now in Dot.h -} // especially as that might be rewritten as cwise product .sum() which would make that automatic. +} void test_mixingtypes() { // check that our operator new is indeed called: - CALL_SUBTEST(mixingtypes<3>(3)); - CALL_SUBTEST(mixingtypes<4>(4)); + CALL_SUBTEST(mixingtypes<3>()); + CALL_SUBTEST(mixingtypes<4>()); CALL_SUBTEST(mixingtypes<Dynamic>(20)); + + CALL_SUBTEST(mixingtypes_small<4>()); + CALL_SUBTEST(mixingtypes_large(20)); } diff --git a/test/product.h b/test/product.h index 157f6262b..40773ad90 100644 --- a/test/product.h +++ b/test/product.h @@ -81,7 +81,7 @@ template<typename MatrixType> void product(const MatrixType& m) m3 = m1; m3 *= m1.transpose() * m2; VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2)); - VERIFY_IS_APPROX(m3, m1.lazy() * (m1.transpose()*m2)); + VERIFY_IS_APPROX(m3, m1 * (m1.transpose()*m2)); // continue testing Product.h: distributivity VERIFY_IS_APPROX(square*(m1 + m2), square*m1+square*m2); @@ -109,26 +109,26 @@ template<typename MatrixType> void product(const MatrixType& m) // test optimized operator+= path res = square; - res += (m1 * m2.transpose()).lazy(); + res.noalias() += m1 * m2.transpose(); VERIFY_IS_APPROX(res, square + m1 * m2.transpose()); if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1) { VERIFY(areNotApprox(res,square + m2 * m1.transpose())); } vcres = vc2; - vcres += (m1.transpose() * v1).lazy(); + vcres.noalias() += m1.transpose() * v1; VERIFY_IS_APPROX(vcres, vc2 + m1.transpose() * v1); // test optimized operator-= path res = square; - res -= (m1 * m2.transpose()).lazy(); + res.noalias() -= m1 * m2.transpose(); VERIFY_IS_APPROX(res, square - (m1 * m2.transpose())); if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1) { VERIFY(areNotApprox(res,square - m2 * m1.transpose())); } vcres = vc2; - vcres -= (m1.transpose() * v1).lazy(); + vcres.noalias() -= m1.transpose() * v1; VERIFY_IS_APPROX(vcres, vc2 - m1.transpose() * v1); tm1 = m1; @@ -145,7 +145,7 @@ template<typename MatrixType> void product(const MatrixType& m) VERIFY_IS_APPROX(res, m1 * m2.transpose()); res2 = square2; - res2 += (m1.transpose() * m2).lazy(); + res2.noalias() += m1.transpose() * m2; VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2); if (NumTraits<Scalar>::HasFloatingPoint && std::min(rows,cols)>1) { diff --git a/test/product_extra.cpp b/test/product_extra.cpp index 213dbced6..3ad99fc7a 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp @@ -59,16 +59,13 @@ template<typename MatrixType> void product_extra(const MatrixType& m) // r0 = ei_random<int>(0,rows/2-1), // r1 = ei_random<int>(rows/2,rows); - // all the expressions in this test should be compiled as a single matrix product - // TODO: add internal checks to verify that - - VERIFY_IS_APPROX(m3 = (m1 * m2.adjoint()).lazy(), m1 * m2.adjoint().eval()); - VERIFY_IS_APPROX(m3 = (m1.adjoint() * square.adjoint()).lazy(), m1.adjoint().eval() * square.adjoint().eval()); - VERIFY_IS_APPROX(m3 = (m1.adjoint() * m2).lazy(), m1.adjoint().eval() * m2); - VERIFY_IS_APPROX(m3 = ((s1 * m1.adjoint()) * m2).lazy(), (s1 * m1.adjoint()).eval() * m2); - VERIFY_IS_APPROX(m3 = ((- m1.adjoint() * s1) * (s3 * m2)).lazy(), (- m1.adjoint() * s1).eval() * (s3 * m2).eval()); - VERIFY_IS_APPROX(m3 = ((s2 * m1.adjoint() * s1) * m2).lazy(), (s2 * m1.adjoint() * s1).eval() * m2); - VERIFY_IS_APPROX(m3 = ((-m1*s2) * s1*m2.adjoint()).lazy(), (-m1*s2).eval() * (s1*m2.adjoint()).eval()); + VERIFY_IS_APPROX(m3.noalias() = m1 * m2.adjoint(), m1 * m2.adjoint().eval()); + VERIFY_IS_APPROX(m3.noalias() = m1.adjoint() * square.adjoint(), m1.adjoint().eval() * square.adjoint().eval()); + VERIFY_IS_APPROX(m3.noalias() = m1.adjoint() * m2, m1.adjoint().eval() * m2); + VERIFY_IS_APPROX(m3.noalias() = (s1 * m1.adjoint()) * m2, (s1 * m1.adjoint()).eval() * m2); + VERIFY_IS_APPROX(m3.noalias() = (- m1.adjoint() * s1) * (s3 * m2), (- m1.adjoint() * s1).eval() * (s3 * m2).eval()); + VERIFY_IS_APPROX(m3.noalias() = (s2 * m1.adjoint() * s1) * m2, (s2 * m1.adjoint() * s1).eval() * m2); + VERIFY_IS_APPROX(m3.noalias() = (-m1*s2) * s1*m2.adjoint(), (-m1*s2).eval() * (s1*m2.adjoint()).eval()); // a very tricky case where a scale factor has to be automatically conjugated: VERIFY_IS_APPROX( m1.adjoint() * (s1*m2).conjugate(), (m1.adjoint()).eval() * ((s1*m2).conjugate()).eval()); @@ -76,7 +73,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) // test all possible conjugate combinations for the four matrix-vector product cases: -// std::cerr << "a\n"; VERIFY_IS_APPROX((-m1.conjugate() * s2) * (s1 * vc2), (-m1.conjugate()*s2).eval() * (s1 * vc2).eval()); VERIFY_IS_APPROX((-m1 * s2) * (s1 * vc2.conjugate()), @@ -84,7 +80,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((-m1.conjugate() * s2) * (s1 * vc2.conjugate()), (-m1.conjugate()*s2).eval() * (s1 * vc2.conjugate()).eval()); -// std::cerr << "b\n"; VERIFY_IS_APPROX((s1 * vc2.transpose()) * (-m1.adjoint() * s2), (s1 * vc2.transpose()).eval() * (-m1.adjoint()*s2).eval()); VERIFY_IS_APPROX((s1 * vc2.adjoint()) * (-m1.transpose() * s2), @@ -92,7 +87,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((s1 * vc2.adjoint()) * (-m1.adjoint() * s2), (s1 * vc2.adjoint()).eval() * (-m1.adjoint()*s2).eval()); -// std::cerr << "c\n"; VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.transpose()), (-m1.adjoint()*s2).eval() * (s1 * v1.transpose()).eval()); VERIFY_IS_APPROX((-m1.transpose() * s2) * (s1 * v1.adjoint()), @@ -100,7 +94,6 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()), (-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval()); -// std::cerr << "d\n"; VERIFY_IS_APPROX((s1 * v1) * (-m1.conjugate() * s2), (s1 * v1).eval() * (-m1.conjugate()*s2).eval()); VERIFY_IS_APPROX((s1 * v1.conjugate()) * (-m1 * s2), @@ -111,13 +104,24 @@ template<typename MatrixType> void product_extra(const MatrixType& m) VERIFY_IS_APPROX((-m1.adjoint() * s2) * (s1 * v1.adjoint()), (-m1.adjoint()*s2).eval() * (s1 * v1.adjoint()).eval()); + // test the vector-matrix product with non aligned starts + int i = ei_random<int>(0,m1.rows()-2); + int j = ei_random<int>(0,m1.cols()-2); + int r = ei_random<int>(1,m1.rows()-i); + int c = ei_random<int>(1,m1.cols()-j); + int i2 = ei_random<int>(0,m1.rows()-1); + int j2 = ei_random<int>(0,m1.cols()-1); + + VERIFY_IS_APPROX(m1.col(j2).adjoint() * m1.block(0,j,m1.rows(),c), m1.col(j2).adjoint().eval() * m1.block(0,j,m1.rows(),c).eval()); + VERIFY_IS_APPROX(m1.block(i,0,r,m1.cols()) * m1.row(i2).adjoint(), m1.block(i,0,r,m1.cols()).eval() * m1.row(i2).adjoint().eval()); + } void test_product_extra() { for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(1,320), ei_random<int>(1,320))) ); + CALL_SUBTEST( product_extra(MatrixXf(ei_random<int>(2,320), ei_random<int>(2,320))) ); CALL_SUBTEST( product_extra(MatrixXcf(ei_random<int>(50,50), ei_random<int>(50,50))) ); - CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(1,50), ei_random<int>(1,50))) ); + CALL_SUBTEST( product_extra(Matrix<std::complex<double>,Dynamic,Dynamic,RowMajor>(ei_random<int>(2,50), ei_random<int>(2,50))) ); } } diff --git a/test/product_notemporary.cpp b/test/product_notemporary.cpp index f6e8f4101..1a3d65291 100644 --- a/test/product_notemporary.cpp +++ b/test/product_notemporary.cpp @@ -69,53 +69,47 @@ template<typename MatrixType> void product_notemporary(const MatrixType& m) r1 = ei_random<int>(8,rows-r0); VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()), 1); - VERIFY_EVALUATION_COUNT( m3 = (m1 * m2.adjoint()).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3.noalias() = m1 * m2.adjoint(), 0); - // NOTE in this case the slow product is used: - // FIXME: VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * (m1 * m2.transpose()), 0); - VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * m2.adjoint()).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3 = (s1 * m1 * s2 * (m1*s3+m2*s2).adjoint()).lazy(), 1); - VERIFY_EVALUATION_COUNT( m3 = ((s1 * m1).adjoint() * s2 * m2).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3 += (s1 * (-m1*s3).adjoint() * (s2 * m2 * s3)).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * m2.adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = s1 * m1 * s2 * (m1*s3+m2*s2).adjoint(), 1); + VERIFY_EVALUATION_COUNT( m3.noalias() = (s1 * m1).adjoint() * s2 * m2, 0); VERIFY_EVALUATION_COUNT( m3.noalias() += s1 * (-m1*s3).adjoint() * (s2 * m2 * s3), 0); - VERIFY_EVALUATION_COUNT( m3 -= (s1 * (m1.transpose() * m2)).lazy(), 0); VERIFY_EVALUATION_COUNT( m3.noalias() -= s1 * (m1.transpose() * m2), 0); - VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) += (-m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint()).lazy() ), 0); - VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) -= (s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1)).lazy() ), 0); + VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() += -m1.block(r0,c0,r1,c1) * (s2*m2.block(r0,c0,r1,c1)).adjoint() ), 0); + VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() -= s1 * m1.block(r0,c0,r1,c1) * m2.block(c0,r0,c1,r1) ), 0); // NOTE this is because the Block expression is not handled yet by our expression analyser - VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1) = (s1 * m1.block(r0,c0,r1,c1) * (s1*m2).block(c0,r0,c1,r1)).lazy() ), 1); + VERIFY_EVALUATION_COUNT(( m3.block(r0,r0,r1,r1).noalias() = s1 * m1.block(r0,c0,r1,c1) * (s1*m2).block(c0,r0,c1,r1) ), 1); - VERIFY_EVALUATION_COUNT( m3 -= ((s1 * m1).template triangularView<LowerTriangular>() * m2).lazy(), 0); - VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template triangularView<UpperTriangular>() * (m2+m2)).lazy(), 1); - VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * m2.adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).template triangularView<LowerTriangular>() * m2, 0); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UpperTriangular>() * (m2+m2), 1); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * m2.adjoint(), 0); - VERIFY_EVALUATION_COUNT( rm3.col(c0) = ((s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * (s2*m2.row(c0)).adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( rm3.col(c0).noalias() = (s1 * m1.adjoint()).template triangularView<UnitUpperTriangular>() * (s2*m2.row(c0)).adjoint(), 0); VERIFY_EVALUATION_COUNT( m1.template triangularView<LowerTriangular>().solveInPlace(m3), 0); VERIFY_EVALUATION_COUNT( m1.adjoint().template triangularView<LowerTriangular>().solveInPlace(m3.transpose()), 0); - VERIFY_EVALUATION_COUNT( m3 -= ((s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2*s3).adjoint()).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3 = (s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<UpperTriangular>()).lazy(), 0); - VERIFY_EVALUATION_COUNT( rm3 = ((s1 * m1.adjoint()).template selfadjointView<LowerTriangular>() * m2.adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() -= (s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2*s3).adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = s2 * m2.adjoint() * (s1 * m1.adjoint()).template selfadjointView<UpperTriangular>(), 0); + VERIFY_EVALUATION_COUNT( rm3.noalias() = (s1 * m1.adjoint()).template selfadjointView<LowerTriangular>() * m2.adjoint(), 0); - VERIFY_EVALUATION_COUNT( m3.col(c0) = ((s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2.row(c0)*s3).adjoint()).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3.col(c0) -= ((s1 * m1).adjoint().template selfadjointView<UpperTriangular>() * (-m2.row(c0)*s3).adjoint()).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() = (s1 * m1).adjoint().template selfadjointView<LowerTriangular>() * (-m2.row(c0)*s3).adjoint(), 0); + VERIFY_EVALUATION_COUNT( m3.col(c0).noalias() -= (s1 * m1).adjoint().template selfadjointView<UpperTriangular>() * (-m2.row(c0)*s3).adjoint(), 0); - VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1) += ((m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * (s1*m2.block(r0,c0,r1,c1)) )).lazy(), 0); - VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1) = ((m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() += m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * (s1*m2.block(r0,c0,r1,c1)), 0); + VERIFY_EVALUATION_COUNT( m3.block(r0,c0,r1,c1).noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<UpperTriangular>() * m2.block(r0,c0,r1,c1), 0); VERIFY_EVALUATION_COUNT( m3.template selfadjointView<LowerTriangular>().rankUpdate(m2.adjoint()), 0); m3.resize(1,1); - VERIFY_EVALUATION_COUNT( m3 = ((m1.block(r0,r0,r1,r1).template selfadjointView<LowerTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template selfadjointView<LowerTriangular>() * m2.block(r0,c0,r1,c1), 0); m3.resize(1,1); - VERIFY_EVALUATION_COUNT( m3 = ((m1.block(r0,r0,r1,r1).template triangularView<UnitUpperTriangular>() * m2.block(r0,c0,r1,c1) )).lazy(), 0); + VERIFY_EVALUATION_COUNT( m3.noalias() = m1.block(r0,r0,r1,r1).template triangularView<UnitUpperTriangular>() * m2.block(r0,c0,r1,c1), 0); } void test_product_notemporary() diff --git a/test/product_symm.cpp b/test/product_symm.cpp index 1300928a2..cf0299c64 100644 --- a/test/product_symm.cpp +++ b/test/product_symm.cpp @@ -96,7 +96,7 @@ template<typename Scalar, int Size, int OtherSize> void symm(int size = Size, in m2 = m1.template triangularView<UpperTriangular>(); rhs13 = rhs12; - VERIFY_IS_APPROX(rhs12 += (s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate())).lazy(), + VERIFY_IS_APPROX(rhs12.noalias() += s1 * ((m2.adjoint()).template selfadjointView<LowerTriangular>() * (s2*rhs3).conjugate()), rhs13 += (s1*m1.adjoint()) * (s2*rhs3).conjugate()); // test matrix * selfadjoint diff --git a/test/product_trsm.cpp b/test/product_trsm.cpp index 9f372df91..756034df9 100644 --- a/test/product_trsm.cpp +++ b/test/product_trsm.cpp @@ -40,8 +40,8 @@ template<typename Scalar> void trsm(int size,int cols) Matrix<Scalar,Dynamic,Dynamic,ColMajor> cmRhs(size,cols), ref(size,cols); Matrix<Scalar,Dynamic,Dynamic,RowMajor> rmRhs(size,cols); - cmLhs.setRandom(); cmLhs *= 0.1; cmLhs.diagonal().cwise() += 1; - rmLhs.setRandom(); rmLhs *= 0.1; rmLhs.diagonal().cwise() += 1; + cmLhs.setRandom(); cmLhs *= static_cast<RealScalar>(0.1); cmLhs.diagonal().cwise() += static_cast<RealScalar>(1); + rmLhs.setRandom(); rmLhs *= static_cast<RealScalar>(0.1); rmLhs.diagonal().cwise() += static_cast<RealScalar>(1); VERIFY_TRSM(cmLhs.conjugate().template triangularView<LowerTriangular>(), cmRhs); VERIFY_TRSM(cmLhs .template triangularView<UpperTriangular>(), cmRhs); diff --git a/test/stable_norm.cpp b/test/stable_norm.cpp new file mode 100644 index 000000000..b8fbf5271 --- /dev/null +++ b/test/stable_norm.cpp @@ -0,0 +1,79 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" + +template<typename MatrixType> void stable_norm(const MatrixType& m) +{ + /* this test covers the following files: + StableNorm.h + */ + + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits<Scalar>::Real RealScalar; + + int rows = m.rows(); + int cols = m.cols(); + + Scalar big = ei_random<Scalar>() * std::numeric_limits<RealScalar>::max() * RealScalar(1e-4); + Scalar small = static_cast<RealScalar>(1)/big; + + MatrixType vzero = MatrixType::Zero(rows, cols), + vrand = MatrixType::Random(rows, cols), + vbig(rows, cols), + vsmall(rows,cols); + + vbig.fill(big); + vsmall.fill(small); + + VERIFY_IS_MUCH_SMALLER_THAN(vzero.norm(), static_cast<RealScalar>(1)); + VERIFY_IS_APPROX(vrand.stableNorm(), vrand.norm()); + VERIFY_IS_APPROX(vrand.blueNorm(), vrand.norm()); + VERIFY_IS_APPROX(vrand.hypotNorm(), vrand.norm()); + + RealScalar size = static_cast<RealScalar>(m.size()); + + // test overflow + VERIFY_IS_NOT_APPROX(static_cast<Scalar>(vbig.norm()), ei_sqrt(size)*big); // here the default norm must fail + VERIFY_IS_APPROX(static_cast<Scalar>(vbig.stableNorm()), ei_sqrt(size)*big); + VERIFY_IS_APPROX(static_cast<Scalar>(vbig.blueNorm()), ei_sqrt(size)*big); + VERIFY_IS_APPROX(static_cast<Scalar>(vbig.hypotNorm()), ei_sqrt(size)*big); + + // test underflow + VERIFY_IS_NOT_APPROX(static_cast<Scalar>(vsmall.norm()), ei_sqrt(size)*small); // here the default norm must fail + VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.stableNorm()), ei_sqrt(size)*small); + VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.blueNorm()), ei_sqrt(size)*small); + VERIFY_IS_APPROX(static_cast<Scalar>(vsmall.hypotNorm()), ei_sqrt(size)*small); +} + +void test_stable_norm() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( stable_norm(Matrix<float, 1, 1>()) ); + CALL_SUBTEST( stable_norm(Vector4d()) ); + CALL_SUBTEST( stable_norm(VectorXd(ei_random<int>(10,2000))) ); + CALL_SUBTEST( stable_norm(VectorXf(ei_random<int>(10,2000))) ); + CALL_SUBTEST( stable_norm(VectorXcd(ei_random<int>(10,2000))) ); + } +} diff --git a/test/submatrices.cpp b/test/submatrices.cpp index a819cadc2..6fe86c281 100644 --- a/test/submatrices.cpp +++ b/test/submatrices.cpp @@ -170,6 +170,48 @@ template<typename MatrixType> void submatrices(const MatrixType& m) VERIFY(ei_real(ones.row(r1).dot(ones.row(r2))) == RealScalar(cols)); } + +template<typename MatrixType> +void compare_using_data_and_stride(const MatrixType& m) +{ + int rows = m.rows(); + int cols = m.cols(); + int size = m.size(); + int stride = m.stride(); + const typename MatrixType::Scalar* data = m.data(); + + for(int j=0;j<cols;++j) + for(int i=0;i<rows;++i) + VERIFY_IS_APPROX(m.coeff(i,j), data[(MatrixType::Flags&RowMajorBit) ? i*stride+j : j*stride + i]); + + if(MatrixType::IsVectorAtCompileTime) + { + VERIFY_IS_APPROX(stride, int((&m.coeff(1))-(&m.coeff(0)))); + for (int i=0;i<size;++i) + VERIFY_IS_APPROX(m.coeff(i), data[i*stride]); + } +} + +template<typename MatrixType> +void data_and_stride(const MatrixType& m) +{ + int rows = m.rows(); + int cols = m.cols(); + + int r1 = ei_random<int>(0,rows-1); + int r2 = ei_random<int>(r1,rows-1); + int c1 = ei_random<int>(0,cols-1); + int c2 = ei_random<int>(c1,cols-1); + + MatrixType m1 = MatrixType::Random(rows, cols); + compare_using_data_and_stride(m1.block(r1, c1, r2-r1+1, c2-c1+1)); + compare_using_data_and_stride(m1.transpose().block(c1, r1, c2-c1+1, r2-r1+1)); + compare_using_data_and_stride(m1.row(r1)); + compare_using_data_and_stride(m1.col(c1)); + compare_using_data_and_stride(m1.row(r1).transpose()); + compare_using_data_and_stride(m1.col(c1).transpose()); +} + void test_submatrices() { for(int i = 0; i < g_repeat; i++) { @@ -179,5 +221,8 @@ void test_submatrices() CALL_SUBTEST( submatrices(MatrixXi(8, 12)) ); CALL_SUBTEST( submatrices(MatrixXcd(20, 20)) ); CALL_SUBTEST( submatrices(MatrixXf(20, 20)) ); + + CALL_SUBTEST( data_and_stride(MatrixXf(ei_random(5,50), ei_random(5,50))) ); + CALL_SUBTEST( data_and_stride(Matrix<int,Dynamic,Dynamic,RowMajor>(ei_random(5,50), ei_random(5,50))) ); } } diff --git a/test/svd.cpp b/test/svd.cpp index 2ccd94764..e6a32bd3f 100644 --- a/test/svd.cpp +++ b/test/svd.cpp @@ -41,15 +41,11 @@ template<typename MatrixType> void svd(const MatrixType& m) Matrix<Scalar, MatrixType::RowsAtCompileTime, 1>::Random(rows,1); Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> x(cols,1), x2(cols,1); - RealScalar largerEps = test_precision<RealScalar>(); - if (ei_is_same_type<RealScalar,float>::ret) - largerEps = 1e-3f; - { SVD<MatrixType> svd(a); MatrixType sigma = MatrixType::Zero(rows,cols); MatrixType matU = MatrixType::Zero(rows,rows); - sigma.block(0,0,cols,cols) = svd.singularValues().asDiagonal(); + sigma.diagonal() = svd.singularValues(); matU = svd.matrixU(); VERIFY_IS_APPROX(a, matU * sigma * svd.matrixV().transpose()); } diff --git a/test/swap.cpp b/test/swap.cpp new file mode 100644 index 000000000..8b325992c --- /dev/null +++ b/test/swap.cpp @@ -0,0 +1,98 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#define EIGEN_NO_STATIC_ASSERT +#include "main.h" + +template<typename T> +struct other_matrix_type +{ + typedef int type; +}; + +template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int _MaxCols> +struct other_matrix_type<Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> > +{ + typedef Matrix<_Scalar, _Rows, _Cols, _Options^RowMajor, _MaxRows, _MaxCols> type; +}; + +template<typename MatrixType> void swap(const MatrixType& m) +{ + typedef typename other_matrix_type<MatrixType>::type OtherMatrixType; + typedef typename MatrixType::Scalar Scalar; + + ei_assert((!ei_is_same_type<MatrixType,OtherMatrixType>::ret)); + int rows = m.rows(); + int cols = m.cols(); + + // construct 3 matrix guaranteed to be distinct + MatrixType m1 = MatrixType::Random(rows,cols); + MatrixType m2 = MatrixType::Random(rows,cols) + Scalar(100) * MatrixType::Identity(rows,cols); + OtherMatrixType m3 = OtherMatrixType::Random(rows,cols) + Scalar(200) * OtherMatrixType::Identity(rows,cols); + + MatrixType m1_copy = m1; + MatrixType m2_copy = m2; + OtherMatrixType m3_copy = m3; + + // test swapping 2 matrices of same type + m1.swap(m2); + VERIFY_IS_APPROX(m1,m2_copy); + VERIFY_IS_APPROX(m2,m1_copy); + m1 = m1_copy; + m2 = m2_copy; + + // test swapping 2 matrices of different types + m1.swap(m3); + VERIFY_IS_APPROX(m1,m3_copy); + VERIFY_IS_APPROX(m3,m1_copy); + m1 = m1_copy; + m3 = m3_copy; + + // test swapping matrix with expression + m1.swap(m2.block(0,0,rows,cols)); + VERIFY_IS_APPROX(m1,m2_copy); + VERIFY_IS_APPROX(m2,m1_copy); + m1 = m1_copy; + m2 = m2_copy; + + // test swapping two expressions of different types + m1.transpose().swap(m3.transpose()); + VERIFY_IS_APPROX(m1,m3_copy); + VERIFY_IS_APPROX(m3,m1_copy); + m1 = m1_copy; + m3 = m3_copy; + + // test assertion on mismatching size -- matrix case + VERIFY_RAISES_ASSERT(m1.swap(m1.row(0))); + // test assertion on mismatching size -- xpr case + VERIFY_RAISES_ASSERT(m1.row(0).swap(m1)); +} + +void test_swap() +{ + CALL_SUBTEST( swap(Matrix3f()) ); // fixed size, no vectorization + CALL_SUBTEST( swap(Matrix4d()) ); // fixed size, possible vectorization + CALL_SUBTEST( swap(MatrixXd(3,3)) ); // dyn size, no vectorization + CALL_SUBTEST( swap(MatrixXf(30,30)) ); // dyn size, possible vectorization +} diff --git a/test/testsuite.cmake b/test/testsuite.cmake index 6a44ce239..37ee87565 100644 --- a/test/testsuite.cmake +++ b/test/testsuite.cmake @@ -27,6 +27,13 @@ # - EIGEN_WORK_DIR: directory used to download the source files and make the builds # default: folder which contains this script # - EIGEN_CMAKE_ARGS: additional arguments passed to cmake +# - EIGEN_GENERATOR_TYPE: allows to overwrite the generator type +# default: nmake (windows +# See http://www.cmake.org/cmake/help/cmake2.6docs.html#section_Generators for a complete +# list of supported generators. +# - EIGEN_NO_UPDATE: allows to submit dash boards from local repositories +# This might be interesting in case you want to submit dashboards +# including local changes. # - CTEST_SOURCE_DIRECTORY: path to eigen's src (use a new and empty folder, not the one you are working on) # default: <EIGEN_WORK_DIR>/src # - CTEST_BINARY_DIRECTORY: build directory @@ -132,11 +139,11 @@ endif(NOT EIGEN_MODE) ## mandatory variables (the default should be ok in most cases): -if(NOT IGNORE_CVS) +if(NOT EIGEN_NO_UPDATE) SET (CTEST_CVS_COMMAND "hg") SET (CTEST_CVS_CHECKOUT "${CTEST_CVS_COMMAND} clone http://bitbucket.org/eigen/eigen2 \"${CTEST_SOURCE_DIRECTORY}\"") SET(CTEST_BACKUP_AND_RESTORE TRUE) # the backup is CVS related ... -endif(NOT IGNORE_CVS) +endif(NOT EIGEN_NO_UPDATE) # which ctest command to use for running the dashboard SET (CTEST_COMMAND "${EIGEN_CMAKE_DIR}ctest -D ${EIGEN_MODE}") @@ -158,15 +165,24 @@ SET(CTEST_START_WITH_EMPTY_BINARY_DIRECTORY TRUE) # any quotes inside of this string if you use it if(WIN32 AND NOT UNIX) #message(SEND_ERROR "win32") - set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"NMake Makefiles\" -DCMAKE_MAKE_PROGRAM=nmake") - SET (CTEST_INITIAL_CACHE " - MAKECOMMAND:STRING=nmake -i - CMAKE_MAKE_PROGRAM:FILEPATH=nmake - CMAKE_GENERATOR:INTERNAL=NMake Makefiles - CMAKE_BUILD_TYPE:STRING=Release - BUILDNAME:STRING=${EIGEN_BUILD_STRING} - SITE:STRING=${EIGEN_SITE} - ") + if(EIGEN_GENERATOR_TYPE) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"${EIGEN_GENERATOR_TYPE}\"") + SET (CTEST_INITIAL_CACHE " + CMAKE_BUILD_TYPE:STRING=Release + BUILDNAME:STRING=${EIGEN_BUILD_STRING} + SITE:STRING=${EIGEN_SITE} + ") + else(EIGEN_GENERATOR_TYPE) + set(CTEST_CMAKE_COMMAND "${CTEST_CMAKE_COMMAND} -G \"NMake Makefiles\" -DCMAKE_MAKE_PROGRAM=nmake") + SET (CTEST_INITIAL_CACHE " + MAKECOMMAND:STRING=nmake -i + CMAKE_MAKE_PROGRAM:FILEPATH=nmake + CMAKE_GENERATOR:INTERNAL=NMake Makefiles + CMAKE_BUILD_TYPE:STRING=Release + BUILDNAME:STRING=${EIGEN_BUILD_STRING} + SITE:STRING=${EIGEN_SITE} + ") + endif(EIGEN_GENERATOR_TYPE) else(WIN32 AND NOT UNIX) SET (CTEST_INITIAL_CACHE " BUILDNAME:STRING=${EIGEN_BUILD_STRING} diff --git a/test/umeyama.cpp b/test/umeyama.cpp index b6d394883..0999c59c9 100644 --- a/test/umeyama.cpp +++ b/test/umeyama.cpp @@ -127,7 +127,7 @@ void run_test(int dim, int num_elements) MatrixX src = MatrixX::Random(dim+1, num_elements); src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1)); - MatrixX dst = (cR_t*src).lazy(); + MatrixX dst = cR_t*src; MatrixX cR_t_umeyama = umeyama(src.block(0,0,dim,num_elements), dst.block(0,0,dim,num_elements)); @@ -160,7 +160,7 @@ void run_fixed_size_test(int num_elements) MatrixX src = MatrixX::Random(dim+1, num_elements); src.row(dim) = Matrix<Scalar, 1, Dynamic>::Constant(num_elements, Scalar(1)); - MatrixX dst = (cR_t*src).lazy(); + MatrixX dst = cR_t*src; Block<MatrixX, Dimension, Dynamic> src_block(src,0,0,dim,num_elements); Block<MatrixX, Dimension, Dynamic> dst_block(dst,0,0,dim,num_elements); diff --git a/unsupported/Eigen/BVH b/unsupported/Eigen/BVH index 7b9c3c7c6..f307da2f7 100644 --- a/unsupported/Eigen/BVH +++ b/unsupported/Eigen/BVH @@ -94,8 +94,6 @@ namespace Eigen { * responsibility of the intersectObject function to keep track of the results in whatever manner is appropriate. * The cartesian product intersection and the BVMinimize queries are similar--see their individual documentation. * - * \addexample BVH_Example \label How to use a BVH to find the closest pair between two point sets - * * The following is a simple but complete example for how to use the BVH to accelerate the search for a closest red-blue point pair: * \include BVH_Example.cpp * Output: \verbinclude BVH_Example.out diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index dd25d7f3d..36d13b7eb 100644 --- a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h +++ b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h @@ -25,6 +25,10 @@ #ifndef EIGEN_MATRIX_EXPONENTIAL #define EIGEN_MATRIX_EXPONENTIAL +#ifdef _MSC_VER + template <typename Scalar> Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); } +#endif + /** \brief Compute the matrix exponential. * * \param M matrix whose exponential is to be computed. @@ -61,257 +65,243 @@ template <typename Derived> EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase<Derived> &M, typename MatrixBase<Derived>::PlainMatrixType* result); +/** \brief Class for computing the matrix exponential.*/ +template <typename MatrixType> +class MatrixExponential { -/** \internal \brief Internal helper functions for computing the - * matrix exponential. - */ -namespace MatrixExponentialInternal { + public: + + /** \brief Compute the matrix exponential. + * + * \param M matrix whose exponential is to be computed. + * \param result pointer to the matrix in which to store the result. + */ + MatrixExponential(const MatrixType &M, MatrixType *result); -#ifdef _MSC_VER - template <typename Scalar> Scalar log2(Scalar v) { return std::log(v)/std::log(Scalar(2)); } -#endif + private: + + // Prevent copying + MatrixExponential(const MatrixExponential&); + MatrixExponential& operator=(const MatrixExponential&); + + /** \brief Compute the (3,3)-Padé approximant to the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. + * + * \param A Argument of matrix exponential + */ + void pade3(const MatrixType &A); + + /** \brief Compute the (5,5)-Padé approximant to the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. + * + * \param A Argument of matrix exponential + */ + void pade5(const MatrixType &A); + + /** \brief Compute the (7,7)-Padé approximant to the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. + * + * \param A Argument of matrix exponential + */ + void pade7(const MatrixType &A); + + /** \brief Compute the (9,9)-Padé approximant to the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. + * + * \param A Argument of matrix exponential + */ + void pade9(const MatrixType &A); + + /** \brief Compute the (13,13)-Padé approximant to the exponential. + * + * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé + * approximant of \f$ \exp(A) \f$ around \f$ A = 0 \f$. + * + * \param A Argument of matrix exponential + */ + void pade13(const MatrixType &A); + + /** \brief Compute Padé approximant to the exponential. + * + * Computes \c m_U, \c m_V and \c m_squarings such that + * \f$ (V+U)(V-U)^{-1} \f$ is a Padé of + * \f$ \exp(2^{-\mbox{squarings}}M) \f$ around \f$ M = 0 \f$. The + * degree of the Padé approximant and the value of + * squarings are chosen such that the approximation error is no + * more than the round-off error. + * + * The argument of this function should correspond with the (real + * part of) the entries of \c m_M. It is used to select the + * correct implementation using overloading. + */ + void computeUV(double); + + /** \brief Compute Padé approximant to the exponential. + * + * \sa computeUV(double); + */ + void computeUV(float); - /** \internal \brief Compute the (3,3)-Padé approximant to - * the exponential. - * - * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé - * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. - * - * \param M Argument of matrix exponential - * \param Id Identity matrix of same size as M - * \param tmp Temporary storage, to be provided by the caller - * \param M2 Temporary storage, to be provided by the caller - * \param U Even-degree terms in numerator of Padé approximant - * \param V Odd-degree terms in numerator of Padé approximant - */ - template <typename MatrixType> - EIGEN_STRONG_INLINE void pade3(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, - MatrixType& M2, MatrixType& U, MatrixType& V) - { - typedef typename ei_traits<MatrixType>::Scalar Scalar; - const Scalar b[] = {120., 60., 12., 1.}; - M2 = (M * M).lazy(); - tmp = b[3]*M2 + b[1]*Id; - U = (M * tmp).lazy(); - V = b[2]*M2 + b[0]*Id; - } - - /** \internal \brief Compute the (5,5)-Padé approximant to - * the exponential. - * - * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé - * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. - * - * \param M Argument of matrix exponential - * \param Id Identity matrix of same size as M - * \param tmp Temporary storage, to be provided by the caller - * \param M2 Temporary storage, to be provided by the caller - * \param U Even-degree terms in numerator of Padé approximant - * \param V Odd-degree terms in numerator of Padé approximant - */ - template <typename MatrixType> - EIGEN_STRONG_INLINE void pade5(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, - MatrixType& M2, MatrixType& U, MatrixType& V) - { - typedef typename ei_traits<MatrixType>::Scalar Scalar; - const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.}; - M2 = (M * M).lazy(); - MatrixType M4 = (M2 * M2).lazy(); - tmp = b[5]*M4 + b[3]*M2 + b[1]*Id; - U = (M * tmp).lazy(); - V = b[4]*M4 + b[2]*M2 + b[0]*Id; - } - - /** \internal \brief Compute the (7,7)-Padé approximant to - * the exponential. - * - * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé - * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. - * - * \param M Argument of matrix exponential - * \param Id Identity matrix of same size as M - * \param tmp Temporary storage, to be provided by the caller - * \param M2 Temporary storage, to be provided by the caller - * \param U Even-degree terms in numerator of Padé approximant - * \param V Odd-degree terms in numerator of Padé approximant - */ - template <typename MatrixType> - EIGEN_STRONG_INLINE void pade7(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, - MatrixType& M2, MatrixType& U, MatrixType& V) - { - typedef typename ei_traits<MatrixType>::Scalar Scalar; - const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.}; - M2 = (M * M).lazy(); - MatrixType M4 = (M2 * M2).lazy(); - MatrixType M6 = (M4 * M2).lazy(); - tmp = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; - U = (M * tmp).lazy(); - V = b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id; - } - - /** \internal \brief Compute the (9,9)-Padé approximant to - * the exponential. - * - * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé - * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. - * - * \param M Argument of matrix exponential - * \param Id Identity matrix of same size as M - * \param tmp Temporary storage, to be provided by the caller - * \param M2 Temporary storage, to be provided by the caller - * \param U Even-degree terms in numerator of Padé approximant - * \param V Odd-degree terms in numerator of Padé approximant - */ - template <typename MatrixType> - EIGEN_STRONG_INLINE void pade9(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, - MatrixType& M2, MatrixType& U, MatrixType& V) - { typedef typename ei_traits<MatrixType>::Scalar Scalar; - const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240., + typedef typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real RealScalar; + + /** \brief Pointer to matrix whose exponential is to be computed. */ + const MatrixType* m_M; + + /** \brief Even-degree terms in numerator of Padé approximant. */ + MatrixType m_U; + + /** \brief Odd-degree terms in numerator of Padé approximant. */ + MatrixType m_V; + + /** \brief Used for temporary storage. */ + MatrixType m_tmp1; + + /** \brief Used for temporary storage. */ + MatrixType m_tmp2; + + /** \brief Identity matrix of the same size as \c m_M. */ + MatrixType m_Id; + + /** \brief Number of squarings required in the last step. */ + int m_squarings; + + /** \brief L1 norm of m_M. */ + float m_l1norm; +}; + +template <typename MatrixType> +MatrixExponential<MatrixType>::MatrixExponential(const MatrixType &M, MatrixType *result) : + m_M(&M), + m_U(M.rows(),M.cols()), + m_V(M.rows(),M.cols()), + m_tmp1(M.rows(),M.cols()), + m_tmp2(M.rows(),M.cols()), + m_Id(MatrixType::Identity(M.rows(), M.cols())), + m_squarings(0), + m_l1norm(static_cast<float>(M.cwise().abs().colwise().sum().maxCoeff())) +{ + computeUV(RealScalar()); + m_tmp1 = m_U + m_V; // numerator of Pade approximant + m_tmp2 = -m_U + m_V; // denominator of Pade approximant + m_tmp2.partialLu().solve(m_tmp1, result); + for (int i=0; i<m_squarings; i++) + *result *= *result; // undo scaling by repeated squaring +} + +template <typename MatrixType> +EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade3(const MatrixType &A) +{ + const Scalar b[] = {120., 60., 12., 1.}; + m_tmp1.noalias() = A * A; + m_tmp2 = b[3]*m_tmp1 + b[1]*m_Id; + m_U.noalias() = A * m_tmp2; + m_V = b[2]*m_tmp1 + b[0]*m_Id; +} + +template <typename MatrixType> +EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade5(const MatrixType &A) +{ + const Scalar b[] = {30240., 15120., 3360., 420., 30., 1.}; + MatrixType A2 = A * A; + m_tmp1.noalias() = A2 * A2; + m_tmp2 = b[5]*m_tmp1 + b[3]*A2 + b[1]*m_Id; + m_U.noalias() = A * m_tmp2; + m_V = b[4]*m_tmp1 + b[2]*A2 + b[0]*m_Id; +} + +template <typename MatrixType> +EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade7(const MatrixType &A) +{ + const Scalar b[] = {17297280., 8648640., 1995840., 277200., 25200., 1512., 56., 1.}; + MatrixType A2 = A * A; + MatrixType A4 = A2 * A2; + m_tmp1.noalias() = A4 * A2; + m_tmp2 = b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; + m_U.noalias() = A * m_tmp2; + m_V = b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; +} + +template <typename MatrixType> +EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade9(const MatrixType &A) +{ + const Scalar b[] = {17643225600., 8821612800., 2075673600., 302702400., 30270240., 2162160., 110880., 3960., 90., 1.}; - M2 = (M * M).lazy(); - MatrixType M4 = (M2 * M2).lazy(); - MatrixType M6 = (M4 * M2).lazy(); - MatrixType M8 = (M6 * M2).lazy(); - tmp = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; - U = (M * tmp).lazy(); - V = b[8]*M8 + b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id; - } - - /** \internal \brief Compute the (13,13)-Padé approximant to - * the exponential. - * - * After exit, \f$ (V+U)(V-U)^{-1} \f$ is the Padé - * approximant of \f$ \exp(M) \f$ around \f$ M = 0 \f$. - * - * \param M Argument of matrix exponential - * \param Id Identity matrix of same size as M - * \param tmp Temporary storage, to be provided by the caller - * \param M2 Temporary storage, to be provided by the caller - * \param U Even-degree terms in numerator of Padé approximant - * \param V Odd-degree terms in numerator of Padé approximant - */ - template <typename MatrixType> - EIGEN_STRONG_INLINE void pade13(const MatrixType &M, const MatrixType& Id, MatrixType& tmp, - MatrixType& M2, MatrixType& U, MatrixType& V) - { - typedef typename ei_traits<MatrixType>::Scalar Scalar; - const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600., + MatrixType A2 = A * A; + MatrixType A4 = A2 * A2; + MatrixType A6 = A4 * A2; + m_tmp1.noalias() = A6 * A2; + m_tmp2 = b[9]*m_tmp1 + b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; + m_U.noalias() = A * m_tmp2; + m_V = b[8]*m_tmp1 + b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; +} + +template <typename MatrixType> +EIGEN_STRONG_INLINE void MatrixExponential<MatrixType>::pade13(const MatrixType &A) +{ + const Scalar b[] = {64764752532480000., 32382376266240000., 7771770303897600., 1187353796428800., 129060195264000., 10559470521600., 670442572800., 33522128640., 1323241920., 40840800., 960960., 16380., 182., 1.}; - M2 = (M * M).lazy(); - MatrixType M4 = (M2 * M2).lazy(); - MatrixType M6 = (M4 * M2).lazy(); - V = b[13]*M6 + b[11]*M4 + b[9]*M2; - tmp = (M6 * V).lazy(); - tmp += b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; - U = (M * tmp).lazy(); - tmp = b[12]*M6 + b[10]*M4 + b[8]*M2; - V = (M6 * tmp).lazy(); - V += b[6]*M6 + b[4]*M4 + b[2]*M2 + b[0]*Id; - } - - /** \internal \brief Helper class for computing Padé - * approximants to the exponential. - */ - template <typename MatrixType, typename RealScalar = typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real> - struct computeUV_selector - { - /** \internal \brief Compute Padé approximant to the exponential. - * - * Computes \p U, \p V and \p squarings such that \f$ (V+U)(V-U)^{-1} \f$ - * is a Padé of \f$ \exp(2^{-\mbox{squarings}}M) \f$ - * around \f$ M = 0 \f$. The degree of the Padé - * approximant and the value of squarings are chosen such that - * the approximation error is no more than the round-off error. - * - * \param M Argument of matrix exponential - * \param Id Identity matrix of same size as M - * \param tmp1 Temporary storage, to be provided by the caller - * \param tmp2 Temporary storage, to be provided by the caller - * \param U Even-degree terms in numerator of Padé approximant - * \param V Odd-degree terms in numerator of Padé approximant - * \param l1norm L<sub>1</sub> norm of M - * \param squarings Pointer to integer containing number of times - * that the result needs to be squared to find the - * matrix exponential - */ - static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, - MatrixType& U, MatrixType& V, float l1norm, int* squarings); - }; - - template <typename MatrixType> - struct computeUV_selector<MatrixType, float> - { - static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, - MatrixType& U, MatrixType& V, float l1norm, int* squarings) - { - *squarings = 0; - if (l1norm < 4.258730016922831e-001) { - pade3(M, Id, tmp1, tmp2, U, V); - } else if (l1norm < 1.880152677804762e+000) { - pade5(M, Id, tmp1, tmp2, U, V); - } else { - const float maxnorm = 3.925724783138660; - *squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm))); - MatrixType A = M / std::pow(typename ei_traits<MatrixType>::Scalar(2), *squarings); - pade7(A, Id, tmp1, tmp2, U, V); - } - } - }; - - template <typename MatrixType> - struct computeUV_selector<MatrixType, double> - { - static void run(const MatrixType &M, const MatrixType& Id, MatrixType& tmp1, MatrixType& tmp2, - MatrixType& U, MatrixType& V, float l1norm, int* squarings) - { - *squarings = 0; - if (l1norm < 1.495585217958292e-002) { - pade3(M, Id, tmp1, tmp2, U, V); - } else if (l1norm < 2.539398330063230e-001) { - pade5(M, Id, tmp1, tmp2, U, V); - } else if (l1norm < 9.504178996162932e-001) { - pade7(M, Id, tmp1, tmp2, U, V); - } else if (l1norm < 2.097847961257068e+000) { - pade9(M, Id, tmp1, tmp2, U, V); - } else { - const double maxnorm = 5.371920351148152; - *squarings = std::max(0, (int)ceil(log2(l1norm / maxnorm))); - MatrixType A = M / std::pow(typename ei_traits<MatrixType>::Scalar(2), *squarings); - pade13(A, Id, tmp1, tmp2, U, V); - } - } - }; - - /** \internal \brief Compute the matrix exponential. - * - * \param M matrix whose exponential is to be computed. - * \param result pointer to the matrix in which to store the result. - */ - template <typename MatrixType> - void compute(const MatrixType &M, MatrixType* result) - { - MatrixType num, den, U, V; - MatrixType Id = MatrixType::Identity(M.rows(), M.cols()); - float l1norm = M.cwise().abs().colwise().sum().maxCoeff(); - int squarings; - computeUV_selector<MatrixType>::run(M, Id, num, den, U, V, l1norm, &squarings); - num = U + V; // numerator of Pade approximant - den = -U + V; // denominator of Pade approximant - den.partialLu().solve(num, result); - for (int i=0; i<squarings; i++) - *result *= *result; // undo scaling by repeated squaring + MatrixType A2 = A * A; + MatrixType A4 = A2 * A2; + m_tmp1.noalias() = A4 * A2; + m_V = b[13]*m_tmp1 + b[11]*A4 + b[9]*A2; // used for temporary storage + m_tmp2.noalias() = m_tmp1 * m_V; + m_tmp2 += b[7]*m_tmp1 + b[5]*A4 + b[3]*A2 + b[1]*m_Id; + m_U.noalias() = A * m_tmp2; + m_tmp2 = b[12]*m_tmp1 + b[10]*A4 + b[8]*A2; + m_V.noalias() = m_tmp1 * m_tmp2; + m_V += b[6]*m_tmp1 + b[4]*A4 + b[2]*A2 + b[0]*m_Id; +} + +template <typename MatrixType> +void MatrixExponential<MatrixType>::computeUV(float) +{ + if (m_l1norm < 4.258730016922831e-001) { + pade3(*m_M); + } else if (m_l1norm < 1.880152677804762e+000) { + pade5(*m_M); + } else { + const float maxnorm = 3.925724783138660f; + m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm))); + MatrixType A = *m_M / std::pow(Scalar(2), m_squarings); + pade7(A); } +} -} // end of namespace MatrixExponentialInternal +template <typename MatrixType> +void MatrixExponential<MatrixType>::computeUV(double) +{ + if (m_l1norm < 1.495585217958292e-002) { + pade3(*m_M); + } else if (m_l1norm < 2.539398330063230e-001) { + pade5(*m_M); + } else if (m_l1norm < 9.504178996162932e-001) { + pade7(*m_M); + } else if (m_l1norm < 2.097847961257068e+000) { + pade9(*m_M); + } else { + const double maxnorm = 5.371920351148152; + m_squarings = std::max(0, (int)ceil(log2(m_l1norm / maxnorm))); + MatrixType A = *m_M / std::pow(Scalar(2), m_squarings); + pade13(A); + } +} template <typename Derived> EIGEN_STRONG_INLINE void ei_matrix_exponential(const MatrixBase<Derived> &M, typename MatrixBase<Derived>::PlainMatrixType* result) { ei_assert(M.rows() == M.cols()); - MatrixExponentialInternal::compute(M.eval(), result); + MatrixExponential<typename MatrixBase<Derived>::PlainMatrixType>(M, result); } #endif // EIGEN_MATRIX_EXPONENTIAL diff --git a/unsupported/doc/Doxyfile.in b/unsupported/doc/Doxyfile.in index c33986224..561f18fe7 100644 --- a/unsupported/doc/Doxyfile.in +++ b/unsupported/doc/Doxyfile.in @@ -211,7 +211,6 @@ ALIASES = "only_for_vectors=This is only for vectors (either row- "svd_module=This is defined in the %SVD module. \code #include <Eigen/SVD> \endcode" \ "geometry_module=This is defined in the %Geometry module. \code #include <Eigen/Geometry> \endcode" \ "leastsquares_module=This is defined in the %LeastSquares module. \code #include <Eigen/LeastSquares> \endcode" \ - "addexample=\anchor" \ "label=\bug" \ "redstar=<a href='#warningarraymodule' style='color:red;text-decoration: none;'>*</a>" \ "nonstableyet=\warning This is not considered to be part of the stable public API yet. Changes may happen in future releases. See \ref Experimental \"Experimental parts of Eigen\"" diff --git a/unsupported/test/matrixExponential.cpp b/unsupported/test/matrixExponential.cpp index 9a6c335d8..7d65a701a 100644 --- a/unsupported/test/matrixExponential.cpp +++ b/unsupported/test/matrixExponential.cpp @@ -43,10 +43,10 @@ void test2dRotation(double tol) A << 0, 1, -1, 0; for (int i=0; i<=20; i++) { - angle = pow(10, i / 5. - 2); + angle = static_cast<T>(pow(10, i / 5. - 2)); B << cos(angle), sin(angle), -sin(angle), cos(angle); ei_matrix_exponential(angle*A, &C); - VERIFY(C.isApprox(B, tol)); + VERIFY(C.isApprox(B, static_cast<T>(tol))); } } @@ -59,13 +59,13 @@ void test2dHyperbolicRotation(double tol) for (int i=0; i<=20; i++) { - angle = (i-10) / 2.0; + angle = static_cast<T>((i-10) / 2.0); ch = std::cosh(angle); sh = std::sinh(angle); A << 0, angle*imagUnit, -angle*imagUnit, 0; B << ch, sh*imagUnit, -sh*imagUnit, ch; ei_matrix_exponential(A, &C); - VERIFY(C.isApprox(B, tol)); + VERIFY(C.isApprox(B, static_cast<T>(tol))); } } @@ -77,13 +77,13 @@ void testPascal(double tol) Matrix<T,Dynamic,Dynamic> A(size,size), B(size,size), C(size,size); A.setZero(); for (int i=0; i<size-1; i++) - A(i+1,i) = i+1; + A(i+1,i) = static_cast<T>(i+1); B.setZero(); for (int i=0; i<size; i++) for (int j=0; j<=i; j++) - B(i,j) = binom(i,j); + B(i,j) = static_cast<T>(binom(i,j)); ei_matrix_exponential(A, &C); - VERIFY(C.isApprox(B, tol)); + VERIFY(C.isApprox(B, static_cast<T>(tol))); } } @@ -98,11 +98,13 @@ void randomTest(const MatrixType& m, double tol) MatrixType m1(rows, cols), m2(rows, cols), m3(rows, cols), identity = MatrixType::Identity(rows, rows); + typedef typename NumTraits<typename ei_traits<MatrixType>::Scalar>::Real RealScalar; + for(int i = 0; i < g_repeat; i++) { m1 = MatrixType::Random(rows, cols); ei_matrix_exponential(m1, &m2); ei_matrix_exponential(-m1, &m3); - VERIFY(identity.isApprox(m2 * m3, tol)); + VERIFY(identity.isApprox(m2 * m3, static_cast<RealScalar>(tol))); } } |