Yoji KURODA / Eigen

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/EigenSolver.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) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
ykuroda 0:13a5d365ba16 5 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
ykuroda 0:13a5d365ba16 6 //
ykuroda 0:13a5d365ba16 7 // This Source Code Form is subject to the terms of the Mozilla
ykuroda 0:13a5d365ba16 8 // Public License v. 2.0. If a copy of the MPL was not distributed
ykuroda 0:13a5d365ba16 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
ykuroda 0:13a5d365ba16 10
ykuroda 0:13a5d365ba16 11 #ifndef EIGEN_EIGENSOLVER_H
ykuroda 0:13a5d365ba16 12 #define EIGEN_EIGENSOLVER_H
ykuroda 0:13a5d365ba16 13
ykuroda 0:13a5d365ba16 14 #include "./RealSchur.h"
ykuroda 0:13a5d365ba16 15
ykuroda 0:13a5d365ba16 16 namespace Eigen {
ykuroda 0:13a5d365ba16 17
ykuroda 0:13a5d365ba16 18 /** \eigenvalues_module \ingroup Eigenvalues_Module
ykuroda 0:13a5d365ba16 19 *
ykuroda 0:13a5d365ba16 20 *
ykuroda 0:13a5d365ba16 21 * \class EigenSolver
ykuroda 0:13a5d365ba16 22 *
ykuroda 0:13a5d365ba16 23 * \brief Computes eigenvalues and eigenvectors of general matrices
ykuroda 0:13a5d365ba16 24 *
ykuroda 0:13a5d365ba16 25 * \tparam _MatrixType the type of the matrix of which we are computing the
ykuroda 0:13a5d365ba16 26 * eigendecomposition; this is expected to be an instantiation of the Matrix
ykuroda 0:13a5d365ba16 27 * class template. Currently, only real matrices are supported.
ykuroda 0:13a5d365ba16 28 *
ykuroda 0:13a5d365ba16 29 * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
ykuroda 0:13a5d365ba16 30 * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v \f$. If
ykuroda 0:13a5d365ba16 31 * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
ykuroda 0:13a5d365ba16 32 * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
ykuroda 0:13a5d365ba16 33 * V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
ykuroda 0:13a5d365ba16 34 * have \f$ A = V D V^{-1} \f$. This is called the eigendecomposition.
ykuroda 0:13a5d365ba16 35 *
ykuroda 0:13a5d365ba16 36 * The eigenvalues and eigenvectors of a matrix may be complex, even when the
ykuroda 0:13a5d365ba16 37 * matrix is real. However, we can choose real matrices \f$ V \f$ and \f$ D
ykuroda 0:13a5d365ba16 38 * \f$ satisfying \f$ A V = V D \f$, just like the eigendecomposition, if the
ykuroda 0:13a5d365ba16 39 * matrix \f$ D \f$ is not required to be diagonal, but if it is allowed to
ykuroda 0:13a5d365ba16 40 * have blocks of the form
ykuroda 0:13a5d365ba16 41 * \f[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f]
ykuroda 0:13a5d365ba16 42 * (where \f$ u \f$ and \f$ v \f$ are real numbers) on the diagonal. These
ykuroda 0:13a5d365ba16 43 * blocks correspond to complex eigenvalue pairs \f$ u \pm iv \f$. We call
ykuroda 0:13a5d365ba16 44 * this variant of the eigendecomposition the pseudo-eigendecomposition.
ykuroda 0:13a5d365ba16 45 *
ykuroda 0:13a5d365ba16 46 * Call the function compute() to compute the eigenvalues and eigenvectors of
ykuroda 0:13a5d365ba16 47 * a given matrix. Alternatively, you can use the
ykuroda 0:13a5d365ba16 48 * EigenSolver(const MatrixType&, bool) constructor which computes the
ykuroda 0:13a5d365ba16 49 * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
ykuroda 0:13a5d365ba16 50 * eigenvectors are computed, they can be retrieved with the eigenvalues() and
ykuroda 0:13a5d365ba16 51 * eigenvectors() functions. The pseudoEigenvalueMatrix() and
ykuroda 0:13a5d365ba16 52 * pseudoEigenvectors() methods allow the construction of the
ykuroda 0:13a5d365ba16 53 * pseudo-eigendecomposition.
ykuroda 0:13a5d365ba16 54 *
ykuroda 0:13a5d365ba16 55 * The documentation for EigenSolver(const MatrixType&, bool) contains an
ykuroda 0:13a5d365ba16 56 * example of the typical use of this class.
ykuroda 0:13a5d365ba16 57 *
ykuroda 0:13a5d365ba16 58 * \note The implementation is adapted from
ykuroda 0:13a5d365ba16 59 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
ykuroda 0:13a5d365ba16 60 * Their code is based on EISPACK.
ykuroda 0:13a5d365ba16 61 *
ykuroda 0:13a5d365ba16 62 * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
ykuroda 0:13a5d365ba16 63 */
ykuroda 0:13a5d365ba16 64 template<typename _MatrixType> class EigenSolver
ykuroda 0:13a5d365ba16 65 {
ykuroda 0:13a5d365ba16 66 public:
ykuroda 0:13a5d365ba16 67
ykuroda 0:13a5d365ba16 68 /** \brief Synonym for the template parameter \p _MatrixType. */
ykuroda 0:13a5d365ba16 69 typedef _MatrixType MatrixType;
ykuroda 0:13a5d365ba16 70
ykuroda 0:13a5d365ba16 71 enum {
ykuroda 0:13a5d365ba16 72 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ykuroda 0:13a5d365ba16 73 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
ykuroda 0:13a5d365ba16 74 Options = MatrixType::Options,
ykuroda 0:13a5d365ba16 75 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
ykuroda 0:13a5d365ba16 76 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
ykuroda 0:13a5d365ba16 77 };
ykuroda 0:13a5d365ba16 78
ykuroda 0:13a5d365ba16 79 /** \brief Scalar type for matrices of type #MatrixType. */
ykuroda 0:13a5d365ba16 80 typedef typename MatrixType::Scalar Scalar;
ykuroda 0:13a5d365ba16 81 typedef typename NumTraits<Scalar>::Real RealScalar;
ykuroda 0:13a5d365ba16 82 typedef typename MatrixType::Index Index;
ykuroda 0:13a5d365ba16 83
ykuroda 0:13a5d365ba16 84 /** \brief Complex scalar type for #MatrixType.
ykuroda 0:13a5d365ba16 85 *
ykuroda 0:13a5d365ba16 86 * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
ykuroda 0:13a5d365ba16 87 * \c float or \c double) and just \c Scalar if #Scalar is
ykuroda 0:13a5d365ba16 88 * complex.
ykuroda 0:13a5d365ba16 89 */
ykuroda 0:13a5d365ba16 90 typedef std::complex<RealScalar> ComplexScalar;
ykuroda 0:13a5d365ba16 91
ykuroda 0:13a5d365ba16 92 /** \brief Type for vector of eigenvalues as returned by eigenvalues().
ykuroda 0:13a5d365ba16 93 *
ykuroda 0:13a5d365ba16 94 * This is a column vector with entries of type #ComplexScalar.
ykuroda 0:13a5d365ba16 95 * The length of the vector is the size of #MatrixType.
ykuroda 0:13a5d365ba16 96 */
ykuroda 0:13a5d365ba16 97 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
ykuroda 0:13a5d365ba16 98
ykuroda 0:13a5d365ba16 99 /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
ykuroda 0:13a5d365ba16 100 *
ykuroda 0:13a5d365ba16 101 * This is a square matrix with entries of type #ComplexScalar.
ykuroda 0:13a5d365ba16 102 * The size is the same as the size of #MatrixType.
ykuroda 0:13a5d365ba16 103 */
ykuroda 0:13a5d365ba16 104 typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
ykuroda 0:13a5d365ba16 105
ykuroda 0:13a5d365ba16 106 /** \brief Default constructor.
ykuroda 0:13a5d365ba16 107 *
ykuroda 0:13a5d365ba16 108 * The default constructor is useful in cases in which the user intends to
ykuroda 0:13a5d365ba16 109 * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
ykuroda 0:13a5d365ba16 110 *
ykuroda 0:13a5d365ba16 111 * \sa compute() for an example.
ykuroda 0:13a5d365ba16 112 */
ykuroda 0:13a5d365ba16 113 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
ykuroda 0:13a5d365ba16 114
ykuroda 0:13a5d365ba16 115 /** \brief Default constructor with memory preallocation
ykuroda 0:13a5d365ba16 116 *
ykuroda 0:13a5d365ba16 117 * Like the default constructor but with preallocation of the internal data
ykuroda 0:13a5d365ba16 118 * according to the specified problem \a size.
ykuroda 0:13a5d365ba16 119 * \sa EigenSolver()
ykuroda 0:13a5d365ba16 120 */
ykuroda 0:13a5d365ba16 121 EigenSolver(Index size)
ykuroda 0:13a5d365ba16 122 : m_eivec(size, size),
ykuroda 0:13a5d365ba16 123 m_eivalues(size),
ykuroda 0:13a5d365ba16 124 m_isInitialized(false),
ykuroda 0:13a5d365ba16 125 m_eigenvectorsOk(false),
ykuroda 0:13a5d365ba16 126 m_realSchur(size),
ykuroda 0:13a5d365ba16 127 m_matT(size, size),
ykuroda 0:13a5d365ba16 128 m_tmp(size)
ykuroda 0:13a5d365ba16 129 {}
ykuroda 0:13a5d365ba16 130
ykuroda 0:13a5d365ba16 131 /** \brief Constructor; computes eigendecomposition of given matrix.
ykuroda 0:13a5d365ba16 132 *
ykuroda 0:13a5d365ba16 133 * \param[in] matrix Square matrix whose eigendecomposition is to be computed.
ykuroda 0:13a5d365ba16 134 * \param[in] computeEigenvectors If true, both the eigenvectors and the
ykuroda 0:13a5d365ba16 135 * eigenvalues are computed; if false, only the eigenvalues are
ykuroda 0:13a5d365ba16 136 * computed.
ykuroda 0:13a5d365ba16 137 *
ykuroda 0:13a5d365ba16 138 * This constructor calls compute() to compute the eigenvalues
ykuroda 0:13a5d365ba16 139 * and eigenvectors.
ykuroda 0:13a5d365ba16 140 *
ykuroda 0:13a5d365ba16 141 * Example: \include EigenSolver_EigenSolver_MatrixType.cpp
ykuroda 0:13a5d365ba16 142 * Output: \verbinclude EigenSolver_EigenSolver_MatrixType.out
ykuroda 0:13a5d365ba16 143 *
ykuroda 0:13a5d365ba16 144 * \sa compute()
ykuroda 0:13a5d365ba16 145 */
ykuroda 0:13a5d365ba16 146 EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
ykuroda 0:13a5d365ba16 147 : m_eivec(matrix.rows(), matrix.cols()),
ykuroda 0:13a5d365ba16 148 m_eivalues(matrix.cols()),
ykuroda 0:13a5d365ba16 149 m_isInitialized(false),
ykuroda 0:13a5d365ba16 150 m_eigenvectorsOk(false),
ykuroda 0:13a5d365ba16 151 m_realSchur(matrix.cols()),
ykuroda 0:13a5d365ba16 152 m_matT(matrix.rows(), matrix.cols()),
ykuroda 0:13a5d365ba16 153 m_tmp(matrix.cols())
ykuroda 0:13a5d365ba16 154 {
ykuroda 0:13a5d365ba16 155 compute(matrix, computeEigenvectors);
ykuroda 0:13a5d365ba16 156 }
ykuroda 0:13a5d365ba16 157
ykuroda 0:13a5d365ba16 158 /** \brief Returns the eigenvectors of given matrix.
ykuroda 0:13a5d365ba16 159 *
ykuroda 0:13a5d365ba16 160 * \returns %Matrix whose columns are the (possibly complex) eigenvectors.
ykuroda 0:13a5d365ba16 161 *
ykuroda 0:13a5d365ba16 162 * \pre Either the constructor
ykuroda 0:13a5d365ba16 163 * EigenSolver(const MatrixType&,bool) or the member function
ykuroda 0:13a5d365ba16 164 * compute(const MatrixType&, bool) has been called before, and
ykuroda 0:13a5d365ba16 165 * \p computeEigenvectors was set to true (the default).
ykuroda 0:13a5d365ba16 166 *
ykuroda 0:13a5d365ba16 167 * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
ykuroda 0:13a5d365ba16 168 * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The
ykuroda 0:13a5d365ba16 169 * eigenvectors are normalized to have (Euclidean) norm equal to one. The
ykuroda 0:13a5d365ba16 170 * matrix returned by this function is the matrix \f$ V \f$ in the
ykuroda 0:13a5d365ba16 171 * eigendecomposition \f$ A = V D V^{-1} \f$, if it exists.
ykuroda 0:13a5d365ba16 172 *
ykuroda 0:13a5d365ba16 173 * Example: \include EigenSolver_eigenvectors.cpp
ykuroda 0:13a5d365ba16 174 * Output: \verbinclude EigenSolver_eigenvectors.out
ykuroda 0:13a5d365ba16 175 *
ykuroda 0:13a5d365ba16 176 * \sa eigenvalues(), pseudoEigenvectors()
ykuroda 0:13a5d365ba16 177 */
ykuroda 0:13a5d365ba16 178 EigenvectorsType eigenvectors() const;
ykuroda 0:13a5d365ba16 179
ykuroda 0:13a5d365ba16 180 /** \brief Returns the pseudo-eigenvectors of given matrix.
ykuroda 0:13a5d365ba16 181 *
ykuroda 0:13a5d365ba16 182 * \returns Const reference to matrix whose columns are the pseudo-eigenvectors.
ykuroda 0:13a5d365ba16 183 *
ykuroda 0:13a5d365ba16 184 * \pre Either the constructor
ykuroda 0:13a5d365ba16 185 * EigenSolver(const MatrixType&,bool) or the member function
ykuroda 0:13a5d365ba16 186 * compute(const MatrixType&, bool) has been called before, and
ykuroda 0:13a5d365ba16 187 * \p computeEigenvectors was set to true (the default).
ykuroda 0:13a5d365ba16 188 *
ykuroda 0:13a5d365ba16 189 * The real matrix \f$ V \f$ returned by this function and the
ykuroda 0:13a5d365ba16 190 * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix()
ykuroda 0:13a5d365ba16 191 * satisfy \f$ AV = VD \f$.
ykuroda 0:13a5d365ba16 192 *
ykuroda 0:13a5d365ba16 193 * Example: \include EigenSolver_pseudoEigenvectors.cpp
ykuroda 0:13a5d365ba16 194 * Output: \verbinclude EigenSolver_pseudoEigenvectors.out
ykuroda 0:13a5d365ba16 195 *
ykuroda 0:13a5d365ba16 196 * \sa pseudoEigenvalueMatrix(), eigenvectors()
ykuroda 0:13a5d365ba16 197 */
ykuroda 0:13a5d365ba16 198 const MatrixType& pseudoEigenvectors() const
ykuroda 0:13a5d365ba16 199 {
ykuroda 0:13a5d365ba16 200 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
ykuroda 0:13a5d365ba16 201 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
ykuroda 0:13a5d365ba16 202 return m_eivec;
ykuroda 0:13a5d365ba16 203 }
ykuroda 0:13a5d365ba16 204
ykuroda 0:13a5d365ba16 205 /** \brief Returns the block-diagonal matrix in the pseudo-eigendecomposition.
ykuroda 0:13a5d365ba16 206 *
ykuroda 0:13a5d365ba16 207 * \returns A block-diagonal matrix.
ykuroda 0:13a5d365ba16 208 *
ykuroda 0:13a5d365ba16 209 * \pre Either the constructor
ykuroda 0:13a5d365ba16 210 * EigenSolver(const MatrixType&,bool) or the member function
ykuroda 0:13a5d365ba16 211 * compute(const MatrixType&, bool) has been called before.
ykuroda 0:13a5d365ba16 212 *
ykuroda 0:13a5d365ba16 213 * The matrix \f$ D \f$ returned by this function is real and
ykuroda 0:13a5d365ba16 214 * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2
ykuroda 0:13a5d365ba16 215 * blocks of the form
ykuroda 0:13a5d365ba16 216 * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$.
ykuroda 0:13a5d365ba16 217 * These blocks are not sorted in any particular order.
ykuroda 0:13a5d365ba16 218 * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by
ykuroda 0:13a5d365ba16 219 * pseudoEigenvectors() satisfy \f$ AV = VD \f$.
ykuroda 0:13a5d365ba16 220 *
ykuroda 0:13a5d365ba16 221 * \sa pseudoEigenvectors() for an example, eigenvalues()
ykuroda 0:13a5d365ba16 222 */
ykuroda 0:13a5d365ba16 223 MatrixType pseudoEigenvalueMatrix() const;
ykuroda 0:13a5d365ba16 224
ykuroda 0:13a5d365ba16 225 /** \brief Returns the eigenvalues of given matrix.
ykuroda 0:13a5d365ba16 226 *
ykuroda 0:13a5d365ba16 227 * \returns A const reference to the column vector containing the eigenvalues.
ykuroda 0:13a5d365ba16 228 *
ykuroda 0:13a5d365ba16 229 * \pre Either the constructor
ykuroda 0:13a5d365ba16 230 * EigenSolver(const MatrixType&,bool) or the member function
ykuroda 0:13a5d365ba16 231 * compute(const MatrixType&, bool) has been called before.
ykuroda 0:13a5d365ba16 232 *
ykuroda 0:13a5d365ba16 233 * The eigenvalues are repeated according to their algebraic multiplicity,
ykuroda 0:13a5d365ba16 234 * so there are as many eigenvalues as rows in the matrix. The eigenvalues
ykuroda 0:13a5d365ba16 235 * are not sorted in any particular order.
ykuroda 0:13a5d365ba16 236 *
ykuroda 0:13a5d365ba16 237 * Example: \include EigenSolver_eigenvalues.cpp
ykuroda 0:13a5d365ba16 238 * Output: \verbinclude EigenSolver_eigenvalues.out
ykuroda 0:13a5d365ba16 239 *
ykuroda 0:13a5d365ba16 240 * \sa eigenvectors(), pseudoEigenvalueMatrix(),
ykuroda 0:13a5d365ba16 241 * MatrixBase::eigenvalues()
ykuroda 0:13a5d365ba16 242 */
ykuroda 0:13a5d365ba16 243 const EigenvalueType& eigenvalues() const
ykuroda 0:13a5d365ba16 244 {
ykuroda 0:13a5d365ba16 245 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
ykuroda 0:13a5d365ba16 246 return m_eivalues;
ykuroda 0:13a5d365ba16 247 }
ykuroda 0:13a5d365ba16 248
ykuroda 0:13a5d365ba16 249 /** \brief Computes eigendecomposition of given matrix.
ykuroda 0:13a5d365ba16 250 *
ykuroda 0:13a5d365ba16 251 * \param[in] matrix Square matrix whose eigendecomposition is to be computed.
ykuroda 0:13a5d365ba16 252 * \param[in] computeEigenvectors If true, both the eigenvectors and the
ykuroda 0:13a5d365ba16 253 * eigenvalues are computed; if false, only the eigenvalues are
ykuroda 0:13a5d365ba16 254 * computed.
ykuroda 0:13a5d365ba16 255 * \returns Reference to \c *this
ykuroda 0:13a5d365ba16 256 *
ykuroda 0:13a5d365ba16 257 * This function computes the eigenvalues of the real matrix \p matrix.
ykuroda 0:13a5d365ba16 258 * The eigenvalues() function can be used to retrieve them. If
ykuroda 0:13a5d365ba16 259 * \p computeEigenvectors is true, then the eigenvectors are also computed
ykuroda 0:13a5d365ba16 260 * and can be retrieved by calling eigenvectors().
ykuroda 0:13a5d365ba16 261 *
ykuroda 0:13a5d365ba16 262 * The matrix is first reduced to real Schur form using the RealSchur
ykuroda 0:13a5d365ba16 263 * class. The Schur decomposition is then used to compute the eigenvalues
ykuroda 0:13a5d365ba16 264 * and eigenvectors.
ykuroda 0:13a5d365ba16 265 *
ykuroda 0:13a5d365ba16 266 * The cost of the computation is dominated by the cost of the
ykuroda 0:13a5d365ba16 267 * Schur decomposition, which is very approximately \f$ 25n^3 \f$
ykuroda 0:13a5d365ba16 268 * (where \f$ n \f$ is the size of the matrix) if \p computeEigenvectors
ykuroda 0:13a5d365ba16 269 * is true, and \f$ 10n^3 \f$ if \p computeEigenvectors is false.
ykuroda 0:13a5d365ba16 270 *
ykuroda 0:13a5d365ba16 271 * This method reuses of the allocated data in the EigenSolver object.
ykuroda 0:13a5d365ba16 272 *
ykuroda 0:13a5d365ba16 273 * Example: \include EigenSolver_compute.cpp
ykuroda 0:13a5d365ba16 274 * Output: \verbinclude EigenSolver_compute.out
ykuroda 0:13a5d365ba16 275 */
ykuroda 0:13a5d365ba16 276 EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
ykuroda 0:13a5d365ba16 277
ykuroda 0:13a5d365ba16 278 ComputationInfo info() const
ykuroda 0:13a5d365ba16 279 {
ykuroda 0:13a5d365ba16 280 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
ykuroda 0:13a5d365ba16 281 return m_realSchur.info();
ykuroda 0:13a5d365ba16 282 }
ykuroda 0:13a5d365ba16 283
ykuroda 0:13a5d365ba16 284 /** \brief Sets the maximum number of iterations allowed. */
ykuroda 0:13a5d365ba16 285 EigenSolver& setMaxIterations(Index maxIters)
ykuroda 0:13a5d365ba16 286 {
ykuroda 0:13a5d365ba16 287 m_realSchur.setMaxIterations(maxIters);
ykuroda 0:13a5d365ba16 288 return *this;
ykuroda 0:13a5d365ba16 289 }
ykuroda 0:13a5d365ba16 290
ykuroda 0:13a5d365ba16 291 /** \brief Returns the maximum number of iterations. */
ykuroda 0:13a5d365ba16 292 Index getMaxIterations()
ykuroda 0:13a5d365ba16 293 {
ykuroda 0:13a5d365ba16 294 return m_realSchur.getMaxIterations();
ykuroda 0:13a5d365ba16 295 }
ykuroda 0:13a5d365ba16 296
ykuroda 0:13a5d365ba16 297 private:
ykuroda 0:13a5d365ba16 298 void doComputeEigenvectors();
ykuroda 0:13a5d365ba16 299
ykuroda 0:13a5d365ba16 300 protected:
ykuroda 0:13a5d365ba16 301
ykuroda 0:13a5d365ba16 302 static void check_template_parameters()
ykuroda 0:13a5d365ba16 303 {
ykuroda 0:13a5d365ba16 304 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
ykuroda 0:13a5d365ba16 305 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
ykuroda 0:13a5d365ba16 306 }
ykuroda 0:13a5d365ba16 307
ykuroda 0:13a5d365ba16 308 MatrixType m_eivec;
ykuroda 0:13a5d365ba16 309 EigenvalueType m_eivalues;
ykuroda 0:13a5d365ba16 310 bool m_isInitialized;
ykuroda 0:13a5d365ba16 311 bool m_eigenvectorsOk;
ykuroda 0:13a5d365ba16 312 RealSchur<MatrixType> m_realSchur;
ykuroda 0:13a5d365ba16 313 MatrixType m_matT;
ykuroda 0:13a5d365ba16 314
ykuroda 0:13a5d365ba16 315 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
ykuroda 0:13a5d365ba16 316 ColumnVectorType m_tmp;
ykuroda 0:13a5d365ba16 317 };
ykuroda 0:13a5d365ba16 318
ykuroda 0:13a5d365ba16 319 template<typename MatrixType>
ykuroda 0:13a5d365ba16 320 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
ykuroda 0:13a5d365ba16 321 {
ykuroda 0:13a5d365ba16 322 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
ykuroda 0:13a5d365ba16 323 Index n = m_eivalues.rows();
ykuroda 0:13a5d365ba16 324 MatrixType matD = MatrixType::Zero(n,n);
ykuroda 0:13a5d365ba16 325 for (Index i=0; i<n; ++i)
ykuroda 0:13a5d365ba16 326 {
ykuroda 0:13a5d365ba16 327 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i))))
ykuroda 0:13a5d365ba16 328 matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
ykuroda 0:13a5d365ba16 329 else
ykuroda 0:13a5d365ba16 330 {
ykuroda 0:13a5d365ba16 331 matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
ykuroda 0:13a5d365ba16 332 -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i));
ykuroda 0:13a5d365ba16 333 ++i;
ykuroda 0:13a5d365ba16 334 }
ykuroda 0:13a5d365ba16 335 }
ykuroda 0:13a5d365ba16 336 return matD;
ykuroda 0:13a5d365ba16 337 }
ykuroda 0:13a5d365ba16 338
ykuroda 0:13a5d365ba16 339 template<typename MatrixType>
ykuroda 0:13a5d365ba16 340 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
ykuroda 0:13a5d365ba16 341 {
ykuroda 0:13a5d365ba16 342 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
ykuroda 0:13a5d365ba16 343 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
ykuroda 0:13a5d365ba16 344 Index n = m_eivec.cols();
ykuroda 0:13a5d365ba16 345 EigenvectorsType matV(n,n);
ykuroda 0:13a5d365ba16 346 for (Index j=0; j<n; ++j)
ykuroda 0:13a5d365ba16 347 {
ykuroda 0:13a5d365ba16 348 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n)
ykuroda 0:13a5d365ba16 349 {
ykuroda 0:13a5d365ba16 350 // we have a real eigen value
ykuroda 0:13a5d365ba16 351 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
ykuroda 0:13a5d365ba16 352 matV.col(j).normalize();
ykuroda 0:13a5d365ba16 353 }
ykuroda 0:13a5d365ba16 354 else
ykuroda 0:13a5d365ba16 355 {
ykuroda 0:13a5d365ba16 356 // we have a pair of complex eigen values
ykuroda 0:13a5d365ba16 357 for (Index i=0; i<n; ++i)
ykuroda 0:13a5d365ba16 358 {
ykuroda 0:13a5d365ba16 359 matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
ykuroda 0:13a5d365ba16 360 matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
ykuroda 0:13a5d365ba16 361 }
ykuroda 0:13a5d365ba16 362 matV.col(j).normalize();
ykuroda 0:13a5d365ba16 363 matV.col(j+1).normalize();
ykuroda 0:13a5d365ba16 364 ++j;
ykuroda 0:13a5d365ba16 365 }
ykuroda 0:13a5d365ba16 366 }
ykuroda 0:13a5d365ba16 367 return matV;
ykuroda 0:13a5d365ba16 368 }
ykuroda 0:13a5d365ba16 369
ykuroda 0:13a5d365ba16 370 template<typename MatrixType>
ykuroda 0:13a5d365ba16 371 EigenSolver<MatrixType>&
ykuroda 0:13a5d365ba16 372 EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
ykuroda 0:13a5d365ba16 373 {
ykuroda 0:13a5d365ba16 374 check_template_parameters();
ykuroda 0:13a5d365ba16 375
ykuroda 0:13a5d365ba16 376 using std::sqrt;
ykuroda 0:13a5d365ba16 377 using std::abs;
ykuroda 0:13a5d365ba16 378 eigen_assert(matrix.cols() == matrix.rows());
ykuroda 0:13a5d365ba16 379
ykuroda 0:13a5d365ba16 380 // Reduce to real Schur form.
ykuroda 0:13a5d365ba16 381 m_realSchur.compute(matrix, computeEigenvectors);
ykuroda 0:13a5d365ba16 382
ykuroda 0:13a5d365ba16 383 if (m_realSchur.info() == Success)
ykuroda 0:13a5d365ba16 384 {
ykuroda 0:13a5d365ba16 385 m_matT = m_realSchur.matrixT();
ykuroda 0:13a5d365ba16 386 if (computeEigenvectors)
ykuroda 0:13a5d365ba16 387 m_eivec = m_realSchur.matrixU();
ykuroda 0:13a5d365ba16 388
ykuroda 0:13a5d365ba16 389 // Compute eigenvalues from matT
ykuroda 0:13a5d365ba16 390 m_eivalues.resize(matrix.cols());
ykuroda 0:13a5d365ba16 391 Index i = 0;
ykuroda 0:13a5d365ba16 392 while (i < matrix.cols())
ykuroda 0:13a5d365ba16 393 {
ykuroda 0:13a5d365ba16 394 if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
ykuroda 0:13a5d365ba16 395 {
ykuroda 0:13a5d365ba16 396 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
ykuroda 0:13a5d365ba16 397 ++i;
ykuroda 0:13a5d365ba16 398 }
ykuroda 0:13a5d365ba16 399 else
ykuroda 0:13a5d365ba16 400 {
ykuroda 0:13a5d365ba16 401 Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
ykuroda 0:13a5d365ba16 402 Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
ykuroda 0:13a5d365ba16 403 m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
ykuroda 0:13a5d365ba16 404 m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
ykuroda 0:13a5d365ba16 405 i += 2;
ykuroda 0:13a5d365ba16 406 }
ykuroda 0:13a5d365ba16 407 }
ykuroda 0:13a5d365ba16 408
ykuroda 0:13a5d365ba16 409 // Compute eigenvectors.
ykuroda 0:13a5d365ba16 410 if (computeEigenvectors)
ykuroda 0:13a5d365ba16 411 doComputeEigenvectors();
ykuroda 0:13a5d365ba16 412 }
ykuroda 0:13a5d365ba16 413
ykuroda 0:13a5d365ba16 414 m_isInitialized = true;
ykuroda 0:13a5d365ba16 415 m_eigenvectorsOk = computeEigenvectors;
ykuroda 0:13a5d365ba16 416
ykuroda 0:13a5d365ba16 417 return *this;
ykuroda 0:13a5d365ba16 418 }
ykuroda 0:13a5d365ba16 419
ykuroda 0:13a5d365ba16 420 // Complex scalar division.
ykuroda 0:13a5d365ba16 421 template<typename Scalar>
ykuroda 0:13a5d365ba16 422 std::complex<Scalar> cdiv(const Scalar& xr, const Scalar& xi, const Scalar& yr, const Scalar& yi)
ykuroda 0:13a5d365ba16 423 {
ykuroda 0:13a5d365ba16 424 using std::abs;
ykuroda 0:13a5d365ba16 425 Scalar r,d;
ykuroda 0:13a5d365ba16 426 if (abs(yr) > abs(yi))
ykuroda 0:13a5d365ba16 427 {
ykuroda 0:13a5d365ba16 428 r = yi/yr;
ykuroda 0:13a5d365ba16 429 d = yr + r*yi;
ykuroda 0:13a5d365ba16 430 return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
ykuroda 0:13a5d365ba16 431 }
ykuroda 0:13a5d365ba16 432 else
ykuroda 0:13a5d365ba16 433 {
ykuroda 0:13a5d365ba16 434 r = yr/yi;
ykuroda 0:13a5d365ba16 435 d = yi + r*yr;
ykuroda 0:13a5d365ba16 436 return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
ykuroda 0:13a5d365ba16 437 }
ykuroda 0:13a5d365ba16 438 }
ykuroda 0:13a5d365ba16 439
ykuroda 0:13a5d365ba16 440
ykuroda 0:13a5d365ba16 441 template<typename MatrixType>
ykuroda 0:13a5d365ba16 442 void EigenSolver<MatrixType>::doComputeEigenvectors()
ykuroda 0:13a5d365ba16 443 {
ykuroda 0:13a5d365ba16 444 using std::abs;
ykuroda 0:13a5d365ba16 445 const Index size = m_eivec.cols();
ykuroda 0:13a5d365ba16 446 const Scalar eps = NumTraits<Scalar>::epsilon();
ykuroda 0:13a5d365ba16 447
ykuroda 0:13a5d365ba16 448 // inefficient! this is already computed in RealSchur
ykuroda 0:13a5d365ba16 449 Scalar norm(0);
ykuroda 0:13a5d365ba16 450 for (Index j = 0; j < size; ++j)
ykuroda 0:13a5d365ba16 451 {
ykuroda 0:13a5d365ba16 452 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
ykuroda 0:13a5d365ba16 453 }
ykuroda 0:13a5d365ba16 454
ykuroda 0:13a5d365ba16 455 // Backsubstitute to find vectors of upper triangular form
ykuroda 0:13a5d365ba16 456 if (norm == 0.0)
ykuroda 0:13a5d365ba16 457 {
ykuroda 0:13a5d365ba16 458 return;
ykuroda 0:13a5d365ba16 459 }
ykuroda 0:13a5d365ba16 460
ykuroda 0:13a5d365ba16 461 for (Index n = size-1; n >= 0; n--)
ykuroda 0:13a5d365ba16 462 {
ykuroda 0:13a5d365ba16 463 Scalar p = m_eivalues.coeff(n).real();
ykuroda 0:13a5d365ba16 464 Scalar q = m_eivalues.coeff(n).imag();
ykuroda 0:13a5d365ba16 465
ykuroda 0:13a5d365ba16 466 // Scalar vector
ykuroda 0:13a5d365ba16 467 if (q == Scalar(0))
ykuroda 0:13a5d365ba16 468 {
ykuroda 0:13a5d365ba16 469 Scalar lastr(0), lastw(0);
ykuroda 0:13a5d365ba16 470 Index l = n;
ykuroda 0:13a5d365ba16 471
ykuroda 0:13a5d365ba16 472 m_matT.coeffRef(n,n) = 1.0;
ykuroda 0:13a5d365ba16 473 for (Index i = n-1; i >= 0; i--)
ykuroda 0:13a5d365ba16 474 {
ykuroda 0:13a5d365ba16 475 Scalar w = m_matT.coeff(i,i) - p;
ykuroda 0:13a5d365ba16 476 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
ykuroda 0:13a5d365ba16 477
ykuroda 0:13a5d365ba16 478 if (m_eivalues.coeff(i).imag() < 0.0)
ykuroda 0:13a5d365ba16 479 {
ykuroda 0:13a5d365ba16 480 lastw = w;
ykuroda 0:13a5d365ba16 481 lastr = r;
ykuroda 0:13a5d365ba16 482 }
ykuroda 0:13a5d365ba16 483 else
ykuroda 0:13a5d365ba16 484 {
ykuroda 0:13a5d365ba16 485 l = i;
ykuroda 0:13a5d365ba16 486 if (m_eivalues.coeff(i).imag() == 0.0)
ykuroda 0:13a5d365ba16 487 {
ykuroda 0:13a5d365ba16 488 if (w != 0.0)
ykuroda 0:13a5d365ba16 489 m_matT.coeffRef(i,n) = -r / w;
ykuroda 0:13a5d365ba16 490 else
ykuroda 0:13a5d365ba16 491 m_matT.coeffRef(i,n) = -r / (eps * norm);
ykuroda 0:13a5d365ba16 492 }
ykuroda 0:13a5d365ba16 493 else // Solve real equations
ykuroda 0:13a5d365ba16 494 {
ykuroda 0:13a5d365ba16 495 Scalar x = m_matT.coeff(i,i+1);
ykuroda 0:13a5d365ba16 496 Scalar y = m_matT.coeff(i+1,i);
ykuroda 0:13a5d365ba16 497 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
ykuroda 0:13a5d365ba16 498 Scalar t = (x * lastr - lastw * r) / denom;
ykuroda 0:13a5d365ba16 499 m_matT.coeffRef(i,n) = t;
ykuroda 0:13a5d365ba16 500 if (abs(x) > abs(lastw))
ykuroda 0:13a5d365ba16 501 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
ykuroda 0:13a5d365ba16 502 else
ykuroda 0:13a5d365ba16 503 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
ykuroda 0:13a5d365ba16 504 }
ykuroda 0:13a5d365ba16 505
ykuroda 0:13a5d365ba16 506 // Overflow control
ykuroda 0:13a5d365ba16 507 Scalar t = abs(m_matT.coeff(i,n));
ykuroda 0:13a5d365ba16 508 if ((eps * t) * t > Scalar(1))
ykuroda 0:13a5d365ba16 509 m_matT.col(n).tail(size-i) /= t;
ykuroda 0:13a5d365ba16 510 }
ykuroda 0:13a5d365ba16 511 }
ykuroda 0:13a5d365ba16 512 }
ykuroda 0:13a5d365ba16 513 else if (q < Scalar(0) && n > 0) // Complex vector
ykuroda 0:13a5d365ba16 514 {
ykuroda 0:13a5d365ba16 515 Scalar lastra(0), lastsa(0), lastw(0);
ykuroda 0:13a5d365ba16 516 Index l = n-1;
ykuroda 0:13a5d365ba16 517
ykuroda 0:13a5d365ba16 518 // Last vector component imaginary so matrix is triangular
ykuroda 0:13a5d365ba16 519 if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
ykuroda 0:13a5d365ba16 520 {
ykuroda 0:13a5d365ba16 521 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
ykuroda 0:13a5d365ba16 522 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
ykuroda 0:13a5d365ba16 523 }
ykuroda 0:13a5d365ba16 524 else
ykuroda 0:13a5d365ba16 525 {
ykuroda 0:13a5d365ba16 526 std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
ykuroda 0:13a5d365ba16 527 m_matT.coeffRef(n-1,n-1) = numext::real(cc);
ykuroda 0:13a5d365ba16 528 m_matT.coeffRef(n-1,n) = numext::imag(cc);
ykuroda 0:13a5d365ba16 529 }
ykuroda 0:13a5d365ba16 530 m_matT.coeffRef(n,n-1) = 0.0;
ykuroda 0:13a5d365ba16 531 m_matT.coeffRef(n,n) = 1.0;
ykuroda 0:13a5d365ba16 532 for (Index i = n-2; i >= 0; i--)
ykuroda 0:13a5d365ba16 533 {
ykuroda 0:13a5d365ba16 534 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
ykuroda 0:13a5d365ba16 535 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
ykuroda 0:13a5d365ba16 536 Scalar w = m_matT.coeff(i,i) - p;
ykuroda 0:13a5d365ba16 537
ykuroda 0:13a5d365ba16 538 if (m_eivalues.coeff(i).imag() < 0.0)
ykuroda 0:13a5d365ba16 539 {
ykuroda 0:13a5d365ba16 540 lastw = w;
ykuroda 0:13a5d365ba16 541 lastra = ra;
ykuroda 0:13a5d365ba16 542 lastsa = sa;
ykuroda 0:13a5d365ba16 543 }
ykuroda 0:13a5d365ba16 544 else
ykuroda 0:13a5d365ba16 545 {
ykuroda 0:13a5d365ba16 546 l = i;
ykuroda 0:13a5d365ba16 547 if (m_eivalues.coeff(i).imag() == RealScalar(0))
ykuroda 0:13a5d365ba16 548 {
ykuroda 0:13a5d365ba16 549 std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
ykuroda 0:13a5d365ba16 550 m_matT.coeffRef(i,n-1) = numext::real(cc);
ykuroda 0:13a5d365ba16 551 m_matT.coeffRef(i,n) = numext::imag(cc);
ykuroda 0:13a5d365ba16 552 }
ykuroda 0:13a5d365ba16 553 else
ykuroda 0:13a5d365ba16 554 {
ykuroda 0:13a5d365ba16 555 // Solve complex equations
ykuroda 0:13a5d365ba16 556 Scalar x = m_matT.coeff(i,i+1);
ykuroda 0:13a5d365ba16 557 Scalar y = m_matT.coeff(i+1,i);
ykuroda 0:13a5d365ba16 558 Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
ykuroda 0:13a5d365ba16 559 Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
ykuroda 0:13a5d365ba16 560 if ((vr == 0.0) && (vi == 0.0))
ykuroda 0:13a5d365ba16 561 vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
ykuroda 0:13a5d365ba16 562
ykuroda 0:13a5d365ba16 563 std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
ykuroda 0:13a5d365ba16 564 m_matT.coeffRef(i,n-1) = numext::real(cc);
ykuroda 0:13a5d365ba16 565 m_matT.coeffRef(i,n) = numext::imag(cc);
ykuroda 0:13a5d365ba16 566 if (abs(x) > (abs(lastw) + abs(q)))
ykuroda 0:13a5d365ba16 567 {
ykuroda 0:13a5d365ba16 568 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
ykuroda 0:13a5d365ba16 569 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
ykuroda 0:13a5d365ba16 570 }
ykuroda 0:13a5d365ba16 571 else
ykuroda 0:13a5d365ba16 572 {
ykuroda 0:13a5d365ba16 573 cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
ykuroda 0:13a5d365ba16 574 m_matT.coeffRef(i+1,n-1) = numext::real(cc);
ykuroda 0:13a5d365ba16 575 m_matT.coeffRef(i+1,n) = numext::imag(cc);
ykuroda 0:13a5d365ba16 576 }
ykuroda 0:13a5d365ba16 577 }
ykuroda 0:13a5d365ba16 578
ykuroda 0:13a5d365ba16 579 // Overflow control
ykuroda 0:13a5d365ba16 580 using std::max;
ykuroda 0:13a5d365ba16 581 Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
ykuroda 0:13a5d365ba16 582 if ((eps * t) * t > Scalar(1))
ykuroda 0:13a5d365ba16 583 m_matT.block(i, n-1, size-i, 2) /= t;
ykuroda 0:13a5d365ba16 584
ykuroda 0:13a5d365ba16 585 }
ykuroda 0:13a5d365ba16 586 }
ykuroda 0:13a5d365ba16 587
ykuroda 0:13a5d365ba16 588 // We handled a pair of complex conjugate eigenvalues, so need to skip them both
ykuroda 0:13a5d365ba16 589 n--;
ykuroda 0:13a5d365ba16 590 }
ykuroda 0:13a5d365ba16 591 else
ykuroda 0:13a5d365ba16 592 {
ykuroda 0:13a5d365ba16 593 eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
ykuroda 0:13a5d365ba16 594 }
ykuroda 0:13a5d365ba16 595 }
ykuroda 0:13a5d365ba16 596
ykuroda 0:13a5d365ba16 597 // Back transformation to get eigenvectors of original matrix
ykuroda 0:13a5d365ba16 598 for (Index j = size-1; j >= 0; j--)
ykuroda 0:13a5d365ba16 599 {
ykuroda 0:13a5d365ba16 600 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
ykuroda 0:13a5d365ba16 601 m_eivec.col(j) = m_tmp;
ykuroda 0:13a5d365ba16 602 }
ykuroda 0:13a5d365ba16 603 }
ykuroda 0:13a5d365ba16 604
ykuroda 0:13a5d365ba16 605 } // end namespace Eigen
ykuroda 0:13a5d365ba16 606
ykuroda 0:13a5d365ba16 607 #endif // EIGEN_EIGENSOLVER_H