70 template<
typename MatrixType,
typename Rhs,
typename Dest,
typename Preconditioner>
71 bool gmres(
const MatrixType & mat,
const Rhs & rhs, Dest & x,
const Preconditioner & precond,
72 int &iters,
const int &restart,
typename Dest::RealScalar & tol_error) {
77 typedef typename Dest::RealScalar RealScalar;
78 typedef typename Dest::Scalar Scalar;
79 typedef Matrix < RealScalar, Dynamic, 1 > RealVectorType;
80 typedef Matrix < Scalar, Dynamic, 1 > VectorType;
81 typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType;
83 RealScalar tol = tol_error;
84 const int maxIters = iters;
87 const int m = mat.rows();
89 VectorType p0 = rhs - mat*x;
90 VectorType r0 = precond.solve(p0);
93 VectorType w = VectorType::Zero(restart + 1);
95 FMatrixType H = FMatrixType::Zero(m, restart + 1);
96 VectorType tau = VectorType::Zero(restart + 1);
97 std::vector < JacobiRotation < Scalar > > G(restart);
102 r0.makeHouseholder(e, tau.coeffRef(0), beta);
104 H.bottomLeftCorner(m - 1, 1) = e;
106 for (
int k = 1; k <= restart; ++k) {
110 VectorType v = VectorType::Unit(m, k - 1), workspace(m);
113 for (
int i = k - 1; i >= 0; --i) {
114 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
122 for (
int i = 0; i < k; ++i) {
123 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
126 if (v.tail(m - k).norm() != 0.0) {
133 v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta);
134 H.col(k).tail(m - k - 1) = e;
137 v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data());
143 for (
int i = 0; i < k - 1; ++i) {
145 v.applyOnTheLeft(i, i + 1, G[i].adjoint());
149 if (k<m && v(k) != (Scalar) 0) {
151 G[k - 1].makeGivens(v(k - 1), v(k));
154 v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
155 w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
160 H.col(k - 1).head(k) = v.head(k);
162 bool stop=(k==m ||
abs(w(k)) < tol || iters == maxIters);
164 if (stop || k == restart) {
167 VectorType y = w.head(k);
168 H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y);
171 VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
174 x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data());
176 for (
int i = k - 2; i >= 0; --i) {
177 x_new += y(i) * VectorType::Unit(m, i);
179 x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
191 VectorType p1=precond.solve(p0);
194 w = VectorType::Zero(restart + 1);
195 H = FMatrixType::Zero(m, restart + 1);
196 tau = VectorType::Zero(restart + 1);
200 r0.makeHouseholder(e, tau.coeffRef(0), beta);
202 H.bottomLeftCorner(m - 1, 1) = e;
218 template<
typename _MatrixType,
219 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
224 template<
typename _MatrixType,
typename _Preconditioner>
225 struct traits<GMRES<_MatrixType,_Preconditioner> >
227 typedef _MatrixType MatrixType;
228 typedef _Preconditioner Preconditioner;
278 template<
typename _MatrixType,
typename _Preconditioner>
279 class GMRES :
public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> >
281 typedef IterativeSolverBase<GMRES> Base;
282 using Base::mp_matrix;
284 using Base::m_iterations;
286 using Base::m_isInitialized;
292 typedef _MatrixType MatrixType;
293 typedef typename MatrixType::Scalar Scalar;
294 typedef typename MatrixType::Index Index;
295 typedef typename MatrixType::RealScalar RealScalar;
296 typedef _Preconditioner Preconditioner;
313 GMRES(
const MatrixType& A) : Base(A), m_restart(30) {}
331 template<
typename Rhs,
typename Guess>
332 inline const internal::solve_retval_with_guess<GMRES, Rhs, Guess>
335 eigen_assert(m_isInitialized &&
"GMRES is not initialized.");
336 eigen_assert(Base::rows()==b.rows()
337 &&
"GMRES::solve(): invalid number of rows of the right hand side matrix b");
338 return internal::solve_retval_with_guess
339 <
GMRES, Rhs, Guess>(*
this, b.derived(), x0);
343 template<
typename Rhs,
typename Dest>
344 void _solveWithGuess(
const Rhs& b, Dest& x)
const
347 for(
int j=0; j<b.cols(); ++j)
349 m_iterations = Base::maxIterations();
350 m_error = Base::m_tolerance;
352 typename Dest::ColXpr xj(x,j);
353 if(!internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
357 : m_error <= Base::m_tolerance ?
Success
359 m_isInitialized =
true;
363 template<
typename Rhs,
typename Dest>
364 void _solve(
const Rhs& b, Dest& x)
const
367 _solveWithGuess(b,x);
377 template<
typename _MatrixType,
typename _Preconditioner,
typename Rhs>
378 struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs>
379 : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs>
381 typedef GMRES<_MatrixType, _Preconditioner> Dec;
382 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
384 template<typename Dest>
void evalTo(Dest& dst)
const
386 dec()._solve(rhs(),dst);
394 #endif // EIGEN_GMRES_H