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; 662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 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 */ 942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit 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 */ 1122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 1132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true) 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_matT(matrix.rows(),matrix.cols()), 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matU(matrix.rows(),matrix.cols()), 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_hess(matrix.rows()), 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized(false), 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matUisUptodate(false), 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_maxIters(-1) 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang compute(matrix.derived(), computeU); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the unitary matrix in the Schur decomposition. 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the matrix U. 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * It is assumed that either the constructor 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexSchur(const MatrixType& matrix, bool computeU) or the 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * member function compute(const MatrixType& matrix, bool computeU) 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * has been called before to compute the Schur decomposition of a 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix, and that \p computeU was set to true (the default 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * value). 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexSchur_matrixU.cpp 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexSchur_matrixU.out 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const ComplexMatrixType& matrixU() const 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition."); 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matU; 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the triangular matrix in the Schur decomposition. 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the matrix T. 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * It is assumed that either the constructor 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexSchur(const MatrixType& matrix, bool computeU) or the 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * member function compute(const MatrixType& matrix, bool computeU) 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * has been called before to compute the Schur decomposition of a 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix. 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Note that this function returns a plain square matrix. If you want to reference 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * only the upper triangular part, use: 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \code schur.matrixT().triangularView<Upper>() \endcode 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexSchur_matrixT.cpp 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexSchur_matrixT.out 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const ComplexMatrixType& matrixT() const 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matT; 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes Schur decomposition of given matrix. 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose Schur decomposition is to be computed. 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeU If true, both T and U are computed; if false, only T is computed. 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The Schur decomposition is computed by first reducing the 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix to Hessenberg form using the class 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * HessenbergDecomposition. The Hessenberg matrix is then reduced 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to triangular form by performing QR iterations with a single 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * shift. The cost of computing the Schur decomposition depends 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * on the number of iterations; as a rough guide, it may be taken 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * on the number of iterations; as a rough guide, it may be taken 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * if \a computeU is false. 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexSchur_compute.cpp 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexSchur_compute.out 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute(const MatrixType&, bool, Index) 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 1912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true); 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Compute Schur decomposition from a given Hessenberg matrix 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] matrixH Matrix in Hessenberg form H 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param computeU Computes the matriX U of the Schur vectors 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \return Reference to \c *this 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This routine assumes that the matrix is already reduced in Hessenberg form matrixH 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * using either the class HessenbergDecomposition or another mean. 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * It computes the upper quasi-triangular matrix T of the Schur decomposition of H 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * When computeU is true, this routine computes the matrix U such that 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * is not available, the user should give an identity matrix (Q.setIdentity()) 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute(const MatrixType&, bool) 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename HessMatrixType, typename OrthMatrixType> 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true); 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Sets the maximum number of iterations allowed. 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * of the matrix. 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComplexSchur& setMaxIterations(Index maxIters) 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_maxIters = maxIters; 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns the maximum number of iterations. */ 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index getMaxIterations() 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_maxIters; 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Maximum number of iterations per row. 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * If not otherwise specified, the maximum number of iterations is this number times the size of the 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * matrix. It is currently set to 30. 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static const int m_maxIterationsPerRow = 30; 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexMatrixType m_matT, m_matU; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HessenbergDecomposition<MatrixType> m_hess; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo m_info; 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_matUisUptodate; 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m_maxIters; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool subdiagonalEntryIsNeglegible(Index i); 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar computeShift(Index iu, Index iter); 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void reduceToTriangularForm(bool computeU); 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** If m_matT(i+1,i) is neglegible in floating point arithmetic 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * return true, else return false. */ 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i) 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1)); 2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar sd = numext::norm1(m_matT.coeff(i+1,i)); 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon())) 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,i) = ComplexScalar(0); 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return false; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Compute the shift in the current QR iteration. */ 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter) 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iter == 10 || iter == 20) 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2))); 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // compute the shift as one of the eigenvalues of t, the 2x2 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // diagonal block on the bottom of the active submatrix 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar normt = t.cwiseAbs().sum(); 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath t /= normt; // the normalization by sf is to avoid under/overflow 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar b = t.coeff(0,1) * t.coeff(1,0); 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar c = t.coeff(0,0) - t.coeff(1,1); 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar disc = sqrt(c*c + RealScalar(4)*b); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b; 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar eival1 = (trace + disc) / RealScalar(2); 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar eival2 = (trace - disc) / RealScalar(2); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(numext::norm1(eival1) > numext::norm1(eival2)) 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eival2 = det / eival1; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eival1 = det / eival2; 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // choose the eigenvalue closest to the bottom entry of the diagonal 3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1))) 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return normt * eival1; 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return normt * eival2; 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 3182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename InputType> 3192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matUisUptodate = false; 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(matrix.cols() == matrix.rows()); 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.cols() == 1) 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matT = matrix.derived().template cast<ComplexScalar>(); 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) m_matU = ComplexMatrixType::Identity(1,1); 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matUisUptodate = computeU; 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU); 3357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez computeFromHessenberg(m_matT, m_matU, computeU); 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename HessMatrixType, typename OrthMatrixType> 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) 3427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 3437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT = matrixH; 3447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(computeU) 3457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matU = matrixQ; 3467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez reduceToTriangularForm(computeU); 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 3487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/* Reduce given matrix to Hessenberg form */ 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, bool IsComplex> 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct complex_schur_reduce_to_hessenberg 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // this is the implementation for the case IsComplex = true 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_hess.compute(matrix); 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_matT = _this.m_hess.matrixH(); 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) _this.m_matU = _this.m_hess.matrixQ(); 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct complex_schur_reduce_to_hessenberg<MatrixType, false> 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar; 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_hess.compute(matrix); 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>(); 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // This may cause an allocation which seems to be avoidable 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType Q = _this.m_hess.matrixQ(); 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _this.m_matU = Q.template cast<ComplexScalar>(); 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU) 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index maxIters = m_maxIters; 3907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (maxIters == -1) 3917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez maxIters = m_maxIterationsPerRow * m_matT.rows(); 3927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // The matrix m_matT is divided in three parts. 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Rows il,...,iu is the part we are working on (the active submatrix). 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Rows iu+1,...,end are already brought in triangular form. 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index iu = m_matT.cols() - 1; 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index il; 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index iter = 0; // number of iterations we are working on the (iu,iu) element 4007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index totalIter = 0; // number of iterations for whole matrix 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(true) 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // find iu, the bottom row of the active submatrix 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(iu > 0) 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(!subdiagonalEntryIsNeglegible(iu-1)) break; 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iter = 0; 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath --iu; 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // if iu is zero then we are done; the whole matrix is triangularized 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(iu==0) break; 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // if we spent too many iterations, we give up 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath iter++; 4177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez totalIter++; 4187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(totalIter > maxIters) break; 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // find il, the top row of the active submatrix 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath il = iu-1; 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while(il > 0 && !subdiagonalEntryIsNeglegible(il-1)) 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath --il; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* perform the QR step using Givens rotations. The first rotation 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath creates a bulge; the (il+2,il) element becomes nonzero. This 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bulge is chased down to the bottom of the active submatrix. */ 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar shift = computeShift(iu, iter); 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<ComplexScalar> rot; 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il)); 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint()); 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot); 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) m_matU.applyOnTheRight(il, il+1, rot); 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=il+1 ; i<iu ; i++) 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1)); 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,i-1) = ComplexScalar(0); 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint()); 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot); 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeU) m_matU.applyOnTheRight(i, i+1, rot); 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(totalIter <= maxIters) 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NoConvergence; 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matUisUptodate = computeU; 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_COMPLEX_SCHUR_H 460