aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
-rw-r--r--Eigen/Core1
-rw-r--r--Eigen/Householder1
-rw-r--r--Eigen/src/Core/AnyMatrixBase.h153
-rw-r--r--Eigen/src/Core/MatrixBase.h58
-rw-r--r--Eigen/src/Core/Product.h14
-rw-r--r--Eigen/src/Core/TriangularMatrix.h12
-rw-r--r--Eigen/src/Core/util/ForwardDeclarations.h1
-rw-r--r--Eigen/src/Householder/HouseholderSequence.h146
-rw-r--r--Eigen/src/QR/HouseholderQR.h37
-rw-r--r--test/householder.cpp9
-rw-r--r--test/jacobisvd.cpp8
-rw-r--r--test/qr.cpp2
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>();