123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462 |
- // This file is part of Eigen, a lightweight C++ template library
- // for linear algebra.
- //
- // Copyright (C) 2009 Claire Maurice
- // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
- // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
- //
- // This Source Code Form is subject to the terms of the Mozilla
- // Public License v. 2.0. If a copy of the MPL was not distributed
- // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
- #ifndef EIGEN_COMPLEX_SCHUR_H
- #define EIGEN_COMPLEX_SCHUR_H
- #include "./HessenbergDecomposition.h"
- namespace Eigen {
- namespace internal {
- template<typename MatrixType, bool IsComplex> struct complex_schur_reduce_to_hessenberg;
- }
- /** \eigenvalues_module \ingroup Eigenvalues_Module
- *
- *
- * \class ComplexSchur
- *
- * \brief Performs a complex Schur decomposition of a real or complex square matrix
- *
- * \tparam _MatrixType the type of the matrix of which we are
- * computing the Schur decomposition; this is expected to be an
- * instantiation of the Matrix class template.
- *
- * Given a real or complex square matrix A, this class computes the
- * Schur decomposition: \f$ A = U T U^*\f$ where U is a unitary
- * complex matrix, and T is a complex upper triangular matrix. The
- * diagonal of the matrix T corresponds to the eigenvalues of the
- * matrix A.
- *
- * Call the function compute() to compute the Schur decomposition of
- * a given matrix. Alternatively, you can use the
- * ComplexSchur(const MatrixType&, bool) constructor which computes
- * the Schur decomposition at construction time. Once the
- * decomposition is computed, you can use the matrixU() and matrixT()
- * functions to retrieve the matrices U and V in the decomposition.
- *
- * \note This code is inspired from Jampack
- *
- * \sa class RealSchur, class EigenSolver, class ComplexEigenSolver
- */
- template<typename _MatrixType> class ComplexSchur
- {
- public:
- typedef _MatrixType MatrixType;
- enum {
- RowsAtCompileTime = MatrixType::RowsAtCompileTime,
- ColsAtCompileTime = MatrixType::ColsAtCompileTime,
- Options = MatrixType::Options,
- MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
- MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
- };
- /** \brief Scalar type for matrices of type \p _MatrixType. */
- typedef typename MatrixType::Scalar Scalar;
- typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
- /** \brief Complex scalar type for \p _MatrixType.
- *
- * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
- * \c float or \c double) and just \c Scalar if #Scalar is
- * complex.
- */
- typedef std::complex<RealScalar> ComplexScalar;
- /** \brief Type for the matrices in the Schur decomposition.
- *
- * This is a square matrix with entries of type #ComplexScalar.
- * The size is the same as the size of \p _MatrixType.
- */
- typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> ComplexMatrixType;
- /** \brief Default constructor.
- *
- * \param [in] size Positive integer, size of the matrix whose Schur decomposition will be computed.
- *
- * The default constructor is useful in cases in which the user
- * intends to perform decompositions via compute(). The \p size
- * parameter is only used as a hint. It is not an error to give a
- * wrong \p size, but it may impair performance.
- *
- * \sa compute() for an example.
- */
- explicit ComplexSchur(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
- : m_matT(size,size),
- m_matU(size,size),
- m_hess(size),
- m_isInitialized(false),
- m_matUisUptodate(false),
- m_maxIters(-1)
- {}
- /** \brief Constructor; computes Schur decomposition of given matrix.
- *
- * \param[in] matrix Square matrix whose Schur decomposition is to be computed.
- * \param[in] computeU If true, both T and U are computed; if false, only T is computed.
- *
- * This constructor calls compute() to compute the Schur decomposition.
- *
- * \sa matrixT() and matrixU() for examples.
- */
- template<typename InputType>
- explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true)
- : m_matT(matrix.rows(),matrix.cols()),
- m_matU(matrix.rows(),matrix.cols()),
- m_hess(matrix.rows()),
- m_isInitialized(false),
- m_matUisUptodate(false),
- m_maxIters(-1)
- {
- compute(matrix.derived(), computeU);
- }
- /** \brief Returns the unitary matrix in the Schur decomposition.
- *
- * \returns A const reference to the matrix U.
- *
- * It is assumed that either the constructor
- * ComplexSchur(const MatrixType& matrix, bool computeU) or the
- * member function compute(const MatrixType& matrix, bool computeU)
- * has been called before to compute the Schur decomposition of a
- * matrix, and that \p computeU was set to true (the default
- * value).
- *
- * Example: \include ComplexSchur_matrixU.cpp
- * Output: \verbinclude ComplexSchur_matrixU.out
- */
- const ComplexMatrixType& matrixU() const
- {
- eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
- eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
- return m_matU;
- }
- /** \brief Returns the triangular matrix in the Schur decomposition.
- *
- * \returns A const reference to the matrix T.
- *
- * It is assumed that either the constructor
- * ComplexSchur(const MatrixType& matrix, bool computeU) or the
- * member function compute(const MatrixType& matrix, bool computeU)
- * has been called before to compute the Schur decomposition of a
- * matrix.
- *
- * Note that this function returns a plain square matrix. If you want to reference
- * only the upper triangular part, use:
- * \code schur.matrixT().triangularView<Upper>() \endcode
- *
- * Example: \include ComplexSchur_matrixT.cpp
- * Output: \verbinclude ComplexSchur_matrixT.out
- */
- const ComplexMatrixType& matrixT() const
- {
- eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
- return m_matT;
- }
- /** \brief Computes Schur decomposition of given matrix.
- *
- * \param[in] matrix Square matrix whose Schur decomposition is to be computed.
- * \param[in] computeU If true, both T and U are computed; if false, only T is computed.
- * \returns Reference to \c *this
- *
- * The Schur decomposition is computed by first reducing the
- * matrix to Hessenberg form using the class
- * HessenbergDecomposition. The Hessenberg matrix is then reduced
- * to triangular form by performing QR iterations with a single
- * shift. The cost of computing the Schur decomposition depends
- * on the number of iterations; as a rough guide, it may be taken
- * on the number of iterations; as a rough guide, it may be taken
- * to be \f$25n^3\f$ complex flops, or \f$10n^3\f$ complex flops
- * if \a computeU is false.
- *
- * Example: \include ComplexSchur_compute.cpp
- * Output: \verbinclude ComplexSchur_compute.out
- *
- * \sa compute(const MatrixType&, bool, Index)
- */
- template<typename InputType>
- ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);
-
- /** \brief Compute Schur decomposition from a given Hessenberg matrix
- * \param[in] matrixH Matrix in Hessenberg form H
- * \param[in] matrixQ orthogonal matrix Q that transform a matrix A to H : A = Q H Q^T
- * \param computeU Computes the matriX U of the Schur vectors
- * \return Reference to \c *this
- *
- * This routine assumes that the matrix is already reduced in Hessenberg form matrixH
- * using either the class HessenbergDecomposition or another mean.
- * It computes the upper quasi-triangular matrix T of the Schur decomposition of H
- * When computeU is true, this routine computes the matrix U such that
- * A = U T U^T = (QZ) T (QZ)^T = Q H Q^T where A is the initial matrix
- *
- * NOTE Q is referenced if computeU is true; so, if the initial orthogonal matrix
- * is not available, the user should give an identity matrix (Q.setIdentity())
- *
- * \sa compute(const MatrixType&, bool)
- */
- template<typename HessMatrixType, typename OrthMatrixType>
- ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU=true);
- /** \brief Reports whether previous computation was successful.
- *
- * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
- */
- ComputationInfo info() const
- {
- eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
- return m_info;
- }
- /** \brief Sets the maximum number of iterations allowed.
- *
- * If not specified by the user, the maximum number of iterations is m_maxIterationsPerRow times the size
- * of the matrix.
- */
- ComplexSchur& setMaxIterations(Index maxIters)
- {
- m_maxIters = maxIters;
- return *this;
- }
- /** \brief Returns the maximum number of iterations. */
- Index getMaxIterations()
- {
- return m_maxIters;
- }
- /** \brief Maximum number of iterations per row.
- *
- * If not otherwise specified, the maximum number of iterations is this number times the size of the
- * matrix. It is currently set to 30.
- */
- static const int m_maxIterationsPerRow = 30;
- protected:
- ComplexMatrixType m_matT, m_matU;
- HessenbergDecomposition<MatrixType> m_hess;
- ComputationInfo m_info;
- bool m_isInitialized;
- bool m_matUisUptodate;
- Index m_maxIters;
- private:
- bool subdiagonalEntryIsNeglegible(Index i);
- ComplexScalar computeShift(Index iu, Index iter);
- void reduceToTriangularForm(bool computeU);
- friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
- };
- /** If m_matT(i+1,i) is neglegible in floating point arithmetic
- * compared to m_matT(i,i) and m_matT(j,j), then set it to zero and
- * return true, else return false. */
- template<typename MatrixType>
- inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i)
- {
- RealScalar d = numext::norm1(m_matT.coeff(i,i)) + numext::norm1(m_matT.coeff(i+1,i+1));
- RealScalar sd = numext::norm1(m_matT.coeff(i+1,i));
- if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon()))
- {
- m_matT.coeffRef(i+1,i) = ComplexScalar(0);
- return true;
- }
- return false;
- }
- /** Compute the shift in the current QR iteration. */
- template<typename MatrixType>
- typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter)
- {
- using std::abs;
- if (iter == 10 || iter == 20)
- {
- // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
- return abs(numext::real(m_matT.coeff(iu,iu-1))) + abs(numext::real(m_matT.coeff(iu-1,iu-2)));
- }
- // compute the shift as one of the eigenvalues of t, the 2x2
- // diagonal block on the bottom of the active submatrix
- Matrix<ComplexScalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
- RealScalar normt = t.cwiseAbs().sum();
- t /= normt; // the normalization by sf is to avoid under/overflow
- ComplexScalar b = t.coeff(0,1) * t.coeff(1,0);
- ComplexScalar c = t.coeff(0,0) - t.coeff(1,1);
- ComplexScalar disc = sqrt(c*c + RealScalar(4)*b);
- ComplexScalar det = t.coeff(0,0) * t.coeff(1,1) - b;
- ComplexScalar trace = t.coeff(0,0) + t.coeff(1,1);
- ComplexScalar eival1 = (trace + disc) / RealScalar(2);
- ComplexScalar eival2 = (trace - disc) / RealScalar(2);
- RealScalar eival1_norm = numext::norm1(eival1);
- RealScalar eival2_norm = numext::norm1(eival2);
- // A division by zero can only occur if eival1==eival2==0.
- // In this case, det==0, and all we have to do is checking that eival2_norm!=0
- if(eival1_norm > eival2_norm)
- eival2 = det / eival1;
- else if(eival2_norm!=RealScalar(0))
- eival1 = det / eival2;
- // choose the eigenvalue closest to the bottom entry of the diagonal
- if(numext::norm1(eival1-t.coeff(1,1)) < numext::norm1(eival2-t.coeff(1,1)))
- return normt * eival1;
- else
- return normt * eival2;
- }
- template<typename MatrixType>
- template<typename InputType>
- ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
- {
- m_matUisUptodate = false;
- eigen_assert(matrix.cols() == matrix.rows());
- if(matrix.cols() == 1)
- {
- m_matT = matrix.derived().template cast<ComplexScalar>();
- if(computeU) m_matU = ComplexMatrixType::Identity(1,1);
- m_info = Success;
- m_isInitialized = true;
- m_matUisUptodate = computeU;
- return *this;
- }
- internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>::run(*this, matrix.derived(), computeU);
- computeFromHessenberg(m_matT, m_matU, computeU);
- return *this;
- }
- template<typename MatrixType>
- template<typename HessMatrixType, typename OrthMatrixType>
- ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ, bool computeU)
- {
- m_matT = matrixH;
- if(computeU)
- m_matU = matrixQ;
- reduceToTriangularForm(computeU);
- return *this;
- }
- namespace internal {
- /* Reduce given matrix to Hessenberg form */
- template<typename MatrixType, bool IsComplex>
- struct complex_schur_reduce_to_hessenberg
- {
- // this is the implementation for the case IsComplex = true
- static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
- {
- _this.m_hess.compute(matrix);
- _this.m_matT = _this.m_hess.matrixH();
- if(computeU) _this.m_matU = _this.m_hess.matrixQ();
- }
- };
- template<typename MatrixType>
- struct complex_schur_reduce_to_hessenberg<MatrixType, false>
- {
- static void run(ComplexSchur<MatrixType>& _this, const MatrixType& matrix, bool computeU)
- {
- typedef typename ComplexSchur<MatrixType>::ComplexScalar ComplexScalar;
- // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
- _this.m_hess.compute(matrix);
- _this.m_matT = _this.m_hess.matrixH().template cast<ComplexScalar>();
- if(computeU)
- {
- // This may cause an allocation which seems to be avoidable
- MatrixType Q = _this.m_hess.matrixQ();
- _this.m_matU = Q.template cast<ComplexScalar>();
- }
- }
- };
- } // end namespace internal
- // Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
- template<typename MatrixType>
- void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU)
- {
- Index maxIters = m_maxIters;
- if (maxIters == -1)
- maxIters = m_maxIterationsPerRow * m_matT.rows();
- // The matrix m_matT is divided in three parts.
- // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
- // Rows il,...,iu is the part we are working on (the active submatrix).
- // Rows iu+1,...,end are already brought in triangular form.
- Index iu = m_matT.cols() - 1;
- Index il;
- Index iter = 0; // number of iterations we are working on the (iu,iu) element
- Index totalIter = 0; // number of iterations for whole matrix
- while(true)
- {
- // find iu, the bottom row of the active submatrix
- while(iu > 0)
- {
- if(!subdiagonalEntryIsNeglegible(iu-1)) break;
- iter = 0;
- --iu;
- }
- // if iu is zero then we are done; the whole matrix is triangularized
- if(iu==0) break;
- // if we spent too many iterations, we give up
- iter++;
- totalIter++;
- if(totalIter > maxIters) break;
- // find il, the top row of the active submatrix
- il = iu-1;
- while(il > 0 && !subdiagonalEntryIsNeglegible(il-1))
- {
- --il;
- }
- /* perform the QR step using Givens rotations. The first rotation
- creates a bulge; the (il+2,il) element becomes nonzero. This
- bulge is chased down to the bottom of the active submatrix. */
- ComplexScalar shift = computeShift(iu, iter);
- JacobiRotation<ComplexScalar> rot;
- rot.makeGivens(m_matT.coeff(il,il) - shift, m_matT.coeff(il+1,il));
- m_matT.rightCols(m_matT.cols()-il).applyOnTheLeft(il, il+1, rot.adjoint());
- m_matT.topRows((std::min)(il+2,iu)+1).applyOnTheRight(il, il+1, rot);
- if(computeU) m_matU.applyOnTheRight(il, il+1, rot);
- for(Index i=il+1 ; i<iu ; i++)
- {
- rot.makeGivens(m_matT.coeffRef(i,i-1), m_matT.coeffRef(i+1,i-1), &m_matT.coeffRef(i,i-1));
- m_matT.coeffRef(i+1,i-1) = ComplexScalar(0);
- m_matT.rightCols(m_matT.cols()-i).applyOnTheLeft(i, i+1, rot.adjoint());
- m_matT.topRows((std::min)(i+2,iu)+1).applyOnTheRight(i, i+1, rot);
- if(computeU) m_matU.applyOnTheRight(i, i+1, rot);
- }
- }
- if(totalIter <= maxIters)
- m_info = Success;
- else
- m_info = NoConvergence;
- m_isInitialized = true;
- m_matUisUptodate = computeU;
- }
- } // end namespace Eigen
- #endif // EIGEN_COMPLEX_SCHUR_H
|