Jaesung Oh / Eigen

Eigne Matrix Class Library

Dependents:   MPC_current_control HydraulicControlBoard_SW AHRS Test_ekf ... more

src/Eigenvalues/HessenbergDecomposition.h@1:3b8049da21b8, 2019-09-24 (annotated)

Committer:
jsoh91
Date:
Tue Sep 24 00:18:23 2019 +0000
Revision:
1:3b8049da21b8
Parent:
0:13a5d365ba16 
ignore and revise some of error parts

Who changed what in which revision?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
ykuroda 0:13a5d365ba16 5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
ykuroda 0:13a5d365ba16 6 //
ykuroda 0:13a5d365ba16 7 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 8 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 10
ykuroda 0:13a5d365ba16 11 #ifndef EIGEN_HESSENBERGDECOMPOSITION_H
ykuroda 0:13a5d365ba16 12 #define EIGEN_HESSENBERGDECOMPOSITION_H
ykuroda 0:13a5d365ba16 13
ykuroda 0:13a5d365ba16 14 namespace Eigen {
ykuroda 0:13a5d365ba16 15
ykuroda 0:13a5d365ba16 16 namespace internal {
ykuroda 0:13a5d365ba16 17
ykuroda 0:13a5d365ba16 18 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
ykuroda 0:13a5d365ba16 19 template<typename MatrixType>
ykuroda 0:13a5d365ba16 20 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
ykuroda 0:13a5d365ba16 21 {
ykuroda 0:13a5d365ba16 22 typedef MatrixType ReturnType;
ykuroda 0:13a5d365ba16 23 };
ykuroda 0:13a5d365ba16 24
ykuroda 0:13a5d365ba16 25 }
ykuroda 0:13a5d365ba16 26
ykuroda 0:13a5d365ba16 27 /** \eigenvalues_module \ingroup Eigenvalues_Module
ykuroda 0:13a5d365ba16 28 *
ykuroda 0:13a5d365ba16 29 *
ykuroda 0:13a5d365ba16 30 * \class HessenbergDecomposition
ykuroda 0:13a5d365ba16 31 *
ykuroda 0:13a5d365ba16 32 * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation
ykuroda 0:13a5d365ba16 33 *
ykuroda 0:13a5d365ba16 34 * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition
ykuroda 0:13a5d365ba16 35 *
ykuroda 0:13a5d365ba16 36 * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In
ykuroda 0:13a5d365ba16 37 * the real case, the Hessenberg decomposition consists of an orthogonal
ykuroda 0:13a5d365ba16 38 * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H
ykuroda 0:13a5d365ba16 39 * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its
ykuroda 0:13a5d365ba16 40 * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the
ykuroda 0:13a5d365ba16 41 * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition
ykuroda 0:13a5d365ba16 42 * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is,
ykuroda 0:13a5d365ba16 43 * \f$ Q^{-1} = Q^* \f$).
ykuroda 0:13a5d365ba16 44 *
ykuroda 0:13a5d365ba16 45 * Call the function compute() to compute the Hessenberg decomposition of a
ykuroda 0:13a5d365ba16 46 * given matrix. Alternatively, you can use the
ykuroda 0:13a5d365ba16 47 * HessenbergDecomposition(const MatrixType&) constructor which computes the
ykuroda 0:13a5d365ba16 48 * Hessenberg decomposition at construction time. Once the decomposition is
ykuroda 0:13a5d365ba16 49 * computed, you can use the matrixH() and matrixQ() functions to construct
ykuroda 0:13a5d365ba16 50 * the matrices H and Q in the decomposition.
ykuroda 0:13a5d365ba16 51 *
ykuroda 0:13a5d365ba16 52 * The documentation for matrixH() contains an example of the typical use of
ykuroda 0:13a5d365ba16 53 * this class.
ykuroda 0:13a5d365ba16 54 *
ykuroda 0:13a5d365ba16 55 * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module"
ykuroda 0:13a5d365ba16 56 */
ykuroda 0:13a5d365ba16 57 template<typename _MatrixType> class HessenbergDecomposition
ykuroda 0:13a5d365ba16 58 {
ykuroda 0:13a5d365ba16 59 public:
ykuroda 0:13a5d365ba16 60
ykuroda 0:13a5d365ba16 61 /** \brief Synonym for the template parameter \p _MatrixType. */
ykuroda 0:13a5d365ba16 62 typedef _MatrixType MatrixType;
ykuroda 0:13a5d365ba16 63
ykuroda 0:13a5d365ba16 64 enum {
ykuroda 0:13a5d365ba16 65 Size = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 66 SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
ykuroda 0:13a5d365ba16 67 Options = MatrixType::Options,
ykuroda 0:13a5d365ba16 68 MaxSize = MatrixType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 69 MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
ykuroda 0:13a5d365ba16 70 };
ykuroda 0:13a5d365ba16 71
ykuroda 0:13a5d365ba16 72 /** \brief Scalar type for matrices of type #MatrixType. */
ykuroda 0:13a5d365ba16 73 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 74 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 75
ykuroda 0:13a5d365ba16 76 /** \brief Type for vector of Householder coefficients.
ykuroda 0:13a5d365ba16 77 *
ykuroda 0:13a5d365ba16 78 * This is column vector with entries of type #Scalar. The length of the
ykuroda 0:13a5d365ba16 79 * vector is one less than the size of #MatrixType, if it is a fixed-side
ykuroda 0:13a5d365ba16 80 * type.
ykuroda 0:13a5d365ba16 81 */
ykuroda 0:13a5d365ba16 82 typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
ykuroda 0:13a5d365ba16 83
ykuroda 0:13a5d365ba16 84 /** \brief Return type of matrixQ() */
ykuroda 0:13a5d365ba16 85 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType;
ykuroda 0:13a5d365ba16 86
ykuroda 0:13a5d365ba16 87 typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
ykuroda 0:13a5d365ba16 88
ykuroda 0:13a5d365ba16 89 /** \brief Default constructor; the decomposition will be computed later.
ykuroda 0:13a5d365ba16 90 *
ykuroda 0:13a5d365ba16 91 * \param [in] size The size of the matrix whose Hessenberg decomposition will be computed.
ykuroda 0:13a5d365ba16 92 *
ykuroda 0:13a5d365ba16 93 * The default constructor is useful in cases in which the user intends to
ykuroda 0:13a5d365ba16 94 * perform decompositions via compute(). The \p size parameter is only
ykuroda 0:13a5d365ba16 95 * used as a hint. It is not an error to give a wrong \p size, but it may
ykuroda 0:13a5d365ba16 96 * impair performance.
ykuroda 0:13a5d365ba16 97 *
ykuroda 0:13a5d365ba16 98 * \sa compute() for an example.
ykuroda 0:13a5d365ba16 99 */
ykuroda 0:13a5d365ba16 100 HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
ykuroda 0:13a5d365ba16 101 : m_matrix(size,size),
ykuroda 0:13a5d365ba16 102 m_temp(size),
ykuroda 0:13a5d365ba16 103 m_isInitialized(false)
ykuroda 0:13a5d365ba16 104 {
ykuroda 0:13a5d365ba16 105 if(size>1)
ykuroda 0:13a5d365ba16 106 m_hCoeffs.resize(size-1);
ykuroda 0:13a5d365ba16 107 }
ykuroda 0:13a5d365ba16 108
ykuroda 0:13a5d365ba16 109 /** \brief Constructor; computes Hessenberg decomposition of given matrix.
ykuroda 0:13a5d365ba16 110 *
ykuroda 0:13a5d365ba16 111 * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed.
ykuroda 0:13a5d365ba16 112 *
ykuroda 0:13a5d365ba16 113 * This constructor calls compute() to compute the Hessenberg
ykuroda 0:13a5d365ba16 114 * decomposition.
ykuroda 0:13a5d365ba16 115 *
ykuroda 0:13a5d365ba16 116 * \sa matrixH() for an example.
ykuroda 0:13a5d365ba16 117 */
ykuroda 0:13a5d365ba16 118 HessenbergDecomposition(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 119 : m_matrix(matrix),
ykuroda 0:13a5d365ba16 120 m_temp(matrix.rows()),
ykuroda 0:13a5d365ba16 121 m_isInitialized(false)
ykuroda 0:13a5d365ba16 122 {
ykuroda 0:13a5d365ba16 123 if(matrix.rows()<2)
ykuroda 0:13a5d365ba16 124 {
ykuroda 0:13a5d365ba16 125 m_isInitialized = true;
ykuroda 0:13a5d365ba16 126 return;
ykuroda 0:13a5d365ba16 127 }
ykuroda 0:13a5d365ba16 128 m_hCoeffs.resize(matrix.rows()-1,1);
ykuroda 0:13a5d365ba16 129 _compute(m_matrix, m_hCoeffs, m_temp);
ykuroda 0:13a5d365ba16 130 m_isInitialized = true;
ykuroda 0:13a5d365ba16 131 }
ykuroda 0:13a5d365ba16 132
ykuroda 0:13a5d365ba16 133 /** \brief Computes Hessenberg decomposition of given matrix.
ykuroda 0:13a5d365ba16 134 *
ykuroda 0:13a5d365ba16 135 * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed.
ykuroda 0:13a5d365ba16 136 * \returns Reference to \c *this
ykuroda 0:13a5d365ba16 137 *
ykuroda 0:13a5d365ba16 138 * The Hessenberg decomposition is computed by bringing the columns of the
ykuroda 0:13a5d365ba16 139 * matrix successively in the required form using Householder reflections
ykuroda 0:13a5d365ba16 140 * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
ykuroda 0:13a5d365ba16 141 * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
ykuroda 0:13a5d365ba16 142 * denotes the size of the given matrix.
ykuroda 0:13a5d365ba16 143 *
ykuroda 0:13a5d365ba16 144 * This method reuses of the allocated data in the HessenbergDecomposition
ykuroda 0:13a5d365ba16 145 * object.
ykuroda 0:13a5d365ba16 146 *
ykuroda 0:13a5d365ba16 147 * Example: \include HessenbergDecomposition_compute.cpp
ykuroda 0:13a5d365ba16 148 * Output: \verbinclude HessenbergDecomposition_compute.out
ykuroda 0:13a5d365ba16 149 */
ykuroda 0:13a5d365ba16 150 HessenbergDecomposition& compute(const MatrixType& matrix)
ykuroda 0:13a5d365ba16 151 {
ykuroda 0:13a5d365ba16 152 m_matrix = matrix;
ykuroda 0:13a5d365ba16 153 if(matrix.rows()<2)
ykuroda 0:13a5d365ba16 154 {
ykuroda 0:13a5d365ba16 155 m_isInitialized = true;
ykuroda 0:13a5d365ba16 156 return *this;
ykuroda 0:13a5d365ba16 157 }
ykuroda 0:13a5d365ba16 158 m_hCoeffs.resize(matrix.rows()-1,1);
ykuroda 0:13a5d365ba16 159 _compute(m_matrix, m_hCoeffs, m_temp);
ykuroda 0:13a5d365ba16 160 m_isInitialized = true;
ykuroda 0:13a5d365ba16 161 return *this;
ykuroda 0:13a5d365ba16 162 }
ykuroda 0:13a5d365ba16 163
ykuroda 0:13a5d365ba16 164 /** \brief Returns the Householder coefficients.
ykuroda 0:13a5d365ba16 165 *
ykuroda 0:13a5d365ba16 166 * \returns a const reference to the vector of Householder coefficients
ykuroda 0:13a5d365ba16 167 *
ykuroda 0:13a5d365ba16 168 * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
ykuroda 0:13a5d365ba16 169 * or the member function compute(const MatrixType&) has been called
ykuroda 0:13a5d365ba16 170 * before to compute the Hessenberg decomposition of a matrix.
ykuroda 0:13a5d365ba16 171 *
ykuroda 0:13a5d365ba16 172 * The Householder coefficients allow the reconstruction of the matrix
ykuroda 0:13a5d365ba16 173 * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
ykuroda 0:13a5d365ba16 174 *
ykuroda 0:13a5d365ba16 175 * \sa packedMatrix(), \ref Householder_Module "Householder module"
ykuroda 0:13a5d365ba16 176 */
ykuroda 0:13a5d365ba16 177 const CoeffVectorType& householderCoefficients() const
ykuroda 0:13a5d365ba16 178 {
ykuroda 0:13a5d365ba16 179 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
ykuroda 0:13a5d365ba16 180 return m_hCoeffs;
ykuroda 0:13a5d365ba16 181 }
ykuroda 0:13a5d365ba16 182
ykuroda 0:13a5d365ba16 183 /** \brief Returns the internal representation of the decomposition
ykuroda 0:13a5d365ba16 184 *
ykuroda 0:13a5d365ba16 185 * \returns a const reference to a matrix with the internal representation
ykuroda 0:13a5d365ba16 186 * of the decomposition.
ykuroda 0:13a5d365ba16 187 *
ykuroda 0:13a5d365ba16 188 * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
ykuroda 0:13a5d365ba16 189 * or the member function compute(const MatrixType&) has been called
ykuroda 0:13a5d365ba16 190 * before to compute the Hessenberg decomposition of a matrix.
ykuroda 0:13a5d365ba16 191 *
ykuroda 0:13a5d365ba16 192 * The returned matrix contains the following information:
ykuroda 0:13a5d365ba16 193 * - the upper part and lower sub-diagonal represent the Hessenberg matrix H
ykuroda 0:13a5d365ba16 194 * - the rest of the lower part contains the Householder vectors that, combined with
ykuroda 0:13a5d365ba16 195 * Householder coefficients returned by householderCoefficients(),
ykuroda 0:13a5d365ba16 196 * allows to reconstruct the matrix Q as
ykuroda 0:13a5d365ba16 197 * \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
ykuroda 0:13a5d365ba16 198 * Here, the matrices \f$ H_i \f$ are the Householder transformations
ykuroda 0:13a5d365ba16 199 * \f$ H_i = (I - h_i v_i v_i^T) \f$
ykuroda 0:13a5d365ba16 200 * where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
ykuroda 0:13a5d365ba16 201 * \f$ v_i \f$ is the Householder vector defined by
ykuroda 0:13a5d365ba16 202 * \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
ykuroda 0:13a5d365ba16 203 * with M the matrix returned by this function.
ykuroda 0:13a5d365ba16 204 *
ykuroda 0:13a5d365ba16 205 * See LAPACK for further details on this packed storage.
ykuroda 0:13a5d365ba16 206 *
ykuroda 0:13a5d365ba16 207 * Example: \include HessenbergDecomposition_packedMatrix.cpp
ykuroda 0:13a5d365ba16 208 * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
ykuroda 0:13a5d365ba16 209 *
ykuroda 0:13a5d365ba16 210 * \sa householderCoefficients()
ykuroda 0:13a5d365ba16 211 */
ykuroda 0:13a5d365ba16 212 const MatrixType& packedMatrix() const
ykuroda 0:13a5d365ba16 213 {
ykuroda 0:13a5d365ba16 214 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
ykuroda 0:13a5d365ba16 215 return m_matrix;
ykuroda 0:13a5d365ba16 216 }
ykuroda 0:13a5d365ba16 217
ykuroda 0:13a5d365ba16 218 /** \brief Reconstructs the orthogonal matrix Q in the decomposition
ykuroda 0:13a5d365ba16 219 *
ykuroda 0:13a5d365ba16 220 * \returns object representing the matrix Q
ykuroda 0:13a5d365ba16 221 *
ykuroda 0:13a5d365ba16 222 * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
ykuroda 0:13a5d365ba16 223 * or the member function compute(const MatrixType&) has been called
ykuroda 0:13a5d365ba16 224 * before to compute the Hessenberg decomposition of a matrix.
ykuroda 0:13a5d365ba16 225 *
ykuroda 0:13a5d365ba16 226 * This function returns a light-weight object of template class
ykuroda 0:13a5d365ba16 227 * HouseholderSequence. You can either apply it directly to a matrix or
ykuroda 0:13a5d365ba16 228 * you can convert it to a matrix of type #MatrixType.
ykuroda 0:13a5d365ba16 229 *
ykuroda 0:13a5d365ba16 230 * \sa matrixH() for an example, class HouseholderSequence
ykuroda 0:13a5d365ba16 231 */
ykuroda 0:13a5d365ba16 232 HouseholderSequenceType matrixQ() const
ykuroda 0:13a5d365ba16 233 {
ykuroda 0:13a5d365ba16 234 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
ykuroda 0:13a5d365ba16 235 return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
ykuroda 0:13a5d365ba16 236 .setLength(m_matrix.rows() - 1)
ykuroda 0:13a5d365ba16 237 .setShift(1);
ykuroda 0:13a5d365ba16 238 }
ykuroda 0:13a5d365ba16 239
ykuroda 0:13a5d365ba16 240 /** \brief Constructs the Hessenberg matrix H in the decomposition
ykuroda 0:13a5d365ba16 241 *
ykuroda 0:13a5d365ba16 242 * \returns expression object representing the matrix H
ykuroda 0:13a5d365ba16 243 *
ykuroda 0:13a5d365ba16 244 * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
ykuroda 0:13a5d365ba16 245 * or the member function compute(const MatrixType&) has been called
ykuroda 0:13a5d365ba16 246 * before to compute the Hessenberg decomposition of a matrix.
ykuroda 0:13a5d365ba16 247 *
ykuroda 0:13a5d365ba16 248 * The object returned by this function constructs the Hessenberg matrix H
ykuroda 0:13a5d365ba16 249 * when it is assigned to a matrix or otherwise evaluated. The matrix H is
ykuroda 0:13a5d365ba16 250 * constructed from the packed matrix as returned by packedMatrix(): The
ykuroda 0:13a5d365ba16 251 * upper part (including the subdiagonal) of the packed matrix contains
ykuroda 0:13a5d365ba16 252 * the matrix H. It may sometimes be better to directly use the packed
ykuroda 0:13a5d365ba16 253 * matrix instead of constructing the matrix H.
ykuroda 0:13a5d365ba16 254 *
ykuroda 0:13a5d365ba16 255 * Example: \include HessenbergDecomposition_matrixH.cpp
ykuroda 0:13a5d365ba16 256 * Output: \verbinclude HessenbergDecomposition_matrixH.out
ykuroda 0:13a5d365ba16 257 *
ykuroda 0:13a5d365ba16 258 * \sa matrixQ(), packedMatrix()
ykuroda 0:13a5d365ba16 259 */
ykuroda 0:13a5d365ba16 260 MatrixHReturnType matrixH() const
ykuroda 0:13a5d365ba16 261 {
ykuroda 0:13a5d365ba16 262 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
ykuroda 0:13a5d365ba16 263 return MatrixHReturnType(*this);
ykuroda 0:13a5d365ba16 264 }
ykuroda 0:13a5d365ba16 265
ykuroda 0:13a5d365ba16 266 private:
ykuroda 0:13a5d365ba16 267
ykuroda 0:13a5d365ba16 268 typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
ykuroda 0:13a5d365ba16 269 typedef typename NumTraits<Scalar>::Real RealScalar;
ykuroda 0:13a5d365ba16 270 static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
ykuroda 0:13a5d365ba16 271
ykuroda 0:13a5d365ba16 272 protected:
ykuroda 0:13a5d365ba16 273 MatrixType m_matrix;
ykuroda 0:13a5d365ba16 274 CoeffVectorType m_hCoeffs;
ykuroda 0:13a5d365ba16 275 VectorType m_temp;
ykuroda 0:13a5d365ba16 276 bool m_isInitialized;
ykuroda 0:13a5d365ba16 277 };
ykuroda 0:13a5d365ba16 278
ykuroda 0:13a5d365ba16 279 /** \internal
ykuroda 0:13a5d365ba16 280 * Performs a tridiagonal decomposition of \a matA in place.
ykuroda 0:13a5d365ba16 281 *
ykuroda 0:13a5d365ba16 282 * \param matA the input selfadjoint matrix
ykuroda 0:13a5d365ba16 283 * \param hCoeffs returned Householder coefficients
ykuroda 0:13a5d365ba16 284 *
ykuroda 0:13a5d365ba16 285 * The result is written in the lower triangular part of \a matA.
ykuroda 0:13a5d365ba16 286 *
ykuroda 0:13a5d365ba16 287 * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1.
ykuroda 0:13a5d365ba16 288 *
ykuroda 0:13a5d365ba16 289 * \sa packedMatrix()
ykuroda 0:13a5d365ba16 290 */
ykuroda 0:13a5d365ba16 291 template<typename MatrixType>
ykuroda 0:13a5d365ba16 292 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
ykuroda 0:13a5d365ba16 293 {
ykuroda 0:13a5d365ba16 294 eigen_assert(matA.rows()==matA.cols());
ykuroda 0:13a5d365ba16 295 Index n = matA.rows();
ykuroda 0:13a5d365ba16 296 temp.resize(n);
ykuroda 0:13a5d365ba16 297 for (Index i = 0; i<n-1; ++i)
ykuroda 0:13a5d365ba16 298 {
ykuroda 0:13a5d365ba16 299 // let's consider the vector v = i-th column starting at position i+1
ykuroda 0:13a5d365ba16 300 Index remainingSize = n-i-1;
ykuroda 0:13a5d365ba16 301 RealScalar beta;
ykuroda 0:13a5d365ba16 302 Scalar h;
ykuroda 0:13a5d365ba16 303 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
ykuroda 0:13a5d365ba16 304 matA.col(i).coeffRef(i+1) = beta;
ykuroda 0:13a5d365ba16 305 hCoeffs.coeffRef(i) = h;
ykuroda 0:13a5d365ba16 306
ykuroda 0:13a5d365ba16 307 // Apply similarity transformation to remaining columns,
ykuroda 0:13a5d365ba16 308 // i.e., compute A = H A H'
ykuroda 0:13a5d365ba16 309
ykuroda 0:13a5d365ba16 310 // A = H A
ykuroda 0:13a5d365ba16 311 matA.bottomRightCorner(remainingSize, remainingSize)
ykuroda 0:13a5d365ba16 312 .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
ykuroda 0:13a5d365ba16 313
ykuroda 0:13a5d365ba16 314 // A = A H'
ykuroda 0:13a5d365ba16 315 matA.rightCols(remainingSize)
ykuroda 0:13a5d365ba16 316 .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0));
ykuroda 0:13a5d365ba16 317 }
ykuroda 0:13a5d365ba16 318 }
ykuroda 0:13a5d365ba16 319
ykuroda 0:13a5d365ba16 320 namespace internal {
ykuroda 0:13a5d365ba16 321
ykuroda 0:13a5d365ba16 322 /** \eigenvalues_module \ingroup Eigenvalues_Module
ykuroda 0:13a5d365ba16 323 *
ykuroda 0:13a5d365ba16 324 *
ykuroda 0:13a5d365ba16 325 * \brief Expression type for return value of HessenbergDecomposition::matrixH()
ykuroda 0:13a5d365ba16 326 *
ykuroda 0:13a5d365ba16 327 * \tparam MatrixType type of matrix in the Hessenberg decomposition
ykuroda 0:13a5d365ba16 328 *
ykuroda 0:13a5d365ba16 329 * Objects of this type represent the Hessenberg matrix in the Hessenberg
ykuroda 0:13a5d365ba16 330 * decomposition of some matrix. The object holds a reference to the
ykuroda 0:13a5d365ba16 331 * HessenbergDecomposition class until the it is assigned or evaluated for
ykuroda 0:13a5d365ba16 332 * some other reason (the reference should remain valid during the life time
ykuroda 0:13a5d365ba16 333 * of this object). This class is the return type of
ykuroda 0:13a5d365ba16 334 * HessenbergDecomposition::matrixH(); there is probably no other use for this
ykuroda 0:13a5d365ba16 335 * class.
ykuroda 0:13a5d365ba16 336 */
ykuroda 0:13a5d365ba16 337 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
ykuroda 0:13a5d365ba16 338 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
ykuroda 0:13a5d365ba16 339 {
ykuroda 0:13a5d365ba16 340 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 341 public:
ykuroda 0:13a5d365ba16 342 /** \brief Constructor.
ykuroda 0:13a5d365ba16 343 *
ykuroda 0:13a5d365ba16 344 * \param[in] hess Hessenberg decomposition
ykuroda 0:13a5d365ba16 345 */
ykuroda 0:13a5d365ba16 346 HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
ykuroda 0:13a5d365ba16 347
ykuroda 0:13a5d365ba16 348 /** \brief Hessenberg matrix in decomposition.
ykuroda 0:13a5d365ba16 349 *
ykuroda 0:13a5d365ba16 350 * \param[out] result Hessenberg matrix in decomposition \p hess which
ykuroda 0:13a5d365ba16 351 * was passed to the constructor
ykuroda 0:13a5d365ba16 352 */
ykuroda 0:13a5d365ba16 353 template <typename ResultType>
ykuroda 0:13a5d365ba16 354 inline void evalTo(ResultType& result) const
ykuroda 0:13a5d365ba16 355 {
ykuroda 0:13a5d365ba16 356 result = m_hess.packedMatrix();
ykuroda 0:13a5d365ba16 357 Index n = result.rows();
ykuroda 0:13a5d365ba16 358 if (n>2)
ykuroda 0:13a5d365ba16 359 result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
ykuroda 0:13a5d365ba16 360 }
ykuroda 0:13a5d365ba16 361
ykuroda 0:13a5d365ba16 362 Index rows() const { return m_hess.packedMatrix().rows(); }
ykuroda 0:13a5d365ba16 363 Index cols() const { return m_hess.packedMatrix().cols(); }
ykuroda 0:13a5d365ba16 364
ykuroda 0:13a5d365ba16 365 protected:
ykuroda 0:13a5d365ba16 366 const HessenbergDecomposition<MatrixType>& m_hess;
ykuroda 0:13a5d365ba16 367 };
ykuroda 0:13a5d365ba16 368
ykuroda 0:13a5d365ba16 369 } // end namespace internal
ykuroda 0:13a5d365ba16 370
ykuroda 0:13a5d365ba16 371 } // end namespace Eigen
ykuroda 0:13a5d365ba16 372
ykuroda 0:13a5d365ba16 373 #endif // EIGEN_HESSENBERGDECOMPOSITION_H