<commit_before>// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008 <NAME> <<EMAIL>>
//
// 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;

    EigenSolver(const MatrixType& matrix)
      : m_eivec(matrix.rows(), matrix.cols()),
        m_eivalues(matrix.cols())
    {
      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 { return m_eivec; }

    MatrixType pseudoEigenvalueMatrix() const;

    /** \returns the eigenvalues as a column vector */
    EigenvalueType eigenvalues() const { 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;
};

/** \returns the real block diagonal matrix D of the eigenvalues.
  *
  * See pseudoEigenvectors() for the details.
  */
template<typename MatrixType>
MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
{
  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
{
  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);
    }
    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);
}

// 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 <NAME> 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 <NAME> 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 = pow(2.0,-52.0);
  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).real() = matH.coeff(j,j);
      m_eivalues.coeffRef(j).imag() = 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).real() = matH.coeff(n,n);
      m_eivalues.coeffRef(n).imag() = 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)) / 2.0;
      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).real() = x + z;
        m_eivalues.coeffRef(n).real() = m_eivalues.coeff(n-1).real();
        if (z != 0.0)
          m_eivalues.coeffRef(n).real() = x - w / z;

        m_eivalues.coeffRef(n-1).imag() = 0.0;
        m_eivalues.coeffRef(n).imag() = 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).real() = x + p;
        m_eivalues.coeffRef(n).real() = x + p;
        m_eivalues.coeffRef(n-1).imag() = z;
        m_eivalues.coeffRef(n).imag() = -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 = 0.75 * s;
        w = -0.4375 * s * s;
      }

      // MATLAB's new ad hoc shift
      if (iter == 30)
      {<commit_msg>fix non standard non-const std::complex::real/imag issue<commit_after>// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// Copyright (C) 2008 <NAME> <<EMAIL>>
//
// 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;

    EigenSolver(const MatrixType& matrix)
      : m_eivec(matrix.rows(), matrix.cols()),
        m_eivalues(matrix.cols())
    {
      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 { return m_eivec; }

    MatrixType pseudoEigenvalueMatrix() const;

    /** \returns the eigenvalues as a column vector */
    EigenvalueType eigenvalues() const { 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;
};

/** \returns the real block diagonal matrix D of the eigenvalues.
  *
  * See pseudoEigenvectors() for the details.
  */
template<typename MatrixType>
MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
{
  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
{
  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);
    }
    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);
}

// 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 <NAME> 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 <NAME> 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 = pow(2.0,-52.0);
  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)) / 2.0;
      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 = 0.75 * s;
        w = -0.4375 * s * s;
      }

      // MATLAB's new ad hoc shift
      if (iter == 30)
      {
{'fraction_non_alphanumeric': 0.15489893595981172, 'fraction_numerical': 0.011550176021518137, 'mean_word_length': 2.17413684871312, 'pattern_counts': {'":': 0, '<': 89, '<?xml version=': 0, '>': 79, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '4792061', 'n_tokens_mistral': 9961, 'n_tokens_neox': 9563, 'n_words': 2322}