diff options
45 files changed, 1128 insertions, 579 deletions
diff --git a/Eigen/Core b/Eigen/Core index 854f737d6..42eb363a9 100644 --- a/Eigen/Core +++ b/Eigen/Core @@ -143,6 +143,7 @@ namespace Eigen { #include "src/Core/Functors.h" #include "src/Core/MatrixBase.h" +#include "src/Core/AnyMatrixBase.h" #include "src/Core/Coeffs.h" #ifndef EIGEN_PARSED_BY_DOXYGEN // work around Doxygen bug triggered by Assign.h r814874 diff --git a/Eigen/Householder b/Eigen/Householder index ba06bd8fb..ef3e61373 100644 --- a/Eigen/Householder +++ b/Eigen/Householder @@ -16,6 +16,7 @@ namespace Eigen { */ #include "src/Householder/Householder.h" +#include "src/Householder/HouseholderSequence.h" } // namespace Eigen diff --git a/Eigen/src/Core/AnyMatrixBase.h b/Eigen/src/Core/AnyMatrixBase.h new file mode 100644 index 000000000..cd354d7b1 --- /dev/null +++ b/Eigen/src/Core/AnyMatrixBase.h @@ -0,0 +1,153 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// 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_ANYMATRIXBASE_H +#define EIGEN_ANYMATRIXBASE_H + + +/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T). + * + * In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase. + * + * Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc. + * + * Notice that this class is trivial, it is only used to disambiguate overloaded functions. + */ +template<typename Derived> struct AnyMatrixBase +{ + typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType; + + Derived& derived() { return *static_cast<Derived*>(this); } + const Derived& derived() const { return *static_cast<const Derived*>(this); } + + /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ + inline int rows() const { return derived().rows(); } + /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ + inline int cols() const { return derived().cols(); } + + /** \internal Don't use it, but do the equivalent: \code dst = *this; \endcode */ + template<typename Dest> inline void evalTo(Dest& dst) const + { derived().evalTo(dst); } + + /** \internal Don't use it, but do the equivalent: \code dst += *this; \endcode */ + template<typename Dest> inline void addToDense(Dest& dst) const + { + // This is the default implementation, + // derived class can reimplement it in a more optimized way. + typename Dest::PlainMatrixType res(rows(),cols()); + evalTo(res); + dst += res; + } + + /** \internal Don't use it, but do the equivalent: \code dst -= *this; \endcode */ + template<typename Dest> inline void subToDense(Dest& dst) const + { + // This is the default implementation, + // derived class can reimplement it in a more optimized way. + typename Dest::PlainMatrixType res(rows(),cols()); + evalTo(res); + dst -= res; + } + + /** \internal Don't use it, but do the equivalent: \code dst.applyOnTheRight(*this); \endcode */ + template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const + { + // This is the default implementation, + // derived class can reimplement it in a more optimized way. + dst = dst * this->derived(); + } + + /** \internal Don't use it, but do the equivalent: \code dst.applyOnTheLeft(*this); \endcode */ + template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const + { + // This is the default implementation, + // derived class can reimplement it in a more optimized way. + dst = this->derived() * dst; + } + +}; + +/*************************************************************************** +* Implementation of matrix base methods +***************************************************************************/ + +/** Copies the generic expression \a other into *this. \returns a reference to *this. + * The expression must provide a (templated) evalToDense(Derived& dst) const function + * which does the actual job. In practice, this allows any user to write its own + * special matrix without having to modify MatrixBase */ +template<typename Derived> +template<typename OtherDerived> +Derived& MatrixBase<Derived>::operator=(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().evalTo(derived()); + return derived(); +} + +template<typename Derived> +template<typename OtherDerived> +Derived& MatrixBase<Derived>::operator+=(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().addToDense(derived()); + return derived(); +} + +template<typename Derived> +template<typename OtherDerived> +Derived& MatrixBase<Derived>::operator-=(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().subToDense(derived()); + return derived(); +} + +/** replaces \c *this by \c *this * \a other. + * + * \returns a reference to \c *this + */ +template<typename Derived> +template<typename OtherDerived> +inline Derived& +MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheRight(derived()); + return derived(); +} + +/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=() */ +template<typename Derived> +template<typename OtherDerived> +inline void MatrixBase<Derived>::applyOnTheRight(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheRight(derived()); +} + +/** replaces \c *this by \c *this * \a other. */ +template<typename Derived> +template<typename OtherDerived> +inline void MatrixBase<Derived>::applyOnTheLeft(const AnyMatrixBase<OtherDerived> &other) +{ + other.derived().applyThisOnTheLeft(derived()); +} + +#endif // EIGEN_ANYMATRIXBASE_H diff --git a/Eigen/src/Core/Matrix.h b/Eigen/src/Core/Matrix.h index 0975b3b77..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 * @@ -308,7 +309,7 @@ class Matrix */ template<typename OtherDerived> EIGEN_STRONG_INLINE void resizeLike(const MatrixBase<OtherDerived>& other) - { + { if(RowsAtCompileTime == 1) { ei_assert(other.isVector()); @@ -324,40 +325,28 @@ class Matrix /** 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, although it is legal to call it on any - * matrix as long as fixed dimensions are left unchanged. If you only want to change the number + * 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 per - * default and set to zero when init_with_zero is set to true. + * of *this. In case values need to be appended to the matrix they will be uninitialized. */ - inline void conservativeResize(int rows, int cols, bool init_with_zero = false) + EIGEN_STRONG_INLINE void conservativeResize(int rows, int cols) { - // 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(Matrix) - PlainMatrixType tmp = init_with_zero ? PlainMatrixType::Zero(rows, cols) : PlainMatrixType(rows,cols); - const int common_rows = std::min(rows, this->rows()); - const int common_cols = std::min(cols, this->cols()); - tmp.block(0,0,common_rows,common_cols) = this->block(0,0,common_rows,common_cols); - this->derived().swap(tmp); + conservativeResizeLike(PlainMatrixType(rows, cols)); } - EIGEN_STRONG_INLINE void conservativeResize(int rows, NoChange_t, bool init_with_zero = false) + EIGEN_STRONG_INLINE void conservativeResize(int rows, NoChange_t) { - // Note: see the comment in conservativeResize(int,int,bool) - conservativeResize(rows, cols(), init_with_zero); + // Note: see the comment in conservativeResize(int,int) + conservativeResize(rows, cols()); } - EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, int cols, bool init_with_zero = false) + EIGEN_STRONG_INLINE void conservativeResize(NoChange_t, int cols) { - // Note: see the comment in conservativeResize(int,int,bool) - conservativeResize(rows(), cols, init_with_zero); + // 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. @@ -366,21 +355,17 @@ class Matrix * 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 per default and set - * to zero when init_with_zero is set to true. + * When values are appended, they will be uninitialized. */ - inline void conservativeResize(int size, bool init_with_zero = false) + EIGEN_STRONG_INLINE void conservativeResize(int size) { - EIGEN_STATIC_ASSERT_VECTOR_ONLY(Matrix) - EIGEN_STATIC_ASSERT_DYNAMIC_SIZE(Matrix) + conservativeResizeLike(PlainMatrixType(size)); + } - if (RowsAtCompileTime == 1 || ColsAtCompileTime == 1) - { - PlainMatrixType tmp = init_with_zero ? PlainMatrixType::Zero(size) : PlainMatrixType(size); - const int common_size = std::min<int>(this->size(),size); - tmp.segment(0,common_size) = this->segment(0,common_size); - this->derived().swap(tmp); - } + 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. @@ -713,13 +698,45 @@ class Matrix m_storage.data()[1] = y; } - template<typename MatrixType, typename OtherDerived, bool IsSameType, bool IsDynamicSize> + template<typename MatrixType, typename OtherDerived, bool SwapPointers> friend struct ei_matrix_swap_impl; }; -template<typename MatrixType, typename OtherDerived, - bool IsSameType = ei_is_same_type<MatrixType, OtherDerived>::ret, - bool IsDynamicSize = MatrixType::SizeAtCompileTime==Dynamic> +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) @@ -729,7 +746,7 @@ struct ei_matrix_swap_impl }; template<typename MatrixType, typename OtherDerived> -struct ei_matrix_swap_impl<MatrixType, OtherDerived, true, true> +struct ei_matrix_swap_impl<MatrixType, OtherDerived, true> { static inline void run(MatrixType& matrix, MatrixBase<OtherDerived>& other) { @@ -741,7 +758,8 @@ template<typename _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows, int template<typename OtherDerived> inline void Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>::swap(const MatrixBase<OtherDerived>& other) { - ei_matrix_swap_impl<Matrix, OtherDerived>::run(*this, *const_cast<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 diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index ad5fde562..4835f167c 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -26,46 +26,6 @@ #ifndef EIGEN_MATRIXBASE_H #define EIGEN_MATRIXBASE_H - -/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T). - * - * In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase. - * - * Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc. - * - * Notice that this class is trivial, it is only used to disambiguate overloaded functions. - */ -template<typename Derived> struct AnyMatrixBase - : public ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar, - typename NumTraits<typename ei_traits<Derived>::Scalar>::Real> -{ - typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType; - - Derived& derived() { return *static_cast<Derived*>(this); } - const Derived& derived() const { return *static_cast<const Derived*>(this); } - /** \returns the number of rows. \sa cols(), RowsAtCompileTime */ - inline int rows() const { return derived().rows(); } - /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/ - inline int cols() const { return derived().cols(); } - - template<typename Dest> inline void evalTo(Dest& dst) const - { derived().evalTo(dst); } - - template<typename Dest> inline void addToDense(Dest& dst) const - { - typename Dest::PlainMatrixType res(rows(),cols()); - evalToDense(res); - dst += res; - } - - template<typename Dest> inline void subToDense(Dest& dst) const - { - typename Dest::PlainMatrixType res(rows(),cols()); - evalToDense(res); - dst -= res; - } -}; - /** \class MatrixBase * * \brief Base class for all matrices, vectors, and expressions @@ -93,11 +53,11 @@ 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> #endif // not EIGEN_PARSED_BY_DOXYGEN { public: - #ifndef EIGEN_PARSED_BY_DOXYGEN using ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar, typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>::operator*; @@ -302,21 +262,14 @@ template<typename Derived> class MatrixBase */ Derived& operator=(const MatrixBase& other); - /** Copies the generic expression \a other into *this. \returns a reference to *this. - * The expression must provide a (templated) evalToDense(Derived& dst) const function - * which does the actual job. In practice, this allows any user to write its own - * special matrix without having to modify MatrixBase */ template<typename OtherDerived> - Derived& operator=(const AnyMatrixBase<OtherDerived> &other) - { other.derived().evalToDense(derived()); return derived(); } + Derived& operator=(const AnyMatrixBase<OtherDerived> &other); template<typename OtherDerived> - Derived& operator+=(const AnyMatrixBase<OtherDerived> &other) - { other.derived().addToDense(derived()); return derived(); } + Derived& operator+=(const AnyMatrixBase<OtherDerived> &other); template<typename OtherDerived> - Derived& operator-=(const AnyMatrixBase<OtherDerived> &other) - { other.derived().subToDense(derived()); return derived(); } + Derived& operator-=(const AnyMatrixBase<OtherDerived> &other); template<typename OtherDerived,typename OtherEvalType> Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func); @@ -437,6 +390,12 @@ template<typename Derived> class MatrixBase template<typename OtherDerived> Derived& operator*=(const AnyMatrixBase<OtherDerived>& other); + template<typename OtherDerived> + void applyOnTheLeft(const AnyMatrixBase<OtherDerived>& other); + + template<typename OtherDerived> + void applyOnTheRight(const AnyMatrixBase<OtherDerived>& other); + template<typename DiagonalDerived> const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight> operator*(const DiagonalBase<DiagonalDerived> &diagonal) const; @@ -676,8 +635,11 @@ template<typename Derived> class MatrixBase typename ei_traits<Derived>::Scalar minCoeff() const; typename ei_traits<Derived>::Scalar maxCoeff() const; - typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col = 0) const; - typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col = 0) const; + typename ei_traits<Derived>::Scalar minCoeff(int* row, int* col) const; + typename ei_traits<Derived>::Scalar maxCoeff(int* row, int* col) const; + + typename ei_traits<Derived>::Scalar minCoeff(int* index) const; + typename ei_traits<Derived>::Scalar maxCoeff(int* index) const; template<typename BinaryOp> typename ei_result_of<BinaryOp(typename ei_traits<Derived>::Scalar)>::type diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h index e7227d4f6..7f0c2df6e 100644 --- a/Eigen/src/Core/Product.h +++ b/Eigen/src/Core/Product.h @@ -434,18 +434,4 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived()); } - - -/** replaces \c *this by \c *this * \a other. - * - * \returns a reference to \c *this - */ -template<typename Derived> -template<typename OtherDerived> -inline Derived & -MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other) -{ - return derived() = derived() * other.derived(); -} - #endif // EIGEN_PRODUCT_H diff --git a/Eigen/src/Core/StableNorm.h b/Eigen/src/Core/StableNorm.h index 77fe79782..facab9dbd 100644 --- a/Eigen/src/Core/StableNorm.h +++ b/Eigen/src/Core/StableNorm.h @@ -56,7 +56,7 @@ MatrixBase<Derived>::stableNorm() const { const int blockSize = 4096; RealScalar scale = 0; - RealScalar invScale; + RealScalar invScale = 1; RealScalar ssq = 0; // sum of square enum { Alignment = (int(Flags)&DirectAccessBit) || (int(Flags)&AlignedBit) ? ForceAligned : AsRequested diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h index b0362f20c..17726bca3 100644 --- a/Eigen/src/Core/TriangularMatrix.h +++ b/Eigen/src/Core/TriangularMatrix.h @@ -91,9 +91,9 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived> #endif // not EIGEN_PARSED_BY_DOXYGEN template<typename DenseDerived> - void evalToDense(MatrixBase<DenseDerived> &other) const; + void evalTo(MatrixBase<DenseDerived> &other) const; template<typename DenseDerived> - void evalToDenseLazy(MatrixBase<DenseDerived> &other) const; + void evalToLazy(MatrixBase<DenseDerived> &other) const; protected: @@ -546,23 +546,23 @@ void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDeri * If the matrix is triangular, the opposite part is set to zero. */ template<typename Derived> template<typename DenseDerived> -void TriangularBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const +void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const { if(ei_traits<Derived>::Flags & EvalBeforeAssigningBit) { typename Derived::PlainMatrixType other_evaluated(rows(), cols()); - evalToDenseLazy(other_evaluated); + evalToLazy(other_evaluated); other.derived().swap(other_evaluated); } else - evalToDenseLazy(other.derived()); + evalToLazy(other.derived()); } /** Assigns a triangular or selfadjoint matrix to a dense matrix. * If the matrix is triangular, the opposite part is set to zero. */ template<typename Derived> template<typename DenseDerived> -void TriangularBase<Derived>::evalToDenseLazy(MatrixBase<DenseDerived> &other) const +void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const { const bool unroll = DenseDerived::SizeAtCompileTime * Derived::CoeffReadCost / 2 <= EIGEN_UNROLLING_LIMIT; diff --git a/Eigen/src/Core/VectorBlock.h b/Eigen/src/Core/VectorBlock.h index b291f7b1a..65268b626 100644 --- a/Eigen/src/Core/VectorBlock.h +++ b/Eigen/src/Core/VectorBlock.h @@ -77,11 +77,12 @@ template<typename VectorType, int Size, int PacketAccess> class VectorBlock typedef Block<VectorType, ei_traits<VectorType>::RowsAtCompileTime==1 ? 1 : Size, ei_traits<VectorType>::ColsAtCompileTime==1 ? 1 : Size, - PacketAccess> Base; + PacketAccess> _Base; enum { IsColVector = ei_traits<VectorType>::ColsAtCompileTime==1 }; public: + _EIGEN_GENERIC_PUBLIC_INTERFACE(VectorBlock, _Base) using Base::operator=; using Base::operator+=; diff --git a/Eigen/src/Core/Visitor.h b/Eigen/src/Core/Visitor.h index 598c2db8d..590efc766 100644 --- a/Eigen/src/Core/Visitor.h +++ b/Eigen/src/Core/Visitor.h @@ -164,7 +164,7 @@ struct ei_functor_traits<ei_max_coeff_visitor<Scalar> > { /** \returns the minimum of all coefficients of *this * and puts in *row and *col its location. * - * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff() + * \sa MatrixBase::minCoeff(int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff() */ template<typename Derived> typename ei_traits<Derived>::Scalar @@ -177,6 +177,22 @@ MatrixBase<Derived>::minCoeff(int* row, int* col) const return minVisitor.res; } +/** \returns the minimum of all coefficients of *this + * and puts in *index its location. + * + * \sa MatrixBase::minCoeff(int*,int*), MatrixBase::maxCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::minCoeff() + */ +template<typename Derived> +typename ei_traits<Derived>::Scalar +MatrixBase<Derived>::minCoeff(int* index) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + ei_min_coeff_visitor<Scalar> minVisitor; + this->visit(minVisitor); + *index = (RowsAtCompileTime==1) ? minVisitor.col : minVisitor.row; + return minVisitor.res; +} + /** \returns the maximum of all coefficients of *this * and puts in *row and *col its location. * @@ -193,5 +209,20 @@ MatrixBase<Derived>::maxCoeff(int* row, int* col) const return maxVisitor.res; } +/** \returns the maximum of all coefficients of *this + * and puts in *index its location. + * + * \sa MatrixBase::maxCoeff(int*,int*), MatrixBase::minCoeff(int*,int*), MatrixBase::visitor(), MatrixBase::maxCoeff() + */ +template<typename Derived> +typename ei_traits<Derived>::Scalar +MatrixBase<Derived>::maxCoeff(int* index) const +{ + EIGEN_STATIC_ASSERT_VECTOR_ONLY(Derived) + ei_max_coeff_visitor<Scalar> maxVisitor; + this->visit(maxVisitor); + *index = (RowsAtCompileTime==1) ? maxVisitor.col : maxVisitor.row; + return maxVisitor.res; +} #endif // EIGEN_VISITOR_H diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index c5f27d80b..3f66738f0 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -123,6 +123,7 @@ 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 VectorsType, typename CoeffsType> class HouseholderSequence; template<typename Scalar> class PlanarRotation; // Geometry module: diff --git a/Eigen/src/Core/util/XprHelper.h b/Eigen/src/Core/util/XprHelper.h index 2f8d35d05..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 diff --git a/Eigen/src/Eigenvalues/ComplexSchur.h b/Eigen/src/Eigenvalues/ComplexSchur.h index 58e2ea440..0534715c4 100644 --- a/Eigen/src/Eigenvalues/ComplexSchur.h +++ b/Eigen/src/Eigenvalues/ComplexSchur.h @@ -31,8 +31,15 @@ * * \class ComplexShur * - * \brief Performs a complex Shur decomposition of a real or complex square matrix + * \brief Performs a complex Schur decomposition of a real or complex square matrix * + * Given a real or complex square matrix A, this class computes the Schur decomposition: + * \f$ A = U T U^*\f$ where U is a unitary complex matrix, and T is a complex upper + * triangular matrix. + * + * The diagonal of the matrix T corresponds to the eigenvalues of the matrix A. + * + * \sa class RealSchur, class EigenSolver */ template<typename _MatrixType> class ComplexSchur { @@ -42,41 +49,56 @@ template<typename _MatrixType> class ComplexSchur typedef typename NumTraits<Scalar>::Real RealScalar; typedef std::complex<RealScalar> Complex; typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType; + enum { + Size = MatrixType::RowsAtCompileTime + }; - /** - * \brief Default Constructor. + /** \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to - * perform decompositions via ComplexSchur::compute(const MatrixType&). + * perform decompositions via ComplexSchur::compute(). */ - ComplexSchur() : m_matT(), m_matU(), m_isInitialized(false) + ComplexSchur(int size = Size==Dynamic ? 0 : Size) + : m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false) {} - ComplexSchur(const MatrixType& matrix) + /** Constructor computing the Schur decomposition of the matrix \a matrix. + * If \a skipU is true, then the matrix U is not computed. */ + ComplexSchur(const MatrixType& matrix, bool skipU = false) : m_matT(matrix.rows(),matrix.cols()), m_matU(matrix.rows(),matrix.cols()), - m_isInitialized(false) + m_isInitialized(false), + m_matUisUptodate(false) { - compute(matrix); + compute(matrix, skipU); } - ComplexMatrixType matrixU() const + /** \returns a const reference to the matrix U of the respective Schur decomposition. */ + const ComplexMatrixType& matrixU() const { ei_assert(m_isInitialized && "ComplexSchur is not initialized."); + ei_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition."); return m_matU; } - ComplexMatrixType matrixT() const + /** \returns a const reference to the matrix T of the respective Schur decomposition. + * Note that this function returns a plain square matrix. If you want to reference + * only the upper triangular part, use: + * \code schur.matrixT().triangularView<Upper>() \endcode. */ + const ComplexMatrixType& matrixT() const { ei_assert(m_isInitialized && "ComplexShur is not initialized."); return m_matT; } - void compute(const MatrixType& matrix); + /** Computes the Schur decomposition of the matrix \a matrix. + * If \a skipU is true, then the matrix U is not computed. */ + void compute(const MatrixType& matrix, bool skipU = false); protected: ComplexMatrixType m_matT, m_matU; bool m_isInitialized; + bool m_matUisUptodate; }; /** Computes the principal value of the square root of the complex \a z. */ @@ -117,17 +139,20 @@ std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z) } template<typename MatrixType> -void ComplexSchur<MatrixType>::compute(const MatrixType& matrix) +void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU) { // this code is inspired from Jampack + + m_matUisUptodate = false; assert(matrix.cols() == matrix.rows()); int n = matrix.cols(); // Reduce to Hessenberg form + // TODO skip Q if skipU = true HessenbergDecomposition<MatrixType> hess(matrix); m_matT = hess.matrixH(); - m_matU = hess.matrixQ(); + if(!skipU) m_matU = hess.matrixQ(); int iu = m_matT.cols() - 1; int il; @@ -206,7 +231,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix) { 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(!skipU) m_matU.applyOnTheRight(i, i+1, rot); if(i != iu-1) { @@ -232,6 +257,7 @@ void ComplexSchur<MatrixType>::compute(const MatrixType& matrix) */ m_isInitialized = true; + m_matUisUptodate = !skipU; } #endif // EIGEN_COMPLEX_SCHUR_H diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index b1e21d4ee..bb7e3fcfc 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -88,14 +88,14 @@ template<typename _MatrixType> class HessenbergDecomposition _compute(m_matrix, m_hCoeffs); } - /** \returns the householder coefficients allowing to + /** \returns a const reference to the householder coefficients allowing to * reconstruct the matrix Q from the packed data. * * \sa packedMatrix() */ - CoeffVectorType householderCoefficients() const { return m_hCoeffs; } + const CoeffVectorType& householderCoefficients() const { return m_hCoeffs; } - /** \returns the internal result of the decomposition. + /** \returns a const reference to the internal representation of the decomposition. * * The returned matrix contains the following information: * - the upper part and lower sub-diagonal represent the Hessenberg matrix H diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 9eb8ed535..dcb41435b 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -395,7 +395,7 @@ public: Transform& fromPositionOrientationScale(const MatrixBase<PositionDerived> &position, const OrientationType& orientation, const MatrixBase<ScaleDerived> &scale); - inline const MatrixType inverse(TransformTraits traits = (TransformTraits)Mode) const; + inline Transform inverse(TransformTraits traits = (TransformTraits)Mode) const; /** \returns a const pointer to the column major internal matrix */ const Scalar* data() const { return m_matrix.data(); } @@ -874,7 +874,7 @@ Transform<Scalar,Dim,Mode>::fromPositionOrientationScale(const MatrixBase<Positi /** \nonstableyet * - * \returns the inverse transformation matrix according to some given knowledge + * \returns the inverse transformation according to some given knowledge * on \c *this. * * \param traits allows to optimize the inversion process when the transformion @@ -892,37 +892,37 @@ Transform<Scalar,Dim,Mode>::fromPositionOrientationScale(const MatrixBase<Positi * \sa MatrixBase::inverse() */ template<typename Scalar, int Dim, int Mode> -const typename Transform<Scalar,Dim,Mode>::MatrixType +Transform<Scalar,Dim,Mode> Transform<Scalar,Dim,Mode>::inverse(TransformTraits hint) const { + Transform res; if (hint == Projective) { - return m_matrix.inverse(); + res.matrix() = m_matrix.inverse(); } else { - MatrixType res; if (hint == Isometry) { - res.template corner<Dim,Dim>(TopLeft) = linear().transpose(); + res.matrix().template corner<Dim,Dim>(TopLeft) = linear().transpose(); } else if(hint&Affine) { - res.template corner<Dim,Dim>(TopLeft) = linear().inverse(); + res.matrix().template corner<Dim,Dim>(TopLeft) = linear().inverse(); } else { ei_assert(false && "Invalid transform traits in Transform::Inverse"); } // translation and remaining parts - res.template corner<Dim,1>(TopRight) = - res.template corner<Dim,Dim>(TopLeft) * translation(); + res.matrix().template corner<Dim,1>(TopRight) = - res.matrix().template corner<Dim,Dim>(TopLeft) * translation(); if(int(Mode)!=int(AffineCompact)) { - res.template block<1,Dim>(Dim,0).setZero(); - res.coeffRef(Dim,Dim) = 1; + res.matrix().template block<1,Dim>(Dim,0).setZero(); + res.matrix().coeffRef(Dim,Dim) = 1; } - return res; } + return res; } /***************************************************** diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h new file mode 100644 index 000000000..16e362814 --- /dev/null +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -0,0 +1,168 @@ +// 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/>. + +#ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H +#define EIGEN_HOUSEHOLDER_SEQUENCE_H + +/** \ingroup Householder_Module + * \householder_module + * \class HouseholderSequence + * \brief Represents a sequence of householder reflections with decreasing size + * + * This class represents a product sequence of householder reflections \f$ H = \Pi_0^{n-1} H_i \f$ + * where \f$ H_i \f$ is the i-th householder transformation \f$ I - h_i v_i v_i^* \f$, + * \f$ v_i \f$ is the i-th householder vector \f$ [ 1, m_vectors(i+1,i), m_vectors(i+2,i), ...] \f$ + * and \f$ h_i \f$ is the i-th householder coefficient \c m_coeffs[i]. + * + * Typical usages are listed below, where H is a HouseholderSequence: + * \code + * A.applyOnTheRight(H); // A = A * H + * A.applyOnTheLeft(H); // A = H * A + * A.applyOnTheRight(H.adjoint()); // A = A * H^* + * A.applyOnTheLeft(H.adjoint()); // A = H^* * A + * MatrixXd Q = H; // conversion to a dense matrix + * \endcode + * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate. + * + * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight() + */ + +template<typename VectorsType, typename CoeffsType> +struct ei_traits<HouseholderSequence<VectorsType,CoeffsType> > +{ + typedef typename VectorsType::Scalar Scalar; + enum { + RowsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime, + ColsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime, + MaxRowsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime, + MaxColsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime, + Flags = 0 + }; +}; + +template<typename VectorsType, typename CoeffsType> class HouseholderSequence + : public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType> > +{ + typedef typename VectorsType::Scalar Scalar; + public: + + typedef HouseholderSequence<VectorsType, + typename ei_meta_if<NumTraits<Scalar>::IsComplex, + NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >, + CoeffsType>::ret> ConjugateReturnType; + + HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false) + : m_vectors(v), m_coeffs(h), m_trans(trans) + {} + + int rows() const { return m_vectors.rows(); } + int cols() const { return m_vectors.rows(); } + + HouseholderSequence transpose() const + { return HouseholderSequence(m_vectors, m_coeffs, !m_trans); } + + ConjugateReturnType conjugate() const + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); } + + ConjugateReturnType adjoint() const + { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); } + + ConjugateReturnType inverse() const { return adjoint(); } + + /** \internal */ + template<typename DestType> void evalTo(DestType& dst) const + { + int vecs = std::min(m_vectors.cols(),m_vectors.rows()); + int length = m_vectors.rows(); + dst.setIdentity(); + Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows()); + for(int k = vecs-1; k >= 0; --k) + { + if(m_trans) + dst.corner(BottomRight, length-k, length-k) + .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0)); + else + dst.corner(BottomRight, length-k, length-k) + .applyHouseholderOnTheLeft(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(k)); + } + } + + /** \internal */ + template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const + { + int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors + int length = m_vectors.rows(); // size of the largest householder vector + Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.rows()); + for(int k = 0; k < vecs; ++k) + { + int actual_k = m_trans ? vecs-k-1 : k; + dst.corner(BottomRight, dst.rows(), length-k) + .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0)); + } + } + + /** \internal */ + template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const + { + int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors + int length = m_vectors.rows(); // size of the largest householder vector + Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols()); + for(int k = 0; k < vecs; ++k) + { + int actual_k = m_trans ? k : vecs-k-1; + dst.corner(BottomRight, length-actual_k, dst.cols()) + .applyHouseholderOnTheLeft(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0)); + } + } + + template<typename OtherDerived> + typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other) const + { + typename OtherDerived::PlainMatrixType res(other); + applyThisOnTheLeft(res); + return res; + } + + template<typename OtherDerived> friend + typename OtherDerived::PlainMatrixType operator*(const MatrixBase<OtherDerived>& other, const HouseholderSequence& h) + { + typename OtherDerived::PlainMatrixType res(other); + h.applyThisOnTheRight(res); + return res; + } + + protected: + + typename VectorsType::Nested m_vectors; + typename CoeffsType::Nested m_coeffs; + bool m_trans; +}; + +template<typename VectorsType, typename CoeffsType> +HouseholderSequence<VectorsType,CoeffsType> makeHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false) +{ + return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans); +} + +#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H diff --git a/Eigen/src/Jacobi/Jacobi.h b/Eigen/src/Jacobi/Jacobi.h index 3905f4d8f..eeb81c178 100644 --- a/Eigen/src/Jacobi/Jacobi.h +++ b/Eigen/src/Jacobi/Jacobi.h @@ -123,7 +123,7 @@ bool PlanarRotation<Scalar>::makeJacobi(RealScalar x, Scalar y, RealScalar z) } /** 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} \\ \overline \text{this}_{pq} & \text{this}_{qq} \end{array} \right )\f$ yields + * \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 diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h index 0ef59bac7..e467c62f0 100644 --- a/Eigen/src/LU/PartialLU.h +++ b/Eigen/src/LU/PartialLU.h @@ -2,6 +2,7 @@ // for linear algebra. // // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com> +// 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 @@ -215,10 +216,10 @@ struct ei_partial_lu_impl typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU; typedef Block<MapLU, Dynamic, Dynamic> MatrixType; typedef Block<MatrixType,Dynamic,Dynamic> BlockType; - + /** \internal performs the LU decomposition in-place of the matrix \a lu * using an unblocked algorithm. - * + * * In addition, this function returns the row transpositions in the * vector \a row_transpositions which must have a size equal to the number * of columns of the matrix \a lu, and an integer \a nb_transpositions @@ -232,7 +233,7 @@ struct ei_partial_lu_impl for(int k = 0; k < size; ++k) { int row_of_biggest_in_col; - lu.block(k,k,rows-k,1).cwise().abs().maxCoeff(&row_of_biggest_in_col); + lu.col(k).end(rows-k).cwise().abs().maxCoeff(&row_of_biggest_in_col); row_of_biggest_in_col += k; row_transpositions[k] = row_of_biggest_in_col; @@ -295,7 +296,7 @@ struct ei_partial_lu_impl int bs = std::min(size-k,blockSize); // actual size of the block int trows = rows - k - bs; // trailing rows int tsize = size - k - bs; // trailing size - + // partition the matrix: // A00 | A01 | A02 // lu = A10 | A11 | A12 @@ -343,7 +344,7 @@ void ei_partial_lu_inplace(MatrixType& lu, IntVector& row_transpositions, int& n { ei_assert(lu.cols() == row_transpositions.size()); ei_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); - + ei_partial_lu_impl <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor> ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.stride(), &row_transpositions.coeffRef(0), nb_transpositions); diff --git a/Eigen/src/QR/ColPivotingHouseholderQR.h b/Eigen/src/QR/ColPivotingHouseholderQR.h index 8024e3b9d..c4c7d2d55 100644 --- a/Eigen/src/QR/ColPivotingHouseholderQR.h +++ b/Eigen/src/QR/ColPivotingHouseholderQR.h @@ -45,14 +45,14 @@ template<typename MatrixType> class ColPivotingHouseholderQR { public: - + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) }; - + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; @@ -62,6 +62,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; typedef Matrix<RealScalar, 1, ColsAtCompileTime> RealRowVectorType; + typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; /** * \brief Default Constructor. @@ -99,7 +100,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR template<typename OtherDerived, typename ResultType> bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; - MatrixQType matrixQ(void) const; + HouseholderSequenceType matrixQ(void) const; /** \returns a reference to the matrix where the Householder QR decomposition is stored */ @@ -110,13 +111,13 @@ template<typename MatrixType> class ColPivotingHouseholderQR } ColPivotingHouseholderQR& compute(const MatrixType& matrix); - + const IntRowVectorType& colsPermutation() const { ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); return m_cols_permutation; } - + /** \returns the absolute value of the determinant of the matrix of which * *this is the QR decomposition. It has only linear complexity * (that is, O(n) where n is the dimension of the square matrix) @@ -145,7 +146,7 @@ template<typename MatrixType> class ColPivotingHouseholderQR * \sa absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; - + /** \returns the rank of the matrix of which *this is the QR decomposition. * * \note This is computed at the time of the construction of the QR decomposition. This @@ -268,7 +269,7 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp int cols = matrix.cols(); int size = std::min(rows,cols); m_rank = size; - + m_qr = matrix; m_hCoeffs.resize(size); @@ -279,18 +280,18 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp IntRowVectorType cols_transpositions(matrix.cols()); m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; - + RealRowVectorType colSqNorms(cols); for(int k = 0; k < cols; ++k) colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm(); RealScalar biggestColSqNorm = colSqNorms.maxCoeff(); - + for (int k = 0; k < size; ++k) { int biggest_col_in_corner; RealScalar biggestColSqNormInCorner = colSqNorms.end(cols-k).maxCoeff(&biggest_col_in_corner); biggest_col_in_corner += k; - + // if the corner is negligible, then we have less than full rank, and we can finish early if(ei_isMuchSmallerThan(biggestColSqNormInCorner, biggestColSqNorm, m_precision)) { @@ -302,10 +303,11 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp } break; } - + cols_transpositions.coeffRef(k) = biggest_col_in_corner; if(k != biggest_col_in_corner) { m_qr.col(k).swap(m_qr.col(biggest_col_in_corner)); + std::swap(colSqNorms.coeffRef(k), colSqNorms.coeffRef(biggest_col_in_corner)); ++number_of_transpositions; } @@ -315,7 +317,7 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp m_qr.corner(BottomRight, rows-k, cols-k-1) .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), m_hCoeffs.coeffRef(k), &temp.coeffRef(k+1)); - + colSqNorms.end(cols-k-1) -= m_qr.row(k).end(cols-k-1).cwise().abs2(); } @@ -325,7 +327,7 @@ ColPivotingHouseholderQR<MatrixType>& ColPivotingHouseholderQR<MatrixType>::comp m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - + return *this; } @@ -351,16 +353,11 @@ bool ColPivotingHouseholderQR<MatrixType>::solve( const int rows = m_qr.rows(); const int cols = b.cols(); ei_assert(b.rows() == rows); - + typename OtherDerived::PlainMatrixType c(b); - - Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols); - for (int k = 0; k < m_rank; ++k) - { - int remainingSize = rows-k; - c.corner(BottomRight, remainingSize, cols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); - } + + // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T + c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,m_rank), m_hCoeffs.start(m_rank)).transpose()); if(!isSurjective()) { @@ -380,25 +377,12 @@ bool ColPivotingHouseholderQR<MatrixType>::solve( return true; } -/** \returns the matrix Q */ +/** \returns the matrix Q as a sequence of householder transformations */ template<typename MatrixType> -typename ColPivotingHouseholderQR<MatrixType>::MatrixQType ColPivotingHouseholderQR<MatrixType>::matrixQ() const +typename ColPivotingHouseholderQR<MatrixType>::HouseholderSequenceType ColPivotingHouseholderQR<MatrixType>::matrixQ() const { ei_assert(m_isInitialized && "ColPivotingHouseholderQR is not initialized."); - // compute the product H'_0 H'_1 ... H'_n-1, - // where H_k is the k-th Householder transformation I - h_k v_k v_k' - // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] - int rows = m_qr.rows(); - int cols = m_qr.cols(); - int size = std::min(rows,cols); - MatrixQType res = MatrixQType::Identity(rows, rows); - Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows); - for (int k = size-1; k >= 0; k--) - { - res.block(k, k, rows-k, rows-k) - .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); - } - return res; + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); } #endif // EIGEN_HIDE_HEAVY_CODE diff --git a/Eigen/src/QR/FullPivotingHouseholderQR.h b/Eigen/src/QR/FullPivotingHouseholderQR.h index 0d542cf7a..9fee77803 100644 --- a/Eigen/src/QR/FullPivotingHouseholderQR.h +++ b/Eigen/src/QR/FullPivotingHouseholderQR.h @@ -45,14 +45,14 @@ template<typename MatrixType> class FullPivotingHouseholderQR { public: - + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) }; - + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; @@ -106,13 +106,13 @@ template<typename MatrixType> class FullPivotingHouseholderQR } FullPivotingHouseholderQR& compute(const MatrixType& matrix); - + const IntRowVectorType& colsPermutation() const { ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); return m_cols_permutation; } - + const IntColVectorType& rowsTranspositions() const { ei_assert(m_isInitialized && "FullPivotingHouseholderQR is not initialized."); @@ -147,7 +147,7 @@ template<typename MatrixType> class FullPivotingHouseholderQR * \sa absDeterminant(), MatrixBase::determinant() */ typename MatrixType::RealScalar logAbsDeterminant() const; - + /** \returns the rank of the matrix of which *this is the QR decomposition. * * \note This is computed at the time of the construction of the QR decomposition. This @@ -271,7 +271,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co int cols = matrix.cols(); int size = std::min(rows,cols); m_rank = size; - + m_qr = matrix; m_hCoeffs.resize(size); @@ -283,9 +283,9 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co IntRowVectorType cols_transpositions(matrix.cols()); m_cols_permutation.resize(matrix.cols()); int number_of_transpositions = 0; - + RealScalar biggest(0); - + for (int k = 0; k < size; ++k) { int row_of_biggest_in_corner, col_of_biggest_in_corner; @@ -297,7 +297,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co row_of_biggest_in_corner += k; col_of_biggest_in_corner += k; if(k==0) biggest = biggest_in_corner; - + // if the corner is negligible, then we have less than full rank, and we can finish early if(ei_isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) { @@ -336,7 +336,7 @@ FullPivotingHouseholderQR<MatrixType>& FullPivotingHouseholderQR<MatrixType>::co m_det_pq = (number_of_transpositions%2) ? -1 : 1; m_isInitialized = true; - + return *this; } @@ -358,13 +358,13 @@ bool FullPivotingHouseholderQR<MatrixType>::solve( } else return false; } - + const int rows = m_qr.rows(); const int cols = b.cols(); ei_assert(b.rows() == rows); - + typename OtherDerived::PlainMatrixType c(b); - + Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols); for (int k = 0; k < m_rank; ++k) { diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h index a89305869..39edda80c 100644 --- a/Eigen/src/QR/HouseholderQR.h +++ b/Eigen/src/QR/HouseholderQR.h @@ -56,12 +56,13 @@ template<typename MatrixType> class HouseholderQR Options = MatrixType::Options, DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime) }; - + typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType; typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType; typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; + typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType; /** * \brief Default Constructor. @@ -97,7 +98,12 @@ template<typename MatrixType> class HouseholderQR template<typename OtherDerived, typename ResultType> void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const; - MatrixQType matrixQ(void) const; + MatrixQType matrixQ() const; + + HouseholderSequenceType matrixQAsHouseholderSequence() const + { + return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); + } /** \returns a reference to the matrix where the Householder QR decomposition is stored * in a LAPACK-compatible way. @@ -169,7 +175,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& int rows = matrix.rows(); int cols = matrix.cols(); int size = std::min(rows,cols); - + m_qr = matrix; m_hCoeffs.resize(size); @@ -206,15 +212,7 @@ void HouseholderQR<MatrixType>::solve( result->resize(rows, cols); *result = b; - - Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols); - for (int k = 0; k < cols; ++k) - { - int remainingSize = rows-k; - - result->corner(BottomRight, remainingSize, cols) - .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0)); - } + result->applyOnTheLeft(matrixQAsHouseholderSequence().inverse()); const int rank = std::min(result->rows(), result->cols()); m_qr.corner(TopLeft, rank, rank) @@ -227,20 +225,7 @@ template<typename MatrixType> typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const { ei_assert(m_isInitialized && "HouseholderQR is not initialized."); - // compute the product H'_0 H'_1 ... H'_n-1, - // where H_k is the k-th Householder transformation I - h_k v_k v_k' - // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] - int rows = m_qr.rows(); - int cols = m_qr.cols(); - int size = std::min(rows,cols); - MatrixQType res = MatrixQType::Identity(rows, rows); - Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows); - for (int k = size-1; k >= 0; k--) - { - res.block(k, k, rows-k, rows-k) - .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k)); - } - return res; + return matrixQAsHouseholderSequence(); } #endif // EIGEN_HIDE_HEAVY_CODE diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index 2801ee077..4b69e67c4 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -25,6 +25,22 @@ #ifndef EIGEN_JACOBISVD_H #define EIGEN_JACOBISVD_H +// forward declarations (needed by ICC) +template<typename MatrixType, unsigned int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex> +struct ei_svd_precondition_2x2_block_to_be_real; + +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; + +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; + /** \ingroup SVD_Module * \nonstableyet * @@ -118,8 +134,8 @@ template<typename MatrixType, unsigned int Options> class JacobiSVD 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 +template<typename MatrixType, unsigned int Options> +struct ei_svd_precondition_2x2_block_to_be_real<MatrixType, Options, false> { typedef JacobiSVD<MatrixType, Options> SVD; static void run(typename SVD::WorkMatrixType&, JacobiSVD<MatrixType, Options>&, int, int) {} @@ -195,10 +211,7 @@ void ei_real_2x2_jacobi_svd(const MatrixType& matrix, int p, int q, *j_left = rot1 * j_right->transpose(); } -template<typename MatrixType, unsigned int Options, - bool PossiblyMoreRowsThanCols = (Options & AtLeastAsManyColsAsRows) == 0 - && (MatrixType::RowsAtCompileTime==Dynamic - || MatrixType::RowsAtCompileTime>MatrixType::ColsAtCompileTime)> +template<typename MatrixType, unsigned int Options, bool PossiblyMoreRowsThanCols> struct ei_svd_precondition_if_more_rows_than_cols { typedef JacobiSVD<MatrixType, Options> SVD; @@ -231,10 +244,7 @@ struct ei_svd_precondition_if_more_rows_than_cols<MatrixType, Options, true> } }; -template<typename MatrixType, unsigned int Options, - bool PossiblyMoreColsThanRows = (Options & AtLeastAsManyRowsAsCols) == 0 - && (MatrixType::ColsAtCompileTime==Dynamic - || MatrixType::ColsAtCompileTime>MatrixType::RowsAtCompileTime)> +template<typename MatrixType, unsigned int Options, bool PossiblyMoreColsThanRows> struct ei_svd_precondition_if_more_cols_than_rows { typedef JacobiSVD<MatrixType, Options> SVD; @@ -256,7 +266,7 @@ struct ei_svd_precondition_if_more_cols_than_rows<MatrixType, Options, true> MaxColsAtCompileTime = SVD::MaxColsAtCompileTime, MatrixOptions = SVD::MatrixOptions }; - + static bool run(const MatrixType& matrix, typename SVD::WorkMatrixType& work_matrix, SVD& svd) { int rows = matrix.rows(); 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 98d598809..708f177e8 100644 --- a/Eigen/src/Sparse/SuperLUSupport.h +++ b/Eigen/src/Sparse/SuperLUSupport.h @@ -161,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; } 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/test/CMakeLists.txt b/test/CMakeLists.txt index 03fbd48fc..f3c15612f 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -94,6 +94,7 @@ ei_add_test(basicstuff) ei_add_test(linearstructure) ei_add_test(cwiseop) ei_add_test(redux) +ei_add_test(visitor) ei_add_test(product_small) ei_add_test(product_large ${EI_OFLAG}) ei_add_test(product_extra ${EI_OFLAG}) @@ -115,6 +116,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}) 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/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 index f0d025283..b92dd5449 100644 --- a/test/conservative_resize.cpp +++ b/test/conservative_resize.cpp @@ -65,7 +65,7 @@ void run_matrix_tests() const int rows = ei_random<int>(50,75); const int cols = ei_random<int>(50,75); m = n = MatrixType::Random(50,50); - m.conservativeResize(rows,cols,true); + 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) ); @@ -102,7 +102,7 @@ void run_vector_tests() { const int size = ei_random<int>(50,100); m = n = MatrixType::Random(50); - m.conservativeResize(size,true); + m.conservativeResizeLike(MatrixType::Zero(size)); VERIFY_IS_APPROX(m.segment(0,50), n); VERIFY( size<=50 || m.segment(50,size-50).sum() == Scalar(0) ); } 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/geo_transformations.cpp b/test/geo_transformations.cpp index 914dbaf74..2b05f2457 100644 --- a/test/geo_transformations.cpp +++ b/test/geo_transformations.cpp @@ -296,10 +296,10 @@ template<typename Scalar, int Mode> void transformations(void) t0.setIdentity(); t0.translate(v0); t0.linear().setRandom(); - VERIFY_IS_APPROX(t0.inverse(Affine), t0.matrix().inverse()); + VERIFY_IS_APPROX(t0.inverse(Affine).matrix(), t0.matrix().inverse()); t0.setIdentity(); t0.translate(v0).rotate(q1); - VERIFY_IS_APPROX(t0.inverse(Isometry), t0.matrix().inverse()); + VERIFY_IS_APPROX(t0.inverse(Isometry).matrix(), t0.matrix().inverse()); } // test extract rotation and aligned scaling diff --git a/test/householder.cpp b/test/householder.cpp index 7d300899f..b27279479 100644 --- a/test/householder.cpp +++ b/test/householder.cpp @@ -43,7 +43,7 @@ template<typename MatrixType> void householder(const MatrixType& m) Matrix<Scalar, EIGEN_ENUM_MAX(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime), 1> _tmp(std::max(rows,cols)); Scalar* tmp = &_tmp.coeffRef(0,0); - + Scalar beta; RealScalar alpha; EssentialVectorType essential; @@ -58,7 +58,7 @@ template<typename MatrixType> void householder(const MatrixType& m) v2 = v1; v1.applyHouseholderOnTheLeft(essential,beta,tmp); VERIFY_IS_APPROX(v1.norm(), v2.norm()); - + MatrixType m1(rows, cols), m2(rows, cols); @@ -72,7 +72,7 @@ template<typename MatrixType> void householder(const MatrixType& m) VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm()); VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m1(0,0)), ei_real(m1(0,0))); VERIFY_IS_APPROX(ei_real(m1(0,0)), alpha); - + v1 = VectorType::Random(rows); if(even) v1.end(rows-1).setZero(); SquareMatrixType m3(rows,rows), m4(rows,rows); @@ -84,6 +84,9 @@ template<typename MatrixType> void householder(const MatrixType& m) VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm()); VERIFY_IS_MUCH_SMALLER_THAN(ei_imag(m3(0,0)), ei_real(m3(0,0))); VERIFY_IS_APPROX(ei_real(m3(0,0)), alpha); + + // test householder sequence + // TODO test HouseholderSequence } void test_householder() diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 5940b8497..2e3f089a0 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -36,14 +36,14 @@ template<typename MatrixType, unsigned int Options> void svd(const MatrixType& m 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; @@ -53,7 +53,7 @@ template<typename MatrixType, unsigned int Options> void svd(const MatrixType& m 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); @@ -98,7 +98,7 @@ void test_jacobisvd() } 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>() )); diff --git a/test/main.h b/test/main.h index 619fc9e06..8c93e856c 100644 --- a/test/main.h +++ b/test/main.h @@ -40,6 +40,11 @@ #define DEFAULT_REPEAT 10 +#ifdef __ICC +// disable warning #279: controlling expression is constant +#pragma warning disable 279 +#endif + namespace Eigen { static std::vector<std::string> g_test_stack; diff --git a/test/mixingtypes.cpp b/test/mixingtypes.cpp index 7dc57e6f7..3e322c7fe 100644 --- a/test/mixingtypes.cpp +++ b/test/mixingtypes.cpp @@ -175,9 +175,9 @@ void test_mixingtypes() { // check that our operator new is indeed called: CALL_SUBTEST(mixingtypes<3>()); -// CALL_SUBTEST(mixingtypes<4>()); -// CALL_SUBTEST(mixingtypes<Dynamic>(20)); -// -// CALL_SUBTEST(mixingtypes_small<4>()); -// CALL_SUBTEST(mixingtypes_large(20)); + CALL_SUBTEST(mixingtypes<4>()); + CALL_SUBTEST(mixingtypes<Dynamic>(20)); + + CALL_SUBTEST(mixingtypes_small<4>()); + CALL_SUBTEST(mixingtypes_large(20)); } diff --git a/test/product_extra.cpp b/test/product_extra.cpp index fcec362a5..3ad99fc7a 100644 --- a/test/product_extra.cpp +++ b/test/product_extra.cpp @@ -104,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/qr.cpp b/test/qr.cpp index f2e2eda61..f185ac86e 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -78,7 +78,7 @@ template<typename MatrixType> void qr_invertible() m3 = MatrixType::Random(size,size); qr.solve(m3, &m2); VERIFY_IS_APPROX(m3, m1*m2); - + // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>(); diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 283855451..588a41e56 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -42,18 +42,18 @@ template<typename MatrixType> void qr() VERIFY(!qr.isInjective()); VERIFY(!qr.isInvertible()); VERIFY(!qr.isSurjective()); - + MatrixType r = qr.matrixQR(); // FIXME need better way to construct trapezoid for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); - + MatrixType b = qr.matrixQ() * r; MatrixType c = MatrixType::Zero(rows,cols); - + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); VERIFY_IS_APPROX(m1, c); - + MatrixType m2 = MatrixType::Random(cols,cols2); MatrixType m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); @@ -116,9 +116,7 @@ template<typename MatrixType> void qr_verify_assert() void test_qr_colpivoting() { - for(int i = 0; i < 1; i++) { - // FIXME : very weird bug here -// CALL_SUBTEST( qr(Matrix2f()) ); + for(int i = 0; i < 1; i++) { CALL_SUBTEST( qr<MatrixXf>() ); CALL_SUBTEST( qr<MatrixXd>() ); CALL_SUBTEST( qr<MatrixXcd>() ); diff --git a/test/qr_fullpivoting.cpp b/test/qr_fullpivoting.cpp index 525c669a5..3a37bcb46 100644 --- a/test/qr_fullpivoting.cpp +++ b/test/qr_fullpivoting.cpp @@ -46,14 +46,14 @@ template<typename MatrixType> void qr() MatrixType r = qr.matrixQR(); // FIXME need better way to construct trapezoid for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0); - + MatrixType b = qr.matrixQ() * r; MatrixType c = MatrixType::Zero(rows,cols); - + for(int i = 0; i < cols; ++i) c.col(qr.colsPermutation().coeff(i)) = b.col(i); VERIFY_IS_APPROX(m1, c); - + MatrixType m2 = MatrixType::Random(cols,cols2); MatrixType m3 = m1*m2; m2 = MatrixType::Random(cols,cols2); @@ -88,7 +88,7 @@ template<typename MatrixType> void qr_invertible() m3 = MatrixType::Random(size,size); VERIFY(qr.solve(m3, &m2)); VERIFY_IS_APPROX(m3, m1*m2); - + // now construct a matrix with prescribed determinant m1.setZero(); for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>(); diff --git a/test/redux.cpp b/test/redux.cpp index 2a0dc97f1..951b34bca 100644 --- a/test/redux.cpp +++ b/test/redux.cpp @@ -120,7 +120,7 @@ void test_redux() CALL_SUBTEST( matrixRedux(MatrixXi(8, 12)) ); } for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST( vectorRedux(VectorXf(5)) ); + CALL_SUBTEST( vectorRedux(Vector4f()) ); CALL_SUBTEST( vectorRedux(VectorXd(10)) ); CALL_SUBTEST( vectorRedux(VectorXf(33)) ); } 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/visitor.cpp b/test/visitor.cpp new file mode 100644 index 000000000..b78782b78 --- /dev/null +++ b/test/visitor.cpp @@ -0,0 +1,131 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2008 Benoit Jacob <jacob.benoit.1@gmail.com> +// +// Eigen is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 3 of the License, or (at your option) any later version. +// +// Alternatively, you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation; either version 2 of +// the License, or (at your option) any later version. +// +// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY +// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License and a copy of the GNU General Public License along with +// Eigen. If not, see <http://www.gnu.org/licenses/>. + +#include "main.h" + +template<typename MatrixType> void matrixVisitor(const MatrixType& p) +{ + typedef typename MatrixType::Scalar Scalar; + + int rows = p.rows(); + int cols = p.cols(); + + // construct a random matrix where all coefficients are different + MatrixType m; + m = MatrixType::Random(rows, cols); + for(int i = 0; i < m.size(); i++) + for(int i2 = 0; i2 < i; i2++) + while(m(i) == m(i2)) // yes, == + m(i) = ei_random<Scalar>(); + + Scalar minc = Scalar(1000), maxc = Scalar(-1000); + int minrow,mincol,maxrow,maxcol; + for(int j = 0; j < cols; j++) + for(int i = 0; i < rows; i++) + { + if(m(i,j) < minc) + { + minc = m(i,j); + minrow = i; + mincol = j; + } + if(m(i,j) > maxc) + { + maxc = m(i,j); + maxrow = i; + maxcol = j; + } + } + int eigen_minrow, eigen_mincol, eigen_maxrow, eigen_maxcol; + Scalar eigen_minc, eigen_maxc; + eigen_minc = m.minCoeff(&eigen_minrow,&eigen_mincol); + eigen_maxc = m.maxCoeff(&eigen_maxrow,&eigen_maxcol); + VERIFY(minrow == eigen_minrow); + VERIFY(maxrow == eigen_maxrow); + VERIFY(mincol == eigen_mincol); + VERIFY(maxcol == eigen_maxcol); + VERIFY_IS_APPROX(minc, eigen_minc); + VERIFY_IS_APPROX(maxc, eigen_maxc); + VERIFY_IS_APPROX(minc, m.minCoeff()); + VERIFY_IS_APPROX(maxc, m.maxCoeff()); +} + +template<typename VectorType> void vectorVisitor(const VectorType& w) +{ + typedef typename VectorType::Scalar Scalar; + + int size = w.size(); + + // construct a random vector where all coefficients are different + VectorType v; + v = VectorType::Random(size); + for(int i = 0; i < size; i++) + for(int i2 = 0; i2 < i; i2++) + while(v(i) == v(i2)) // yes, == + v(i) = ei_random<Scalar>(); + + Scalar minc = Scalar(1000), maxc = Scalar(-1000); + int minidx,maxidx; + for(int i = 0; i < size; i++) + { + if(v(i) < minc) + { + minc = v(i); + minidx = i; + } + if(v(i) > maxc) + { + maxc = v(i); + maxidx = i; + } + } + int eigen_minidx, eigen_maxidx; + Scalar eigen_minc, eigen_maxc; + eigen_minc = v.minCoeff(&eigen_minidx); + eigen_maxc = v.maxCoeff(&eigen_maxidx); + VERIFY(minidx == eigen_minidx); + VERIFY(maxidx == eigen_maxidx); + VERIFY_IS_APPROX(minc, eigen_minc); + VERIFY_IS_APPROX(maxc, eigen_maxc); + VERIFY_IS_APPROX(minc, v.minCoeff()); + VERIFY_IS_APPROX(maxc, v.maxCoeff()); +} + +void test_visitor() +{ + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( matrixVisitor(Matrix<float, 1, 1>()) ); + CALL_SUBTEST( matrixVisitor(Matrix2f()) ); + CALL_SUBTEST( matrixVisitor(Matrix4d()) ); + CALL_SUBTEST( matrixVisitor(MatrixXd(8, 12)) ); + CALL_SUBTEST( matrixVisitor(Matrix<double,Dynamic,Dynamic,RowMajor>(20, 20)) ); + CALL_SUBTEST( matrixVisitor(MatrixXi(8, 12)) ); + } + for(int i = 0; i < g_repeat; i++) { + CALL_SUBTEST( vectorVisitor(Vector4f()) ); + CALL_SUBTEST( vectorVisitor(VectorXd(10)) ); + CALL_SUBTEST( vectorVisitor(RowVectorXd(10)) ); + CALL_SUBTEST( vectorVisitor(VectorXf(33)) ); + } +} diff --git a/unsupported/Eigen/AlignedVector3 b/unsupported/Eigen/AlignedVector3 index 854f0601f..f20fad6d1 100644 --- a/unsupported/Eigen/AlignedVector3 +++ b/unsupported/Eigen/AlignedVector3 @@ -63,37 +63,37 @@ template<typename _Scalar> class AlignedVector3 typedef Matrix<_Scalar,4,1> CoeffType; CoeffType m_coeffs; public: - + EIGEN_GENERIC_PUBLIC_INTERFACE(AlignedVector3) using Base::operator*; - + inline int rows() const { return 3; } inline int cols() const { return 1; } - + inline const Scalar& coeff(int row, int col) const { return m_coeffs.coeff(row, col); } - + inline Scalar& coeffRef(int row, int col) { return m_coeffs.coeffRef(row, col); } - + inline const Scalar& coeff(int index) const { return m_coeffs.coeff(index); } inline Scalar& coeffRef(int index) { return m_coeffs.coeffRef(index);} - - + + inline AlignedVector3(const Scalar& x, const Scalar& y, const Scalar& z) : m_coeffs(x, y, z, Scalar(0)) {} - + inline AlignedVector3(const AlignedVector3& other) - : m_coeffs(other.m_coeffs) + : Base(), m_coeffs(other.m_coeffs) {} - + template<typename XprType, int Size=XprType::SizeAtCompileTime> struct generic_assign_selector {}; - + template<typename XprType> struct generic_assign_selector<XprType,4> { inline static void run(AlignedVector3& dest, const XprType& src) @@ -101,7 +101,7 @@ template<typename _Scalar> class AlignedVector3 dest.m_coeffs = src; } }; - + template<typename XprType> struct generic_assign_selector<XprType,3> { inline static void run(AlignedVector3& dest, const XprType& src) @@ -110,49 +110,48 @@ template<typename _Scalar> class AlignedVector3 dest.m_coeffs.w() = Scalar(0); } }; - + template<typename Derived> inline explicit AlignedVector3(const MatrixBase<Derived>& other) { generic_assign_selector<Derived>::run(*this,other.derived()); } - + inline AlignedVector3& operator=(const AlignedVector3& other) { m_coeffs = other.m_coeffs; return *this; } - - + + inline AlignedVector3 operator+(const AlignedVector3& other) const { return AlignedVector3(m_coeffs + other.m_coeffs); } - + inline AlignedVector3& operator+=(const AlignedVector3& other) { m_coeffs += other.m_coeffs; return *this; } - + inline AlignedVector3 operator-(const AlignedVector3& other) const { return AlignedVector3(m_coeffs - other.m_coeffs); } - + inline AlignedVector3 operator-=(const AlignedVector3& other) { m_coeffs -= other.m_coeffs; return *this; } - + inline AlignedVector3 operator*(const Scalar& s) const { return AlignedVector3(m_coeffs * s); } - + inline friend AlignedVector3 operator*(const Scalar& s,const AlignedVector3& vec) { return AlignedVector3(s * vec.m_coeffs); } - + inline AlignedVector3& operator*=(const Scalar& s) { m_coeffs *= s; return *this; } - + inline AlignedVector3 operator/(const Scalar& s) const { return AlignedVector3(m_coeffs / s); } - + inline AlignedVector3& operator/=(const Scalar& s) { m_coeffs /= s; return *this; } - + inline Scalar dot(const AlignedVector3& other) const { ei_assert(m_coeffs.w()==Scalar(0)); ei_assert(other.m_coeffs.w()==Scalar(0)); - Scalar r = m_coeffs.dot(other.m_coeffs); return m_coeffs.dot(other.m_coeffs); } @@ -165,29 +164,29 @@ template<typename _Scalar> class AlignedVector3 { return AlignedVector3(m_coeffs / norm()); } - + inline Scalar sum() const { ei_assert(m_coeffs.w()==Scalar(0)); return m_coeffs.sum(); } - + inline Scalar squaredNorm() const { ei_assert(m_coeffs.w()==Scalar(0)); return m_coeffs.squaredNorm(); } - + inline Scalar norm() const { return ei_sqrt(squaredNorm()); } - + inline AlignedVector3 cross(const AlignedVector3& other) const { return AlignedVector3(m_coeffs.cross3(other.m_coeffs)); } - + template<typename Derived> inline bool isApprox(const MatrixBase<Derived>& other, RealScalar eps=precision<Scalar>()) const { diff --git a/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h b/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h index bf5b79955..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,260 +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.noalias() = M * M; - tmp = b[3]*M2 + b[1]*Id; - U.noalias() = M * tmp; - 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.noalias() = M * M; - MatrixType M4 = M2 * M2; - tmp = b[5]*M4 + b[3]*M2 + b[1]*Id; - U.noalias() = M * tmp; - 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.noalias() = M * M; - MatrixType M4 = M2 * M2; - MatrixType M6 = M4 * M2; - tmp = b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; - U.noalias() = M * tmp; - 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.noalias() = M * M; - MatrixType M4 = M2 * M2; - MatrixType M6 = M4 * M2; - MatrixType M8 = M6 * M2; - tmp = b[9]*M8 + b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; - U.noalias() = M * tmp; - 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.noalias() = M * M; - MatrixType M4 = M2 * M2; - MatrixType M6 = M4 * M2; - V = b[13]*M6 + b[11]*M4 + b[9]*M2; - tmp.noalias() = M6 * V; - tmp += b[7]*M6 + b[5]*M4 + b[3]*M2 + b[1]*Id; - U.noalias() = M * tmp; - tmp = b[12]*M6 + b[10]*M4 + b[8]*M2; - V.noalias() = M6 * tmp; - 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.925724783138660f; - *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(M.rows(), M.cols()); - MatrixType den(M.rows(), M.cols()); - MatrixType U(M.rows(), M.cols()); - MatrixType V(M.rows(), M.cols()); - MatrixType Id = MatrixType::Identity(M.rows(), M.cols()); - float l1norm = static_cast<float>(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 |