Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
 
Loading...
Searching...
No Matches
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
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
38template <typename VectorType, typename IndexType>
39void 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}
100template< typename _MatrixType, typename _Preconditioner>
101class 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;
116 typedef typename MatrixType::RealScalar RealScalar;
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};
226template< typename _MatrixType, typename _Preconditioner>
227template<typename Rhs, typename Dest>
228void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
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
284template< typename _MatrixType, typename _Preconditioner>
285template<typename Dest>
286Index 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
377template< 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
387template< typename _MatrixType, typename _Preconditioner>
388inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const ComplexSchur<DenseMatrix>& schurofH) const
389{
390 return schurofH.matrixT().diagonal();
391}
392
393template< typename _MatrixType, typename _Preconditioner>
394inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(const RealSchur<DenseMatrix>& schurofH) const
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
418template< typename _MatrixType, typename _Preconditioner>
419Index DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(const MatrixType& mat, const Preconditioner& precond, const Index& it, StorageIndex& neig) const
420{
421 // First, find the Schur form of the Hessenberg matrix H
422 typename internal::conditional<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> >::type schurofH;
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);
429 Matrix<StorageIndex,Dynamic,1>perm(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}
501template<typename _MatrixType, typename _Preconditioner>
502template<typename RhsType, typename DestType>
503Index DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(const RhsType &x, DestType &y) const
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
A Restarted GMRES with deflation. This class implements a modification of the GMRES solver for sparse...
Definition: DGMRES.h:102
DGMRES()
Definition: DGMRES.h:126
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
void setMaxEigenv(const Index maxNeig)
Definition: DGMRES.h:182
DGMRES(const EigenBase< MatrixDerived > &A)
Definition: DGMRES.h:139
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
Index deflSize()
Definition: DGMRES.h:177
void set_restart(const Index restart)
Definition: DGMRES.h:163
void setEigenv(const Index neig)
Definition: DGMRES.h:168
Index restart()
Definition: DGMRES.h:158
Derived & setZero(Index size)
Namespace containing all symbols from the Eigen library.
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
const int Dynamic