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