diff options
Diffstat (limited to 'Eigen')
24 files changed, 582 insertions, 247 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; } |