DGMRES.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// 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_DGMRES_H
#define EIGEN_DGMRES_H
 
#include <Eigen/Eigenvalues>
 
namespace Eigen { 
  
template< typename _MatrixType,
          typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
class DGMRES;
 
namespace internal {
 
template< typename _MatrixType, typename _Preconditioner>
struct traits<DGMRES<_MatrixType,_Preconditioner> >
{
  typedef _MatrixType MatrixType;
  typedef _Preconditioner Preconditioner;
};

template <typename VectorType, typename IndexType>
void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut)
{
  eigen_assert(vec.size() == perm.size());
  bool flag; 
  for (Index k  = 0; k < ncut; k++)
  {
    flag = false;
    for (Index j = 0; j < vec.size()-1; j++)
    {
      if ( vec(perm(j)) < vec(perm(j+1)) )
      {
        std::swap(perm(j),perm(j+1)); 
        flag = true;
      }
      if (!flag) break; // The vector is in sorted order
    }
  }
}
 
}
template< typename _MatrixType, typename _Preconditioner>
class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
{
    typedef IterativeSolverBase<DGMRES> Base;
    using Base::matrix;
    using Base::m_error;
    using Base::m_iterations;
    using Base::m_info;
    using Base::m_isInitialized;
    using Base::m_tolerance; 
  public:
    using Base::_solve_impl;
    typedef _MatrixType MatrixType;
    typedef typename MatrixType::Scalar Scalar;
    typedef typename MatrixType::StorageIndex StorageIndex;
    typedef typename MatrixType::RealScalar RealScalar;
    typedef _Preconditioner Preconditioner;
    typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 
    typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix; 
    typedef Matrix<Scalar,Dynamic,1> DenseVector;
    typedef Matrix<RealScalar,Dynamic,1> DenseRealVector; 
    typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector;
 
    
  DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}

  template<typename MatrixDerived>
  explicit DGMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
 
  ~DGMRES() {}
  
  template<typename Rhs,typename Dest>
  void _solve_with_guess_impl(const Rhs& b, Dest& x) const
  {    
    bool failed = false;
    for(Index j=0; j<b.cols(); ++j)
    {
      m_iterations = Base::maxIterations();
      m_error = Base::m_tolerance;
      
      typename Dest::ColXpr xj(x,j);
      dgmres(matrix(), b.col(j), xj, Base::m_preconditioner);
    }
    m_info = failed ? NumericalIssue
           : m_error <= Base::m_tolerance ? Success
           : NoConvergence;
    m_isInitialized = true;
  }

  template<typename Rhs,typename Dest>
  void _solve_impl(const Rhs& b, MatrixBase<Dest>& x) const
  {
    x = b;
    _solve_with_guess_impl(b,x.derived());
  }
  Index restart() { return m_restart; }
  
  void set_restart(const Index restart) { m_restart=restart; }
  
  void setEigenv(const Index neig) 
  {
    m_neig = neig;
    if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
  }
   
  Index deflSize() {return m_r; }
  
  void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; }
  
  protected:
    // DGMRES algorithm 
    template<typename Rhs, typename Dest>
    void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
    // Perform one cycle of GMRES
    template<typename Dest>
    Index dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const; 
    // Compute data to use for deflation 
    Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const;
    // Apply deflation to a vector
    template<typename RhsType, typename DestType>
    Index dgmresApplyDeflation(const RhsType& In, DestType& Out) const; 
    ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
    ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
    // Init data for deflation
    void dgmresInitDeflation(Index& rows) const; 
    mutable DenseMatrix m_V; // Krylov basis vectors
    mutable DenseMatrix m_H; // Hessenberg matrix 
    mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied
    mutable Index m_restart; // Maximum size of the Krylov subspace
    mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace 
    mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
    mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
    mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
    mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart
    mutable Index m_r; // Current number of deflated eigenvalues, size of m_U
    mutable Index m_maxNeig; // Maximum number of eigenvalues to deflate
    mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
    mutable bool m_isDeflAllocated;
    mutable bool m_isDeflInitialized;
    
    //Adaptive strategy 
    mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed
    mutable bool m_force; // Force the use of deflation at each restart
    
}; 

template< typename _MatrixType, typename _Preconditioner>
template<typename Rhs, typename Dest>
void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
              const Preconditioner& precond) const
{
  //Initialization
  Index n = mat.rows(); 
  DenseVector r0(n); 
  Index nbIts = 0; 
  m_H.resize(m_restart+1, m_restart);
  m_Hes.resize(m_restart, m_restart);
  m_V.resize(n,m_restart+1);
  //Initial residual vector and intial norm
  x = precond.solve(x);
  r0 = rhs - mat * x; 
  RealScalar beta = r0.norm(); 
  RealScalar normRhs = rhs.norm();
  m_error = beta/normRhs; 
  if(m_error < m_tolerance)
    m_info = Success; 
  else
    m_info = NoConvergence;
  
  // Iterative process
  while (nbIts < m_iterations && m_info == NoConvergence)
  {
    dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts); 
    
    // Compute the new residual vector for the restart 
    if (nbIts < m_iterations && m_info == NoConvergence)
      r0 = rhs - mat * x; 
  }
} 

template< typename _MatrixType, typename _Preconditioner>
template<typename Dest>
Index DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const
{
  //Initialization 
  DenseVector g(m_restart+1); // Right hand side of the least square problem
  g.setZero();  
  g(0) = Scalar(beta); 
  m_V.col(0) = r0/beta; 
  m_info = NoConvergence; 
  std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
  Index it = 0; // Number of inner iterations 
  Index n = mat.rows();
  DenseVector tv1(n), tv2(n);  //Temporary vectors
  while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
  {    
    // Apply preconditioner(s) at right
    if (m_isDeflInitialized )
    {
      dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
      tv2 = precond.solve(tv1); 
    }
    else
    {
      tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
    }
    tv1 = mat * tv2; 
   
    // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
    Scalar coef; 
    for (Index i = 0; i <= it; ++i)
    { 
      coef = tv1.dot(m_V.col(i));
      tv1 = tv1 - coef * m_V.col(i); 
      m_H(i,it) = coef; 
      m_Hes(i,it) = coef; 
    }
    // Normalize the vector 
    coef = tv1.norm(); 
    m_V.col(it+1) = tv1/coef;
    m_H(it+1, it) = coef;
//     m_Hes(it+1,it) = coef; 
    
    // FIXME Check for happy breakdown 
    
    // Update Hessenberg matrix with Givens rotations
    for (Index i = 1; i <= it; ++i) 
    {
      m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
    }
    // Compute the new plane rotation 
    gr[it].makeGivens(m_H(it, it), m_H(it+1,it)); 
    // Apply the new rotation
    m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
    g.applyOnTheLeft(it,it+1, gr[it].adjoint()); 
    
    beta = std::abs(g(it+1));
    m_error = beta/normRhs; 
    // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
    it++; nbIts++; 
    
    if (m_error < m_tolerance)
    {
      // The method has converged
      m_info = Success;
      break;
    }
  }
  
  // Compute the new coefficients by solving the least square problem
//   it++;
  //FIXME  Check first if the matrix is singular ... zero diagonal
  DenseVector nrs(m_restart); 
  nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it)); 
  
  // Form the new solution
  if (m_isDeflInitialized)
  {
    tv1 = m_V.leftCols(it) * nrs; 
    dgmresApplyDeflation(tv1, tv2); 
    x = x + precond.solve(tv2);
  }
  else
    x = x + precond.solve(m_V.leftCols(it) * nrs); 
  
  // Go for a new cycle and compute data for deflation
  if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
    dgmresComputeDeflationData(mat, precond, it, m_neig); 
  return 0; 
  
}
 
 
template< typename _MatrixType, typename _Preconditioner>
void DGMRES<_MatrixType, _Preconditioner>::dgmresInitDeflation(Index& rows) const
{
  m_U.resize(rows, m_maxNeig);
  m_MU.resize(rows, m_maxNeig); 
  m_T.resize(m_maxNeig, m_maxNeig);
  m_lambdaN = 0.0; 
  m_isDeflAllocated = true; 
}
 
template< typename _MatrixType, typename _Preconditioner>
inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const
{
  return schurofH.matrixT().diagonal();
}
 
template< typename _MatrixType, typename _Preconditioner>
inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const
{
  const DenseMatrix& T = schurofH.matrixT();
  Index it = T.rows();
  ComplexVector eig(it);
  Index j = 0;
  while (j < it-1)
  {
    if (T(j+1,j) ==Scalar(0))
    {
      eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0)); 
      j++; 
    }
    else
    {
      eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j)); 
      eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1));
      j++;
    }
  }
  if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
  return eig;
}
 
template< typename _MatrixType, typename _Preconditioner>
Index DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const
{
  // First, find the Schur form of the Hessenberg matrix H
  typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH; 
  bool computeU = true;
  DenseMatrix matrixQ(it,it); 
  matrixQ.setIdentity();
  schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU); 
  
  ComplexVector eig(it);
  Matrix<StorageIndex,Dynamic,1>perm(it);
  eig = this->schurValues(schurofH);
  
  // Reorder the absolute values of Schur values
  DenseRealVector modulEig(it); 
  for (Index j=0; j<it; ++j) modulEig(j) = std::abs(eig(j)); 
  perm.setLinSpaced(it,0,internal::convert_index<StorageIndex>(it-1));
  internal::sortWithPermutation(modulEig, perm, neig);
  
  if (!m_lambdaN)
  {
    m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
  }
  //Count the real number of extracted eigenvalues (with complex conjugates)
  Index nbrEig = 0; 
  while (nbrEig < neig)
  {
    if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++; 
    else nbrEig += 2; 
  }
  // Extract the  Schur vectors corresponding to the smallest Ritz values
  DenseMatrix Sr(it, nbrEig); 
  Sr.setZero();
  for (Index j = 0; j < nbrEig; j++)
  {
    Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
  }
  
  // Form the Schur vectors of the initial matrix using the Krylov basis
  DenseMatrix X; 
  X = m_V.leftCols(it) * Sr;
  if (m_r)
  {
   // Orthogonalize X against m_U using modified Gram-Schmidt
   for (Index j = 0; j < nbrEig; j++)
     for (Index k =0; k < m_r; k++)
      X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k); 
  }
  
  // Compute m_MX = A * M^-1 * X
  Index m = m_V.rows();
  if (!m_isDeflAllocated) 
    dgmresInitDeflation(m); 
  DenseMatrix MX(m, nbrEig);
  DenseVector tv1(m);
  for (Index j = 0; j < nbrEig; j++)
  {
    tv1 = mat * X.col(j);
    MX.col(j) = precond.solve(tv1);
  }
  
  //Update m_T = [U'MU U'MX; X'MU X'MX]
  m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX; 
  if(m_r)
  {
    m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX; 
    m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r);
  }
  
  // Save X into m_U and m_MX in m_MU
  for (Index j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
  for (Index j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
  // Increase the size of the invariant subspace
  m_r += nbrEig; 
  
  // Factorize m_T into m_luT
  m_luT.compute(m_T.topLeftCorner(m_r, m_r));
  
  //FIXME CHeck if the factorization was correctly done (nonsingular matrix)
  m_isDeflInitialized = true;
  return 0; 
}
template<typename _MatrixType, typename _Preconditioner>
template<typename RhsType, typename DestType>
Index DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
{
  DenseVector x1 = m_U.leftCols(m_r).transpose() * x; 
  y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1);
  return 0; 
}
 
} // end namespace Eigen
#endif