diff options
| -rw-r--r-- | Eigen/Core | 1 | ||||
| -rw-r--r-- | Eigen/Householder | 1 | ||||
| -rw-r--r-- | Eigen/src/Core/AnyMatrixBase.h | 153 | ||||
| -rw-r--r-- | Eigen/src/Core/MatrixBase.h | 58 | ||||
| -rw-r--r-- | Eigen/src/Core/Product.h | 14 | ||||
| -rw-r--r-- | Eigen/src/Core/TriangularMatrix.h | 12 | ||||
| -rw-r--r-- | Eigen/src/Core/util/ForwardDeclarations.h | 1 | ||||
| -rw-r--r-- | Eigen/src/Householder/HouseholderSequence.h | 146 | ||||
| -rw-r--r-- | Eigen/src/QR/HouseholderQR.h | 37 | ||||
| -rw-r--r-- | test/householder.cpp | 9 | ||||
| -rw-r--r-- | test/jacobisvd.cpp | 8 | ||||
| -rw-r--r-- | test/qr.cpp | 2 |
12 files changed, 339 insertions, 103 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/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index f267aa34d..4835f167c 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -26,44 +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 -{ - 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 @@ -96,7 +58,6 @@ template<typename Derived> class MatrixBase #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*; @@ -301,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); @@ -436,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; 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/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/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/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h new file mode 100644 index 000000000..86395213b --- /dev/null +++ b/Eigen/src/Householder/HouseholderSequence.h @@ -0,0 +1,146 @@ +// 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)); + } + } + + protected: + + typename VectorsType::Nested m_vectors; + typename CoeffsType::Nested m_coeffs; + bool m_trans; +}; + +#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H 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/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/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>(); |