1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Claire Maurice 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_COMPLEX_SCHUR_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_COMPLEX_SCHUR_H 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "./HessenbergDecomposition.h" 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg; 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \eigenvalues_module \ingroup Eigenvalues_Module 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class ComplexSchur 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Performs a complex Schur decomposition of a real or complex square matrix 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the matrix of which we are 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computing the Schur decomposition; this is expected to be an 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * instantiation of the Matrix class template. 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Given a real or complex square matrix A, this class computes the 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * complex matrix, and T is a complex upper triangular matrix. The 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * diagonal of the matrix T corresponds to the eigenvalues of the 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix A. 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Call the function compute() to compute the Schur decomposition of 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * a given matrix. Alternatively, you can use the 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexSchur(const MatrixType&, bool) constructor which computes 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the Schur decomposition at construction time. Once the 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition is computed, you can use the matrixU() and matrixT() 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * functions to retrieve the matrices U and V in the decomposition. 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \note This code is inspired from Jampack 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> class ComplexSchur 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options, 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Scalar type for matrices of type \p _MatrixType. */ 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Complex scalar type for \p _MatrixType. 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is \c std::complex<Scalar> if #Scalar is real (e.g., 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c float or \c double) and just \c Scalar if #Scalar is 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * complex. 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef std::complex<RealScalar> ComplexScalar; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Type for the matrices in the Schur decomposition. 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is a square matrix with entries of type #ComplexScalar. 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The size is the same as the size of \p _MatrixType. 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType; 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default constructor. 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed. 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * intends to perform decompositions via compute(). The \p size 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * parameter is only used as a hint. It is not an error to give a 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * wrong \p size, but it may impair performance. 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() for an example. 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matT(size,size), 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matU(size,size), 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_hess(size), 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false), 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matUisUptodate(false), 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_maxIters(-1) 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor; computes Schur decomposition of given matrix. 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor calls compute() to compute the Schur decomposition. 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa matrixT() and matrixU() for examples. 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexSchur(const MatrixType& matrix, bool computeU = true) 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_matT(matrix.rows(),matrix.cols()), 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matU(matrix.rows(),matrix.cols()), 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_hess(matrix.rows()), 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized(false), 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matUisUptodate(false), 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_maxIters(-1) 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix, computeU); 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the unitary matrix in the Schur decomposition. 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the matrix U. 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * It is assumed that either the constructor 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexSchur(const MatrixType& matrix, bool computeU) or the 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * member function compute(const MatrixType& matrix, bool computeU) 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * has been called before to compute the Schur decomposition of a 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix, and that \p computeU was set to true (the default 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * value). 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexSchur_matrixU.cpp 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexSchur_matrixU.out 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const ComplexMatrixType& matrixU() const 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition."); 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matU; 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the triangular matrix in the Schur decomposition. 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the matrix T. 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * It is assumed that either the constructor 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexSchur(const MatrixType& matrix, bool computeU) or the 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * member function compute(const MatrixType& matrix, bool computeU) 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * has been called before to compute the Schur decomposition of a 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix. 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Note that this function returns a plain square matrix. If you want to reference 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * only the upper triangular part, use: 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \code schur.matrixT().triangularView<Upper>() \endcode 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexSchur_matrixT.cpp 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexSchur_matrixT.out 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const ComplexMatrixType& matrixT() const 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matT; 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes Schur decomposition of given matrix. 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The Schur decomposition is computed by first reducing the 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix to Hessenberg form using the class 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * HessenbergDecomposition. The Hessenberg matrix is then reduced 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to triangular form by performing QR iterations with a single 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * shift. The cost of computing the Schur decomposition depends 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * on the number of iterations; as a rough guide, it may be taken 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * on the number of iterations; as a rough guide, it may be taken 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * if \a computeU is false. 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexSchur_compute.cpp 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexSchur_compute.out 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute(const MatrixType&, bool, Index) 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexSchur& compute(const MatrixType& matrix, bool computeU = true); 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Compute Schur decomposition from a given Hessenberg matrix 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] matrixH Matrix in Hessenberg form H 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param computeU Computes the matriX U of the Schur vectors 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \return Reference to \c *this 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This routine assumes that the matrix is already reduced in Hessenberg form matrixH 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * using either the class HessenbergDecomposition or another mean. 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * It computes the upper quasi-triangular matrix T of the Schur decomposition of H 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * When computeU is true, this routine computes the matrix U such that 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * is not available, the user should give an identity matrix (Q.setIdentity()) 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute(const MatrixType&, bool) 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename HessMatrixType, typename OrthMatrixType> 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true); 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Sets the maximum number of iterations allowed. 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * of the matrix. 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComplexSchur& setMaxIterations(Index maxIters) 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_maxIters = maxIters; 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns the maximum number of iterations. */ 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index getMaxIterations() 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_maxIters; 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Maximum number of iterations per row. 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * If not otherwise specified, the maximum number of iterations is this number times the size of the 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * matrix. It is currently set to 30. 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const int m_maxIterationsPerRow = 30; 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexMatrixType m_matT, m_matU; 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HessenbergDecomposition<MatrixType> m_hess; 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo m_info; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_matUisUptodate; 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m_maxIters; 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool subdiagonalEntryIsNeglegible(Index i); 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar computeShift(Index iu, Index iter); 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void reduceToTriangularForm(bool computeU); 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** If m_matT(i+1,i) is neglegible in floating point arithmetic 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * return true, else return false. */ 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i) 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1)); 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar sd = numext::norm1(m_matT.coeff(i+1,i)); 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon())) 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,i) = ComplexScalar(0); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Compute the shift in the current QR iteration. */ 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter) 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iter == 10 || iter == 20) 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2))); 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute the shift as one of the eigenvalues of t, the 2x2 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // diagonal block on the bottom of the active submatrix 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar normt = t.cwiseAbs().sum(); 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t /= normt; // the normalization by sf is to avoid under/overflow 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar b = t.coeff(0,1) * t.coeff(1,0); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar c = t.coeff(0,0) - t.coeff(1,1); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar disc = sqrt(c*c + RealScalar(4)*b); 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b; 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar eival1 = (trace + disc) / RealScalar(2); 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar eival2 = (trace - disc) / RealScalar(2); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(numext::norm1(eival1) > numext::norm1(eival2)) 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eival2 = det / eival1; 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eival1 = det / eival2; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // choose the eigenvalue closest to the bottom entry of the diagonal 3087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1))) 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return normt * eival1; 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return normt * eival2; 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matUisUptodate = false; 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(matrix.cols() == matrix.rows()); 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.cols() == 1) 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT = matrix.template cast<ComplexScalar>(); 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) m_matU = ComplexMatrixType::Identity(1,1); 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matUisUptodate = computeU; 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU); 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez computeFromHessenberg(m_matT, m_matU, computeU); 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 3377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename HessMatrixType, typename OrthMatrixType> 3387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) 3397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT = matrixH; 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(computeU) 3427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matU = matrixQ; 3437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez reduceToTriangularForm(computeU); 3447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 3457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/* Reduce given matrix to Hessenberg form */ 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, bool IsComplex> 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct complex_schur_reduce_to_hessenberg 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // this is the implementation for the case IsComplex = true 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_hess.compute(matrix); 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_matT = _this.m_hess.matrixH(); 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) _this.m_matU = _this.m_hess.matrixQ(); 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct complex_schur_reduce_to_hessenberg<MatrixType, false> 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar; 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_hess.compute(matrix); 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>(); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // This may cause an allocation which seems to be avoidable 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType Q = _this.m_hess.matrixQ(); 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_matU = Q.template cast<ComplexScalar>(); 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU) 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 3867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index maxIters = m_maxIters; 3877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (maxIters == -1) 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez maxIters = m_maxIterationsPerRow * m_matT.rows(); 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // The matrix m_matT is divided in three parts. 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Rows il,...,iu is the part we are working on (the active submatrix). 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Rows iu+1,...,end are already brought in triangular form. 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index iu = m_matT.cols() - 1; 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index il; 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index iter = 0; // number of iterations we are working on the (iu,iu) element 3977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index totalIter = 0; // number of iterations for whole matrix 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(true) 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // find iu, the bottom row of the active submatrix 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(iu > 0) 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!subdiagonalEntryIsNeglegible(iu-1)) break; 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iter = 0; 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath --iu; 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // if iu is zero then we are done; the whole matrix is triangularized 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(iu==0) break; 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // if we spent too many iterations, we give up 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iter++; 4147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez totalIter++; 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(totalIter > maxIters) break; 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // find il, the top row of the active submatrix 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath il = iu-1; 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(il > 0 && !subdiagonalEntryIsNeglegible(il-1)) 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath --il; 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* perform the QR step using Givens rotations. The first rotation 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath creates a bulge; the (il+2,il) element becomes nonzero. This 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bulge is chased down to the bottom of the active submatrix. */ 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar shift = computeShift(iu, iter); 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<ComplexScalar> rot; 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il)); 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint()); 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot); 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) m_matU.applyOnTheRight(il, il+1, rot); 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=il+1 ; i<iu ; i++) 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1)); 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,i-1) = ComplexScalar(0); 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint()); 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot); 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) m_matU.applyOnTheRight(i, i+1, rot); 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(totalIter <= maxIters) 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NoConvergence; 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matUisUptodate = computeU; 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_COMPLEX_SCHUR_H 457