10 #ifndef EIGEN_DGMRES_H
11 #define EIGEN_DGMRES_H
13 #include <Eigen/Eigenvalues>
17 template<
typename _MatrixType,
18 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
23 template<
typename _MatrixType,
typename _Preconditioner>
24 struct traits<
DGMRES<_MatrixType,_Preconditioner> >
26 typedef _MatrixType MatrixType;
27 typedef _Preconditioner Preconditioner;
38 template <
typename VectorType,
typename IndexType>
39 void sortWithPermutation (VectorType& vec, IndexType& perm,
typename IndexType::Scalar& ncut)
41 eigen_assert(vec.size() == perm.size());
43 for (
Index k = 0; k < ncut; k++)
46 for (
Index j = 0; j < vec.size()-1; j++)
48 if ( vec(perm(j)) < vec(perm(j+1)) )
50 std::swap(perm(j),perm(j+1));
100 template<
typename _MatrixType,
typename _Preconditioner>
101 class DGMRES :
public IterativeSolverBase<DGMRES<_MatrixType,_Preconditioner> >
103 typedef IterativeSolverBase<DGMRES> Base;
106 using Base::m_iterations;
108 using Base::m_isInitialized;
109 using Base::m_tolerance;
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;
125 DGMRES() :
Base(),m_restart(30),m_neig(0),m_r(0),m_maxNeig(5),m_isDeflAllocated(false),m_isDeflInitialized(false) {}
137 template<
typename MatrixDerived>
143 template<
typename Rhs,
typename Dest>
144 void _solve_with_guess_impl(
const Rhs& b, Dest& x)
const
147 for(
Index j=0; j<b.cols(); ++j)
150 m_error = Base::m_tolerance;
152 typename Dest::ColXpr xj(x,j);
153 dgmres(matrix(), b.col(j), xj, Base::m_preconditioner);
156 : m_error <= Base::m_tolerance ?
Success
158 m_isInitialized =
true;
162 template<
typename Rhs,
typename Dest>
163 void _solve_impl(
const Rhs& b, MatrixBase<Dest>& x)
const
166 _solve_with_guess_impl(b,x.derived());
184 if (neig+1 > m_maxNeig) m_maxNeig = neig+1;
199 template<
typename Rhs,
typename Dest>
200 void dgmres(
const MatrixType& mat,
const Rhs& rhs, Dest& x,
const Preconditioner& precond)
const;
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;
205 Index dgmresComputeDeflationData(
const MatrixType& mat,
const Preconditioner& precond,
const Index& it, StorageIndex& neig)
const;
207 template<
typename RhsType,
typename DestType>
208 Index dgmresApplyDeflation(
const RhsType& In, DestType& Out)
const;
212 void dgmresInitDeflation(
Index& rows)
const;
213 mutable DenseMatrix m_V;
214 mutable DenseMatrix m_H;
215 mutable DenseMatrix m_Hes;
216 mutable Index m_restart;
217 mutable DenseMatrix m_U;
218 mutable DenseMatrix m_MU;
219 mutable DenseMatrix m_T;
221 mutable StorageIndex m_neig;
223 mutable Index m_maxNeig;
224 mutable RealScalar m_lambdaN;
225 mutable bool m_isDeflAllocated;
226 mutable bool m_isDeflInitialized;
229 mutable RealScalar m_smv;
230 mutable bool m_force;
239 template<
typename _MatrixType,
typename _Preconditioner>
240 template<
typename Rhs,
typename Dest>
242 const Preconditioner& precond)
const
245 Index n = mat.rows();
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);
252 x = precond.solve(x);
254 RealScalar beta = r0.
norm();
255 RealScalar normRhs = rhs.norm();
256 m_error = beta/normRhs;
257 if(m_error < m_tolerance)
265 dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts);
283 template<
typename _MatrixType,
typename _Preconditioner>
284 template<
typename Dest>
291 m_V.col(0) = r0/beta;
293 std::vector<JacobiRotation<Scalar> >gr(m_restart);
295 Index n = mat.rows();
297 while (m_info ==
NoConvergence && it < m_restart && nbIts < m_iterations)
300 if (m_isDeflInitialized )
302 dgmresApplyDeflation(m_V.col(it), tv1);
303 tv2 = precond.solve(tv1);
307 tv2 = precond.solve(m_V.col(it));
313 for (
Index i = 0; i <= it; ++i)
315 coef = tv1.
dot(m_V.col(i));
316 tv1 = tv1 - coef * m_V.col(i);
322 m_V.col(it+1) = tv1/coef;
323 m_H(it+1, it) = coef;
329 for (
Index i = 1; i <= it; ++i)
331 m_H.col(it).applyOnTheLeft(i-1,i,gr[i-1].adjoint());
334 gr[it].makeGivens(m_H(it, it), m_H(it+1,it));
336 m_H.col(it).applyOnTheLeft(it,it+1,gr[it].adjoint());
339 beta = std::abs(g(it+1));
340 m_error = beta/normRhs;
344 if (m_error < m_tolerance)
356 nrs = m_H.topLeftCorner(it,it).template triangularView<Upper>().solve(g.head(it));
359 if (m_isDeflInitialized)
361 tv1 = m_V.leftCols(it) * nrs;
362 dgmresApplyDeflation(tv1, tv2);
363 x = x + precond.solve(tv2);
366 x = x + precond.solve(m_V.leftCols(it) * nrs);
369 if(nbIts < m_iterations && m_info == NoConvergence && m_neig > 0 && (m_r+m_neig) < m_maxNeig)
370 dgmresComputeDeflationData(mat, precond, it, m_neig);
376 template<
typename _MatrixType,
typename _Preconditioner>
379 m_U.resize(rows, m_maxNeig);
380 m_MU.resize(rows, m_maxNeig);
381 m_T.resize(m_maxNeig, m_maxNeig);
383 m_isDeflAllocated =
true;
386 template<
typename _MatrixType,
typename _Preconditioner>
387 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(
const ComplexSchur<DenseMatrix>& schurofH)
const
389 return schurofH.matrixT().diagonal();
392 template<
typename _MatrixType,
typename _Preconditioner>
393 inline typename DGMRES<_MatrixType, _Preconditioner>::ComplexVector DGMRES<_MatrixType, _Preconditioner>::schurValues(
const RealSchur<DenseMatrix>& schurofH)
const
395 const DenseMatrix& T = schurofH.matrixT();
397 ComplexVector eig(it);
401 if (T(j+1,j) ==Scalar(0))
403 eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
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));
413 if (j < it-1) eig(j) = std::complex<RealScalar>(T(j,j),RealScalar(0));
417 template<
typename _MatrixType,
typename _Preconditioner>
418 Index DGMRES<_MatrixType, _Preconditioner>::dgmresComputeDeflationData(
const MatrixType& mat,
const Preconditioner& precond,
const Index& it, StorageIndex& neig)
const
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);
427 ComplexVector eig(it);
428 Matrix<StorageIndex,Dynamic,1>perm(it);
429 eig = this->schurValues(schurofH);
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);
439 m_lambdaN = (std::max)(modulEig.maxCoeff(), m_lambdaN);
443 while (nbrEig < neig)
445 if(eig(perm(it-nbrEig-1)).imag() == RealScalar(0)) nbrEig++;
449 DenseMatrix Sr(it, nbrEig);
451 for (
Index j = 0; j < nbrEig; j++)
453 Sr.col(j) = schurofH.matrixU().col(perm(it-j-1));
458 X = m_V.leftCols(it) * Sr;
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);
468 Index m = m_V.rows();
469 if (!m_isDeflAllocated)
470 dgmresInitDeflation(m);
471 DenseMatrix MX(m, nbrEig);
473 for (
Index j = 0; j < nbrEig; j++)
475 tv1 = mat * X.col(j);
476 MX.col(j) = precond.solve(tv1);
480 m_T.block(m_r, m_r, nbrEig, nbrEig) = X.transpose() * MX;
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);
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);
494 m_luT.compute(m_T.topLeftCorner(m_r, m_r));
497 m_isDeflInitialized =
true;
500 template<
typename _MatrixType,
typename _Preconditioner>
501 template<
typename RhsType,
typename DestType>
502 Index DGMRES<_MatrixType, _Preconditioner>::dgmresApplyDeflation(
const RhsType &x, DestType &y)
const
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);
A Restarted GMRES with deflation. This class implements a modification of the GMRES solver for sparse...
Definition: DGMRES.h:19
void applyOnTheLeft(const EigenBase< OtherDerived > &other)
void setEigenv(const Index neig)
Definition: DGMRES.h:181
Index deflSize()
Definition: DGMRES.h:190
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 maxIterations() const
Derived & setZero(Index size)
void set_restart(const Index restart)
Definition: DGMRES.h:176
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
Index restart()
Definition: DGMRES.h:171
ScalarBinaryOpTraits< typename internal::traits< Derived >::Scalar, typename internal::traits< OtherDerived >::Scalar >::ReturnType dot(const MatrixBase< OtherDerived > &other) const
DGMRES()
Definition: DGMRES.h:125
void setMaxEigenv(const Index maxNeig)
Definition: DGMRES.h:195
DGMRES(const EigenBase< MatrixDerived > &A)
Definition: DGMRES.h:138
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