Eigne Matrix Class Library
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
Eigen Matrix Class Library for mbed.
Finally, you can use Eigen on your mbed!!!
src/Eigenvalues/ComplexSchur.h@0:13a5d365ba16, 2016-10-13 (annotated)
- Committer:
- ykuroda
- Date:
- Thu Oct 13 04:07:23 2016 +0000
- Revision:
- 0:13a5d365ba16
First commint, Eigne Matrix Class Library
Who changed what in which revision?
User | Revision | Line number | New 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) 2009 Claire Maurice |
ykuroda | 0:13a5d365ba16 | 5 | // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> |
ykuroda | 0:13a5d365ba16 | 6 | // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> |
ykuroda | 0:13a5d365ba16 | 7 | // |
ykuroda | 0:13a5d365ba16 | 8 | // This Source Code Form is subject to the terms of the Mozilla |
ykuroda | 0:13a5d365ba16 | 9 | // Public License v. 2.0. If a copy of the MPL was not distributed |
ykuroda | 0:13a5d365ba16 | 10 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
ykuroda | 0:13a5d365ba16 | 11 | |
ykuroda | 0:13a5d365ba16 | 12 | #ifndef EIGEN_COMPLEX_SCHUR_H |
ykuroda | 0:13a5d365ba16 | 13 | #define EIGEN_COMPLEX_SCHUR_H |
ykuroda | 0:13a5d365ba16 | 14 | |
ykuroda | 0:13a5d365ba16 | 15 | #include "./HessenbergDecomposition.h" |
ykuroda | 0:13a5d365ba16 | 16 | |
ykuroda | 0:13a5d365ba16 | 17 | namespace Eigen { |
ykuroda | 0:13a5d365ba16 | 18 | |
ykuroda | 0:13a5d365ba16 | 19 | namespace internal { |
ykuroda | 0:13a5d365ba16 | 20 | template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg; |
ykuroda | 0:13a5d365ba16 | 21 | } |
ykuroda | 0:13a5d365ba16 | 22 | |
ykuroda | 0:13a5d365ba16 | 23 | /** \eigenvalues_module \ingroup Eigenvalues_Module |
ykuroda | 0:13a5d365ba16 | 24 | * |
ykuroda | 0:13a5d365ba16 | 25 | * |
ykuroda | 0:13a5d365ba16 | 26 | * \class ComplexSchur |
ykuroda | 0:13a5d365ba16 | 27 | * |
ykuroda | 0:13a5d365ba16 | 28 | * \brief Performs a complex Schur decomposition of a real or complex square matrix |
ykuroda | 0:13a5d365ba16 | 29 | * |
ykuroda | 0:13a5d365ba16 | 30 | * \tparam _MatrixType the type of the matrix of which we are |
ykuroda | 0:13a5d365ba16 | 31 | * computing the Schur decomposition; this is expected to be an |
ykuroda | 0:13a5d365ba16 | 32 | * instantiation of the Matrix class template. |
ykuroda | 0:13a5d365ba16 | 33 | * |
ykuroda | 0:13a5d365ba16 | 34 | * Given a real or complex square matrix A, this class computes the |
ykuroda | 0:13a5d365ba16 | 35 | * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary |
ykuroda | 0:13a5d365ba16 | 36 | * complex matrix, and T is a complex upper triangular matrix. The |
ykuroda | 0:13a5d365ba16 | 37 | * diagonal of the matrix T corresponds to the eigenvalues of the |
ykuroda | 0:13a5d365ba16 | 38 | * matrix A. |
ykuroda | 0:13a5d365ba16 | 39 | * |
ykuroda | 0:13a5d365ba16 | 40 | * Call the function compute() to compute the Schur decomposition of |
ykuroda | 0:13a5d365ba16 | 41 | * a given matrix. Alternatively, you can use the |
ykuroda | 0:13a5d365ba16 | 42 | * ComplexSchur(const MatrixType&, bool) constructor which computes |
ykuroda | 0:13a5d365ba16 | 43 | * the Schur decomposition at construction time. Once the |
ykuroda | 0:13a5d365ba16 | 44 | * decomposition is computed, you can use the matrixU() and matrixT() |
ykuroda | 0:13a5d365ba16 | 45 | * functions to retrieve the matrices U and V in the decomposition. |
ykuroda | 0:13a5d365ba16 | 46 | * |
ykuroda | 0:13a5d365ba16 | 47 | * \note This code is inspired from Jampack |
ykuroda | 0:13a5d365ba16 | 48 | * |
ykuroda | 0:13a5d365ba16 | 49 | * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver |
ykuroda | 0:13a5d365ba16 | 50 | */ |
ykuroda | 0:13a5d365ba16 | 51 | template<typename _MatrixType> class ComplexSchur |
ykuroda | 0:13a5d365ba16 | 52 | { |
ykuroda | 0:13a5d365ba16 | 53 | public: |
ykuroda | 0:13a5d365ba16 | 54 | typedef _MatrixType MatrixType; |
ykuroda | 0:13a5d365ba16 | 55 | enum { |
ykuroda | 0:13a5d365ba16 | 56 | RowsAtCompileTime = MatrixType::RowsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 57 | ColsAtCompileTime = MatrixType::ColsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 58 | Options = MatrixType::Options, |
ykuroda | 0:13a5d365ba16 | 59 | MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, |
ykuroda | 0:13a5d365ba16 | 60 | MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime |
ykuroda | 0:13a5d365ba16 | 61 | }; |
ykuroda | 0:13a5d365ba16 | 62 | |
ykuroda | 0:13a5d365ba16 | 63 | /** \brief Scalar type for matrices of type \p _MatrixType. */ |
ykuroda | 0:13a5d365ba16 | 64 | typedef typename MatrixType::Scalar Scalar; |
ykuroda | 0:13a5d365ba16 | 65 | typedef typename NumTraits<Scalar>::Real RealScalar; |
ykuroda | 0:13a5d365ba16 | 66 | typedef typename MatrixType::Index Index; |
ykuroda | 0:13a5d365ba16 | 67 | |
ykuroda | 0:13a5d365ba16 | 68 | /** \brief Complex scalar type for \p _MatrixType. |
ykuroda | 0:13a5d365ba16 | 69 | * |
ykuroda | 0:13a5d365ba16 | 70 | * This is \c std::complex<Scalar> if #Scalar is real (e.g., |
ykuroda | 0:13a5d365ba16 | 71 | * \c float or \c double) and just \c Scalar if #Scalar is |
ykuroda | 0:13a5d365ba16 | 72 | * complex. |
ykuroda | 0:13a5d365ba16 | 73 | */ |
ykuroda | 0:13a5d365ba16 | 74 | typedef std::complex<RealScalar> ComplexScalar; |
ykuroda | 0:13a5d365ba16 | 75 | |
ykuroda | 0:13a5d365ba16 | 76 | /** \brief Type for the matrices in the Schur decomposition. |
ykuroda | 0:13a5d365ba16 | 77 | * |
ykuroda | 0:13a5d365ba16 | 78 | * This is a square matrix with entries of type #ComplexScalar. |
ykuroda | 0:13a5d365ba16 | 79 | * The size is the same as the size of \p _MatrixType. |
ykuroda | 0:13a5d365ba16 | 80 | */ |
ykuroda | 0:13a5d365ba16 | 81 | typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType; |
ykuroda | 0:13a5d365ba16 | 82 | |
ykuroda | 0:13a5d365ba16 | 83 | /** \brief Default constructor. |
ykuroda | 0:13a5d365ba16 | 84 | * |
ykuroda | 0:13a5d365ba16 | 85 | * \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed. |
ykuroda | 0:13a5d365ba16 | 86 | * |
ykuroda | 0:13a5d365ba16 | 87 | * The default constructor is useful in cases in which the user |
ykuroda | 0:13a5d365ba16 | 88 | * intends to perform decompositions via compute(). The \p size |
ykuroda | 0:13a5d365ba16 | 89 | * parameter is only used as a hint. It is not an error to give a |
ykuroda | 0:13a5d365ba16 | 90 | * wrong \p size, but it may impair performance. |
ykuroda | 0:13a5d365ba16 | 91 | * |
ykuroda | 0:13a5d365ba16 | 92 | * \sa compute() for an example. |
ykuroda | 0:13a5d365ba16 | 93 | */ |
ykuroda | 0:13a5d365ba16 | 94 | ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) |
ykuroda | 0:13a5d365ba16 | 95 | : m_matT(size,size), |
ykuroda | 0:13a5d365ba16 | 96 | m_matU(size,size), |
ykuroda | 0:13a5d365ba16 | 97 | m_hess(size), |
ykuroda | 0:13a5d365ba16 | 98 | m_isInitialized(false), |
ykuroda | 0:13a5d365ba16 | 99 | m_matUisUptodate(false), |
ykuroda | 0:13a5d365ba16 | 100 | m_maxIters(-1) |
ykuroda | 0:13a5d365ba16 | 101 | {} |
ykuroda | 0:13a5d365ba16 | 102 | |
ykuroda | 0:13a5d365ba16 | 103 | /** \brief Constructor; computes Schur decomposition of given matrix. |
ykuroda | 0:13a5d365ba16 | 104 | * |
ykuroda | 0:13a5d365ba16 | 105 | * \param[in] matrix Square matrix whose Schur decomposition is to be computed. |
ykuroda | 0:13a5d365ba16 | 106 | * \param[in] computeU If true, both T and U are computed; if false, only T is computed. |
ykuroda | 0:13a5d365ba16 | 107 | * |
ykuroda | 0:13a5d365ba16 | 108 | * This constructor calls compute() to compute the Schur decomposition. |
ykuroda | 0:13a5d365ba16 | 109 | * |
ykuroda | 0:13a5d365ba16 | 110 | * \sa matrixT() and matrixU() for examples. |
ykuroda | 0:13a5d365ba16 | 111 | */ |
ykuroda | 0:13a5d365ba16 | 112 | ComplexSchur(const MatrixType& matrix, bool computeU = true) |
ykuroda | 0:13a5d365ba16 | 113 | : m_matT(matrix.rows(),matrix.cols()), |
ykuroda | 0:13a5d365ba16 | 114 | m_matU(matrix.rows(),matrix.cols()), |
ykuroda | 0:13a5d365ba16 | 115 | m_hess(matrix.rows()), |
ykuroda | 0:13a5d365ba16 | 116 | m_isInitialized(false), |
ykuroda | 0:13a5d365ba16 | 117 | m_matUisUptodate(false), |
ykuroda | 0:13a5d365ba16 | 118 | m_maxIters(-1) |
ykuroda | 0:13a5d365ba16 | 119 | { |
ykuroda | 0:13a5d365ba16 | 120 | compute(matrix, computeU); |
ykuroda | 0:13a5d365ba16 | 121 | } |
ykuroda | 0:13a5d365ba16 | 122 | |
ykuroda | 0:13a5d365ba16 | 123 | /** \brief Returns the unitary matrix in the Schur decomposition. |
ykuroda | 0:13a5d365ba16 | 124 | * |
ykuroda | 0:13a5d365ba16 | 125 | * \returns A const reference to the matrix U. |
ykuroda | 0:13a5d365ba16 | 126 | * |
ykuroda | 0:13a5d365ba16 | 127 | * It is assumed that either the constructor |
ykuroda | 0:13a5d365ba16 | 128 | * ComplexSchur(const MatrixType& matrix, bool computeU) or the |
ykuroda | 0:13a5d365ba16 | 129 | * member function compute(const MatrixType& matrix, bool computeU) |
ykuroda | 0:13a5d365ba16 | 130 | * has been called before to compute the Schur decomposition of a |
ykuroda | 0:13a5d365ba16 | 131 | * matrix, and that \p computeU was set to true (the default |
ykuroda | 0:13a5d365ba16 | 132 | * value). |
ykuroda | 0:13a5d365ba16 | 133 | * |
ykuroda | 0:13a5d365ba16 | 134 | * Example: \include ComplexSchur_matrixU.cpp |
ykuroda | 0:13a5d365ba16 | 135 | * Output: \verbinclude ComplexSchur_matrixU.out |
ykuroda | 0:13a5d365ba16 | 136 | */ |
ykuroda | 0:13a5d365ba16 | 137 | const ComplexMatrixType& matrixU() const |
ykuroda | 0:13a5d365ba16 | 138 | { |
ykuroda | 0:13a5d365ba16 | 139 | eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); |
ykuroda | 0:13a5d365ba16 | 140 | eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition."); |
ykuroda | 0:13a5d365ba16 | 141 | return m_matU; |
ykuroda | 0:13a5d365ba16 | 142 | } |
ykuroda | 0:13a5d365ba16 | 143 | |
ykuroda | 0:13a5d365ba16 | 144 | /** \brief Returns the triangular matrix in the Schur decomposition. |
ykuroda | 0:13a5d365ba16 | 145 | * |
ykuroda | 0:13a5d365ba16 | 146 | * \returns A const reference to the matrix T. |
ykuroda | 0:13a5d365ba16 | 147 | * |
ykuroda | 0:13a5d365ba16 | 148 | * It is assumed that either the constructor |
ykuroda | 0:13a5d365ba16 | 149 | * ComplexSchur(const MatrixType& matrix, bool computeU) or the |
ykuroda | 0:13a5d365ba16 | 150 | * member function compute(const MatrixType& matrix, bool computeU) |
ykuroda | 0:13a5d365ba16 | 151 | * has been called before to compute the Schur decomposition of a |
ykuroda | 0:13a5d365ba16 | 152 | * matrix. |
ykuroda | 0:13a5d365ba16 | 153 | * |
ykuroda | 0:13a5d365ba16 | 154 | * Note that this function returns a plain square matrix. If you want to reference |
ykuroda | 0:13a5d365ba16 | 155 | * only the upper triangular part, use: |
ykuroda | 0:13a5d365ba16 | 156 | * \code schur.matrixT().triangularView<Upper>() \endcode |
ykuroda | 0:13a5d365ba16 | 157 | * |
ykuroda | 0:13a5d365ba16 | 158 | * Example: \include ComplexSchur_matrixT.cpp |
ykuroda | 0:13a5d365ba16 | 159 | * Output: \verbinclude ComplexSchur_matrixT.out |
ykuroda | 0:13a5d365ba16 | 160 | */ |
ykuroda | 0:13a5d365ba16 | 161 | const ComplexMatrixType& matrixT() const |
ykuroda | 0:13a5d365ba16 | 162 | { |
ykuroda | 0:13a5d365ba16 | 163 | eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); |
ykuroda | 0:13a5d365ba16 | 164 | return m_matT; |
ykuroda | 0:13a5d365ba16 | 165 | } |
ykuroda | 0:13a5d365ba16 | 166 | |
ykuroda | 0:13a5d365ba16 | 167 | /** \brief Computes Schur decomposition of given matrix. |
ykuroda | 0:13a5d365ba16 | 168 | * |
ykuroda | 0:13a5d365ba16 | 169 | * \param[in] matrix Square matrix whose Schur decomposition is to be computed. |
ykuroda | 0:13a5d365ba16 | 170 | * \param[in] computeU If true, both T and U are computed; if false, only T is computed. |
ykuroda | 0:13a5d365ba16 | 171 | |
ykuroda | 0:13a5d365ba16 | 172 | * \returns Reference to \c *this |
ykuroda | 0:13a5d365ba16 | 173 | * |
ykuroda | 0:13a5d365ba16 | 174 | * The Schur decomposition is computed by first reducing the |
ykuroda | 0:13a5d365ba16 | 175 | * matrix to Hessenberg form using the class |
ykuroda | 0:13a5d365ba16 | 176 | * HessenbergDecomposition. The Hessenberg matrix is then reduced |
ykuroda | 0:13a5d365ba16 | 177 | * to triangular form by performing QR iterations with a single |
ykuroda | 0:13a5d365ba16 | 178 | * shift. The cost of computing the Schur decomposition depends |
ykuroda | 0:13a5d365ba16 | 179 | * on the number of iterations; as a rough guide, it may be taken |
ykuroda | 0:13a5d365ba16 | 180 | * on the number of iterations; as a rough guide, it may be taken |
ykuroda | 0:13a5d365ba16 | 181 | * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops |
ykuroda | 0:13a5d365ba16 | 182 | * if \a computeU is false. |
ykuroda | 0:13a5d365ba16 | 183 | * |
ykuroda | 0:13a5d365ba16 | 184 | * Example: \include ComplexSchur_compute.cpp |
ykuroda | 0:13a5d365ba16 | 185 | * Output: \verbinclude ComplexSchur_compute.out |
ykuroda | 0:13a5d365ba16 | 186 | * |
ykuroda | 0:13a5d365ba16 | 187 | * \sa compute(const MatrixType&, bool, Index) |
ykuroda | 0:13a5d365ba16 | 188 | */ |
ykuroda | 0:13a5d365ba16 | 189 | ComplexSchur& compute(const MatrixType& matrix, bool computeU = true); |
ykuroda | 0:13a5d365ba16 | 190 | |
ykuroda | 0:13a5d365ba16 | 191 | /** \brief Compute Schur decomposition from a given Hessenberg matrix |
ykuroda | 0:13a5d365ba16 | 192 | * \param[in] matrixH Matrix in Hessenberg form H |
ykuroda | 0:13a5d365ba16 | 193 | * \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T |
ykuroda | 0:13a5d365ba16 | 194 | * \param computeU Computes the matriX U of the Schur vectors |
ykuroda | 0:13a5d365ba16 | 195 | * \return Reference to \c *this |
ykuroda | 0:13a5d365ba16 | 196 | * |
ykuroda | 0:13a5d365ba16 | 197 | * This routine assumes that the matrix is already reduced in Hessenberg form matrixH |
ykuroda | 0:13a5d365ba16 | 198 | * using either the class HessenbergDecomposition or another mean. |
ykuroda | 0:13a5d365ba16 | 199 | * It computes the upper quasi-triangular matrix T of the Schur decomposition of H |
ykuroda | 0:13a5d365ba16 | 200 | * When computeU is true, this routine computes the matrix U such that |
ykuroda | 0:13a5d365ba16 | 201 | * A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix |
ykuroda | 0:13a5d365ba16 | 202 | * |
ykuroda | 0:13a5d365ba16 | 203 | * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix |
ykuroda | 0:13a5d365ba16 | 204 | * is not available, the user should give an identity matrix (Q.setIdentity()) |
ykuroda | 0:13a5d365ba16 | 205 | * |
ykuroda | 0:13a5d365ba16 | 206 | * \sa compute(const MatrixType&, bool) |
ykuroda | 0:13a5d365ba16 | 207 | */ |
ykuroda | 0:13a5d365ba16 | 208 | template<typename HessMatrixType, typename OrthMatrixType> |
ykuroda | 0:13a5d365ba16 | 209 | ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true); |
ykuroda | 0:13a5d365ba16 | 210 | |
ykuroda | 0:13a5d365ba16 | 211 | /** \brief Reports whether previous computation was successful. |
ykuroda | 0:13a5d365ba16 | 212 | * |
ykuroda | 0:13a5d365ba16 | 213 | * \returns \c Success if computation was succesful, \c NoConvergence otherwise. |
ykuroda | 0:13a5d365ba16 | 214 | */ |
ykuroda | 0:13a5d365ba16 | 215 | ComputationInfo info() const |
ykuroda | 0:13a5d365ba16 | 216 | { |
ykuroda | 0:13a5d365ba16 | 217 | eigen_assert(m_isInitialized && "ComplexSchur is not initialized."); |
ykuroda | 0:13a5d365ba16 | 218 | return m_info; |
ykuroda | 0:13a5d365ba16 | 219 | } |
ykuroda | 0:13a5d365ba16 | 220 | |
ykuroda | 0:13a5d365ba16 | 221 | /** \brief Sets the maximum number of iterations allowed. |
ykuroda | 0:13a5d365ba16 | 222 | * |
ykuroda | 0:13a5d365ba16 | 223 | * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size |
ykuroda | 0:13a5d365ba16 | 224 | * of the matrix. |
ykuroda | 0:13a5d365ba16 | 225 | */ |
ykuroda | 0:13a5d365ba16 | 226 | ComplexSchur& setMaxIterations(Index maxIters) |
ykuroda | 0:13a5d365ba16 | 227 | { |
ykuroda | 0:13a5d365ba16 | 228 | m_maxIters = maxIters; |
ykuroda | 0:13a5d365ba16 | 229 | return *this; |
ykuroda | 0:13a5d365ba16 | 230 | } |
ykuroda | 0:13a5d365ba16 | 231 | |
ykuroda | 0:13a5d365ba16 | 232 | /** \brief Returns the maximum number of iterations. */ |
ykuroda | 0:13a5d365ba16 | 233 | Index getMaxIterations() |
ykuroda | 0:13a5d365ba16 | 234 | { |
ykuroda | 0:13a5d365ba16 | 235 | return m_maxIters; |
ykuroda | 0:13a5d365ba16 | 236 | } |
ykuroda | 0:13a5d365ba16 | 237 | |
ykuroda | 0:13a5d365ba16 | 238 | /** \brief Maximum number of iterations per row. |
ykuroda | 0:13a5d365ba16 | 239 | * |
ykuroda | 0:13a5d365ba16 | 240 | * If not otherwise specified, the maximum number of iterations is this number times the size of the |
ykuroda | 0:13a5d365ba16 | 241 | * matrix. It is currently set to 30. |
ykuroda | 0:13a5d365ba16 | 242 | */ |
ykuroda | 0:13a5d365ba16 | 243 | static const int m_maxIterationsPerRow = 30; |
ykuroda | 0:13a5d365ba16 | 244 | |
ykuroda | 0:13a5d365ba16 | 245 | protected: |
ykuroda | 0:13a5d365ba16 | 246 | ComplexMatrixType m_matT, m_matU; |
ykuroda | 0:13a5d365ba16 | 247 | HessenbergDecomposition<MatrixType> m_hess; |
ykuroda | 0:13a5d365ba16 | 248 | ComputationInfo m_info; |
ykuroda | 0:13a5d365ba16 | 249 | bool m_isInitialized; |
ykuroda | 0:13a5d365ba16 | 250 | bool m_matUisUptodate; |
ykuroda | 0:13a5d365ba16 | 251 | Index m_maxIters; |
ykuroda | 0:13a5d365ba16 | 252 | |
ykuroda | 0:13a5d365ba16 | 253 | private: |
ykuroda | 0:13a5d365ba16 | 254 | bool subdiagonalEntryIsNeglegible(Index i); |
ykuroda | 0:13a5d365ba16 | 255 | ComplexScalar computeShift(Index iu, Index iter); |
ykuroda | 0:13a5d365ba16 | 256 | void reduceToTriangularForm(bool computeU); |
ykuroda | 0:13a5d365ba16 | 257 | friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>; |
ykuroda | 0:13a5d365ba16 | 258 | }; |
ykuroda | 0:13a5d365ba16 | 259 | |
ykuroda | 0:13a5d365ba16 | 260 | /** If m_matT(i+1,i) is neglegible in floating point arithmetic |
ykuroda | 0:13a5d365ba16 | 261 | * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and |
ykuroda | 0:13a5d365ba16 | 262 | * return true, else return false. */ |
ykuroda | 0:13a5d365ba16 | 263 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 264 | inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i) |
ykuroda | 0:13a5d365ba16 | 265 | { |
ykuroda | 0:13a5d365ba16 | 266 | RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1)); |
ykuroda | 0:13a5d365ba16 | 267 | RealScalar sd = numext::norm1(m_matT.coeff(i+1,i)); |
ykuroda | 0:13a5d365ba16 | 268 | if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon())) |
ykuroda | 0:13a5d365ba16 | 269 | { |
ykuroda | 0:13a5d365ba16 | 270 | m_matT.coeffRef(i+1,i) = ComplexScalar(0); |
ykuroda | 0:13a5d365ba16 | 271 | return true; |
ykuroda | 0:13a5d365ba16 | 272 | } |
ykuroda | 0:13a5d365ba16 | 273 | return false; |
ykuroda | 0:13a5d365ba16 | 274 | } |
ykuroda | 0:13a5d365ba16 | 275 | |
ykuroda | 0:13a5d365ba16 | 276 | |
ykuroda | 0:13a5d365ba16 | 277 | /** Compute the shift in the current QR iteration. */ |
ykuroda | 0:13a5d365ba16 | 278 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 279 | typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter) |
ykuroda | 0:13a5d365ba16 | 280 | { |
ykuroda | 0:13a5d365ba16 | 281 | using std::abs; |
ykuroda | 0:13a5d365ba16 | 282 | if (iter == 10 || iter == 20) |
ykuroda | 0:13a5d365ba16 | 283 | { |
ykuroda | 0:13a5d365ba16 | 284 | // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f |
ykuroda | 0:13a5d365ba16 | 285 | return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2))); |
ykuroda | 0:13a5d365ba16 | 286 | } |
ykuroda | 0:13a5d365ba16 | 287 | |
ykuroda | 0:13a5d365ba16 | 288 | // compute the shift as one of the eigenvalues of t, the 2x2 |
ykuroda | 0:13a5d365ba16 | 289 | // diagonal block on the bottom of the active submatrix |
ykuroda | 0:13a5d365ba16 | 290 | Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1); |
ykuroda | 0:13a5d365ba16 | 291 | RealScalar normt = t.cwiseAbs().sum(); |
ykuroda | 0:13a5d365ba16 | 292 | t /= normt; // the normalization by sf is to avoid under/overflow |
ykuroda | 0:13a5d365ba16 | 293 | |
ykuroda | 0:13a5d365ba16 | 294 | ComplexScalar b = t.coeff(0,1) * t.coeff(1,0); |
ykuroda | 0:13a5d365ba16 | 295 | ComplexScalar c = t.coeff(0,0) - t.coeff(1,1); |
ykuroda | 0:13a5d365ba16 | 296 | ComplexScalar disc = sqrt(c*c + RealScalar(4)*b); |
ykuroda | 0:13a5d365ba16 | 297 | ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b; |
ykuroda | 0:13a5d365ba16 | 298 | ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1); |
ykuroda | 0:13a5d365ba16 | 299 | ComplexScalar eival1 = (trace + disc) / RealScalar(2); |
ykuroda | 0:13a5d365ba16 | 300 | ComplexScalar eival2 = (trace - disc) / RealScalar(2); |
ykuroda | 0:13a5d365ba16 | 301 | |
ykuroda | 0:13a5d365ba16 | 302 | if(numext::norm1(eival1) > numext::norm1(eival2)) |
ykuroda | 0:13a5d365ba16 | 303 | eival2 = det / eival1; |
ykuroda | 0:13a5d365ba16 | 304 | else |
ykuroda | 0:13a5d365ba16 | 305 | eival1 = det / eival2; |
ykuroda | 0:13a5d365ba16 | 306 | |
ykuroda | 0:13a5d365ba16 | 307 | // choose the eigenvalue closest to the bottom entry of the diagonal |
ykuroda | 0:13a5d365ba16 | 308 | if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1))) |
ykuroda | 0:13a5d365ba16 | 309 | return normt * eival1; |
ykuroda | 0:13a5d365ba16 | 310 | else |
ykuroda | 0:13a5d365ba16 | 311 | return normt * eival2; |
ykuroda | 0:13a5d365ba16 | 312 | } |
ykuroda | 0:13a5d365ba16 | 313 | |
ykuroda | 0:13a5d365ba16 | 314 | |
ykuroda | 0:13a5d365ba16 | 315 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 316 | ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool computeU) |
ykuroda | 0:13a5d365ba16 | 317 | { |
ykuroda | 0:13a5d365ba16 | 318 | m_matUisUptodate = false; |
ykuroda | 0:13a5d365ba16 | 319 | eigen_assert(matrix.cols() == matrix.rows()); |
ykuroda | 0:13a5d365ba16 | 320 | |
ykuroda | 0:13a5d365ba16 | 321 | if(matrix.cols() == 1) |
ykuroda | 0:13a5d365ba16 | 322 | { |
ykuroda | 0:13a5d365ba16 | 323 | m_matT = matrix.template cast<ComplexScalar>(); |
ykuroda | 0:13a5d365ba16 | 324 | if(computeU) m_matU = ComplexMatrixType::Identity(1,1); |
ykuroda | 0:13a5d365ba16 | 325 | m_info = Success; |
ykuroda | 0:13a5d365ba16 | 326 | m_isInitialized = true; |
ykuroda | 0:13a5d365ba16 | 327 | m_matUisUptodate = computeU; |
ykuroda | 0:13a5d365ba16 | 328 | return *this; |
ykuroda | 0:13a5d365ba16 | 329 | } |
ykuroda | 0:13a5d365ba16 | 330 | |
ykuroda | 0:13a5d365ba16 | 331 | internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix, computeU); |
ykuroda | 0:13a5d365ba16 | 332 | computeFromHessenberg(m_matT, m_matU, computeU); |
ykuroda | 0:13a5d365ba16 | 333 | return *this; |
ykuroda | 0:13a5d365ba16 | 334 | } |
ykuroda | 0:13a5d365ba16 | 335 | |
ykuroda | 0:13a5d365ba16 | 336 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 337 | template<typename HessMatrixType, typename OrthMatrixType> |
ykuroda | 0:13a5d365ba16 | 338 | ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU) |
ykuroda | 0:13a5d365ba16 | 339 | { |
ykuroda | 0:13a5d365ba16 | 340 | m_matT = matrixH; |
ykuroda | 0:13a5d365ba16 | 341 | if(computeU) |
ykuroda | 0:13a5d365ba16 | 342 | m_matU = matrixQ; |
ykuroda | 0:13a5d365ba16 | 343 | reduceToTriangularForm(computeU); |
ykuroda | 0:13a5d365ba16 | 344 | return *this; |
ykuroda | 0:13a5d365ba16 | 345 | } |
ykuroda | 0:13a5d365ba16 | 346 | namespace internal { |
ykuroda | 0:13a5d365ba16 | 347 | |
ykuroda | 0:13a5d365ba16 | 348 | /* Reduce given matrix to Hessenberg form */ |
ykuroda | 0:13a5d365ba16 | 349 | template<typename MatrixType, bool IsComplex> |
ykuroda | 0:13a5d365ba16 | 350 | struct complex_schur_reduce_to_hessenberg |
ykuroda | 0:13a5d365ba16 | 351 | { |
ykuroda | 0:13a5d365ba16 | 352 | // this is the implementation for the case IsComplex = true |
ykuroda | 0:13a5d365ba16 | 353 | static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) |
ykuroda | 0:13a5d365ba16 | 354 | { |
ykuroda | 0:13a5d365ba16 | 355 | _this.m_hess.compute(matrix); |
ykuroda | 0:13a5d365ba16 | 356 | _this.m_matT = _this.m_hess.matrixH(); |
ykuroda | 0:13a5d365ba16 | 357 | if(computeU) _this.m_matU = _this.m_hess.matrixQ(); |
ykuroda | 0:13a5d365ba16 | 358 | } |
ykuroda | 0:13a5d365ba16 | 359 | }; |
ykuroda | 0:13a5d365ba16 | 360 | |
ykuroda | 0:13a5d365ba16 | 361 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 362 | struct complex_schur_reduce_to_hessenberg<MatrixType, false> |
ykuroda | 0:13a5d365ba16 | 363 | { |
ykuroda | 0:13a5d365ba16 | 364 | static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU) |
ykuroda | 0:13a5d365ba16 | 365 | { |
ykuroda | 0:13a5d365ba16 | 366 | typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar; |
ykuroda | 0:13a5d365ba16 | 367 | |
ykuroda | 0:13a5d365ba16 | 368 | // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar |
ykuroda | 0:13a5d365ba16 | 369 | _this.m_hess.compute(matrix); |
ykuroda | 0:13a5d365ba16 | 370 | _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>(); |
ykuroda | 0:13a5d365ba16 | 371 | if(computeU) |
ykuroda | 0:13a5d365ba16 | 372 | { |
ykuroda | 0:13a5d365ba16 | 373 | // This may cause an allocation which seems to be avoidable |
ykuroda | 0:13a5d365ba16 | 374 | MatrixType Q = _this.m_hess.matrixQ(); |
ykuroda | 0:13a5d365ba16 | 375 | _this.m_matU = Q.template cast<ComplexScalar>(); |
ykuroda | 0:13a5d365ba16 | 376 | } |
ykuroda | 0:13a5d365ba16 | 377 | } |
ykuroda | 0:13a5d365ba16 | 378 | }; |
ykuroda | 0:13a5d365ba16 | 379 | |
ykuroda | 0:13a5d365ba16 | 380 | } // end namespace internal |
ykuroda | 0:13a5d365ba16 | 381 | |
ykuroda | 0:13a5d365ba16 | 382 | // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration. |
ykuroda | 0:13a5d365ba16 | 383 | template<typename MatrixType> |
ykuroda | 0:13a5d365ba16 | 384 | void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU) |
ykuroda | 0:13a5d365ba16 | 385 | { |
ykuroda | 0:13a5d365ba16 | 386 | Index maxIters = m_maxIters; |
ykuroda | 0:13a5d365ba16 | 387 | if (maxIters == -1) |
ykuroda | 0:13a5d365ba16 | 388 | maxIters = m_maxIterationsPerRow * m_matT.rows(); |
ykuroda | 0:13a5d365ba16 | 389 | |
ykuroda | 0:13a5d365ba16 | 390 | // The matrix m_matT is divided in three parts. |
ykuroda | 0:13a5d365ba16 | 391 | // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero. |
ykuroda | 0:13a5d365ba16 | 392 | // Rows il,...,iu is the part we are working on (the active submatrix). |
ykuroda | 0:13a5d365ba16 | 393 | // Rows iu+1,...,end are already brought in triangular form. |
ykuroda | 0:13a5d365ba16 | 394 | Index iu = m_matT.cols() - 1; |
ykuroda | 0:13a5d365ba16 | 395 | Index il; |
ykuroda | 0:13a5d365ba16 | 396 | Index iter = 0; // number of iterations we are working on the (iu,iu) element |
ykuroda | 0:13a5d365ba16 | 397 | Index totalIter = 0; // number of iterations for whole matrix |
ykuroda | 0:13a5d365ba16 | 398 | |
ykuroda | 0:13a5d365ba16 | 399 | while(true) |
ykuroda | 0:13a5d365ba16 | 400 | { |
ykuroda | 0:13a5d365ba16 | 401 | // find iu, the bottom row of the active submatrix |
ykuroda | 0:13a5d365ba16 | 402 | while(iu > 0) |
ykuroda | 0:13a5d365ba16 | 403 | { |
ykuroda | 0:13a5d365ba16 | 404 | if(!subdiagonalEntryIsNeglegible(iu-1)) break; |
ykuroda | 0:13a5d365ba16 | 405 | iter = 0; |
ykuroda | 0:13a5d365ba16 | 406 | --iu; |
ykuroda | 0:13a5d365ba16 | 407 | } |
ykuroda | 0:13a5d365ba16 | 408 | |
ykuroda | 0:13a5d365ba16 | 409 | // if iu is zero then we are done; the whole matrix is triangularized |
ykuroda | 0:13a5d365ba16 | 410 | if(iu==0) break; |
ykuroda | 0:13a5d365ba16 | 411 | |
ykuroda | 0:13a5d365ba16 | 412 | // if we spent too many iterations, we give up |
ykuroda | 0:13a5d365ba16 | 413 | iter++; |
ykuroda | 0:13a5d365ba16 | 414 | totalIter++; |
ykuroda | 0:13a5d365ba16 | 415 | if(totalIter > maxIters) break; |
ykuroda | 0:13a5d365ba16 | 416 | |
ykuroda | 0:13a5d365ba16 | 417 | // find il, the top row of the active submatrix |
ykuroda | 0:13a5d365ba16 | 418 | il = iu-1; |
ykuroda | 0:13a5d365ba16 | 419 | while(il > 0 && !subdiagonalEntryIsNeglegible(il-1)) |
ykuroda | 0:13a5d365ba16 | 420 | { |
ykuroda | 0:13a5d365ba16 | 421 | --il; |
ykuroda | 0:13a5d365ba16 | 422 | } |
ykuroda | 0:13a5d365ba16 | 423 | |
ykuroda | 0:13a5d365ba16 | 424 | /* perform the QR step using Givens rotations. The first rotation |
ykuroda | 0:13a5d365ba16 | 425 | creates a bulge; the (il+2,il) element becomes nonzero. This |
ykuroda | 0:13a5d365ba16 | 426 | bulge is chased down to the bottom of the active submatrix. */ |
ykuroda | 0:13a5d365ba16 | 427 | |
ykuroda | 0:13a5d365ba16 | 428 | ComplexScalar shift = computeShift(iu, iter); |
ykuroda | 0:13a5d365ba16 | 429 | JacobiRotation<ComplexScalar> rot; |
ykuroda | 0:13a5d365ba16 | 430 | rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il)); |
ykuroda | 0:13a5d365ba16 | 431 | m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint()); |
ykuroda | 0:13a5d365ba16 | 432 | m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot); |
ykuroda | 0:13a5d365ba16 | 433 | if(computeU) m_matU.applyOnTheRight(il, il+1, rot); |
ykuroda | 0:13a5d365ba16 | 434 | |
ykuroda | 0:13a5d365ba16 | 435 | for(Index i=il+1 ; i<iu ; i++) |
ykuroda | 0:13a5d365ba16 | 436 | { |
ykuroda | 0:13a5d365ba16 | 437 | rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1)); |
ykuroda | 0:13a5d365ba16 | 438 | m_matT.coeffRef(i+1,i-1) = ComplexScalar(0); |
ykuroda | 0:13a5d365ba16 | 439 | m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint()); |
ykuroda | 0:13a5d365ba16 | 440 | m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot); |
ykuroda | 0:13a5d365ba16 | 441 | if(computeU) m_matU.applyOnTheRight(i, i+1, rot); |
ykuroda | 0:13a5d365ba16 | 442 | } |
ykuroda | 0:13a5d365ba16 | 443 | } |
ykuroda | 0:13a5d365ba16 | 444 | |
ykuroda | 0:13a5d365ba16 | 445 | if(totalIter <= maxIters) |
ykuroda | 0:13a5d365ba16 | 446 | m_info = Success; |
ykuroda | 0:13a5d365ba16 | 447 | else |
ykuroda | 0:13a5d365ba16 | 448 | m_info = NoConvergence; |
ykuroda | 0:13a5d365ba16 | 449 | |
ykuroda | 0:13a5d365ba16 | 450 | m_isInitialized = true; |
ykuroda | 0:13a5d365ba16 | 451 | m_matUisUptodate = computeU; |
ykuroda | 0:13a5d365ba16 | 452 | } |
ykuroda | 0:13a5d365ba16 | 453 | |
ykuroda | 0:13a5d365ba16 | 454 | } // end namespace Eigen |
ykuroda | 0:13a5d365ba16 | 455 | |
ykuroda | 0:13a5d365ba16 | 456 | #endif // EIGEN_COMPLEX_SCHUR_H |