Squashed 'third_party/eigen/' content from commit 61d72f6
Change-Id: Iccc90fa0b55ab44037f018046d2fcffd90d9d025
git-subtree-dir: third_party/eigen
git-subtree-split: 61d72f6383cfa842868c53e30e087b0258177257
diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
new file mode 100644
index 0000000..3db0c01
--- /dev/null
+++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h
@@ -0,0 +1,373 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#ifndef EIGEN_HESSENBERGDECOMPOSITION_H
+#define EIGEN_HESSENBERGDECOMPOSITION_H
+
+namespace Eigen {
+
+namespace internal {
+
+template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
+template<typename MatrixType>
+struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
+{
+ typedef MatrixType ReturnType;
+};
+
+}
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+ *
+ *
+ * \class HessenbergDecomposition
+ *
+ * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation
+ *
+ * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
+ *
+ * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In
+ * the real case, the Hessenberg decomposition consists of an orthogonal
+ * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H
+ * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its
+ * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the
+ * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition
+ * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is,
+ * \f$ Q^{-1} = Q^* \f$).
+ *
+ * Call the function compute() to compute the Hessenberg decomposition of a
+ * given matrix. Alternatively, you can use the
+ * HessenbergDecomposition(const MatrixType&) constructor which computes the
+ * Hessenberg decomposition at construction time. Once the decomposition is
+ * computed, you can use the matrixH() and matrixQ() functions to construct
+ * the matrices H and Q in the decomposition.
+ *
+ * The documentation for matrixH() contains an example of the typical use of
+ * this class.
+ *
+ * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module"
+ */
+template<typename _MatrixType> class HessenbergDecomposition
+{
+ public:
+
+ /** \brief Synonym for the template parameter \p _MatrixType. */
+ typedef _MatrixType MatrixType;
+
+ enum {
+ Size = MatrixType::RowsAtCompileTime,
+ SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
+ Options = MatrixType::Options,
+ MaxSize = MatrixType::MaxRowsAtCompileTime,
+ MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
+ };
+
+ /** \brief Scalar type for matrices of type #MatrixType. */
+ typedef typename MatrixType::Scalar Scalar;
+ typedef typename MatrixType::Index Index;
+
+ /** \brief Type for vector of Householder coefficients.
+ *
+ * This is column vector with entries of type #Scalar. The length of the
+ * vector is one less than the size of #MatrixType, if it is a fixed-side
+ * type.
+ */
+ typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
+
+ /** \brief Return type of matrixQ() */
+ typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType;
+
+ typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
+
+ /** \brief Default constructor; the decomposition will be computed later.
+ *
+ * \param [in] size The size of the matrix whose Hessenberg decomposition will be computed.
+ *
+ * The default constructor is useful in cases in which the user intends to
+ * perform decompositions via compute(). The \p size parameter is only
+ * used as a hint. It is not an error to give a wrong \p size, but it may
+ * impair performance.
+ *
+ * \sa compute() for an example.
+ */
+ HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
+ : m_matrix(size,size),
+ m_temp(size),
+ m_isInitialized(false)
+ {
+ if(size>1)
+ m_hCoeffs.resize(size-1);
+ }
+
+ /** \brief Constructor; computes Hessenberg decomposition of given matrix.
+ *
+ * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed.
+ *
+ * This constructor calls compute() to compute the Hessenberg
+ * decomposition.
+ *
+ * \sa matrixH() for an example.
+ */
+ HessenbergDecomposition(const MatrixType& matrix)
+ : m_matrix(matrix),
+ m_temp(matrix.rows()),
+ m_isInitialized(false)
+ {
+ if(matrix.rows()<2)
+ {
+ m_isInitialized = true;
+ return;
+ }
+ m_hCoeffs.resize(matrix.rows()-1,1);
+ _compute(m_matrix, m_hCoeffs, m_temp);
+ m_isInitialized = true;
+ }
+
+ /** \brief Computes Hessenberg decomposition of given matrix.
+ *
+ * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed.
+ * \returns Reference to \c *this
+ *
+ * The Hessenberg decomposition is computed by bringing the columns of the
+ * matrix successively in the required form using Householder reflections
+ * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
+ * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
+ * denotes the size of the given matrix.
+ *
+ * This method reuses of the allocated data in the HessenbergDecomposition
+ * object.
+ *
+ * Example: \include HessenbergDecomposition_compute.cpp
+ * Output: \verbinclude HessenbergDecomposition_compute.out
+ */
+ HessenbergDecomposition& compute(const MatrixType& matrix)
+ {
+ m_matrix = matrix;
+ if(matrix.rows()<2)
+ {
+ m_isInitialized = true;
+ return *this;
+ }
+ m_hCoeffs.resize(matrix.rows()-1,1);
+ _compute(m_matrix, m_hCoeffs, m_temp);
+ m_isInitialized = true;
+ return *this;
+ }
+
+ /** \brief Returns the Householder coefficients.
+ *
+ * \returns a const reference to the vector of Householder coefficients
+ *
+ * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
+ * or the member function compute(const MatrixType&) has been called
+ * before to compute the Hessenberg decomposition of a matrix.
+ *
+ * The Householder coefficients allow the reconstruction of the matrix
+ * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
+ *
+ * \sa packedMatrix(), \ref Householder_Module "Householder module"
+ */
+ const CoeffVectorType& householderCoefficients() const
+ {
+ eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
+ return m_hCoeffs;
+ }
+
+ /** \brief Returns the internal representation of the decomposition
+ *
+ * \returns a const reference to a matrix with the internal representation
+ * of the decomposition.
+ *
+ * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
+ * or the member function compute(const MatrixType&) has been called
+ * before to compute the Hessenberg decomposition of a matrix.
+ *
+ * The returned matrix contains the following information:
+ * - the upper part and lower sub-diagonal represent the Hessenberg matrix H
+ * - the rest of the lower part contains the Householder vectors that, combined with
+ * Householder coefficients returned by householderCoefficients(),
+ * allows to reconstruct the matrix Q as
+ * \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
+ * Here, the matrices \f$ H_i \f$ are the Householder transformations
+ * \f$ H_i = (I - h_i v_i v_i^T) \f$
+ * where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
+ * \f$ v_i \f$ is the Householder vector defined by
+ * \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
+ * with M the matrix returned by this function.
+ *
+ * See LAPACK for further details on this packed storage.
+ *
+ * Example: \include HessenbergDecomposition_packedMatrix.cpp
+ * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
+ *
+ * \sa householderCoefficients()
+ */
+ const MatrixType& packedMatrix() const
+ {
+ eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
+ return m_matrix;
+ }
+
+ /** \brief Reconstructs the orthogonal matrix Q in the decomposition
+ *
+ * \returns object representing the matrix Q
+ *
+ * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
+ * or the member function compute(const MatrixType&) has been called
+ * before to compute the Hessenberg decomposition of a matrix.
+ *
+ * This function returns a light-weight object of template class
+ * HouseholderSequence. You can either apply it directly to a matrix or
+ * you can convert it to a matrix of type #MatrixType.
+ *
+ * \sa matrixH() for an example, class HouseholderSequence
+ */
+ HouseholderSequenceType matrixQ() const
+ {
+ eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
+ return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
+ .setLength(m_matrix.rows() - 1)
+ .setShift(1);
+ }
+
+ /** \brief Constructs the Hessenberg matrix H in the decomposition
+ *
+ * \returns expression object representing the matrix H
+ *
+ * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
+ * or the member function compute(const MatrixType&) has been called
+ * before to compute the Hessenberg decomposition of a matrix.
+ *
+ * The object returned by this function constructs the Hessenberg matrix H
+ * when it is assigned to a matrix or otherwise evaluated. The matrix H is
+ * constructed from the packed matrix as returned by packedMatrix(): The
+ * upper part (including the subdiagonal) of the packed matrix contains
+ * the matrix H. It may sometimes be better to directly use the packed
+ * matrix instead of constructing the matrix H.
+ *
+ * Example: \include HessenbergDecomposition_matrixH.cpp
+ * Output: \verbinclude HessenbergDecomposition_matrixH.out
+ *
+ * \sa matrixQ(), packedMatrix()
+ */
+ MatrixHReturnType matrixH() const
+ {
+ eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
+ return MatrixHReturnType(*this);
+ }
+
+ private:
+
+ typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
+ typedef typename NumTraits<Scalar>::Real RealScalar;
+ static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
+
+ protected:
+ MatrixType m_matrix;
+ CoeffVectorType m_hCoeffs;
+ VectorType m_temp;
+ bool m_isInitialized;
+};
+
+/** \internal
+ * Performs a tridiagonal decomposition of \a matA in place.
+ *
+ * \param matA the input selfadjoint matrix
+ * \param hCoeffs returned Householder coefficients
+ *
+ * The result is written in the lower triangular part of \a matA.
+ *
+ * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1.
+ *
+ * \sa packedMatrix()
+ */
+template<typename MatrixType>
+void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
+{
+ eigen_assert(matA.rows()==matA.cols());
+ Index n = matA.rows();
+ temp.resize(n);
+ for (Index i = 0; i<n-1; ++i)
+ {
+ // let's consider the vector v = i-th column starting at position i+1
+ Index remainingSize = n-i-1;
+ RealScalar beta;
+ Scalar h;
+ matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
+ matA.col(i).coeffRef(i+1) = beta;
+ hCoeffs.coeffRef(i) = h;
+
+ // Apply similarity transformation to remaining columns,
+ // i.e., compute A = H A H'
+
+ // A = H A
+ matA.bottomRightCorner(remainingSize, remainingSize)
+ .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
+
+ // A = A H'
+ matA.rightCols(remainingSize)
+ .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0));
+ }
+}
+
+namespace internal {
+
+/** \eigenvalues_module \ingroup Eigenvalues_Module
+ *
+ *
+ * \brief Expression type for return value of HessenbergDecomposition::matrixH()
+ *
+ * \tparam MatrixType type of matrix in the Hessenberg decomposition
+ *
+ * Objects of this type represent the Hessenberg matrix in the Hessenberg
+ * decomposition of some matrix. The object holds a reference to the
+ * HessenbergDecomposition class until the it is assigned or evaluated for
+ * some other reason (the reference should remain valid during the life time
+ * of this object). This class is the return type of
+ * HessenbergDecomposition::matrixH(); there is probably no other use for this
+ * class.
+ */
+template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
+: public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
+{
+ typedef typename MatrixType::Index Index;
+ public:
+ /** \brief Constructor.
+ *
+ * \param[in] hess Hessenberg decomposition
+ */
+ HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
+
+ /** \brief Hessenberg matrix in decomposition.
+ *
+ * \param[out] result Hessenberg matrix in decomposition \p hess which
+ * was passed to the constructor
+ */
+ template <typename ResultType>
+ inline void evalTo(ResultType& result) const
+ {
+ result = m_hess.packedMatrix();
+ Index n = result.rows();
+ if (n>2)
+ result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
+ }
+
+ Index rows() const { return m_hess.packedMatrix().rows(); }
+ Index cols() const { return m_hess.packedMatrix().cols(); }
+
+ protected:
+ const HessenbergDecomposition<MatrixType>& m_hess;
+};
+
+} // end namespace internal
+
+} // end namespace Eigen
+
+#endif // EIGEN_HESSENBERGDECOMPOSITION_H