Eigen::DGMRES< MatrixType_, Preconditioner_ > Class Template Reference

A Restarted GMRES with deflation. This class implements a modification of the GMRES solver for sparse linear systems. The basis is built with modified Gram-Schmidt. At each restart, a few approximated eigenvectors corresponding to the smallest eigenvalues are used to build a preconditioner for the next cycle. This preconditioner for deflation can be combined with any other preconditioner, the IncompleteLUT for instance. The preconditioner is applied at right of the matrix and the combination is multiplicative. More...

+ Inheritance diagram for Eigen::DGMRES< MatrixType_, Preconditioner_ >:

Public Types

typedef Matrix< std::complex< RealScalar >, Dynamic, 1 > ComplexVector
 
typedef Matrix< Scalar, Dynamic, DynamicDenseMatrix
 
typedef Matrix< RealScalar, Dynamic, DynamicDenseRealMatrix
 
typedef Matrix< RealScalar, Dynamic, 1 > DenseRealVector
 
typedef Matrix< Scalar, Dynamic, 1 > DenseVector
 
typedef MatrixType_ MatrixType
 
typedef Preconditioner_ Preconditioner
 
typedef MatrixType::RealScalar RealScalar
 
typedef MatrixType::Scalar Scalar
 
typedef MatrixType::StorageIndex StorageIndex
 
- Public Types inherited from Eigen::IterativeSolverBase< DGMRES< MatrixType_, Preconditioner_ > >
typedef internal::traits< Derived >::MatrixType MatrixType
 
typedef internal::traits< Derived >::Preconditioner Preconditioner
 
typedef MatrixType::RealScalar RealScalar
 
typedef MatrixType::Scalar Scalar
 
typedef MatrixType::StorageIndex StorageIndex
 

Public Member Functions

void _solve_impl (const Rhs &b, Dest &x) const
 
template<typename Rhs , typename Dest >
void _solve_vector_with_guess_impl (const Rhs &b, Dest &x) const
 
std::enable_if_t< Rhs::ColsAtCompileTime!=1 &&DestDerived::ColsAtCompileTime!=1 > _solve_with_guess_impl (const Rhs &b, MatrixBase< DestDerived > &aDest) const
 
std::enable_if_t< Rhs::ColsAtCompileTime==1||DestDerived::ColsAtCompileTime==1 > _solve_with_guess_impl (const Rhs &b, MatrixBase< DestDerived > &dest) const
 
void _solve_with_guess_impl (const Rhs &b, SparseMatrixBase< DestDerived > &aDest) const
 
Index deflSize ()
 
 DGMRES ()
 
template<typename MatrixDerived >
 DGMRES (const EigenBase< MatrixDerived > &A)
 
Index restart ()
 
void set_restart (const Index restart)
 
void setEigenv (const Index neig)
 
void setMaxEigenv (const Index maxNeig)
 
 ~DGMRES ()
 
- Public Member Functions inherited from Eigen::IterativeSolverBase< DGMRES< MatrixType_, Preconditioner_ > >
void _solve_impl (const Rhs &b, Dest &x) const
 
std::enable_if_t< Rhs::ColsAtCompileTime!=1 &&DestDerived::ColsAtCompileTime!=1 > _solve_with_guess_impl (const Rhs &b, MatrixBase< DestDerived > &aDest) const
 
std::enable_if_t< Rhs::ColsAtCompileTime==1||DestDerived::ColsAtCompileTime==1 > _solve_with_guess_impl (const Rhs &b, MatrixBase< DestDerived > &dest) const
 
void _solve_with_guess_impl (const Rhs &b, SparseMatrixBase< DestDerived > &aDest) const
 
Derived & analyzePattern (const EigenBase< MatrixDerived > &A)
 
EIGEN_CONSTEXPR Index cols () const EIGEN_NOEXCEPT
 
Derived & compute (const EigenBase< MatrixDerived > &A)
 
Derived & derived ()
 
const Derived & derived () const
 
RealScalar error () const
 
Derived & factorize (const EigenBase< MatrixDerived > &A)
 
ComputationInfo info () const
 
Index iterations () const
 
 IterativeSolverBase ()
 
 IterativeSolverBase (const EigenBase< MatrixDerived > &A)
 
 IterativeSolverBase (IterativeSolverBase &&)=default
 
Index maxIterations () const
 
Preconditionerpreconditioner ()
 
const Preconditionerpreconditioner () const
 
EIGEN_CONSTEXPR Index rows () const EIGEN_NOEXCEPT
 
Derived & setMaxIterations (Index maxIters)
 
Derived & setTolerance (const RealScalar &tolerance)
 
const SolveWithGuess< Derived, Rhs, Guess > solveWithGuess (const MatrixBase< Rhs > &b, const Guess &x0) const
 
RealScalar tolerance () const
 
 ~IterativeSolverBase ()
 
- Public Member Functions inherited from Eigen::SparseSolverBase< class >
Derived & derived ()
 
const Derived & derived () const
 
const Solve< Derived, Rhs > solve (const MatrixBase< Rhs > &b) const
 
const Solve< Derived, Rhs > solve (const SparseMatrixBase< Rhs > &b) const
 
 SparseSolverBase ()
 
 SparseSolverBase (SparseSolverBase &&other)
 
 ~SparseSolverBase ()
 

Protected Member Functions

template<typename Rhs , typename Dest >
void dgmres (const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond) const
 Perform several cycles of restarted GMRES with modified Gram Schmidt,. More...
 
template<typename RhsType , typename DestType >
Index dgmresApplyDeflation (const RhsType &In, DestType &Out) const
 
Index dgmresComputeDeflationData (const MatrixType &mat, const Preconditioner &precond, const Index &it, StorageIndex &neig) const
 
template<typename Dest >
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. More...
 
void dgmresInitDeflation (Index &rows) const
 
ComplexVector schurValues (const ComplexSchur< DenseMatrix > &schurofH) const
 
ComplexVector schurValues (const RealSchur< DenseMatrix > &schurofH) const
 
- Protected Member Functions inherited from Eigen::IterativeSolverBase< DGMRES< MatrixType_, Preconditioner_ > >
void grab (const InputType &A)
 
void init ()
 
const ActualMatrixTypematrix () const
 

Protected Attributes

bool m_force
 
DenseMatrix m_H
 
DenseMatrix m_Hes
 
bool m_isDeflAllocated
 
bool m_isDeflInitialized
 
RealScalar m_lambdaN
 
PartialPivLU< DenseMatrixm_luT
 
Index m_maxNeig
 
DenseMatrix m_MU
 
StorageIndex m_neig
 
Index m_r
 
Index m_restart
 
RealScalar m_smv
 
DenseMatrix m_T
 
DenseMatrix m_U
 
DenseMatrix m_V
 
- Protected Attributes inherited from Eigen::IterativeSolverBase< DGMRES< MatrixType_, Preconditioner_ > >
bool m_analysisIsOk
 
RealScalar m_error
 
bool m_factorizationIsOk
 
ComputationInfo m_info
 
bool m_isInitialized
 
Index m_iterations
 
MatrixWrapper m_matrixWrapper
 
Index m_maxIterations
 
Preconditioner m_preconditioner
 
RealScalar m_tolerance
 
- Protected Attributes inherited from Eigen::SparseSolverBase< class >
bool m_isInitialized
 

Private Types

typedef IterativeSolverBase< DGMRESBase
 

Private Member Functions

const ActualMatrixTypematrix () const
 

Private Attributes

RealScalar m_error
 
ComputationInfo m_info
 
bool m_isInitialized
 
Index m_iterations
 
RealScalar m_tolerance
 

Additional Inherited Members

- Public Attributes inherited from Eigen::IterativeSolverBase< DGMRES< MatrixType_, Preconditioner_ > >
 ColsAtCompileTime
 
 MaxColsAtCompileTime
 
- Protected Types inherited from Eigen::IterativeSolverBase< DGMRES< MatrixType_, Preconditioner_ > >
typedef MatrixWrapper::ActualMatrixType ActualMatrixType
 
typedef SparseSolverBase< Derived > Base
 
typedef internal::generic_matrix_wrapper< MatrixTypeMatrixWrapper
 

Detailed Description

template<typename MatrixType_, typename Preconditioner_>
class Eigen::DGMRES< MatrixType_, Preconditioner_ >

A Restarted GMRES with deflation. This class implements a modification of the GMRES solver for sparse linear systems. The basis is built with modified Gram-Schmidt. At each restart, a few approximated eigenvectors corresponding to the smallest eigenvalues are used to build a preconditioner for the next cycle. This preconditioner for deflation can be combined with any other preconditioner, the IncompleteLUT for instance. The preconditioner is applied at right of the matrix and the combination is multiplicative.

Template Parameters
MatrixType_the type of the sparse matrix A, can be a dense or a sparse matrix.
Preconditioner_the type of the preconditioner. Default is DiagonalPreconditioner Typical usage :
//Fill A and b ...
DGMRES<SparseMatrix<double> > solver;
solver.set_restart(30); // Set restarting value
solver.setEigenv(1); // Set the number of eigenvalues to deflate
solver.compute(A);
x = solver.solve(b);
BiCGSTAB< SparseMatrix< double > > solver
SparseMatrix< double > A(n, n)
Matrix< double, Dynamic, 1 > VectorXd

DGMRES can also be used in a matrix-free context, see the following example .

References : [1] D. NUENTSA WAKAM and F. PACULL, Memory Efficient Hybrid Algebraic Solvers for Linear Systems Arising from Compressible Flows, Computers and Fluids, In Press, https://doi.org/10.1016/j.compfluid.2012.03.023
[2] K. Burrage and J. Erhel, On the performance of various adaptive preconditioned GMRES strategies, 5(1998), 101-121. [3] J. Erhel, K. Burrage and B. Pohl, Restarted GMRES preconditioned by deflation,J. Computational and Applied Mathematics, 69(1996), 303-318.

Definition at line 103 of file DGMRES.h.

Member Typedef Documentation

◆ Base

template<typename MatrixType_ , typename Preconditioner_ >
typedef IterativeSolverBase<DGMRES> Eigen::DGMRES< MatrixType_, Preconditioner_ >::Base
private

Definition at line 105 of file DGMRES.h.

◆ ComplexVector

template<typename MatrixType_ , typename Preconditioner_ >
typedef Matrix<std::complex<RealScalar>, Dynamic, 1> Eigen::DGMRES< MatrixType_, Preconditioner_ >::ComplexVector

Definition at line 124 of file DGMRES.h.

◆ DenseMatrix

template<typename MatrixType_ , typename Preconditioner_ >
typedef Matrix<Scalar,Dynamic,Dynamic> Eigen::DGMRES< MatrixType_, Preconditioner_ >::DenseMatrix

Definition at line 120 of file DGMRES.h.

◆ DenseRealMatrix

template<typename MatrixType_ , typename Preconditioner_ >
typedef Matrix<RealScalar,Dynamic,Dynamic> Eigen::DGMRES< MatrixType_, Preconditioner_ >::DenseRealMatrix

Definition at line 121 of file DGMRES.h.

◆ DenseRealVector

template<typename MatrixType_ , typename Preconditioner_ >
typedef Matrix<RealScalar,Dynamic,1> Eigen::DGMRES< MatrixType_, Preconditioner_ >::DenseRealVector

Definition at line 123 of file DGMRES.h.

◆ DenseVector

template<typename MatrixType_ , typename Preconditioner_ >
typedef Matrix<Scalar,Dynamic,1> Eigen::DGMRES< MatrixType_, Preconditioner_ >::DenseVector

Definition at line 122 of file DGMRES.h.

◆ MatrixType

template<typename MatrixType_ , typename Preconditioner_ >
typedef MatrixType_ Eigen::DGMRES< MatrixType_, Preconditioner_ >::MatrixType

Definition at line 115 of file DGMRES.h.

◆ Preconditioner

template<typename MatrixType_ , typename Preconditioner_ >
typedef Preconditioner_ Eigen::DGMRES< MatrixType_, Preconditioner_ >::Preconditioner

Definition at line 119 of file DGMRES.h.

◆ RealScalar

template<typename MatrixType_ , typename Preconditioner_ >
typedef MatrixType::RealScalar Eigen::DGMRES< MatrixType_, Preconditioner_ >::RealScalar

Definition at line 118 of file DGMRES.h.

◆ Scalar

template<typename MatrixType_ , typename Preconditioner_ >
typedef MatrixType::Scalar Eigen::DGMRES< MatrixType_, Preconditioner_ >::Scalar

Definition at line 116 of file DGMRES.h.

◆ StorageIndex

template<typename MatrixType_ , typename Preconditioner_ >
typedef MatrixType::StorageIndex Eigen::DGMRES< MatrixType_, Preconditioner_ >::StorageIndex

Definition at line 117 of file DGMRES.h.

Constructor & Destructor Documentation

◆ DGMRES() [1/2]

template<typename MatrixType_ , typename Preconditioner_ >
Eigen::DGMRES< MatrixType_, Preconditioner_ >::DGMRES ( )
inline

Default constructor.

Definition at line 128 of file DGMRES.h.

128 : Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
IterativeSolverBase< DGMRES > Base
Definition: DGMRES.h:105
StorageIndex m_neig
Definition: DGMRES.h:210
bool m_isDeflAllocated
Definition: DGMRES.h:214
Index m_restart
Definition: DGMRES.h:205
Index m_maxNeig
Definition: DGMRES.h:212
Index m_r
Definition: DGMRES.h:211
bool m_isDeflInitialized
Definition: DGMRES.h:215

◆ DGMRES() [2/2]

template<typename MatrixType_ , typename Preconditioner_ >
template<typename MatrixDerived >
Eigen::DGMRES< MatrixType_, Preconditioner_ >::DGMRES ( const EigenBase< MatrixDerived > &  A)
inlineexplicit

Initialize the solver with matrix A for further Ax=b solving.

This constructor is a shortcut for the default constructor followed by a call to compute().

Warning
this class stores a reference to the matrix A as well as some precomputed values that depend on it. Therefore, if A is changed this class becomes invalid. Call compute() to update it with the new matrix A, or modify a copy of A.

Definition at line 141 of file DGMRES.h.

141 : Base(A.derived()), m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
SparseMatrix< double, Options_, StorageIndex_ > & derived()

◆ ~DGMRES()

template<typename MatrixType_ , typename Preconditioner_ >
Eigen::DGMRES< MatrixType_, Preconditioner_ >::~DGMRES ( )
inline

Definition at line 143 of file DGMRES.h.

143 {}

Member Function Documentation

◆ _solve_impl()

template<typename MatrixType_ , typename Preconditioner_ >
void Eigen::IterativeSolverBase< class >::_solve_impl

◆ _solve_vector_with_guess_impl()

template<typename MatrixType_ , typename Preconditioner_ >
template<typename Rhs , typename Dest >
void Eigen::DGMRES< MatrixType_, Preconditioner_ >::_solve_vector_with_guess_impl ( const Rhs &  b,
Dest &  x 
) const
inline

Definition at line 147 of file DGMRES.h.

148  {
149  EIGEN_STATIC_ASSERT(Rhs::ColsAtCompileTime==1 || Dest::ColsAtCompileTime==1, YOU_TRIED_CALLING_A_VECTOR_METHOD_ON_A_MATRIX);
150 
153 
155  }
#define EIGEN_STATIC_ASSERT(X, MSG)
RealScalar m_error
const ActualMatrixType & matrix() const
Index m_iterations
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:230
Preconditioner m_preconditioner

◆ _solve_with_guess_impl() [1/3]

template<typename MatrixType_ , typename Preconditioner_ >
std::enable_if_t< Rhs::ColsAtCompileTime!=1 &&DestDerived::ColsAtCompileTime!=1 > Eigen::IterativeSolverBase< class >::_solve_with_guess_impl

◆ _solve_with_guess_impl() [2/3]

template<typename MatrixType_ , typename Preconditioner_ >
std::enable_if_t< Rhs::ColsAtCompileTime==1||DestDerived::ColsAtCompileTime==1 > Eigen::IterativeSolverBase< class >::_solve_with_guess_impl

◆ _solve_with_guess_impl() [3/3]

template<typename MatrixType_ , typename Preconditioner_ >
void Eigen::IterativeSolverBase< class >::_solve_with_guess_impl

◆ deflSize()

template<typename MatrixType_ , typename Preconditioner_ >
Index Eigen::DGMRES< MatrixType_, Preconditioner_ >::deflSize ( )
inline

Get the size of the deflation subspace size

Definition at line 179 of file DGMRES.h.

179 {return m_r; }

◆ dgmres()

template<typename MatrixType_ , typename Preconditioner_ >
template<typename Rhs , typename Dest >
void Eigen::DGMRES< MatrixType_, Preconditioner_ >::dgmres ( const MatrixType mat,
const Rhs &  rhs,
Dest &  x,
const Preconditioner precond 
) const
protected

Perform several cycles of restarted GMRES with modified Gram Schmidt,.

A right preconditioner is used combined with deflation.

Definition at line 230 of file DGMRES.h.

232 {
233  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
234 
235  RealScalar normRhs = rhs.norm();
236  if(normRhs <= considerAsZero)
237  {
238  x.setZero();
239  m_error = 0;
240  return;
241  }
242 
243  //Initialization
244  m_isDeflInitialized = false;
245  Index n = mat.rows();
246  DenseVector r0(n);
247  Index nbIts = 0;
250  m_V.resize(n,m_restart+1);
251  //Initial residual vector and initial norm
252  if(x.squaredNorm()==0)
253  x = precond.solve(rhs);
254  r0 = rhs - mat * x;
255  RealScalar beta = r0.norm();
256 
257  m_error = beta/normRhs;
258  if(m_error < m_tolerance)
259  m_info = Success;
260  else
262 
263  // Iterative process
264  while (nbIts < m_iterations && m_info == NoConvergence)
265  {
266  dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts);
267 
268  // Compute the new residual vector for the restart
269  if (nbIts < m_iterations && m_info == NoConvergence) {
270  r0 = rhs - mat * x;
271  beta = r0.norm();
272  }
273  }
274 }
int n
MatrixXf mat
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:288
RealScalar m_tolerance
ComputationInfo m_info
Matrix< Scalar, Dynamic, 1 > DenseVector
Definition: DGMRES.h:122
DenseMatrix m_H
Definition: DGMRES.h:203
DenseMatrix m_V
Definition: DGMRES.h:202
DenseMatrix m_Hes
Definition: DGMRES.h:204
MatrixType::RealScalar RealScalar
Definition: DGMRES.h:118
constexpr void resize(Index rows, Index cols)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
CleanedUpDerType< DerType >::type() min(const AutoDiffScalar< DerType > &x, const T &y)

◆ dgmresApplyDeflation()

template<typename MatrixType_ , typename Preconditioner_ >
template<typename RhsType , typename DestType >
Index Eigen::DGMRES< MatrixType_, Preconditioner_ >::dgmresApplyDeflation ( const RhsType &  In,
DestType &  Out 
) const
protected

Definition at line 505 of file DGMRES.h.

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 }
Matrix3f y
PartialPivLU< DenseMatrix > m_luT
Definition: DGMRES.h:209
DenseMatrix m_U
Definition: DGMRES.h:206
RealScalar m_lambdaN
Definition: DGMRES.h:213

◆ dgmresComputeDeflationData()

template<typename MatrixType_ , typename Preconditioner_ >
Index Eigen::DGMRES< MatrixType_, Preconditioner_ >::dgmresComputeDeflationData ( const MatrixType mat,
const Preconditioner precond,
const Index it,
StorageIndex neig 
) const
protected

Definition at line 421 of file DGMRES.h.

422 {
423  // First, find the Schur form of the Hessenberg matrix H
424  std::conditional_t<NumTraits<Scalar>::IsComplex, ComplexSchur<DenseMatrix>, RealSchur<DenseMatrix> > 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);
432  eig = this->schurValues(schurofH);
433 
434  // Reorder the absolute values of Schur values
435  DenseRealVector modulEig(it);
436  for (Index j=0; j<it; ++j) modulEig(j) = std::abs(eig(j));
437  perm.setLinSpaced(it,0,internal::convert_index<StorageIndex>(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  Index 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 (Index 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 (Index j = 0; j < nbrEig; j++)
466  for (Index 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)
474  DenseMatrix MX(m, nbrEig);
475  DenseVector tv1(m);
476  for (Index 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 (Index j = 0; j < nbrEig; j++) m_U.col(m_r+j) = X.col(j);
492  for (Index 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 }
Matrix3f m
MatrixXf X
void dgmresInitDeflation(Index &rows) const
Definition: DGMRES.h:380
Matrix< std::complex< RealScalar >, Dynamic, 1 > ComplexVector
Definition: DGMRES.h:124
Matrix< RealScalar, Dynamic, 1 > DenseRealVector
Definition: DGMRES.h:123
DenseMatrix m_T
Definition: DGMRES.h:208
DenseMatrix m_MU
Definition: DGMRES.h:207
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Definition: DGMRES.h:120
ComplexVector schurValues(const ComplexSchur< DenseMatrix > &schurofH) const
Definition: DGMRES.h:390
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
void sortWithPermutation(VectorType &vec, IndexType &perm, typename IndexType::Scalar &ncut)
Computes a permutation vector to have a sorted sequence.
Definition: DGMRES.h:41
CleanedUpDerType< DerType >::type() max(const AutoDiffScalar< DerType > &x, const T &y)
adouble abs(const adouble &x)
Definition: AdolcForward:74
std::ptrdiff_t j

◆ dgmresCycle()

template<typename MatrixType_ , typename Preconditioner_ >
template<typename Dest >
Index Eigen::DGMRES< MatrixType_, Preconditioner_ >::dgmresCycle ( const MatrixType mat,
const Preconditioner precond,
Dest &  x,
DenseVector r0,
RealScalar beta,
const RealScalar normRhs,
Index nbIts 
) const
protected

Perform one restart cycle of DGMRES.

Parameters
matThe coefficient matrix
precondThe preconditioner
xthe new approximated solution
r0The initial residual vector
betaThe norm of the residual computed so far
normRhsThe norm of the right hand side vector
nbItsThe number of iterations

Definition at line 288 of file DGMRES.h.

289 {
290  //Initialization
291  DenseVector g(m_restart+1); // Right hand side of the least square problem
292  g.setZero();
293  g(0) = Scalar(beta);
294  m_V.col(0) = r0/beta;
296  std::vector<JacobiRotation<Scalar> >gr(m_restart); // Givens rotations
297  Index it = 0; // Number of inner iterations
298  Index n = mat.rows();
299  DenseVector tv1(n), tv2(n); //Temporary vectors
300  while (m_info == NoConvergence && it < m_restart && nbIts < m_iterations)
301  {
302  // Apply preconditioner(s) at right
303  if (m_isDeflInitialized )
304  {
305  dgmresApplyDeflation(m_V.col(it), tv1); // Deflation
306  tv2 = precond.solve(tv1);
307  }
308  else
309  {
310  tv2 = precond.solve(m_V.col(it)); // User's selected preconditioner
311  }
312  tv1 = mat * tv2;
313 
314  // Orthogonalize it with the previous basis in the basis using modified Gram-Schmidt
315  Scalar coef;
316  for (Index i = 0; i <= it; ++i)
317  {
318  coef = tv1.dot(m_V.col(i));
319  tv1 = tv1 - coef * m_V.col(i);
320  m_H(i,it) = coef;
321  m_Hes(i,it) = coef;
322  }
323  // Normalize the vector
324  coef = tv1.norm();
325  m_V.col(it+1) = tv1/coef;
326  m_H(it+1, it) = coef;
327 // m_Hes(it+1,it) = coef;
328 
329  // FIXME Check for happy breakdown
330 
331  // Update Hessenberg matrix with Givens rotations
332  for (Index i = 1; i <= it; ++i)
333  {
334  m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
335  }
336  // Compute the new plane rotation
337  gr[it].makeGivens(m_H(it, it), m_H(it+1,it));
338  // Apply the new rotation
339  m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
340  g.applyOnTheLeft(it,it+1, gr[it].adjoint());
341 
342  beta = std::abs(g(it+1));
343  m_error = beta/normRhs;
344  // std::cerr << nbIts << " Relative Residual Norm " << m_error << std::endl;
345  it++; nbIts++;
346 
347  if (m_error < m_tolerance)
348  {
349  // The method has converged
350  m_info = Success;
351  break;
352  }
353  }
354 
355  // Compute the new coefficients by solving the least square problem
356 // it++;
357  //FIXME Check first if the matrix is singular ... zero diagonal
358  DenseVector nrs(m_restart);
359  nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it));
360 
361  // Form the new solution
363  {
364  tv1 = m_V.leftCols(it) * nrs;
365  dgmresApplyDeflation(tv1, tv2);
366  x = x + precond.solve(tv2);
367  }
368  else
369  x = x + precond.solve(m_V.leftCols(it) * nrs);
370 
371  // Go for a new cycle and compute data for deflation
372  if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
373  dgmresComputeDeflationData(mat, precond, it, m_neig);
374  return 0;
375 
376 }
int i
Index dgmresApplyDeflation(const RhsType &In, DestType &Out) const
Definition: DGMRES.h:505
MatrixType::Scalar Scalar
Definition: DGMRES.h:116
Index dgmresComputeDeflationData(const MatrixType &mat, const Preconditioner &precond, const Index &it, StorageIndex &neig) const
Definition: DGMRES.h:421

◆ dgmresInitDeflation()

template<typename MatrixType_ , typename Preconditioner_ >
void Eigen::DGMRES< MatrixType_, Preconditioner_ >::dgmresInitDeflation ( Index rows) const
protected

Definition at line 380 of file DGMRES.h.

381 {
385  m_lambdaN = 0.0;
386  m_isDeflAllocated = true;
387 }

◆ matrix()

template<typename MatrixType_ , typename Preconditioner_ >
const ActualMatrixType & Eigen::IterativeSolverBase< class >::matrix
private

◆ restart()

template<typename MatrixType_ , typename Preconditioner_ >
Index Eigen::DGMRES< MatrixType_, Preconditioner_ >::restart ( )
inline

Get the restart value

Definition at line 160 of file DGMRES.h.

160 { return m_restart; }

◆ schurValues() [1/2]

template<typename MatrixType_ , typename Preconditioner_ >
DGMRES< MatrixType_, Preconditioner_ >::ComplexVector Eigen::DGMRES< MatrixType_, Preconditioner_ >::schurValues ( const ComplexSchur< DenseMatrix > &  schurofH) const
inlineprotected

Definition at line 390 of file DGMRES.h.

391 {
392  return schurofH.matrixT().diagonal();
393 }

◆ schurValues() [2/2]

template<typename MatrixType_ , typename Preconditioner_ >
DGMRES< MatrixType_, Preconditioner_ >::ComplexVector Eigen::DGMRES< MatrixType_, Preconditioner_ >::schurValues ( const RealSchur< DenseMatrix > &  schurofH) const
inlineprotected

Definition at line 396 of file DGMRES.h.

397 {
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 }

◆ set_restart()

template<typename MatrixType_ , typename Preconditioner_ >
void Eigen::DGMRES< MatrixType_, Preconditioner_ >::set_restart ( const Index  restart)
inline

Set the restart value (default is 30)

Definition at line 165 of file DGMRES.h.

165 { m_restart=restart; }
Index restart()
Definition: DGMRES.h:160

◆ setEigenv()

template<typename MatrixType_ , typename Preconditioner_ >
void Eigen::DGMRES< MatrixType_, Preconditioner_ >::setEigenv ( const Index  neig)
inline

Set the number of eigenvalues to deflate at each restart

Definition at line 170 of file DGMRES.h.

171  {
172  m_neig = neig;
173  if (neig+1 > m_maxNeig) m_maxNeig = neig+1; // To allow for complex conjugates
174  }

◆ setMaxEigenv()

template<typename MatrixType_ , typename Preconditioner_ >
void Eigen::DGMRES< MatrixType_, Preconditioner_ >::setMaxEigenv ( const Index  maxNeig)
inline

Set the maximum size of the deflation subspace

Definition at line 184 of file DGMRES.h.

184 { m_maxNeig = maxNeig; }

Member Data Documentation

◆ m_error

template<typename MatrixType_ , typename Preconditioner_ >
RealScalar Eigen::IterativeSolverBase< class >::m_error
private

◆ m_force

template<typename MatrixType_ , typename Preconditioner_ >
bool Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_force
mutableprotected

Definition at line 219 of file DGMRES.h.

◆ m_H

template<typename MatrixType_ , typename Preconditioner_ >
DenseMatrix Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_H
mutableprotected

Definition at line 203 of file DGMRES.h.

◆ m_Hes

template<typename MatrixType_ , typename Preconditioner_ >
DenseMatrix Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_Hes
mutableprotected

Definition at line 204 of file DGMRES.h.

◆ m_info

template<typename MatrixType_ , typename Preconditioner_ >
ComputationInfo Eigen::IterativeSolverBase< class >::m_info
private

◆ m_isDeflAllocated

template<typename MatrixType_ , typename Preconditioner_ >
bool Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_isDeflAllocated
mutableprotected

Definition at line 214 of file DGMRES.h.

◆ m_isDeflInitialized

template<typename MatrixType_ , typename Preconditioner_ >
bool Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_isDeflInitialized
mutableprotected

Definition at line 215 of file DGMRES.h.

◆ m_isInitialized

template<typename MatrixType_ , typename Preconditioner_ >
bool Eigen::IterativeSolverBase< class >::m_isInitialized
private

◆ m_iterations

template<typename MatrixType_ , typename Preconditioner_ >
Index Eigen::IterativeSolverBase< class >::m_iterations
private

◆ m_lambdaN

template<typename MatrixType_ , typename Preconditioner_ >
RealScalar Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_lambdaN
mutableprotected

Definition at line 213 of file DGMRES.h.

◆ m_luT

template<typename MatrixType_ , typename Preconditioner_ >
PartialPivLU<DenseMatrix> Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_luT
mutableprotected

Definition at line 209 of file DGMRES.h.

◆ m_maxNeig

template<typename MatrixType_ , typename Preconditioner_ >
Index Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_maxNeig
mutableprotected

Definition at line 212 of file DGMRES.h.

◆ m_MU

template<typename MatrixType_ , typename Preconditioner_ >
DenseMatrix Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_MU
mutableprotected

Definition at line 207 of file DGMRES.h.

◆ m_neig

template<typename MatrixType_ , typename Preconditioner_ >
StorageIndex Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_neig
mutableprotected

Definition at line 210 of file DGMRES.h.

◆ m_r

template<typename MatrixType_ , typename Preconditioner_ >
Index Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_r
mutableprotected

Definition at line 211 of file DGMRES.h.

◆ m_restart

template<typename MatrixType_ , typename Preconditioner_ >
Index Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_restart
mutableprotected

Definition at line 205 of file DGMRES.h.

◆ m_smv

template<typename MatrixType_ , typename Preconditioner_ >
RealScalar Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_smv
mutableprotected

Definition at line 218 of file DGMRES.h.

◆ m_T

template<typename MatrixType_ , typename Preconditioner_ >
DenseMatrix Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_T
mutableprotected

Definition at line 208 of file DGMRES.h.

◆ m_tolerance

template<typename MatrixType_ , typename Preconditioner_ >
RealScalar Eigen::IterativeSolverBase< class >::m_tolerance
private

◆ m_U

template<typename MatrixType_ , typename Preconditioner_ >
DenseMatrix Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_U
mutableprotected

Definition at line 206 of file DGMRES.h.

◆ m_V

template<typename MatrixType_ , typename Preconditioner_ >
DenseMatrix Eigen::DGMRES< MatrixType_, Preconditioner_ >::m_V
mutableprotected

Definition at line 202 of file DGMRES.h.


The documentation for this class was generated from the following file: