Home | History | Annotate | Download | only in IterativeSolvers
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_DGMRES_H
     11 #define EIGEN_DGMRES_H
     12 
     13 #include <Eigen/Eigenvalues>
     14 
     15 namespace Eigen {
     16 
     17 template< typename _MatrixType,
     18           typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
     19 class DGMRES;
     20 
     21 namespace internal {
     22 
     23 template< typename _MatrixType, typename _Preconditioner>
     24 struct traits<DGMRES<_MatrixType,_Preconditioner> >
     25 {
     26   typedef _MatrixType MatrixType;
     27   typedef _Preconditioner Preconditioner;
     28 };
     29 
     30 /** \brief Computes a permutation vector to have a sorted sequence
     31   * \param vec The vector to reorder.
     32   * \param perm gives the sorted sequence on output. Must be initialized with 0..n-1
     33   * \param ncut Put  the ncut smallest elements at the end of the vector
     34   * WARNING This is an expensive sort, so should be used only
     35   * for small size vectors
     36   * TODO Use modified QuickSplit or std::nth_element to get the smallest values
     37   */
     38 template <typename VectorType, typename IndexType>
     39 void sortWithPermutation (VectorType& vec, IndexType& perm, typename IndexType::Scalar& ncut)
     40 {
     41   eigen_assert(vec.size() == perm.size());
     42   typedef typename IndexType::Scalar Index;
     43   typedef typename VectorType::Scalar Scalar;
     44   bool flag;
     45   for (Index k  = 0; k < ncut; k++)
     46   {
     47     flag = false;
     48     for (Index j = 0; j < vec.size()-1; j++)
     49     {
     50       if ( vec(perm(j)) < vec(perm(j+1)) )
     51       {
     52         std::swap(perm(j),perm(j+1));
     53         flag = true;
     54       }
     55       if (!flag) break; // The vector is in sorted order
     56     }
     57   }
     58 }
     59 
     60 }
     61 /**
     62  * \ingroup IterativeLInearSolvers_Module
     63  * \brief A Restarted GMRES with deflation.
     64  * This class implements a modification of the GMRES solver for
     65  * sparse linear systems. The basis is built with modified
     66  * Gram-Schmidt. At each restart, a few approximated eigenvectors
     67  * corresponding to the smallest eigenvalues are used to build a
     68  * preconditioner for the next cycle. This preconditioner
     69  * for deflation can be combined with any other preconditioner,
     70  * the IncompleteLUT for instance. The preconditioner is applied
     71  * at right of the matrix and the combination is multiplicative.
     72  *
     73  * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
     74  * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
     75  * Typical usage :
     76  * \code
     77  * SparseMatrix<double> A;
     78  * VectorXd x, b;
     79  * //Fill A and b ...
     80  * DGMRES<SparseMatrix<double> > solver;
     81  * solver.set_restart(30); // Set restarting value
     82  * solver.setEigenv(1); // Set the number of eigenvalues to deflate
     83  * solver.compute(A);
     84  * x = solver.solve(b);
     85  * \endcode
     86  *
     87  * References :
     88  * [1] D. NUENTSA WAKAM and F. PACULL, Memory Efficient Hybrid
     89  *  Algebraic Solvers for Linear Systems Arising from Compressible
     90  *  Flows, Computers and Fluids, In Press,
     91  *  http://dx.doi.org/10.1016/j.compfluid.2012.03.023
     92  * [2] K. Burrage and J. Erhel, On the performance of various
     93  * adaptive preconditioned GMRES strategies, 5(1998), 101-121.
     94  * [3] J. Erhel, K. Burrage and B. Pohl, Restarted GMRES
     95  *  preconditioned by deflation,J. Computational and Applied
     96  *  Mathematics, 69(1996), 303-318.
     97 
     98  *
     99  */
    100 template< typename _MatrixType, typename _Preconditioner>
    101 class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
    102 {
    103     typedef IterativeSolverBase<DGMRES> Base;
    104     using Base::mp_matrix;
    105     using Base::m_error;
    106     using Base::m_iterations;
    107     using Base::m_info;
    108     using Base::m_isInitialized;
    109     using Base::m_tolerance;
    110   public:
    111     typedef _MatrixType MatrixType;
    112     typedef typename MatrixType::Scalar Scalar;
    113     typedef typename MatrixType::Index Index;
    114     typedef typename MatrixType::RealScalar RealScalar;
    115     typedef _Preconditioner Preconditioner;
    116     typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
    117     typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix;
    118     typedef Matrix<Scalar,Dynamic,1> DenseVector;
    119     typedef Matrix<RealScalar,Dynamic,1> DenseRealVector;
    120     typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector;
    121 
    122 
    123   /** Default constructor. */
    124   DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
    125 
    126   /** Initialize the solver with matrix \a A for further \c Ax=b solving.
    127     *
    128     * This constructor is a shortcut for the default constructor followed
    129     * by a call to compute().
    130     *
    131     * \warning this class stores a reference to the matrix A as well as some
    132     * precomputed values that depend on it. Therefore, if \a A is changed
    133     * this class becomes invalid. Call compute() to update it with the new
    134     * matrix A, or modify a copy of A.
    135     */
    136   DGMRES(const MatrixType& A) : Base(A),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false)
    137   {}
    138 
    139   ~DGMRES() {}
    140 
    141   /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
    142     * \a x0 as an initial solution.
    143     *
    144     * \sa compute()
    145     */
    146   template<typename Rhs,typename Guess>
    147   inline const internal::solve_retval_with_guess<DGMRES, Rhs, Guess>
    148   solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
    149   {
    150     eigen_assert(m_isInitialized && "DGMRES is not initialized.");
    151     eigen_assert(Base::rows()==b.rows()
    152               && "DGMRES::solve(): invalid number of rows of the right hand side matrix b");
    153     return internal::solve_retval_with_guess
    154             <DGMRES, Rhs, Guess>(*this, b.derived(), x0);
    155   }
    156 
    157   /** \internal */
    158   template<typename Rhs,typename Dest>
    159   void _solveWithGuess(const Rhs& b, Dest& x) const
    160   {
    161     bool failed = false;
    162     for(int j=0; j<b.cols(); ++j)
    163     {
    164       m_iterations = Base::maxIterations();
    165       m_error = Base::m_tolerance;
    166 
    167       typename Dest::ColXpr xj(x,j);
    168       dgmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner);
    169     }
    170     m_info = failed ? NumericalIssue
    171            : m_error <= Base::m_tolerance ? Success
    172            : NoConvergence;
    173     m_isInitialized = true;
    174   }
    175 
    176   /** \internal */
    177   template<typename Rhs,typename Dest>
    178   void _solve(const Rhs& b, Dest& x) const
    179   {
    180     x = b;
    181     _solveWithGuess(b,x);
    182   }
    183   /**
    184    * Get the restart value
    185     */
    186   int restart() { return m_restart; }
    187 
    188   /**
    189    * Set the restart value (default is 30)
    190    */
    191   void set_restart(const int restart) { m_restart=restart; }
    192 
    193   /**
    194    * Set the number of eigenvalues to deflate at each restart
    195    */
    196   void setEigenv(const int neig)
    197   {
    198     m_neig = neig;
    199     if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
    200   }
    201 
    202   /**
    203    * Get the size of the deflation subspace size
    204    */
    205   int deflSize() {return m_r; }
    206 
    207   /**
    208    * Set the maximum size of the deflation subspace
    209    */
    210   void setMaxEigenv(const int maxNeig) { m_maxNeig = maxNeig; }
    211 
    212   protected:
    213     // DGMRES algorithm
    214     template<typename Rhs, typename Dest>
    215     void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
    216     // Perform one cycle of GMRES
    217     template<typename Dest>
    218     int dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const;
    219     // Compute data to use for deflation
    220     int dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const;
    221     // Apply deflation to a vector
    222     template<typename RhsType, typename DestType>
    223     int dgmresApplyDeflation(const RhsType& In, DestType& Out) const;
    224     ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
    225     ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
    226     // Init data for deflation
    227     void dgmresInitDeflation(Index& rows) const;
    228     mutable DenseMatrix m_V; // Krylov basis vectors
    229     mutable DenseMatrix m_H; // Hessenberg matrix
    230     mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied
    231     mutable Index m_restart; // Maximum size of the Krylov subspace
    232     mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace
    233     mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
    234     mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
    235     mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
    236     mutable int m_neig; //Number of eigenvalues to extract at each restart
    237     mutable int m_r; // Current number of deflated eigenvalues, size of m_U
    238     mutable int m_maxNeig; // Maximum number of eigenvalues to deflate
    239     mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
    240     mutable bool m_isDeflAllocated;
    241     mutable bool m_isDeflInitialized;
    242 
    243     //Adaptive strategy
    244     mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed
    245     mutable bool m_force; // Force the use of deflation at each restart
    246 
    247 };
    248 /**
    249  * \brief Perform several cycles of restarted GMRES with modified Gram Schmidt,
    250  *
    251  * A right preconditioner is used combined with deflation.
    252  *
    253  */
    254 template< typename _MatrixType, typename _Preconditioner>
    255 template<typename Rhs, typename Dest>
    256 void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
    257               const Preconditioner& precond) const
    258 {
    259   //Initialization
    260   int n = mat.rows();
    261   DenseVector r0(n);
    262   int nbIts = 0;
    263   m_H.resize(m_restart+1, m_restart);
    264   m_Hes.resize(m_restart, m_restart);
    265   m_V.resize(n,m_restart+1);
    266   //Initial residual vector and intial norm
    267   x = precond.solve(x);
    268   r0 = rhs - mat * x;
    269   RealScalar beta = r0.norm();
    270   RealScalar normRhs = rhs.norm();
    271   m_error = beta/normRhs;
    272   if(m_error < m_tolerance)
    273     m_info = Success;
    274   else
    275     m_info = NoConvergence;
    276 
    277   // Iterative process
    278   while (nbIts < m_iterations && m_info == NoConvergence)
    279   {
    280     dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts);
    281 
    282     // Compute the new residual vector for the restart
    283     if (nbIts < m_iterations && m_info == NoConvergence)
    284       r0 = rhs - mat * x;
    285   }
    286 }
    287 
    288 /**
    289  * \brief Perform one restart cycle of DGMRES
    290  * \param mat The coefficient matrix
    291  * \param precond The preconditioner
    292  * \param x the new approximated solution
    293  * \param r0 The initial residual vector
    294  * \param beta The norm of the residual computed so far
    295  * \param normRhs The norm of the right hand side vector
    296  * \param nbIts The number of iterations
    297  */
    298 template< typename _MatrixType, typename _Preconditioner>
    299 template<typename Dest>
    300 int DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, int& nbIts) const
    301 {
    302   //Initialization
    303   DenseVector g(m_restart+1); // Right hand side of the least square problem
    304   g.setZero();
    305   g(0) = Scalar(beta);
    306   m_V.col(0) = r0/beta;
    307   m_info = NoConvergence;
    308   std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
    309   int it = 0; // Number of inner iterations
    310   int n = mat.rows();
    311   DenseVector tv1(n), tv2(n);  //Temporary vectors
    312   while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
    313   {
    314     // Apply preconditioner(s) at right
    315     if (m_isDeflInitialized )
    316     {
    317       dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
    318       tv2 = precond.solve(tv1);
    319     }
    320     else
    321     {
    322       tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
    323     }
    324     tv1 = mat * tv2;
    325 
    326     // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
    327     Scalar coef;
    328     for (int i = 0; i <= it; ++i)
    329     {
    330       coef = tv1.dot(m_V.col(i));
    331       tv1 = tv1 - coef * m_V.col(i);
    332       m_H(i,it) = coef;
    333       m_Hes(i,it) = coef;
    334     }
    335     // Normalize the vector
    336     coef = tv1.norm();
    337     m_V.col(it+1) = tv1/coef;
    338     m_H(it+1, it) = coef;
    339 //     m_Hes(it+1,it) = coef;
    340 
    341     // FIXME Check for happy breakdown
    342 
    343     // Update Hessenberg matrix with Givens rotations
    344     for (int i = 1; i <= it; ++i)
    345     {
    346       m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
    347     }
    348     // Compute the new plane rotation
    349     gr[it].makeGivens(m_H(it, it), m_H(it+1,it));
    350     // Apply the new rotation
    351     m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
    352     g.applyOnTheLeft(it,it+1, gr[it].adjoint());
    353 
    354     beta = std::abs(g(it+1));
    355     m_error = beta/normRhs;
    356     std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
    357     it++; nbIts++;
    358 
    359     if (m_error < m_tolerance)
    360     {
    361       // The method has converged
    362       m_info = Success;
    363       break;
    364     }
    365   }
    366 
    367   // Compute the new coefficients by solving the least square problem
    368 //   it++;
    369   //FIXME  Check first if the matrix is singular ... zero diagonal
    370   DenseVector nrs(m_restart);
    371   nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it));
    372 
    373   // Form the new solution
    374   if (m_isDeflInitialized)
    375   {
    376     tv1 = m_V.leftCols(it) * nrs;
    377     dgmresApplyDeflation(tv1, tv2);
    378     x = x + precond.solve(tv2);
    379   }
    380   else
    381     x = x + precond.solve(m_V.leftCols(it) * nrs);
    382 
    383   // Go for a new cycle and compute data for deflation
    384   if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
    385     dgmresComputeDeflationData(mat, precond, it, m_neig);
    386   return 0;
    387 
    388 }
    389 
    390 
    391 template< typename _MatrixType, typename _Preconditioner>
    392 void DGMRES<_MatrixType, _Preconditioner>::dgmresInitDeflation(Index& rows) const
    393 {
    394   m_U.resize(rows, m_maxNeig);
    395   m_MU.resize(rows, m_maxNeig);
    396   m_T.resize(m_maxNeig, m_maxNeig);
    397   m_lambdaN = 0.0;
    398   m_isDeflAllocated = true;
    399 }
    400 
    401 template< typename _MatrixType, typename _Preconditioner>
    402 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const
    403 {
    404   return schurofH.matrixT().diagonal();
    405 }
    406 
    407 template< typename _MatrixType, typename _Preconditioner>
    408 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const
    409 {
    410   typedef typename MatrixType::Index Index;
    411   const DenseMatrix& T = schurofH.matrixT();
    412   Index it = T.rows();
    413   ComplexVector eig(it);
    414   Index j = 0;
    415   while (j < it-1)
    416   {
    417     if (T(j+1,j) ==Scalar(0))
    418     {
    419       eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
    420       j++;
    421     }
    422     else
    423     {
    424       eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j));
    425       eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1));
    426       j++;
    427     }
    428   }
    429   if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
    430   return eig;
    431 }
    432 
    433 template< typename _MatrixType, typename _Preconditioner>
    434 int DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, Index& neig) const
    435 {
    436   // First, find the Schur form of the Hessenberg matrix H
    437   typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH;
    438   bool computeU = true;
    439   DenseMatrix matrixQ(it,it);
    440   matrixQ.setIdentity();
    441   schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU);
    442 
    443   ComplexVector eig(it);
    444   Matrix<Index,Dynamic,1>perm(it);
    445   eig = this->schurValues(schurofH);
    446 
    447   // Reorder the absolute values of Schur values
    448   DenseRealVector modulEig(it);
    449   for (int j=0; j<it; ++j) modulEig(j) = std::abs(eig(j));
    450   perm.setLinSpaced(it,0,it-1);
    451   internal::sortWithPermutation(modulEig, perm, neig);
    452 
    453   if (!m_lambdaN)
    454   {
    455     m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
    456   }
    457   //Count the real number of extracted eigenvalues (with complex conjugates)
    458   int nbrEig = 0;
    459   while (nbrEig < neig)
    460   {
    461     if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++;
    462     else nbrEig += 2;
    463   }
    464   // Extract the  Schur vectors corresponding to the smallest Ritz values
    465   DenseMatrix Sr(it, nbrEig);
    466   Sr.setZero();
    467   for (int j = 0; j < nbrEig; j++)
    468   {
    469     Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
    470   }
    471 
    472   // Form the Schur vectors of the initial matrix using the Krylov basis
    473   DenseMatrix X;
    474   X = m_V.leftCols(it) * Sr;
    475   if (m_r)
    476   {
    477    // Orthogonalize X against m_U using modified Gram-Schmidt
    478    for (int j = 0; j < nbrEig; j++)
    479      for (int k =0; k < m_r; k++)
    480       X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k);
    481   }
    482 
    483   // Compute m_MX = A * M^-1 * X
    484   Index m = m_V.rows();
    485   if (!m_isDeflAllocated)
    486     dgmresInitDeflation(m);
    487   DenseMatrix MX(m, nbrEig);
    488   DenseVector tv1(m);
    489   for (int j = 0; j < nbrEig; j++)
    490   {
    491     tv1 = mat * X.col(j);
    492     MX.col(j) = precond.solve(tv1);
    493   }
    494 
    495   //Update m_T = [U'MU U'MX; X'MU X'MX]
    496   m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX;
    497   if(m_r)
    498   {
    499     m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX;
    500     m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r);
    501   }
    502 
    503   // Save X into m_U and m_MX in m_MU
    504   for (int j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
    505   for (int j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
    506   // Increase the size of the invariant subspace
    507   m_r += nbrEig;
    508 
    509   // Factorize m_T into m_luT
    510   m_luT.compute(m_T.topLeftCorner(m_r, m_r));
    511 
    512   //FIXME CHeck if the factorization was correctly done (nonsingular matrix)
    513   m_isDeflInitialized = true;
    514   return 0;
    515 }
    516 template<typename _MatrixType, typename _Preconditioner>
    517 template<typename RhsType, typename DestType>
    518 int DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
    519 {
    520   DenseVector x1 = m_U.leftCols(m_r).transpose() * x;
    521   y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1);
    522   return 0;
    523 }
    524 
    525 namespace internal {
    526 
    527   template<typename _MatrixType, typename _Preconditioner, typename Rhs>
    528 struct solve_retval<DGMRES<_MatrixType, _Preconditioner>, Rhs>
    529   : solve_retval_base<DGMRES<_MatrixType, _Preconditioner>, Rhs>
    530 {
    531   typedef DGMRES<_MatrixType, _Preconditioner> Dec;
    532   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
    533 
    534   template<typename Dest> void evalTo(Dest& dst) const
    535   {
    536     dec()._solve(rhs(),dst);
    537   }
    538 };
    539 } // end namespace internal
    540 
    541 } // end namespace Eigen
    542 #endif
    543