DGMRES.h
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 
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 
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  bool flag;
43  for (Index k = 0; k < ncut; k++)
44  {
45  flag = false;
46  for (Index j = 0; j < vec.size()-1; j++)
47  {
48  if ( vec(perm(j)) < vec(perm(j+1)) )
49  {
50  std::swap(perm(j),perm(j+1));
51  flag = true;
52  }
53  if (!flag) break; // The vector is in sorted order
54  }
55  }
56 }
57 
58 }
100 template< typename _MatrixType, typename _Preconditioner>
101 class DGMRES : public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
102 {
103  typedef IterativeSolverBase<DGMRES> Base;
104  using Base::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  using Base::_solve_impl;
112  typedef _MatrixType MatrixType;
113  typedef typename MatrixType::Scalar Scalar;
114  typedef typename MatrixType::StorageIndex StorageIndex;
115  typedef typename MatrixType::RealScalar RealScalar;
116  typedef _Preconditioner Preconditioner;
117  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
118  typedef Matrix<RealScalar,Dynamic,Dynamic> DenseRealMatrix;
119  typedef Matrix<Scalar,Dynamic,1> DenseVector;
120  typedef Matrix<RealScalar,Dynamic,1> DenseRealVector;
121  typedef Matrix<std::complex<RealScalar>, Dynamic, 1> ComplexVector;
122 
123 
125  DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
126 
137  template<typename MatrixDerived>
138  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) {}
139 
140  ~DGMRES() {}
141 
143  template<typename Rhs,typename Dest>
144  void _solve_with_guess_impl(const Rhs& b, Dest& x) const
145  {
146  bool failed = false;
147  for(Index j=0; j<b.cols(); ++j)
148  {
149  m_iterations = Base::maxIterations();
150  m_error = Base::m_tolerance;
151 
152  typename Dest::ColXpr xj(x,j);
153  dgmres(matrix(), b.col(j), xj, Base::m_preconditioner);
154  }
155  m_info = failed ? NumericalIssue
156  : m_error <= Base::m_tolerance ? Success
157  : NoConvergence;
158  m_isInitialized = true;
159  }
160 
162  template<typename Rhs,typename Dest>
163  void _solve_impl(const Rhs& b, MatrixBase<Dest>& x) const
164  {
165  x = b;
166  _solve_with_guess_impl(b,x.derived());
167  }
171  Index restart() { return m_restart; }
172 
176  void set_restart(const Index restart) { m_restart=restart; }
177 
181  void setEigenv(const Index neig)
182  {
183  m_neig = neig;
184  if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
185  }
186 
190  Index deflSize() {return m_r; }
191 
195  void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; }
196 
197  protected:
198  // DGMRES algorithm
199  template<typename Rhs, typename Dest>
200  void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
201  // Perform one cycle of GMRES
202  template<typename Dest>
203  Index dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const;
204  // Compute data to use for deflation
205  Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const;
206  // Apply deflation to a vector
207  template<typename RhsType, typename DestType>
208  Index dgmresApplyDeflation(const RhsType& In, DestType& Out) const;
209  ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
210  ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
211  // Init data for deflation
212  void dgmresInitDeflation(Index& rows) const;
213  mutable DenseMatrix m_V; // Krylov basis vectors
214  mutable DenseMatrix m_H; // Hessenberg matrix
215  mutable DenseMatrix m_Hes; // Initial hessenberg matrix wihout Givens rotations applied
216  mutable Index m_restart; // Maximum size of the Krylov subspace
217  mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace
218  mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
219  mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
220  mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
221  mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart
222  mutable Index m_r; // Current number of deflated eigenvalues, size of m_U
223  mutable Index m_maxNeig; // Maximum number of eigenvalues to deflate
224  mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
225  mutable bool m_isDeflAllocated;
226  mutable bool m_isDeflInitialized;
227 
228  //Adaptive strategy
229  mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed
230  mutable bool m_force; // Force the use of deflation at each restart
231 
232 };
239 template< typename _MatrixType, typename _Preconditioner>
240 template<typename Rhs, typename Dest>
241 void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
242  const Preconditioner& precond) const
243 {
244  //Initialization
245  Index n = mat.rows();
246  DenseVector r0(n);
247  Index nbIts = 0;
248  m_H.resize(m_restart+1, m_restart);
249  m_Hes.resize(m_restart, m_restart);
250  m_V.resize(n,m_restart+1);
251  //Initial residual vector and intial norm
252  x = precond.solve(x);
253  r0 = rhs - mat * x;
254  RealScalar beta = r0.norm();
255  RealScalar normRhs = rhs.norm();
256  m_error = beta/normRhs;
257  if(m_error < m_tolerance)
258  m_info = Success;
259  else
260  m_info = NoConvergence;
261 
262  // Iterative process
263  while (nbIts < m_iterations && m_info == NoConvergence)
264  {
265  dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts);
266 
267  // Compute the new residual vector for the restart
268  if (nbIts < m_iterations && m_info == NoConvergence)
269  r0 = rhs - mat * x;
270  }
271 }
272 
283 template< typename _MatrixType, typename _Preconditioner>
284 template<typename Dest>
285 Index DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const
286 {
287  //Initialization
288  DenseVector g(m_restart+1); // Right hand side of the least square problem
289  g.setZero();
290  g(0) = Scalar(beta);
291  m_V.col(0) = r0/beta;
292  m_info = NoConvergence;
293  std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
294  Index it = 0; // Number of inner iterations
295  Index n = mat.rows();
296  DenseVector tv1(n), tv2(n); //Temporary vectors
297  while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
298  {
299  // Apply preconditioner(s) at right
300  if (m_isDeflInitialized )
301  {
302  dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
303  tv2 = precond.solve(tv1);
304  }
305  else
306  {
307  tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
308  }
309  tv1 = mat * tv2;
310 
311  // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
312  Scalar coef;
313  for (Index i = 0; i <= it; ++i)
314  {
315  coef = tv1.dot(m_V.col(i));
316  tv1 = tv1 - coef * m_V.col(i);
317  m_H(i,it) = coef;
318  m_Hes(i,it) = coef;
319  }
320  // Normalize the vector
321  coef = tv1.norm();
322  m_V.col(it+1) = tv1/coef;
323  m_H(it+1, it) = coef;
324 // m_Hes(it+1,it) = coef;
325 
326  // FIXME Check for happy breakdown
327 
328  // Update Hessenberg matrix with Givens rotations
329  for (Index i = 1; i <= it; ++i)
330  {
331  m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
332  }
333  // Compute the new plane rotation
334  gr[it].makeGivens(m_H(it, it), m_H(it+1,it));
335  // Apply the new rotation
336  m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
337  g.applyOnTheLeft(it,it+1, gr[it].adjoint());
338 
339  beta = std::abs(g(it+1));
340  m_error = beta/normRhs;
341  // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
342  it++; nbIts++;
343 
344  if (m_error < m_tolerance)
345  {
346  // The method has converged
347  m_info = Success;
348  break;
349  }
350  }
351 
352  // Compute the new coefficients by solving the least square problem
353 // it++;
354  //FIXME Check first if the matrix is singular ... zero diagonal
355  DenseVector nrs(m_restart);
356  nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it));
357 
358  // Form the new solution
359  if (m_isDeflInitialized)
360  {
361  tv1 = m_V.leftCols(it) * nrs;
362  dgmresApplyDeflation(tv1, tv2);
363  x = x + precond.solve(tv2);
364  }
365  else
366  x = x + precond.solve(m_V.leftCols(it) * nrs);
367 
368  // Go for a new cycle and compute data for deflation
369  if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
370  dgmresComputeDeflationData(mat, precond, it, m_neig);
371  return 0;
372 
373 }
374 
375 
376 template< typename _MatrixType, typename _Preconditioner>
378 {
379  m_U.resize(rows, m_maxNeig);
380  m_MU.resize(rows, m_maxNeig);
381  m_T.resize(m_maxNeig, m_maxNeig);
382  m_lambdaN = 0.0;
383  m_isDeflAllocated = true;
384 }
385 
386 template< typename _MatrixType, typename _Preconditioner>
387 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const
388 {
389  return schurofH.matrixT().diagonal();
390 }
391 
392 template< typename _MatrixType, typename _Preconditioner>
393 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const
394 {
395  const DenseMatrix& T = schurofH.matrixT();
396  Index it = T.rows();
397  ComplexVector eig(it);
398  Index j = 0;
399  while (j < it-1)
400  {
401  if (T(j+1,j) ==Scalar(0))
402  {
403  eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
404  j++;
405  }
406  else
407  {
408  eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j));
409  eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1));
410  j++;
411  }
412  }
413  if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
414  return eig;
415 }
416 
417 template< typename _MatrixType, typename _Preconditioner>
418 Index DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const
419 {
420  // First, find the Schur form of the Hessenberg matrix H
421  typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH;
422  bool computeU = true;
423  DenseMatrix matrixQ(it,it);
424  matrixQ.setIdentity();
425  schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU);
426 
427  ComplexVector eig(it);
429  eig = this->schurValues(schurofH);
430 
431  // Reorder the absolute values of Schur values
432  DenseRealVector modulEig(it);
433  for (Index j=0; j<it; ++j) modulEig(j) = std::abs(eig(j));
434  perm.setLinSpaced(it,0,internal::convert_index<StorageIndex>(it-1));
435  internal::sortWithPermutation(modulEig, perm, neig);
436 
437  if (!m_lambdaN)
438  {
439  m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
440  }
441  //Count the real number of extracted eigenvalues (with complex conjugates)
442  Index nbrEig = 0;
443  while (nbrEig < neig)
444  {
445  if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++;
446  else nbrEig += 2;
447  }
448  // Extract the Schur vectors corresponding to the smallest Ritz values
449  DenseMatrix Sr(it, nbrEig);
450  Sr.setZero();
451  for (Index j = 0; j < nbrEig; j++)
452  {
453  Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
454  }
455 
456  // Form the Schur vectors of the initial matrix using the Krylov basis
457  DenseMatrix X;
458  X = m_V.leftCols(it) * Sr;
459  if (m_r)
460  {
461  // Orthogonalize X against m_U using modified Gram-Schmidt
462  for (Index j = 0; j < nbrEig; j++)
463  for (Index k =0; k < m_r; k++)
464  X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k);
465  }
466 
467  // Compute m_MX = A * M^-1 * X
468  Index m = m_V.rows();
469  if (!m_isDeflAllocated)
470  dgmresInitDeflation(m);
471  DenseMatrix MX(m, nbrEig);
472  DenseVector tv1(m);
473  for (Index j = 0; j < nbrEig; j++)
474  {
475  tv1 = mat * X.col(j);
476  MX.col(j) = precond.solve(tv1);
477  }
478 
479  //Update m_T = [U'MU U'MX; X'MU X'MX]
480  m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX;
481  if(m_r)
482  {
483  m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX;
484  m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r);
485  }
486 
487  // Save X into m_U and m_MX in m_MU
488  for (Index j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
489  for (Index j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
490  // Increase the size of the invariant subspace
491  m_r += nbrEig;
492 
493  // Factorize m_T into m_luT
494  m_luT.compute(m_T.topLeftCorner(m_r, m_r));
495 
496  //FIXME CHeck if the factorization was correctly done (nonsingular matrix)
497  m_isDeflInitialized = true;
498  return 0;
499 }
500 template<typename _MatrixType, typename _Preconditioner>
501 template<typename RhsType, typename DestType>
502 Index DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
503 {
504  DenseVector x1 = m_U.leftCols(m_r).transpose() * x;
505  y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1);
506  return 0;
507 }
508 
509 } // end namespace Eigen
510 #endif
A Restarted GMRES with deflation. This class implements a modification of the GMRES solver for sparse...
Definition: DGMRES.h:19
void dgmres(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond) const
Perform several cycles of restarted GMRES with modified Gram Schmidt,.
Definition: DGMRES.h:241
void setEigenv(const Index neig)
Definition: DGMRES.h:181
Index deflSize()
Definition: DGMRES.h:190
Namespace containing all symbols from the Eigen library.
Derived & setZero(Index size)
void set_restart(const Index restart)
Definition: DGMRES.h:176
NumericalIssue
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
Index dgmresCycle(const MatrixType &mat, const Preconditioner &precond, Dest &x, DenseVector &r0, RealScalar &beta, const RealScalar &normRhs, Index &nbIts) const
Perform one restart cycle of DGMRES.
Definition: DGMRES.h:285
Index restart()
Definition: DGMRES.h:171
DGMRES()
Definition: DGMRES.h:125
void setMaxEigenv(const Index maxNeig)
Definition: DGMRES.h:195
const int Dynamic
DGMRES(const EigenBase< MatrixDerived > &A)
Definition: DGMRES.h:138