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, 2014 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>
|
---|
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 < Scalar, Dynamic, 1 > VectorType;
|
---|
65 | typedef Matrix < Scalar, Dynamic, Dynamic > FMatrixType;
|
---|
66 |
|
---|
67 | RealScalar tol = tol_error;
|
---|
68 | const int maxIters = iters;
|
---|
69 | iters = 0;
|
---|
70 |
|
---|
71 | const int m = mat.rows();
|
---|
72 |
|
---|
73 | VectorType p0 = rhs - mat*x;
|
---|
74 | VectorType r0 = precond.solve(p0);
|
---|
75 |
|
---|
76 | // is initial guess already good enough?
|
---|
77 | if(abs(r0.norm()) < tol) {
|
---|
78 | return true;
|
---|
79 | }
|
---|
80 |
|
---|
81 | VectorType w = VectorType::Zero(restart + 1);
|
---|
82 |
|
---|
83 | FMatrixType H = FMatrixType::Zero(m, restart + 1); // Hessenberg matrix
|
---|
84 | VectorType tau = VectorType::Zero(restart + 1);
|
---|
85 | std::vector < JacobiRotation < Scalar > > G(restart);
|
---|
86 |
|
---|
87 | // generate first Householder vector
|
---|
88 | VectorType e(m-1);
|
---|
89 | RealScalar beta;
|
---|
90 | r0.makeHouseholder(e, tau.coeffRef(0), beta);
|
---|
91 | w(0)=(Scalar) beta;
|
---|
92 | H.bottomLeftCorner(m - 1, 1) = e;
|
---|
93 |
|
---|
94 | for (int k = 1; k <= restart; ++k) {
|
---|
95 |
|
---|
96 | ++iters;
|
---|
97 |
|
---|
98 | VectorType v = VectorType::Unit(m, k - 1), workspace(m);
|
---|
99 |
|
---|
100 | // apply Householder reflections H_{1} ... H_{k-1} to v
|
---|
101 | for (int i = k - 1; i >= 0; --i) {
|
---|
102 | v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
|
---|
103 | }
|
---|
104 |
|
---|
105 | // apply matrix M to v: v = mat * v;
|
---|
106 | VectorType t=mat*v;
|
---|
107 | v=precond.solve(t);
|
---|
108 |
|
---|
109 | // apply Householder reflections H_{k-1} ... H_{1} to v
|
---|
110 | for (int i = 0; i < k; ++i) {
|
---|
111 | v.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
|
---|
112 | }
|
---|
113 |
|
---|
114 | if (v.tail(m - k).norm() != 0.0) {
|
---|
115 |
|
---|
116 | if (k <= restart) {
|
---|
117 |
|
---|
118 | // generate new Householder vector
|
---|
119 | VectorType e(m - k - 1);
|
---|
120 | RealScalar beta;
|
---|
121 | v.tail(m - k).makeHouseholder(e, tau.coeffRef(k), beta);
|
---|
122 | H.col(k).tail(m - k - 1) = e;
|
---|
123 |
|
---|
124 | // apply Householder reflection H_{k} to v
|
---|
125 | v.tail(m - k).applyHouseholderOnTheLeft(H.col(k).tail(m - k - 1), tau.coeffRef(k), workspace.data());
|
---|
126 |
|
---|
127 | }
|
---|
128 | }
|
---|
129 |
|
---|
130 | if (k > 1) {
|
---|
131 | for (int i = 0; i < k - 1; ++i) {
|
---|
132 | // apply old Givens rotations to v
|
---|
133 | v.applyOnTheLeft(i, i + 1, G[i].adjoint());
|
---|
134 | }
|
---|
135 | }
|
---|
136 |
|
---|
137 | if (k<m && v(k) != (Scalar) 0) {
|
---|
138 | // determine next Givens rotation
|
---|
139 | G[k - 1].makeGivens(v(k - 1), v(k));
|
---|
140 |
|
---|
141 | // apply Givens rotation to v and w
|
---|
142 | v.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
|
---|
143 | w.applyOnTheLeft(k - 1, k, G[k - 1].adjoint());
|
---|
144 |
|
---|
145 | }
|
---|
146 |
|
---|
147 | // insert coefficients into upper matrix triangle
|
---|
148 | H.col(k - 1).head(k) = v.head(k);
|
---|
149 |
|
---|
150 | bool stop=(k==m || abs(w(k)) < tol || iters == maxIters);
|
---|
151 |
|
---|
152 | if (stop || k == restart) {
|
---|
153 |
|
---|
154 | // solve upper triangular system
|
---|
155 | VectorType y = w.head(k);
|
---|
156 | H.topLeftCorner(k, k).template triangularView < Eigen::Upper > ().solveInPlace(y);
|
---|
157 |
|
---|
158 | // use Horner-like scheme to calculate solution vector
|
---|
159 | VectorType x_new = y(k - 1) * VectorType::Unit(m, k - 1);
|
---|
160 |
|
---|
161 | // apply Householder reflection H_{k} to x_new
|
---|
162 | x_new.tail(m - k + 1).applyHouseholderOnTheLeft(H.col(k - 1).tail(m - k), tau.coeffRef(k - 1), workspace.data());
|
---|
163 |
|
---|
164 | for (int i = k - 2; i >= 0; --i) {
|
---|
165 | x_new += y(i) * VectorType::Unit(m, i);
|
---|
166 | // apply Householder reflection H_{i} to x_new
|
---|
167 | x_new.tail(m - i).applyHouseholderOnTheLeft(H.col(i).tail(m - i - 1), tau.coeffRef(i), workspace.data());
|
---|
168 | }
|
---|
169 |
|
---|
170 | x += x_new;
|
---|
171 |
|
---|
172 | if (stop) {
|
---|
173 | return true;
|
---|
174 | } else {
|
---|
175 | k=0;
|
---|
176 |
|
---|
177 | // reset data for a restart r0 = rhs - mat * x;
|
---|
178 | VectorType p0=mat*x;
|
---|
179 | VectorType p1=precond.solve(p0);
|
---|
180 | r0 = rhs - p1;
|
---|
181 | // r0_sqnorm = r0.squaredNorm();
|
---|
182 | w = VectorType::Zero(restart + 1);
|
---|
183 | H = FMatrixType::Zero(m, restart + 1);
|
---|
184 | tau = VectorType::Zero(restart + 1);
|
---|
185 |
|
---|
186 | // generate first Householder vector
|
---|
187 | RealScalar beta;
|
---|
188 | r0.makeHouseholder(e, tau.coeffRef(0), beta);
|
---|
189 | w(0)=(Scalar) beta;
|
---|
190 | H.bottomLeftCorner(m - 1, 1) = e;
|
---|
191 |
|
---|
192 | }
|
---|
193 |
|
---|
194 | }
|
---|
195 |
|
---|
196 |
|
---|
197 |
|
---|
198 | }
|
---|
199 |
|
---|
200 | return false;
|
---|
201 |
|
---|
202 | }
|
---|
203 |
|
---|
204 | }
|
---|
205 |
|
---|
206 | template< typename _MatrixType,
|
---|
207 | typename _Preconditioner = DiagonalPreconditioner<typename _MatrixType::Scalar> >
|
---|
208 | class GMRES;
|
---|
209 |
|
---|
210 | namespace internal {
|
---|
211 |
|
---|
212 | template< typename _MatrixType, typename _Preconditioner>
|
---|
213 | struct traits<GMRES<_MatrixType,_Preconditioner> >
|
---|
214 | {
|
---|
215 | typedef _MatrixType MatrixType;
|
---|
216 | typedef _Preconditioner Preconditioner;
|
---|
217 | };
|
---|
218 |
|
---|
219 | }
|
---|
220 |
|
---|
221 | /** \ingroup IterativeLinearSolvers_Module
|
---|
222 | * \brief A GMRES solver for sparse square problems
|
---|
223 | *
|
---|
224 | * This class allows to solve for A.x = b sparse linear problems using a generalized minimal
|
---|
225 | * residual method. The vectors x and b can be either dense or sparse.
|
---|
226 | *
|
---|
227 | * \tparam _MatrixType the type of the sparse matrix A, can be a dense or a sparse matrix.
|
---|
228 | * \tparam _Preconditioner the type of the preconditioner. Default is DiagonalPreconditioner
|
---|
229 | *
|
---|
230 | * The maximal number of iterations and tolerance value can be controlled via the setMaxIterations()
|
---|
231 | * and setTolerance() methods. The defaults are the size of the problem for the maximal number of iterations
|
---|
232 | * and NumTraits<Scalar>::epsilon() for the tolerance.
|
---|
233 | *
|
---|
234 | * This class can be used as the direct solver classes. Here is a typical usage example:
|
---|
235 | * \code
|
---|
236 | * int n = 10000;
|
---|
237 | * VectorXd x(n), b(n);
|
---|
238 | * SparseMatrix<double> A(n,n);
|
---|
239 | * // fill A and b
|
---|
240 | * GMRES<SparseMatrix<double> > solver(A);
|
---|
241 | * x = solver.solve(b);
|
---|
242 | * std::cout << "#iterations: " << solver.iterations() << std::endl;
|
---|
243 | * std::cout << "estimated error: " << solver.error() << std::endl;
|
---|
244 | * // update b, and solve again
|
---|
245 | * x = solver.solve(b);
|
---|
246 | * \endcode
|
---|
247 | *
|
---|
248 | * By default the iterations start with x=0 as an initial guess of the solution.
|
---|
249 | * One can control the start using the solveWithGuess() method.
|
---|
250 | *
|
---|
251 | * \sa class SimplicialCholesky, DiagonalPreconditioner, IdentityPreconditioner
|
---|
252 | */
|
---|
253 | template< typename _MatrixType, typename _Preconditioner>
|
---|
254 | class GMRES : public IterativeSolverBase<GMRES<_MatrixType,_Preconditioner> >
|
---|
255 | {
|
---|
256 | typedef IterativeSolverBase<GMRES> Base;
|
---|
257 | using Base::mp_matrix;
|
---|
258 | using Base::m_error;
|
---|
259 | using Base::m_iterations;
|
---|
260 | using Base::m_info;
|
---|
261 | using Base::m_isInitialized;
|
---|
262 |
|
---|
263 | private:
|
---|
264 | int m_restart;
|
---|
265 |
|
---|
266 | public:
|
---|
267 | typedef _MatrixType MatrixType;
|
---|
268 | typedef typename MatrixType::Scalar Scalar;
|
---|
269 | typedef typename MatrixType::Index Index;
|
---|
270 | typedef typename MatrixType::RealScalar RealScalar;
|
---|
271 | typedef _Preconditioner Preconditioner;
|
---|
272 |
|
---|
273 | public:
|
---|
274 |
|
---|
275 | /** Default constructor. */
|
---|
276 | GMRES() : Base(), m_restart(30) {}
|
---|
277 |
|
---|
278 | /** Initialize the solver with matrix \a A for further \c Ax=b solving.
|
---|
279 | *
|
---|
280 | * This constructor is a shortcut for the default constructor followed
|
---|
281 | * by a call to compute().
|
---|
282 | *
|
---|
283 | * \warning this class stores a reference to the matrix A as well as some
|
---|
284 | * precomputed values that depend on it. Therefore, if \a A is changed
|
---|
285 | * this class becomes invalid. Call compute() to update it with the new
|
---|
286 | * matrix A, or modify a copy of A.
|
---|
287 | */
|
---|
288 | template<typename MatrixDerived>
|
---|
289 | explicit GMRES(const EigenBase<MatrixDerived>& A) : Base(A.derived()), m_restart(30) {}
|
---|
290 |
|
---|
291 | ~GMRES() {}
|
---|
292 |
|
---|
293 | /** Get the number of iterations after that a restart is performed.
|
---|
294 | */
|
---|
295 | int get_restart() { return m_restart; }
|
---|
296 |
|
---|
297 | /** Set the number of iterations after that a restart is performed.
|
---|
298 | * \param restart number of iterations for a restarti, default is 30.
|
---|
299 | */
|
---|
300 | void set_restart(const int restart) { m_restart=restart; }
|
---|
301 |
|
---|
302 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A
|
---|
303 | * \a x0 as an initial solution.
|
---|
304 | *
|
---|
305 | * \sa compute()
|
---|
306 | */
|
---|
307 | template<typename Rhs,typename Guess>
|
---|
308 | inline const internal::solve_retval_with_guess<GMRES, Rhs, Guess>
|
---|
309 | solveWithGuess(const MatrixBase<Rhs>& b, const Guess& x0) const
|
---|
310 | {
|
---|
311 | eigen_assert(m_isInitialized && "GMRES is not initialized.");
|
---|
312 | eigen_assert(Base::rows()==b.rows()
|
---|
313 | && "GMRES::solve(): invalid number of rows of the right hand side matrix b");
|
---|
314 | return internal::solve_retval_with_guess
|
---|
315 | <GMRES, Rhs, Guess>(*this, b.derived(), x0);
|
---|
316 | }
|
---|
317 |
|
---|
318 | /** \internal */
|
---|
319 | template<typename Rhs,typename Dest>
|
---|
320 | void _solveWithGuess(const Rhs& b, Dest& x) const
|
---|
321 | {
|
---|
322 | bool failed = false;
|
---|
323 | for(int j=0; j<b.cols(); ++j)
|
---|
324 | {
|
---|
325 | m_iterations = Base::maxIterations();
|
---|
326 | m_error = Base::m_tolerance;
|
---|
327 |
|
---|
328 | typename Dest::ColXpr xj(x,j);
|
---|
329 | if(!internal::gmres(*mp_matrix, b.col(j), xj, Base::m_preconditioner, m_iterations, m_restart, m_error))
|
---|
330 | failed = true;
|
---|
331 | }
|
---|
332 | m_info = failed ? NumericalIssue
|
---|
333 | : m_error <= Base::m_tolerance ? Success
|
---|
334 | : NoConvergence;
|
---|
335 | m_isInitialized = true;
|
---|
336 | }
|
---|
337 |
|
---|
338 | /** \internal */
|
---|
339 | template<typename Rhs,typename Dest>
|
---|
340 | void _solve(const Rhs& b, Dest& x) const
|
---|
341 | {
|
---|
342 | x = b;
|
---|
343 | if(x.squaredNorm() == 0) return; // Check Zero right hand side
|
---|
344 | _solveWithGuess(b,x);
|
---|
345 | }
|
---|
346 |
|
---|
347 | protected:
|
---|
348 |
|
---|
349 | };
|
---|
350 |
|
---|
351 |
|
---|
352 | namespace internal {
|
---|
353 |
|
---|
354 | template<typename _MatrixType, typename _Preconditioner, typename Rhs>
|
---|
355 | struct solve_retval<GMRES<_MatrixType, _Preconditioner>, Rhs>
|
---|
356 | : solve_retval_base<GMRES<_MatrixType, _Preconditioner>, Rhs>
|
---|
357 | {
|
---|
358 | typedef GMRES<_MatrixType, _Preconditioner> Dec;
|
---|
359 | EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
|
---|
360 |
|
---|
361 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
362 | {
|
---|
363 | dec()._solve(rhs(),dst);
|
---|
364 | }
|
---|
365 | };
|
---|
366 |
|
---|
367 | } // end namespace internal
|
---|
368 |
|
---|
369 | } // end namespace Eigen
|
---|
370 |
|
---|
371 | #endif // EIGEN_GMRES_H
|
---|