// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 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_COMPLEX_SCHUR_H
#define EIGEN_COMPLEX_SCHUR_H

/** \eigenvalues_module \ingroup Eigenvalues_Module
  * \nonstableyet
  *
  * \class ComplexShur
  *
  * \brief Performs a complex Schur decomposition of a real or complex square matrix
  *
  * 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.
  *
  * \sa class RealSchur, class EigenSolver
  */
template<typename _MatrixType> class ComplexSchur
{
  public:
    typedef _MatrixType MatrixType;
    typedef typename MatrixType::Scalar Scalar;
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef std::complex<RealScalar> Complex;
    typedef Matrix<Complex, MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime> ComplexMatrixType;
    enum {
      Size = MatrixType::RowsAtCompileTime
    };

    /** \brief Default Constructor.
      *
      * The default constructor is useful in cases in which the user intends to
      * perform decompositions via ComplexSchur::compute().
      */
    ComplexSchur(int size = Size==Dynamic ? 0 : Size)
      : m_matT(size,size), m_matU(size,size), m_isInitialized(false), m_matUisUptodate(false)
    {}

    /** Constructor computing the Schur decomposition of the matrix \a matrix.
      * If \a skipU is true, then the matrix U is not computed. */
    ComplexSchur(const MatrixType& matrix, bool skipU = false)
            : m_matT(matrix.rows(),matrix.cols()),
              m_matU(matrix.rows(),matrix.cols()),
              m_isInitialized(false),
              m_matUisUptodate(false)
    {
      compute(matrix, skipU);
    }

    /** \returns a const reference to the matrix U of the respective Schur decomposition. */
    const ComplexMatrixType& matrixU() const
    {
      ei_assert(m_isInitialized && "ComplexSchur is not initialized.");
      ei_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
      return m_matU;
    }

    /** \returns a const reference to the matrix T of the respective Schur decomposition.
      * 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. */
    const ComplexMatrixType& matrixT() const
    {
      ei_assert(m_isInitialized && "ComplexShur is not initialized.");
      return m_matT;
    }

    /** Computes the Schur decomposition of the matrix \a matrix.
      * If \a skipU is true, then the matrix U is not computed. */
    void compute(const MatrixType& matrix, bool skipU = false);

  protected:
    ComplexMatrixType m_matT, m_matU;
    bool m_isInitialized;
    bool m_matUisUptodate;
};

/** Computes the principal value of the square root of the complex \a z. */
template<typename RealScalar>
std::complex<RealScalar> ei_sqrt(const std::complex<RealScalar> &z)
{
  RealScalar t, tre, tim;

  t = ei_abs(z);

  if (ei_abs(ei_real(z)) <= ei_abs(ei_imag(z)))
  {
    // No cancellation in these formulas
    tre = ei_sqrt(0.5*(t + ei_real(z)));
    tim = ei_sqrt(0.5*(t - ei_real(z)));
  }
  else
  {
    // Stable computation of the above formulas
    if (z.real() > 0)
    {
      tre = t + z.real();
      tim = ei_abs(ei_imag(z))*ei_sqrt(0.5/tre);
      tre = ei_sqrt(0.5*tre);
    }
    else
    {
      tim = t - z.real();
      tre = ei_abs(ei_imag(z))*ei_sqrt(0.5/tim);
      tim = ei_sqrt(0.5*tim);
    }
  }
  if(z.imag() < 0)
    tim = -tim;

  return (std::complex<RealScalar>(tre,tim));

}

template<typename MatrixType>
void ComplexSchur<MatrixType>::compute(const MatrixType& matrix, bool skipU)
{
  // this code is inspired from Jampack

  m_matUisUptodate = false;
  assert(matrix.cols() == matrix.rows());
  int n = matrix.cols();

  // Reduce to Hessenberg form
  // TODO skip Q if skipU = true
  HessenbergDecomposition<MatrixType> hess(matrix);

  m_matT = hess.matrixH();
  if(!skipU)  m_matU = hess.matrixQ();

  int iu = m_matT.cols() - 1;
  int il;
  RealScalar d,sd,sf;
  Complex c,b,disc,r1,r2,kappa;

  RealScalar eps = epsilon<RealScalar>();

  int iter = 0;
  while(true)
  {
    //locate the range in which to iterate
    while(iu > 0)
    {
      d = ei_norm1(m_matT.coeffRef(iu,iu)) + ei_norm1(m_matT.coeffRef(iu-1,iu-1));
      sd = ei_norm1(m_matT.coeffRef(iu,iu-1));

      if(sd >= eps * d) break; // FIXME : precision criterion ??

      m_matT.coeffRef(iu,iu-1) = Complex(0);
      iter = 0;
      --iu;
    }
    if(iu==0) break;
    iter++;

    if(iter >= 30)
    {
      // FIXME : what to do when iter==MAXITER ??
      std::cerr << "MAXITER" << std::endl;
      return;
    }

    il = iu-1;
    while( il > 0 )
    {
      // check if the current 2x2 block on the diagonal is upper triangular
      d = ei_norm1(m_matT.coeffRef(il,il)) + ei_norm1(m_matT.coeffRef(il-1,il-1));
      sd = ei_norm1(m_matT.coeffRef(il,il-1));

      if(sd < eps * d) break; // FIXME : precision criterion ??

      --il;
    }

    if( il != 0 ) m_matT.coeffRef(il,il-1) = Complex(0);

    // compute the shift (the normalization by sf is to avoid under/overflow)
    Matrix<Scalar,2,2> t = m_matT.template block<2,2>(iu-1,iu-1);
    sf = t.cwise().abs().sum();
    t /= sf;

    c = t.determinant();
    b = t.diagonal().sum();

    disc = ei_sqrt(b*b - RealScalar(4)*c);

    r1 = (b+disc)/RealScalar(2);
    r2 = (b-disc)/RealScalar(2);

    if(ei_norm1(r1) > ei_norm1(r2))
      r2 = c/r1;
    else
      r1 = c/r2;

    if(ei_norm1(r1-t.coeff(1,1)) < ei_norm1(r2-t.coeff(1,1)))
      kappa = sf * r1;
    else
      kappa = sf * r2;

    // perform the QR step using Givens rotations
    PlanarRotation<Complex> rot;
    rot.makeGivens(m_matT.coeff(il,il) - kappa, m_matT.coeff(il+1,il));

    for(int i=il ; i<iu ; i++)
    {
      m_matT.block(0,i,n,n-i).applyOnTheLeft(i, i+1, rot.adjoint());
      m_matT.block(0,0,std::min(i+2,iu)+1,n).applyOnTheRight(i, i+1, rot);
      if(!skipU) m_matU.applyOnTheRight(i, i+1, rot);

      if(i != iu-1)
      {
        int i1 = i+1;
        int i2 = i+2;

        rot.makeGivens(m_matT.coeffRef(i1,i), m_matT.coeffRef(i2,i), &m_matT.coeffRef(i1,i));
        m_matT.coeffRef(i2,i) = Complex(0);
      }
    }
  }

  // FIXME : is it necessary ?
  /*
  for(int i=0 ; i<n ; i++)
    for(int j=0 ; j<n ; j++)
    {
      if(ei_abs(ei_real(m_matT.coeff(i,j))) < eps)
        ei_real_ref(m_matT.coeffRef(i,j)) = 0;
      if(ei_imag(ei_abs(m_matT.coeff(i,j))) < eps)
        ei_imag_ref(m_matT.coeffRef(i,j)) = 0;
    }
  */

  m_isInitialized = true;
  m_matUisUptodate = !skipU;
}

#endif // EIGEN_COMPLEX_SCHUR_H