Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat (limited to 'extern/Eigen2/Eigen/src/QR/EigenSolver.h')
-rw-r--r--extern/Eigen2/Eigen/src/QR/EigenSolver.h722
1 files changed, 0 insertions, 722 deletions
diff --git a/extern/Eigen2/Eigen/src/QR/EigenSolver.h b/extern/Eigen2/Eigen/src/QR/EigenSolver.h
deleted file mode 100644
index 70f21cebcdb..00000000000
--- a/extern/Eigen2/Eigen/src/QR/EigenSolver.h
+++ /dev/null
@@ -1,722 +0,0 @@
-// This file is part of Eigen, a lightweight C++ template library
-// for linear algebra.
-//
-// Copyright (C) 2008 Gael Guennebaud <g.gael@free.fr>
-//
-// Eigen is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 3 of the License, or (at your option) any later version.
-//
-// Alternatively, you can redistribute it and/or
-// modify it under the terms of the GNU General Public License as
-// published by the Free Software Foundation; either version 2 of
-// the License, or (at your option) any later version.
-//
-// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
-// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
-// GNU General Public License for more details.
-//
-// You should have received a copy of the GNU Lesser General Public
-// License and a copy of the GNU General Public License along with
-// Eigen. If not, see <http://www.gnu.org/licenses/>.
-
-#ifndef EIGEN_EIGENSOLVER_H
-#define EIGEN_EIGENSOLVER_H
-
-/** \ingroup QR_Module
- * \nonstableyet
- *
- * \class EigenSolver
- *
- * \brief Eigen values/vectors solver for non selfadjoint matrices
- *
- * \param MatrixType the type of the matrix of which we are computing the eigen decomposition
- *
- * Currently it only support real matrices.
- *
- * \note this code was adapted from JAMA (public domain)
- *
- * \sa MatrixBase::eigenvalues(), SelfAdjointEigenSolver
- */
-template<typename _MatrixType> class EigenSolver
-{
- public:
-
- typedef _MatrixType MatrixType;
- typedef typename MatrixType::Scalar Scalar;
- typedef typename NumTraits<Scalar>::Real RealScalar;
- typedef std::complex<RealScalar> Complex;
- typedef Matrix<Complex, MatrixType::ColsAtCompileTime, 1> EigenvalueType;
- typedef Matrix<Complex, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime> EigenvectorType;
- typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
- typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
-
- /**
- * \brief Default Constructor.
- *
- * The default constructor is useful in cases in which the user intends to
- * perform decompositions via EigenSolver::compute(const MatrixType&).
- */
- EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false) {}
-
- EigenSolver(const MatrixType& matrix)
- : m_eivec(matrix.rows(), matrix.cols()),
- m_eivalues(matrix.cols()),
- m_isInitialized(false)
- {
- compute(matrix);
- }
-
-
- EigenvectorType eigenvectors(void) const;
-
- /** \returns a real matrix V of pseudo eigenvectors.
- *
- * Let D be the block diagonal matrix with the real eigenvalues in 1x1 blocks,
- * and any complex values u+iv in 2x2 blocks [u v ; -v u]. Then, the matrices D
- * and V satisfy A*V = V*D.
- *
- * More precisely, if the diagonal matrix of the eigen values is:\n
- * \f$
- * \left[ \begin{array}{cccccc}
- * u+iv & & & & & \\
- * & u-iv & & & & \\
- * & & a+ib & & & \\
- * & & & a-ib & & \\
- * & & & & x & \\
- * & & & & & y \\
- * \end{array} \right]
- * \f$ \n
- * then, we have:\n
- * \f$
- * D =\left[ \begin{array}{cccccc}
- * u & v & & & & \\
- * -v & u & & & & \\
- * & & a & b & & \\
- * & & -b & a & & \\
- * & & & & x & \\
- * & & & & & y \\
- * \end{array} \right]
- * \f$
- *
- * \sa pseudoEigenvalueMatrix()
- */
- const MatrixType& pseudoEigenvectors() const
- {
- ei_assert(m_isInitialized && "EigenSolver is not initialized.");
- return m_eivec;
- }
-
- MatrixType pseudoEigenvalueMatrix() const;
-
- /** \returns the eigenvalues as a column vector */
- EigenvalueType eigenvalues() const
- {
- ei_assert(m_isInitialized && "EigenSolver is not initialized.");
- return m_eivalues;
- }
-
- void compute(const MatrixType& matrix);
-
- private:
-
- void orthes(MatrixType& matH, RealVectorType& ort);
- void hqr2(MatrixType& matH);
-
- protected:
- MatrixType m_eivec;
- EigenvalueType m_eivalues;
- bool m_isInitialized;
-};
-
-/** \returns the real block diagonal matrix D of the eigenvalues.
- *
- * See pseudoEigenvectors() for the details.
- */
-template<typename MatrixType>
-MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
-{
- ei_assert(m_isInitialized && "EigenSolver is not initialized.");
- int n = m_eivec.cols();
- MatrixType matD = MatrixType::Zero(n,n);
- for (int i=0; i<n; ++i)
- {
- if (ei_isMuchSmallerThan(ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i))))
- matD.coeffRef(i,i) = ei_real(m_eivalues.coeff(i));
- else
- {
- matD.template block<2,2>(i,i) << ei_real(m_eivalues.coeff(i)), ei_imag(m_eivalues.coeff(i)),
- -ei_imag(m_eivalues.coeff(i)), ei_real(m_eivalues.coeff(i));
- ++i;
- }
- }
- return matD;
-}
-
-/** \returns the normalized complex eigenvectors as a matrix of column vectors.
- *
- * \sa eigenvalues(), pseudoEigenvectors()
- */
-template<typename MatrixType>
-typename EigenSolver<MatrixType>::EigenvectorType EigenSolver<MatrixType>::eigenvectors(void) const
-{
- ei_assert(m_isInitialized && "EigenSolver is not initialized.");
- int n = m_eivec.cols();
- EigenvectorType matV(n,n);
- for (int j=0; j<n; ++j)
- {
- if (ei_isMuchSmallerThan(ei_abs(ei_imag(m_eivalues.coeff(j))), ei_abs(ei_real(m_eivalues.coeff(j)))))
- {
- // we have a real eigen value
- matV.col(j) = m_eivec.col(j).template cast<Complex>();
- }
- else
- {
- // we have a pair of complex eigen values
- for (int i=0; i<n; ++i)
- {
- matV.coeffRef(i,j) = Complex(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
- matV.coeffRef(i,j+1) = Complex(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
- }
- matV.col(j).normalize();
- matV.col(j+1).normalize();
- ++j;
- }
- }
- return matV;
-}
-
-template<typename MatrixType>
-void EigenSolver<MatrixType>::compute(const MatrixType& matrix)
-{
- assert(matrix.cols() == matrix.rows());
- int n = matrix.cols();
- m_eivalues.resize(n,1);
-
- MatrixType matH = matrix;
- RealVectorType ort(n);
-
- // Reduce to Hessenberg form.
- orthes(matH, ort);
-
- // Reduce Hessenberg to real Schur form.
- hqr2(matH);
-
- m_isInitialized = true;
-}
-
-// Nonsymmetric reduction to Hessenberg form.
-template<typename MatrixType>
-void EigenSolver<MatrixType>::orthes(MatrixType& matH, RealVectorType& ort)
-{
- // This is derived from the Algol procedures orthes and ortran,
- // by Martin and Wilkinson, Handbook for Auto. Comp.,
- // Vol.ii-Linear Algebra, and the corresponding
- // Fortran subroutines in EISPACK.
-
- int n = m_eivec.cols();
- int low = 0;
- int high = n-1;
-
- for (int m = low+1; m <= high-1; ++m)
- {
- // Scale column.
- RealScalar scale = matH.block(m, m-1, high-m+1, 1).cwise().abs().sum();
- if (scale != 0.0)
- {
- // Compute Householder transformation.
- RealScalar h = 0.0;
- // FIXME could be rewritten, but this one looks better wrt cache
- for (int i = high; i >= m; i--)
- {
- ort.coeffRef(i) = matH.coeff(i,m-1)/scale;
- h += ort.coeff(i) * ort.coeff(i);
- }
- RealScalar g = ei_sqrt(h);
- if (ort.coeff(m) > 0)
- g = -g;
- h = h - ort.coeff(m) * g;
- ort.coeffRef(m) = ort.coeff(m) - g;
-
- // Apply Householder similarity transformation
- // H = (I-u*u'/h)*H*(I-u*u')/h)
- int bSize = high-m+1;
- matH.block(m, m, bSize, n-m) -= ((ort.segment(m, bSize)/h)
- * (ort.segment(m, bSize).transpose() * matH.block(m, m, bSize, n-m)).lazy()).lazy();
-
- matH.block(0, m, high+1, bSize) -= ((matH.block(0, m, high+1, bSize) * ort.segment(m, bSize)).lazy()
- * (ort.segment(m, bSize)/h).transpose()).lazy();
-
- ort.coeffRef(m) = scale*ort.coeff(m);
- matH.coeffRef(m,m-1) = scale*g;
- }
- }
-
- // Accumulate transformations (Algol's ortran).
- m_eivec.setIdentity();
-
- for (int m = high-1; m >= low+1; m--)
- {
- if (matH.coeff(m,m-1) != 0.0)
- {
- ort.segment(m+1, high-m) = matH.col(m-1).segment(m+1, high-m);
-
- int bSize = high-m+1;
- m_eivec.block(m, m, bSize, bSize) += ( (ort.segment(m, bSize) / (matH.coeff(m,m-1) * ort.coeff(m) ) )
- * (ort.segment(m, bSize).transpose() * m_eivec.block(m, m, bSize, bSize)).lazy());
- }
- }
-}
-
-// Complex scalar division.
-template<typename Scalar>
-std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
-{
- Scalar r,d;
- if (ei_abs(yr) > ei_abs(yi))
- {
- r = yi/yr;
- d = yr + r*yi;
- return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
- }
- else
- {
- r = yr/yi;
- d = yi + r*yr;
- return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
- }
-}
-
-
-// Nonsymmetric reduction from Hessenberg to real Schur form.
-template<typename MatrixType>
-void EigenSolver<MatrixType>::hqr2(MatrixType& matH)
-{
- // This is derived from the Algol procedure hqr2,
- // by Martin and Wilkinson, Handbook for Auto. Comp.,
- // Vol.ii-Linear Algebra, and the corresponding
- // Fortran subroutine in EISPACK.
-
- // Initialize
- int nn = m_eivec.cols();
- int n = nn-1;
- int low = 0;
- int high = nn-1;
- Scalar eps = ei_pow(Scalar(2),ei_is_same_type<Scalar,float>::ret ? Scalar(-23) : Scalar(-52));
- Scalar exshift = 0.0;
- Scalar p=0,q=0,r=0,s=0,z=0,t,w,x,y;
-
- // Store roots isolated by balanc and compute matrix norm
- // FIXME to be efficient the following would requires a triangular reduxion code
- // Scalar norm = matH.upper().cwise().abs().sum() + matH.corner(BottomLeft,n,n).diagonal().cwise().abs().sum();
- Scalar norm = 0.0;
- for (int j = 0; j < nn; ++j)
- {
- // FIXME what's the purpose of the following since the condition is always false
- if ((j < low) || (j > high))
- {
- m_eivalues.coeffRef(j) = Complex(matH.coeff(j,j), 0.0);
- }
- norm += matH.row(j).segment(std::max(j-1,0), nn-std::max(j-1,0)).cwise().abs().sum();
- }
-
- // Outer loop over eigenvalue index
- int iter = 0;
- while (n >= low)
- {
- // Look for single small sub-diagonal element
- int l = n;
- while (l > low)
- {
- s = ei_abs(matH.coeff(l-1,l-1)) + ei_abs(matH.coeff(l,l));
- if (s == 0.0)
- s = norm;
- if (ei_abs(matH.coeff(l,l-1)) < eps * s)
- break;
- l--;
- }
-
- // Check for convergence
- // One root found
- if (l == n)
- {
- matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
- m_eivalues.coeffRef(n) = Complex(matH.coeff(n,n), 0.0);
- n--;
- iter = 0;
- }
- else if (l == n-1) // Two roots found
- {
- w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
- p = (matH.coeff(n-1,n-1) - matH.coeff(n,n)) * Scalar(0.5);
- q = p * p + w;
- z = ei_sqrt(ei_abs(q));
- matH.coeffRef(n,n) = matH.coeff(n,n) + exshift;
- matH.coeffRef(n-1,n-1) = matH.coeff(n-1,n-1) + exshift;
- x = matH.coeff(n,n);
-
- // Scalar pair
- if (q >= 0)
- {
- if (p >= 0)
- z = p + z;
- else
- z = p - z;
-
- m_eivalues.coeffRef(n-1) = Complex(x + z, 0.0);
- m_eivalues.coeffRef(n) = Complex(z!=0.0 ? x - w / z : m_eivalues.coeff(n-1).real(), 0.0);
-
- x = matH.coeff(n,n-1);
- s = ei_abs(x) + ei_abs(z);
- p = x / s;
- q = z / s;
- r = ei_sqrt(p * p+q * q);
- p = p / r;
- q = q / r;
-
- // Row modification
- for (int j = n-1; j < nn; ++j)
- {
- z = matH.coeff(n-1,j);
- matH.coeffRef(n-1,j) = q * z + p * matH.coeff(n,j);
- matH.coeffRef(n,j) = q * matH.coeff(n,j) - p * z;
- }
-
- // Column modification
- for (int i = 0; i <= n; ++i)
- {
- z = matH.coeff(i,n-1);
- matH.coeffRef(i,n-1) = q * z + p * matH.coeff(i,n);
- matH.coeffRef(i,n) = q * matH.coeff(i,n) - p * z;
- }
-
- // Accumulate transformations
- for (int i = low; i <= high; ++i)
- {
- z = m_eivec.coeff(i,n-1);
- m_eivec.coeffRef(i,n-1) = q * z + p * m_eivec.coeff(i,n);
- m_eivec.coeffRef(i,n) = q * m_eivec.coeff(i,n) - p * z;
- }
- }
- else // Complex pair
- {
- m_eivalues.coeffRef(n-1) = Complex(x + p, z);
- m_eivalues.coeffRef(n) = Complex(x + p, -z);
- }
- n = n - 2;
- iter = 0;
- }
- else // No convergence yet
- {
- // Form shift
- x = matH.coeff(n,n);
- y = 0.0;
- w = 0.0;
- if (l < n)
- {
- y = matH.coeff(n-1,n-1);
- w = matH.coeff(n,n-1) * matH.coeff(n-1,n);
- }
-
- // Wilkinson's original ad hoc shift
- if (iter == 10)
- {
- exshift += x;
- for (int i = low; i <= n; ++i)
- matH.coeffRef(i,i) -= x;
- s = ei_abs(matH.coeff(n,n-1)) + ei_abs(matH.coeff(n-1,n-2));
- x = y = Scalar(0.75) * s;
- w = Scalar(-0.4375) * s * s;
- }
-
- // MATLAB's new ad hoc shift
- if (iter == 30)
- {
- s = Scalar((y - x) / 2.0);
- s = s * s + w;
- if (s > 0)
- {
- s = ei_sqrt(s);
- if (y < x)
- s = -s;
- s = Scalar(x - w / ((y - x) / 2.0 + s));
- for (int i = low; i <= n; ++i)
- matH.coeffRef(i,i) -= s;
- exshift += s;
- x = y = w = Scalar(0.964);
- }
- }
-
- iter = iter + 1; // (Could check iteration count here.)
-
- // Look for two consecutive small sub-diagonal elements
- int m = n-2;
- while (m >= l)
- {
- z = matH.coeff(m,m);
- r = x - z;
- s = y - z;
- p = (r * s - w) / matH.coeff(m+1,m) + matH.coeff(m,m+1);
- q = matH.coeff(m+1,m+1) - z - r - s;
- r = matH.coeff(m+2,m+1);
- s = ei_abs(p) + ei_abs(q) + ei_abs(r);
- p = p / s;
- q = q / s;
- r = r / s;
- if (m == l) {
- break;
- }
- if (ei_abs(matH.coeff(m,m-1)) * (ei_abs(q) + ei_abs(r)) <
- eps * (ei_abs(p) * (ei_abs(matH.coeff(m-1,m-1)) + ei_abs(z) +
- ei_abs(matH.coeff(m+1,m+1)))))
- {
- break;
- }
- m--;
- }
-
- for (int i = m+2; i <= n; ++i)
- {
- matH.coeffRef(i,i-2) = 0.0;
- if (i > m+2)
- matH.coeffRef(i,i-3) = 0.0;
- }
-
- // Double QR step involving rows l:n and columns m:n
- for (int k = m; k <= n-1; ++k)
- {
- int notlast = (k != n-1);
- if (k != m) {
- p = matH.coeff(k,k-1);
- q = matH.coeff(k+1,k-1);
- r = notlast ? matH.coeff(k+2,k-1) : Scalar(0);
- x = ei_abs(p) + ei_abs(q) + ei_abs(r);
- if (x != 0.0)
- {
- p = p / x;
- q = q / x;
- r = r / x;
- }
- }
-
- if (x == 0.0)
- break;
-
- s = ei_sqrt(p * p + q * q + r * r);
-
- if (p < 0)
- s = -s;
-
- if (s != 0)
- {
- if (k != m)
- matH.coeffRef(k,k-1) = -s * x;
- else if (l != m)
- matH.coeffRef(k,k-1) = -matH.coeff(k,k-1);
-
- p = p + s;
- x = p / s;
- y = q / s;
- z = r / s;
- q = q / p;
- r = r / p;
-
- // Row modification
- for (int j = k; j < nn; ++j)
- {
- p = matH.coeff(k,j) + q * matH.coeff(k+1,j);
- if (notlast)
- {
- p = p + r * matH.coeff(k+2,j);
- matH.coeffRef(k+2,j) = matH.coeff(k+2,j) - p * z;
- }
- matH.coeffRef(k,j) = matH.coeff(k,j) - p * x;
- matH.coeffRef(k+1,j) = matH.coeff(k+1,j) - p * y;
- }
-
- // Column modification
- for (int i = 0; i <= std::min(n,k+3); ++i)
- {
- p = x * matH.coeff(i,k) + y * matH.coeff(i,k+1);
- if (notlast)
- {
- p = p + z * matH.coeff(i,k+2);
- matH.coeffRef(i,k+2) = matH.coeff(i,k+2) - p * r;
- }
- matH.coeffRef(i,k) = matH.coeff(i,k) - p;
- matH.coeffRef(i,k+1) = matH.coeff(i,k+1) - p * q;
- }
-
- // Accumulate transformations
- for (int i = low; i <= high; ++i)
- {
- p = x * m_eivec.coeff(i,k) + y * m_eivec.coeff(i,k+1);
- if (notlast)
- {
- p = p + z * m_eivec.coeff(i,k+2);
- m_eivec.coeffRef(i,k+2) = m_eivec.coeff(i,k+2) - p * r;
- }
- m_eivec.coeffRef(i,k) = m_eivec.coeff(i,k) - p;
- m_eivec.coeffRef(i,k+1) = m_eivec.coeff(i,k+1) - p * q;
- }
- } // (s != 0)
- } // k loop
- } // check convergence
- } // while (n >= low)
-
- // Backsubstitute to find vectors of upper triangular form
- if (norm == 0.0)
- {
- return;
- }
-
- for (n = nn-1; n >= 0; n--)
- {
- p = m_eivalues.coeff(n).real();
- q = m_eivalues.coeff(n).imag();
-
- // Scalar vector
- if (q == 0)
- {
- int l = n;
- matH.coeffRef(n,n) = 1.0;
- for (int i = n-1; i >= 0; i--)
- {
- w = matH.coeff(i,i) - p;
- r = (matH.row(i).segment(l,n-l+1) * matH.col(n).segment(l, n-l+1))(0,0);
-
- if (m_eivalues.coeff(i).imag() < 0.0)
- {
- z = w;
- s = r;
- }
- else
- {
- l = i;
- if (m_eivalues.coeff(i).imag() == 0.0)
- {
- if (w != 0.0)
- matH.coeffRef(i,n) = -r / w;
- else
- matH.coeffRef(i,n) = -r / (eps * norm);
- }
- else // Solve real equations
- {
- x = matH.coeff(i,i+1);
- y = matH.coeff(i+1,i);
- q = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
- t = (x * s - z * r) / q;
- matH.coeffRef(i,n) = t;
- if (ei_abs(x) > ei_abs(z))
- matH.coeffRef(i+1,n) = (-r - w * t) / x;
- else
- matH.coeffRef(i+1,n) = (-s - y * t) / z;
- }
-
- // Overflow control
- t = ei_abs(matH.coeff(i,n));
- if ((eps * t) * t > 1)
- matH.col(n).end(nn-i) /= t;
- }
- }
- }
- else if (q < 0) // Complex vector
- {
- std::complex<Scalar> cc;
- int l = n-1;
-
- // Last vector component imaginary so matrix is triangular
- if (ei_abs(matH.coeff(n,n-1)) > ei_abs(matH.coeff(n-1,n)))
- {
- matH.coeffRef(n-1,n-1) = q / matH.coeff(n,n-1);
- matH.coeffRef(n-1,n) = -(matH.coeff(n,n) - p) / matH.coeff(n,n-1);
- }
- else
- {
- cc = cdiv<Scalar>(0.0,-matH.coeff(n-1,n),matH.coeff(n-1,n-1)-p,q);
- matH.coeffRef(n-1,n-1) = ei_real(cc);
- matH.coeffRef(n-1,n) = ei_imag(cc);
- }
- matH.coeffRef(n,n-1) = 0.0;
- matH.coeffRef(n,n) = 1.0;
- for (int i = n-2; i >= 0; i--)
- {
- Scalar ra,sa,vr,vi;
- ra = (matH.block(i,l, 1, n-l+1) * matH.block(l,n-1, n-l+1, 1)).lazy()(0,0);
- sa = (matH.block(i,l, 1, n-l+1) * matH.block(l,n, n-l+1, 1)).lazy()(0,0);
- w = matH.coeff(i,i) - p;
-
- if (m_eivalues.coeff(i).imag() < 0.0)
- {
- z = w;
- r = ra;
- s = sa;
- }
- else
- {
- l = i;
- if (m_eivalues.coeff(i).imag() == 0)
- {
- cc = cdiv(-ra,-sa,w,q);
- matH.coeffRef(i,n-1) = ei_real(cc);
- matH.coeffRef(i,n) = ei_imag(cc);
- }
- else
- {
- // Solve complex equations
- x = matH.coeff(i,i+1);
- y = matH.coeff(i+1,i);
- 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;
- vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
- if ((vr == 0.0) && (vi == 0.0))
- vr = eps * norm * (ei_abs(w) + ei_abs(q) + ei_abs(x) + ei_abs(y) + ei_abs(z));
-
- cc= cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
- matH.coeffRef(i,n-1) = ei_real(cc);
- matH.coeffRef(i,n) = ei_imag(cc);
- if (ei_abs(x) > (ei_abs(z) + ei_abs(q)))
- {
- matH.coeffRef(i+1,n-1) = (-ra - w * matH.coeff(i,n-1) + q * matH.coeff(i,n)) / x;
- matH.coeffRef(i+1,n) = (-sa - w * matH.coeff(i,n) - q * matH.coeff(i,n-1)) / x;
- }
- else
- {
- cc = cdiv(-r-y*matH.coeff(i,n-1),-s-y*matH.coeff(i,n),z,q);
- matH.coeffRef(i+1,n-1) = ei_real(cc);
- matH.coeffRef(i+1,n) = ei_imag(cc);
- }
- }
-
- // Overflow control
- t = std::max(ei_abs(matH.coeff(i,n-1)),ei_abs(matH.coeff(i,n)));
- if ((eps * t) * t > 1)
- matH.block(i, n-1, nn-i, 2) /= t;
-
- }
- }
- }
- }
-
- // Vectors of isolated roots
- for (int i = 0; i < nn; ++i)
- {
- // FIXME again what's the purpose of this test ?
- // in this algo low==0 and high==nn-1 !!
- if (i < low || i > high)
- {
- m_eivec.row(i).end(nn-i) = matH.row(i).end(nn-i);
- }
- }
-
- // Back transformation to get eigenvectors of original matrix
- int bRows = high-low+1;
- for (int j = nn-1; j >= low; j--)
- {
- int bSize = std::min(j,high)-low+1;
- m_eivec.col(j).segment(low, bRows) = (m_eivec.block(low, low, bRows, bSize) * matH.col(j).segment(low, bSize));
- }
-}
-
-#endif // EIGEN_EIGENSOLVER_H