diff options
Diffstat (limited to 'runtimes/nn/depend/external/eigen/Eigen/src/Eigenvalues/ComplexSchur.h')
| -rw-r--r-- | runtimes/nn/depend/external/eigen/Eigen/src/Eigenvalues/ComplexSchur.h | 459 |
1 files changed, 0 insertions, 459 deletions
diff --git a/runtimes/nn/depend/external/eigen/Eigen/src/Eigenvalues/ComplexSchur.h b/runtimes/nn/depend/external/eigen/Eigen/src/Eigenvalues/ComplexSchur.h deleted file mode 100644 index 7f38919f7..000000000 --- a/runtimes/nn/depend/external/eigen/Eigen/src/Eigenvalues/ComplexSchur.h +++ /dev/null @@ -1,459 +0,0 @@ -// 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); - - if(numext::norm1(eival1) > numext::norm1(eival2)) - eival2 = det / eival1; - else - 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 |