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 < RealScalar, Dynamic, 1 > RealVectorType;
65 typedef Matrix < Scalar, Dynamic, 1 > VectorType;
66 typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType;
68 RealScalar tol = tol_error;
69 const int maxIters = iters;
72 const int m = mat.rows();
74 VectorType p0 = rhs - mat*x;
75 VectorType r0 = precond.solve(p0);
78 VectorType w = VectorType::Zero(restart + 1);
80 FMatrixType H = FMatrixType::Zero(m, restart + 1);
81 VectorType tau = VectorType::Zero(restart + 1);
82 std::vector < JacobiRotation < Scalar > > G(restart);
87 r0.makeHouseholder(e, tau.coeffRef(0), beta);
89 H.bottomLeftCorner(m - 1, 1) = e;
91 for (
int k = 1; k <= restart; ++k) {
95 VectorType v = VectorType::Unit(m, k - 1), workspace(m);
98 for (
int i = k - 1; i >= 0; --i) {
99 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
107 for (
int i = 0; i < k; ++i) {
108 v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
111 if (v.tail(m - k).norm() != 0.0) {
116 VectorType e(m - k - 1);
118 v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta);
119 H.col(k).tail(m - k - 1) = e;
122 v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data());
128 for (
int i = 0; i < k - 1; ++i) {
130 v.applyOnTheLeft(i, i + 1, G[i].adjoint());
134 if (k<m && v(k) != (Scalar) 0) {
136 G[k - 1].makeGivens(v(k - 1), v(k));
139 v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
140 w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
145 H.col(k - 1).head(k) = v.head(k);
147 bool stop=(k==m || abs(w(k)) < tol || iters == maxIters);
149 if (stop || k == restart) {
152 VectorType y = w.head(k);
153 H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y);
156 VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
159 x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data());
161 for (
int i = k - 2; i >= 0; --i) {
162 x_new += y(i) * VectorType::Unit(m, i);
164 x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
176 VectorType p1=precond.solve(p0);
179 w = VectorType::Zero(restart + 1);
180 H = FMatrixType::Zero(m, restart + 1);
181 tau = VectorType::Zero(restart + 1);
185 r0.makeHouseholder(e, tau.coeffRef(0), beta);
187 H.bottomLeftCorner(m - 1, 1) = e;
203 template<
typename _MatrixType,
204 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
209 template<
typename _MatrixType,
typename _Preconditioner>
210 struct traits<GMRES<_MatrixType,_Preconditioner> >
212 typedef _MatrixType MatrixType;
213 typedef _Preconditioner Preconditioner;
263 template<
typename _MatrixType,
typename _Preconditioner>
267 using Base::mp_matrix;
269 using Base::m_iterations;
271 using Base::m_isInitialized;
277 typedef _MatrixType MatrixType;
278 typedef typename MatrixType::Scalar Scalar;
279 typedef typename MatrixType::Index Index;
280 typedef typename MatrixType::RealScalar RealScalar;
281 typedef _Preconditioner Preconditioner;
316 template<
typename Rhs,
typename Guess>
317 inline const internal::solve_retval_with_guess<GMRES, Rhs, Guess>
320 eigen_assert(m_isInitialized &&
"GMRES is not initialized.");
321 eigen_assert(Base::rows()==b.rows()
322 &&
"GMRES::solve(): invalid number of rows of the right hand side matrix b");
323 return internal::solve_retval_with_guess
324 <
GMRES, Rhs, Guess>(*
this, b.derived(), x0);
328 template<
typename Rhs,
typename Dest>
329 void _solveWithGuess(
const Rhs& b, Dest& x)
const
332 for(
int j=0; j<b.cols(); ++j)
335 m_error = Base::m_tolerance;
337 typename Dest::ColXpr xj(x,j);
338 if(!internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
342 : m_error <= Base::m_tolerance ?
Success
344 m_isInitialized =
true;
348 template<
typename Rhs,
typename Dest>
349 void _solve(
const Rhs& b, Dest& x)
const
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