1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2012 Kolja Brix <brix@igpm.rwth-aaachen.de>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 #ifndef EIGEN_GMRES_H
12 #define EIGEN_GMRES_H
13
14 namespace Eigen {
15
16 namespace internal {
17
18 /**
19 * Generalized Minimal Residual Algorithm based on the
20 * Arnoldi algorithm implemented with Householder reflections.
21 *
22 * Parameters:
23 * \param mat matrix of linear system of equations
24 * \param Rhs right hand side vector of linear system of equations
25 * \param x on input: initial guess, on output: solution
26 * \param precond preconditioner used
27 * \param iters on input: maximum number of iterations to perform
28 * on output: number of iterations performed
29 * \param restart number of iterations for a restart
30 * \param tol_error on input: residual tolerance
31 * on output: residuum achieved
32 *
33 * \sa IterativeMethods::bicgstab()
34 *
35 *
36 * For references, please see:
37 *
38 * Saad, Y. and Schultz, M. H.
39 * GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems.
40 * SIAM J.Sci.Stat.Comp. 7, 1986, pp. 856 - 869.
41 *
42 * Saad, Y.
43 * Iterative Methods for Sparse Linear Systems.
44 * Society for Industrial and Applied Mathematics, Philadelphia, 2003.
45 *
46 * Walker, H. F.
47 * Implementations of the GMRES method.
48 * Comput.Phys.Comm. 53, 1989, pp. 311 - 320.
49 *
50 * Walker, H. F.
51 * Implementation of the GMRES Method using Householder Transformations.
52 * SIAM J.Sci.Stat.Comp. 9, 1988, pp. 152 - 163.
53 *
54 */
55 template<typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
gmres(const MatrixType & mat,const Rhs & rhs,Dest & x,const Preconditioner & precond,int & iters,const int & restart,typename Dest::RealScalar & tol_error)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) {
58
59 using std::sqrt;
60 using std::abs;
61
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;
67
68 RealScalar tol = tol_error;
69 const int maxIters = iters;
70 iters = 0;
71
72 const int m = mat.rows();
73
74 VectorType p0 = rhs - mat*x;
75 VectorType r0 = precond.solve(p0);
76 // RealScalar r0_sqnorm = r0.squaredNorm();
77
78 VectorType w = VectorType::Zero(restart + 1);
79
80 FMatrixType H = FMatrixType::Zero(m, restart + 1);
81 VectorType tau = VectorType::Zero(restart + 1);
82 std::vector < JacobiRotation < Scalar > > G(restart);
83
84 // generate first Householder vector
85 VectorType e;
86 RealScalar beta;
87 r0.makeHouseholder(e, tau.coeffRef(0), beta);
88 w(0)=(Scalar) beta;
89 H.bottomLeftCorner(m - 1, 1) = e;
90
91 for (int k = 1; k <= restart; ++k) {
92
93 ++iters;
94
95 VectorType v = VectorType::Unit(m, k - 1), workspace(m);
96
97 // apply Householder reflections H_{1} ... H_{k-1} to v
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());
100 }
101
102 // apply matrix M to v: v = mat * v;
103 VectorType t=mat*v;
104 v=precond.solve(t);
105
106 // apply Householder reflections H_{k-1} ... H_{1} to v
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());
109 }
110
111 if (v.tail(m - k).norm() != 0.0) {
112
113 if (k <= restart) {
114
115 // generate new Householder vector
116 VectorType e(m - k - 1);
117 RealScalar beta;
118 v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta);
119 H.col(k).tail(m - k - 1) = e;
120
121 // apply Householder reflection H_{k} to v
122 v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data());
123
124 }
125 }
126
127 if (k > 1) {
128 for (int i = 0; i < k - 1; ++i) {
129 // apply old Givens rotations to v
130 v.applyOnTheLeft(i, i + 1, G[i].adjoint());
131 }
132 }
133
134 if (k<m && v(k) != (Scalar) 0) {
135 // determine next Givens rotation
136 G[k - 1].makeGivens(v(k - 1), v(k));
137
138 // apply Givens rotation to v and w
139 v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
140 w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
141
142 }
143
144 // insert coefficients into upper matrix triangle
145 H.col(k - 1).head(k) = v.head(k);
146
147 bool stop=(k==m || abs(w(k)) < tol || iters == maxIters);
148
149 if (stop || k == restart) {
150
151 // solve upper triangular system
152 VectorType y = w.head(k);
153 H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y);
154
155 // use Horner-like scheme to calculate solution vector
156 VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
157
158 // apply Householder reflection H_{k} to x_new
159 x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data());
160
161 for (int i = k - 2; i >= 0; --i) {
162 x_new += y(i) * VectorType::Unit(m, i);
163 // apply Householder reflection H_{i} to x_new
164 x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
165 }
166
167 x += x_new;
168
169 if (stop) {
170 return true;
171 } else {
172 k=0;
173
174 // reset data for a restart r0 = rhs - mat * x;
175 VectorType p0=mat*x;
176 VectorType p1=precond.solve(p0);
177 r0 = rhs - p1;
178 // r0_sqnorm = r0.squaredNorm();
179 w = VectorType::Zero(restart + 1);
180 H = FMatrixType::Zero(m, restart + 1);
181 tau = VectorType::Zero(restart + 1);
182
183 // generate first Householder vector
184 RealScalar beta;
185 r0.makeHouseholder(e, tau.coeffRef(0), beta);
186 w(0)=(Scalar) beta;
187 H.bottomLeftCorner(m - 1, 1) = e;
188
189 }
190
191 }
192
193
194
195 }
196
197 return false;
198
199 }
200
201 }
202
203 template< typename _MatrixType,
204 typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
205 class GMRES;
206
207 namespace internal {
208
209 template< typename _MatrixType, typename _Preconditioner>
210 struct traits<GMRES<_MatrixType,_Preconditioner> >
211 {
212 typedef _MatrixType MatrixType;
213 typedef _Preconditioner Preconditioner;
214 };
215
216 }
217
218 /** \ingroup IterativeLinearSolvers_Module
219 * \brief A GMRES solver for sparse square problems
220 *
221 * This class allows to solve for A.x = b sparse linear problems using a generalized minimal
222 * residual method. The vectors x and b can be either dense or sparse.
223 *
224 * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
225 * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
226 *
227 * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
228 * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
229 * and NumTraits<Scalar>::epsilon() for the tolerance.
230 *
231 * This class can be used as the direct solver classes. Here is a typical usage example:
232 * \code
233 * int n = 10000;
234 * VectorXd x(n), b(n);
235 * SparseMatrix<double> A(n,n);
236 * // fill A and b
237 * GMRES<SparseMatrix<double> > solver(A);
238 * x = solver.solve(b);
239 * std::cout << "#iterations: " << solver.iterations() << std::endl;
240 * std::cout << "estimated error: " << solver.error() << std::endl;
241 * // update b, and solve again
242 * x = solver.solve(b);
243 * \endcode
244 *
245 * By default the iterations start with x=0 as an initial guess of the solution.
246 * One can control the start using the solveWithGuess() method. Here is a step by
247 * step execution example starting with a random guess and printing the evolution
248 * of the estimated error:
249 * * \code
250 * x = VectorXd::Random(n);
251 * solver.setMaxIterations(1);
252 * int i = 0;
253 * do {
254 * x = solver.solveWithGuess(b,x);
255 * std::cout << i << " : " << solver.error() << std::endl;
256 * ++i;
257 * } while (solver.info()!=Success && i<100);
258 * \endcode
259 * Note that such a step by step excution is slightly slower.
260 *
261 * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
262 */
263 template< typename _MatrixType, typename _Preconditioner>
264 class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> >
265 {
266 typedef IterativeSolverBase<GMRES> Base;
267 using Base::mp_matrix;
268 using Base::m_error;
269 using Base::m_iterations;
270 using Base::m_info;
271 using Base::m_isInitialized;
272
273 private:
274 int m_restart;
275
276 public:
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;
282
283 public:
284
285 /** Default constructor. */
286 GMRES() : Base(), m_restart(30) {}
287
288 /** Initialize the solver with matrix \a A for further \c Ax=b solving.
289 *
290 * This constructor is a shortcut for the default constructor followed
291 * by a call to compute().
292 *
293 * \warning this class stores a reference to the matrix A as well as some
294 * precomputed values that depend on it. Therefore, if \a A is changed
295 * this class becomes invalid. Call compute() to update it with the new
296 * matrix A, or modify a copy of A.
297 */
298 GMRES(const MatrixType& A) : Base(A), m_restart(30) {}
299
300 ~GMRES() {}
301
302 /** Get the number of iterations after that a restart is performed.
303 */
304 int get_restart() { return m_restart; }
305
306 /** Set the number of iterations after that a restart is performed.
307 * \param restart number of iterations for a restarti, default is 30.
308 */
309 void set_restart(const int restart) { m_restart=restart; }
310
311 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
312 * \a x0 as an initial solution.
313 *
314 * \sa compute()
315 */
316 template<typename Rhs,typename Guess>
317 inline const internal::solve_retval_with_guess<GMRES, Rhs, Guess>
318 solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
319 {
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);
325 }
326
327 /** \internal */
328 template<typename Rhs,typename Dest>
329 void _solveWithGuess(const Rhs& b, Dest& x) const
330 {
331 bool failed = false;
332 for(int j=0; j<b.cols(); ++j)
333 {
334 m_iterations = Base::maxIterations();
335 m_error = Base::m_tolerance;
336
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))
339 failed = true;
340 }
341 m_info = failed ? NumericalIssue
342 : m_error <= Base::m_tolerance ? Success
343 : NoConvergence;
344 m_isInitialized = true;
345 }
346
347 /** \internal */
348 template<typename Rhs,typename Dest>
349 void _solve(const Rhs& b, Dest& x) const
350 {
351 x.setZero();
352 _solveWithGuess(b,x);
353 }
354
355 protected:
356
357 };
358
359
360 namespace internal {
361
362 template<typename _MatrixType, typename _Preconditioner, typename Rhs>
363 struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs>
364 : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs>
365 {
366 typedef GMRES<_MatrixType, _Preconditioner> Dec;
367 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
368
369 template<typename Dest> void evalTo(Dest& dst) const
370 {
371 dec()._solve(rhs(),dst);
372 }
373 };
374
375 } // end namespace internal
376
377 } // end namespace Eigen
378
379 #endif // EIGEN_GMRES_H
380