DGMRES.h
Go to the documentation of this file.
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 {
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  using Base::_solve_with_guess_impl;
113  typedef _MatrixType MatrixType;
114  typedef typename MatrixType::Scalar Scalar;
115  typedef typename MatrixType::StorageIndex StorageIndex;
117  typedef _Preconditioner Preconditioner;
123 
124 
126  DGMRES() : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
127 
138  template<typename MatrixDerived>
139  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) {}
140 
141  ~DGMRES() {}
142 
144  template<typename Rhs,typename Dest>
145  void _solve_vector_with_guess_impl(const Rhs& b, Dest& x) const
146  {
147  EIGEN_STATIC_ASSERT(Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1, YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
148 
149  m_iterations = Base::maxIterations();
150  m_error = Base::m_tolerance;
151 
152  dgmres(matrix(), b, x, Base::m_preconditioner);
153  }
154 
158  Index restart() { return m_restart; }
159 
163  void set_restart(const Index restart) { m_restart=restart; }
164 
168  void setEigenv(const Index neig)
169  {
170  m_neig = neig;
171  if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
172  }
173 
177  Index deflSize() {return m_r; }
178 
182  void setMaxEigenv(const Index maxNeig) { m_maxNeig = maxNeig; }
183 
184  protected:
185  // DGMRES algorithm
186  template<typename Rhs, typename Dest>
187  void dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x, const Preconditioner& precond) const;
188  // Perform one cycle of GMRES
189  template<typename Dest>
190  Index dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const;
191  // Compute data to use for deflation
192  Index dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const;
193  // Apply deflation to a vector
194  template<typename RhsType, typename DestType>
195  Index dgmresApplyDeflation(const RhsType& In, DestType& Out) const;
196  ComplexVector schurValues(const ComplexSchur<DenseMatrix>& schurofH) const;
197  ComplexVector schurValues(const RealSchur<DenseMatrix>& schurofH) const;
198  // Init data for deflation
199  void dgmresInitDeflation(Index& rows) const;
200  mutable DenseMatrix m_V; // Krylov basis vectors
201  mutable DenseMatrix m_H; // Hessenberg matrix
202  mutable DenseMatrix m_Hes; // Initial hessenberg matrix without Givens rotations applied
203  mutable Index m_restart; // Maximum size of the Krylov subspace
204  mutable DenseMatrix m_U; // Vectors that form the basis of the invariant subspace
205  mutable DenseMatrix m_MU; // matrix operator applied to m_U (for next cycles)
206  mutable DenseMatrix m_T; /* T=U^T*M^{-1}*A*U */
207  mutable PartialPivLU<DenseMatrix> m_luT; // LU factorization of m_T
208  mutable StorageIndex m_neig; //Number of eigenvalues to extract at each restart
209  mutable Index m_r; // Current number of deflated eigenvalues, size of m_U
210  mutable Index m_maxNeig; // Maximum number of eigenvalues to deflate
211  mutable RealScalar m_lambdaN; //Modulus of the largest eigenvalue of A
212  mutable bool m_isDeflAllocated;
213  mutable bool m_isDeflInitialized;
214 
215  //Adaptive strategy
216  mutable RealScalar m_smv; // Smaller multiple of the remaining number of steps allowed
217  mutable bool m_force; // Force the use of deflation at each restart
218 
219 };
226 template< typename _MatrixType, typename _Preconditioner>
227 template<typename Rhs, typename Dest>
229  const Preconditioner& precond) const
230 {
231  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
232 
233  RealScalar normRhs = rhs.norm();
234  if(normRhs <= considerAsZero)
235  {
236  x.setZero();
237  m_error = 0;
238  return;
239  }
240 
241  //Initialization
242  m_isDeflInitialized = false;
243  Index n = mat.rows();
244  DenseVector r0(n);
245  Index nbIts = 0;
246  m_H.resize(m_restart+1, m_restart);
247  m_Hes.resize(m_restart, m_restart);
248  m_V.resize(n,m_restart+1);
249  //Initial residual vector and initial norm
250  if(x.squaredNorm()==0)
251  x = precond.solve(rhs);
252  r0 = rhs - mat * x;
253  RealScalar beta = r0.norm();
254 
255  m_error = beta/normRhs;
256  if(m_error < m_tolerance)
257  m_info = Success;
258  else
259  m_info = NoConvergence;
260 
261  // Iterative process
262  while (nbIts < m_iterations && m_info == NoConvergence)
263  {
264  dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts);
265 
266  // Compute the new residual vector for the restart
267  if (nbIts < m_iterations && m_info == NoConvergence) {
268  r0 = rhs - mat * x;
269  beta = r0.norm();
270  }
271  }
272 }
273 
284 template< typename _MatrixType, typename _Preconditioner>
285 template<typename Dest>
286 Index DGMRES<_MatrixType, _Preconditioner>::dgmresCycle(const MatrixType& mat, const Preconditioner& precond, Dest& x, DenseVector& r0, RealScalar& beta, const RealScalar& normRhs, Index& nbIts) const
287 {
288  //Initialization
289  DenseVector g(m_restart+1); // Right hand side of the least square problem
290  g.setZero();
291  g(0) = Scalar(beta);
292  m_V.col(0) = r0/beta;
293  m_info = NoConvergence;
294  std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
295  Index it = 0; // Number of inner iterations
296  Index n = mat.rows();
297  DenseVector tv1(n), tv2(n); //Temporary vectors
298  while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
299  {
300  // Apply preconditioner(s) at right
301  if (m_isDeflInitialized )
302  {
303  dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
304  tv2 = precond.solve(tv1);
305  }
306  else
307  {
308  tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
309  }
310  tv1 = mat * tv2;
311 
312  // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
313  Scalar coef;
314  for (Index i = 0; i <= it; ++i)
315  {
316  coef = tv1.dot(m_V.col(i));
317  tv1 = tv1 - coef * m_V.col(i);
318  m_H(i,it) = coef;
319  m_Hes(i,it) = coef;
320  }
321  // Normalize the vector
322  coef = tv1.norm();
323  m_V.col(it+1) = tv1/coef;
324  m_H(it+1, it) = coef;
325 // m_Hes(it+1,it) = coef;
326 
327  // FIXME Check for happy breakdown
328 
329  // Update Hessenberg matrix with Givens rotations
330  for (Index i = 1; i <= it; ++i)
331  {
332  m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
333  }
334  // Compute the new plane rotation
335  gr[it].makeGivens(m_H(it, it), m_H(it+1,it));
336  // Apply the new rotation
337  m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
338  g.applyOnTheLeft(it,it+1, gr[it].adjoint());
339 
340  beta = std::abs(g(it+1));
341  m_error = beta/normRhs;
342  // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
343  it++; nbIts++;
344 
345  if (m_error < m_tolerance)
346  {
347  // The method has converged
348  m_info = Success;
349  break;
350  }
351  }
352 
353  // Compute the new coefficients by solving the least square problem
354 // it++;
355  //FIXME Check first if the matrix is singular ... zero diagonal
356  DenseVector nrs(m_restart);
357  nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it));
358 
359  // Form the new solution
360  if (m_isDeflInitialized)
361  {
362  tv1 = m_V.leftCols(it) * nrs;
363  dgmresApplyDeflation(tv1, tv2);
364  x = x + precond.solve(tv2);
365  }
366  else
367  x = x + precond.solve(m_V.leftCols(it) * nrs);
368 
369  // Go for a new cycle and compute data for deflation
370  if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
371  dgmresComputeDeflationData(mat, precond, it, m_neig);
372  return 0;
373 
374 }
375 
376 
377 template< typename _MatrixType, typename _Preconditioner>
379 {
380  m_U.resize(rows, m_maxNeig);
381  m_MU.resize(rows, m_maxNeig);
382  m_T.resize(m_maxNeig, m_maxNeig);
383  m_lambdaN = 0.0;
384  m_isDeflAllocated = true;
385 }
386 
387 template< typename _MatrixType, typename _Preconditioner>
389 {
390  return schurofH.matrixT().diagonal();
391 }
392 
393 template< typename _MatrixType, typename _Preconditioner>
395 {
396  const DenseMatrix& T = schurofH.matrixT();
397  Index it = T.rows();
398  ComplexVector eig(it);
399  Index j = 0;
400  while (j < it-1)
401  {
402  if (T(j+1,j) ==Scalar(0))
403  {
404  eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
405  j++;
406  }
407  else
408  {
409  eig(j) = std::complex<RealScalar>(T(j,j),T(j+1,j));
410  eig(j+1) = std::complex<RealScalar>(T(j,j+1),T(j+1,j+1));
411  j++;
412  }
413  }
414  if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
415  return eig;
416 }
417 
418 template< typename _MatrixType, typename _Preconditioner>
420 {
421  // First, find the Schur form of the Hessenberg matrix H
423  bool computeU = true;
424  DenseMatrix matrixQ(it,it);
425  matrixQ.setIdentity();
426  schurofH.computeFromHessenberg(m_Hes.topLeftCorner(it,it), matrixQ, computeU);
427 
428  ComplexVector eig(it);
430  eig = this->schurValues(schurofH);
431 
432  // Reorder the absolute values of Schur values
433  DenseRealVector modulEig(it);
434  for (Index j=0; j<it; ++j) modulEig(j) = std::abs(eig(j));
435  perm.setLinSpaced(it,0,internal::convert_index<StorageIndex>(it-1));
436  internal::sortWithPermutation(modulEig, perm, neig);
437 
438  if (!m_lambdaN)
439  {
440  m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
441  }
442  //Count the real number of extracted eigenvalues (with complex conjugates)
443  Index nbrEig = 0;
444  while (nbrEig < neig)
445  {
446  if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++;
447  else nbrEig += 2;
448  }
449  // Extract the Schur vectors corresponding to the smallest Ritz values
450  DenseMatrix Sr(it, nbrEig);
451  Sr.setZero();
452  for (Index j = 0; j < nbrEig; j++)
453  {
454  Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
455  }
456 
457  // Form the Schur vectors of the initial matrix using the Krylov basis
458  DenseMatrix X;
459  X = m_V.leftCols(it) * Sr;
460  if (m_r)
461  {
462  // Orthogonalize X against m_U using modified Gram-Schmidt
463  for (Index j = 0; j < nbrEig; j++)
464  for (Index k =0; k < m_r; k++)
465  X.col(j) = X.col(j) - (m_U.col(k).dot(X.col(j)))*m_U.col(k);
466  }
467 
468  // Compute m_MX = A * M^-1 * X
469  Index m = m_V.rows();
470  if (!m_isDeflAllocated)
471  dgmresInitDeflation(m);
472  DenseMatrix MX(m, nbrEig);
473  DenseVector tv1(m);
474  for (Index j = 0; j < nbrEig; j++)
475  {
476  tv1 = mat * X.col(j);
477  MX.col(j) = precond.solve(tv1);
478  }
479 
480  //Update m_T = [U'MU U'MX; X'MU X'MX]
481  m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX;
482  if(m_r)
483  {
484  m_T.block(0, m_r, m_r, nbrEig) = m_U.leftCols(m_r).transpose() * MX;
485  m_T.block(m_r, 0, nbrEig, m_r) = X.transpose() * m_MU.leftCols(m_r);
486  }
487 
488  // Save X into m_U and m_MX in m_MU
489  for (Index j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
490  for (Index j = 0; j < nbrEig; j++) m_MU.col(m_r+j) = MX.col(j);
491  // Increase the size of the invariant subspace
492  m_r += nbrEig;
493 
494  // Factorize m_T into m_luT
495  m_luT.compute(m_T.topLeftCorner(m_r, m_r));
496 
497  //FIXME CHeck if the factorization was correctly done (nonsingular matrix)
498  m_isDeflInitialized = true;
499  return 0;
500 }
501 template<typename _MatrixType, typename _Preconditioner>
502 template<typename RhsType, typename DestType>
504 {
505  DenseVector x1 = m_U.leftCols(m_r).transpose() * x;
506  y = x + m_U.leftCols(m_r) * ( m_lambdaN * m_luT.solve(x1) - x1);
507  return 0;
508 }
509 
510 } // end namespace Eigen
511 #endif
Matrix3f m
Index dgmresApplyDeflation(const RhsType &In, DestType &Out) const
Definition: DGMRES.h:503
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define max(a, b)
Definition: datatypes.h:20
StorageIndex m_neig
Definition: DGMRES.h:208
EIGEN_DEVICE_FUNC Derived & setZero(Index size)
ComplexVector schurValues(const ComplexSchur< DenseMatrix > &schurofH) const
Definition: DGMRES.h:388
Scalar * b
Definition: benchVecAdd.cpp:17
void adjoint(const MatrixType &m)
Definition: adjoint.cpp:67
A Restarted GMRES with deflation. This class implements a modification of the GMRES solver for sparse...
Definition: DGMRES.h:19
DenseMatrix m_H
Definition: DGMRES.h:201
Performs a real Schur decomposition of a square matrix.
Definition: RealSchur.h:54
bool m_isDeflAllocated
Definition: DGMRES.h:212
#define min(a, b)
Definition: datatypes.h:19
DenseMatrix m_U
Definition: DGMRES.h:204
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:228
int n
void setEigenv(const Index neig)
Definition: DGMRES.h:168
Index deflSize()
Definition: DGMRES.h:177
LU decomposition of a matrix with partial pivoting, and related features.
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
MatrixType::RealScalar RealScalar
Definition: DGMRES.h:116
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
void g(const string &key, int i)
Definition: testBTree.cpp:41
Index dgmresComputeDeflationData(const MatrixType &mat, const Preconditioner &precond, const Index &it, StorageIndex &neig) const
Definition: DGMRES.h:419
void sortWithPermutation(VectorType &vec, IndexType &perm, typename IndexType::Scalar &ncut)
Computes a permutation vector to have a sorted sequence.
Definition: DGMRES.h:39
Index m_r
Definition: DGMRES.h:209
MatrixType::Scalar Scalar
Definition: DGMRES.h:114
bool m_isDeflInitialized
Definition: DGMRES.h:213
void set_restart(const Index restart)
Definition: DGMRES.h:163
DenseMatrix m_T
Definition: DGMRES.h:206
Matrix< RealScalar, Dynamic, 1 > DenseRealVector
Definition: DGMRES.h:121
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
PartialPivLU< DenseMatrix > m_luT
Definition: DGMRES.h:207
#define eigen_assert(x)
Definition: Macros.h:1037
Eigen::Triplet< double > T
SelfAdjointEigenSolver< PlainMatrixType > eig(mat, computeVectors?ComputeEigenvectors:EigenvaluesOnly)
idx_t idx_t idx_t idx_t idx_t * perm
Matrix< RealScalar, Dynamic, Dynamic > DenseRealMatrix
Definition: DGMRES.h:119
DenseMatrix m_V
Definition: DGMRES.h:200
IterativeSolverBase< DGMRES > Base
Definition: DGMRES.h:103
_Preconditioner Preconditioner
Definition: DGMRES.h:117
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
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:286
Matrix< std::complex< RealScalar >, Dynamic, 1 > ComplexVector
Definition: DGMRES.h:122
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Index restart()
Definition: DGMRES.h:158
Index m_restart
Definition: DGMRES.h:203
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition: RealSchur.h:144
_MatrixType MatrixType
Definition: DGMRES.h:113
bool m_force
Definition: DGMRES.h:217
MatrixType::StorageIndex StorageIndex
Definition: DGMRES.h:115
Pose3 x1
Definition: testPose3.cpp:663
Index m_maxNeig
Definition: DGMRES.h:210
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: DGMRES.h:118
void _solve_vector_with_guess_impl(const Rhs &b, Dest &x) const
Definition: DGMRES.h:145
void setMaxEigenv(const Index maxNeig)
Definition: DGMRES.h:182
DenseMatrix m_Hes
Definition: DGMRES.h:202
RealScalar m_lambdaN
Definition: DGMRES.h:211
RealScalar m_smv
Definition: DGMRES.h:216
const int Dynamic
Definition: Constants.h:22
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Performs a complex Schur decomposition of a real or complex square matrix.
Definition: ComplexSchur.h:51
#define X
Definition: icosphere.cpp:20
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define abs(x)
Definition: datatypes.h:17
Base class for linear iterative solvers.
DenseMatrix m_MU
Definition: DGMRES.h:205
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:162
std::ptrdiff_t j
void dgmresInitDeflation(Index &rows) const
Definition: DGMRES.h:378
Definition: pytypes.h:1370
Matrix< Scalar, Dynamic, 1 > DenseVector
Definition: DGMRES.h:120
DGMRES(const EigenBase< MatrixDerived > &A)
Definition: DGMRES.h:139


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:10