55 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
56 bool gmres(
const MatrixType & mat,
const Rhs & rhs, Dest & x,
const Preconditioner & precond,
57 int &iters,
const int &restart,
typename Dest::RealScalar & tol_error) {
62 typedef typename Dest::RealScalar RealScalar;
63 typedef typename Dest::Scalar Scalar;
64 typedef Matrix < Scalar, Dynamic, 1 > VectorType;
65 typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType;
67 RealScalar tol = tol_error;
68 const int maxIters = iters;
71 const int m = mat.rows();
73 VectorType p0 = rhs - mat*x;
74 VectorType r0 = precond.solve(p0);
77 VectorType w = VectorType::Zero(restart + 1);
79 FMatrixType H = FMatrixType::Zero(m, restart + 1);
80 VectorType tau = VectorType::Zero(restart + 1);
81 std::vector < JacobiRotation < Scalar > > G(restart);
86 r0.makeHouseholder(e, tau.coeffRef(0), beta);
88 H.bottomLeftCorner(m - 1, 1) = e;
90 for (
int k = 1; k <= restart; ++k) {
94 VectorType v = VectorType::Unit(m, k - 1), workspace(m);
97 for (
int i = k - 1; i >= 0; --i) {
98 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
106 for (
int i = 0; i < k; ++i) {
107 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
110 if (v.tail(m - k).norm() != 0.0) {
115 VectorType e(m - k - 1);
117 v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta);
118 H.col(k).tail(m - k - 1) = e;
121 v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data());
127 for (
int i = 0; i < k - 1; ++i) {
129 v.applyOnTheLeft(i, i + 1, G[i].adjoint());
133 if (k<m && v(k) != (Scalar) 0) {
135 G[k - 1].makeGivens(v(k - 1), v(k));
138 v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
139 w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
144 H.col(k - 1).head(k) = v.head(k);
146 bool stop=(k==m || abs(w(k)) < tol || iters == maxIters);
148 if (stop || k == restart) {
151 VectorType y = w.head(k);
152 H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y);
155 VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
158 x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data());
160 for (
int i = k - 2; i >= 0; --i) {
161 x_new += y(i) * VectorType::Unit(m, i);
163 x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
175 VectorType p1=precond.solve(p0);
178 w = VectorType::Zero(restart + 1);
179 H = FMatrixType::Zero(m, restart + 1);
180 tau = VectorType::Zero(restart + 1);
184 r0.makeHouseholder(e, tau.coeffRef(0), beta);
186 H.bottomLeftCorner(m - 1, 1) = e;
202 template<
typename _MatrixType,
203 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
208 template<
typename _MatrixType,
typename _Preconditioner>
209 struct traits<
GMRES<_MatrixType,_Preconditioner> >
211 typedef _MatrixType MatrixType;
212 typedef _Preconditioner Preconditioner;
262 template<
typename _MatrixType,
typename _Preconditioner>
263 class GMRES :
public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> >
265 typedef IterativeSolverBase<GMRES> Base;
266 using Base::mp_matrix;
268 using Base::m_iterations;
270 using Base::m_isInitialized;
276 typedef _MatrixType MatrixType;
277 typedef typename MatrixType::Scalar Scalar;
278 typedef typename MatrixType::Index Index;
279 typedef typename MatrixType::RealScalar RealScalar;
280 typedef _Preconditioner Preconditioner;
297 GMRES(
const MatrixType& A) : Base(A), m_restart(30) {}
315 template<
typename Rhs,
typename Guess>
316 inline const internal::solve_retval_with_guess<GMRES, Rhs, Guess>
319 eigen_assert(m_isInitialized &&
"GMRES is not initialized.");
320 eigen_assert(Base::rows()==b.rows()
321 &&
"GMRES::solve(): invalid number of rows of the right hand side matrix b");
322 return internal::solve_retval_with_guess
323 <
GMRES, Rhs, Guess>(*
this, b.derived(), x0);
327 template<
typename Rhs,
typename Dest>
328 void _solveWithGuess(
const Rhs& b, Dest& x)
const
331 for(
int j=0; j<b.cols(); ++j)
333 m_iterations = Base::maxIterations();
334 m_error = Base::m_tolerance;
336 typename Dest::ColXpr xj(x,j);
337 if(!internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
340 m_info = failed ? NumericalIssue
341 : m_error <= Base::m_tolerance ? Success
343 m_isInitialized =
true;
347 template<
typename Rhs,
typename Dest>
348 void _solve(
const Rhs& b, Dest& x)
const
351 if(x.squaredNorm() == 0)
return;
352 _solveWithGuess(b,x);
362 template<
typename _MatrixType,
typename _Preconditioner,
typename Rhs>
363 struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs>
364 : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs>
366 typedef GMRES<_MatrixType, _Preconditioner> Dec;
367 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
369 template<typename Dest>
void evalTo(Dest& dst)
const
371 dec()._solve(rhs(),dst);
379 #endif // EIGEN_GMRES_H