Michael Ernst Peter / Eigen

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?

UserRevisionLine numberNew contents of line
ykuroda 0:13a5d365ba16 1 // This file is part of Eigen, a lightweight C++ template library
ykuroda 0:13a5d365ba16 2 // for linear algebra.
ykuroda 0:13a5d365ba16 3 //
ykuroda 0:13a5d365ba16 4 // Copyright (C) 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