Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GMRES__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GMRES__
35 :
36 : #include <iostream>
37 : #include <string>
38 :
39 : #include "lib_algebra/operator/interface/operator.h"
40 : #include "common/profiler/profiler.h"
41 : #include "lib_algebra/operator/interface/pprocess.h"
42 : #ifdef UG_PARALLEL
43 : #include "lib_algebra/parallelization/parallelization.h"
44 : #endif
45 :
46 : namespace ug{
47 :
48 : /// the GMREs method as a solver for linear operators
49 : /**
50 : * This class implements the GMRES - method for the solution of linear
51 : * operator problems like A*x = b, where the solution x = A^{-1} b is computed.
52 : *
53 : * For detailed description of the algorithm, please refer to:
54 : *
55 : * - Barrett, Berry, Chan, Demmel, Donatom Dongarra, Eijkhout, Pozo, Romine,
56 : * Van der Vorst, "Templates for the Solution of Linear Systems: Building
57 : * Blocks for Iterative Methods"
58 : *
59 : * - Saad, "Iterative Methods For Sparse Linear Systems"
60 : *
61 : * \tparam TVector vector type
62 : */
63 : template <typename TVector>
64 : class GMRES
65 : : public IPreconditionedLinearOperatorInverse<TVector>
66 : {
67 : public:
68 : /// Vector type
69 : typedef TVector vector_type;
70 :
71 : /// Base type
72 : typedef IPreconditionedLinearOperatorInverse<vector_type> base_type;
73 :
74 : protected:
75 : using base_type::convergence_check;
76 : using base_type::linear_operator;
77 : using base_type::preconditioner;
78 : using base_type::write_debug;
79 :
80 : public:
81 : /// default constructor
82 0 : GMRES(size_t restart) : m_restart(restart) {};
83 :
84 : /// constructor setting the preconditioner and the convergence check
85 : GMRES( size_t restart,
86 : SmartPtr<ILinearIterator<vector_type> > spPrecond,
87 : SmartPtr<IConvergenceCheck<vector_type> > spConvCheck)
88 : : base_type(spPrecond, spConvCheck), m_restart(restart)
89 : {};
90 :
91 : /// name of solver
92 0 : virtual const char* name() const {return "GMRES";}
93 :
94 : /// returns if parallel solving is supported
95 0 : virtual bool supports_parallel() const
96 : {
97 0 : if(preconditioner().valid())
98 0 : return preconditioner()->supports_parallel();
99 : return true;
100 : }
101 :
102 : // Solve J(u)*x = b, such that x = J(u)^{-1} b
103 0 : virtual bool apply_return_defect(vector_type& x, vector_type& b)
104 : {
105 : // check correct storage type in parallel
106 : #ifdef UG_PARALLEL
107 : if(!b.has_storage_type(PST_ADDITIVE) || !x.has_storage_type(PST_CONSISTENT))
108 : UG_THROW("GMRES: Inadequate storage format of Vectors.");
109 : #endif
110 :
111 : // copy rhs
112 0 : SmartPtr<vector_type> spR = b.clone();
113 :
114 : // build defect: b := b - A*x
115 0 : linear_operator()->apply_sub(*spR, x);
116 :
117 : // prepare convergence check
118 0 : prepare_conv_check();
119 :
120 : // compute start defect norm
121 0 : convergence_check()->start(*spR);
122 :
123 : // storage for v, h, gamma
124 0 : std::vector<SmartPtr<vector_type> > v(m_restart+1);
125 0 : std::vector<std::vector<number> > h(m_restart+1);
126 0 : for(size_t i = 0; i < h.size(); ++i) h[i].resize(m_restart+1);
127 0 : std::vector<number> gamma(m_restart+1);
128 0 : std::vector<number> c(m_restart+1);
129 0 : std::vector<number> s(m_restart+1);
130 :
131 : // old norm
132 : number oldNorm;
133 :
134 : // Iteration loop
135 0 : while(!convergence_check()->iteration_ended())
136 : {
137 : // get storage for first vector v[0]
138 0 : if(v[0].invalid()) v[0] = x.clone_without_values();
139 :
140 : // apply v[0] = M^-1 * (b-A*x)
141 0 : if(preconditioner().valid()){
142 0 : if(!preconditioner()->apply(*v[0], *spR)){
143 : UG_LOG("GMRES: Cannot apply preconditioner to b-A*x0.\n");
144 : return false;
145 : }
146 : }
147 : // ... or reuse v[0] = (b-A*x)
148 : else{
149 0 : SmartPtr<vector_type> tmp = v[0]; v[0] = spR; spR = tmp;
150 : }
151 :
152 : // make v[0] unique
153 : #ifdef UG_PARALLEL
154 : if(!v[0]->change_storage_type(PST_UNIQUE))
155 : UG_THROW("GMRES: Cannot convert v0 to consistent vector.");
156 : #endif
157 :
158 : // post-process the correction
159 0 : m_corr_post_process.apply (*v[0]);
160 :
161 : // Compute norm of inital residuum:
162 0 : oldNorm = gamma[0] = v[0]->norm();
163 :
164 : // normalize v[0] := v[0] / ||v[0]||
165 0 : *v[0] *= 1./gamma[0];
166 :
167 : // loop gmres iterations
168 : size_t numIter = 0;
169 0 : for(size_t j = 0; j < m_restart; ++j)
170 : {
171 : numIter = j;
172 :
173 : // get storage for v[j+1]
174 0 : if(v[j+1].invalid()) v[j+1] = x.clone_without_values();
175 :
176 : #ifdef UG_PARALLEL
177 : if(!v[j]->change_storage_type(PST_CONSISTENT))
178 : UG_THROW("GMRES: Cannot convert v["<<j+1<<"] to consistent vector.");
179 : #endif
180 :
181 : // compute r = A*v[j]
182 0 : linear_operator()->apply(*spR, *v[j]);
183 :
184 : // apply v[j+1] = M^-1 * A * v[j]
185 0 : if(preconditioner().valid()){
186 0 : if(!preconditioner()->apply(*v[j+1], *spR)){
187 : UG_LOG("GMRES: Cannot apply preconditioner to A*v["<<j<<"].\n");
188 : return false;
189 : }
190 : }
191 : // ... or reuse v[j+1] = A * v[j]
192 : else{
193 0 : SmartPtr<vector_type> tmp = v[j+1]; v[j+1] = spR; spR = tmp;
194 : }
195 :
196 : // make v[j], v[j+1] unique
197 : #ifdef UG_PARALLEL
198 : if(!v[j]->change_storage_type(PST_UNIQUE))
199 : UG_THROW("GMRES: Cannot convert v0 to consistent vector.");
200 : if(!v[j+1]->change_storage_type(PST_UNIQUE))
201 : UG_THROW("GMRES: Cannot convert v["<<j<<"] to consistent vector.");
202 : #endif
203 :
204 : // post-process the correction
205 0 : m_corr_post_process.apply (*v[j+1]);
206 :
207 : // loop previous steps
208 0 : for(size_t i = 0; i <= j; ++i)
209 : {
210 : // h_ij := (r, v[j])
211 0 : h[i][j] = VecProd(*v[j+1], *v[i]);
212 :
213 : // v[j+1] -= h_ij * v[i]
214 0 : VecScaleAppend(*v[j+1], *v[i], (-1)*h[i][j]);
215 : }
216 :
217 : // compute h_{j+1,j}
218 0 : h[j+1][j] = v[j+1]->norm();
219 :
220 : // update h
221 0 : for(size_t i = 0; i < j; ++i)
222 : {
223 0 : const number hij = h[i][j];
224 0 : const number hi1j = h[i+1][j];
225 :
226 0 : h[i][j] = c[i+1]*hij + s[i+1]*hi1j;
227 0 : h[i+1][j] = s[i+1]*hij - c[i+1]*hi1j;
228 : }
229 :
230 : // alpha := sqrt(h_jj ^2 + h_{j+1,j}^2)
231 0 : const number alpha = sqrt(h[j][j]*h[j][j] + h[j+1][j]*h[j+1][j]);
232 :
233 : // update s, c
234 0 : s[j+1] = h[j+1][j] / alpha;
235 0 : c[j+1] = h[j][j] / alpha;
236 0 : h[j][j] = alpha;
237 :
238 : // compute new norm
239 0 : gamma[j+1] = s[j+1]*gamma[j];
240 0 : gamma[j] = c[j+1]*gamma[j];
241 :
242 0 : if(preconditioner().valid()) {
243 0 : UG_LOG(std::string(convergence_check()->get_offset(),' '));
244 0 : UG_LOG("% GMRES "<<std::setw(4) <<j+1<<": "
245 : << gamma[j+1] << " " << gamma[j+1] / oldNorm);
246 : UG_LOG(" (in Precond-Norm) \n");
247 0 : oldNorm = gamma[j+1];
248 : }
249 : else{
250 0 : convergence_check()->update_defect(gamma[j+1]);
251 : }
252 :
253 : // normalize v[j+1]
254 0 : *v[j+1] *= 1./(h[j+1][j]);
255 : }
256 :
257 : // compute current x
258 0 : for(size_t i = numIter; ; --i){
259 0 : for(size_t j = i+1; j <= numIter; ++j)
260 0 : gamma[i] -= h[i][j] * gamma[j];
261 :
262 0 : gamma[i] /= h[i][i];
263 :
264 : // x = x + gamma[i] * v[i]
265 : VecScaleAppend(x, *v[i], gamma[i]);
266 :
267 0 : if(i == 0) break;
268 : }
269 :
270 : // compute fresh defect: b := b - A*x
271 0 : *spR = b;
272 0 : linear_operator()->apply_sub(*spR, x);
273 :
274 0 : if(preconditioner().valid())
275 0 : convergence_check()->update(*spR);
276 : }
277 :
278 : // print ending output
279 0 : return convergence_check()->post();
280 0 : }
281 :
282 : public:
283 0 : virtual std::string config_string() const
284 : {
285 0 : std::stringstream ss;
286 0 : ss << "GMRes ( restart = " << m_restart << ")\n";
287 0 : ss << base_type::config_string_preconditioner_convergence_check();
288 0 : return ss.str();
289 0 : }
290 :
291 : /// adds a post-process for the iterates
292 0 : void add_postprocess_corr (SmartPtr<IPProcessVector<vector_type> > p)
293 : {
294 0 : m_corr_post_process.add (p);
295 0 : }
296 :
297 : /// removes a post-process for the iterates
298 0 : void remove_postprocess_corr (SmartPtr<IPProcessVector<vector_type> > p)
299 : {
300 0 : m_corr_post_process.remove (p);
301 0 : }
302 :
303 : protected:
304 : /// prepares the output of the convergence check
305 0 : void prepare_conv_check()
306 : {
307 : // set iteration symbol and name
308 0 : convergence_check()->set_name(name());
309 0 : convergence_check()->set_symbol('%');
310 :
311 : // set preconditioner string
312 : std::string s;
313 0 : if(preconditioner().valid())
314 0 : s = std::string(" (Precond: ") + preconditioner()->name() + ")";
315 : else
316 : s = " (No Preconditioner) ";
317 0 : convergence_check()->set_info(s);
318 0 : }
319 :
320 : protected:
321 : /// restart parameter
322 : size_t m_restart;
323 :
324 : /// postprocessor for the correction in the iterations
325 : /**
326 : * These postprocess operations are applied to the preconditioned
327 : * defect before the orthogonalization. The goal is to prevent the
328 : * useless kernel parts to prevail in the (floating point) arithmetics.
329 : */
330 : PProcessChain<vector_type> m_corr_post_process;
331 :
332 : /// adds a scaled vector to a second one
333 : bool VecScaleAppend(vector_type& a, vector_type& b, number s)
334 : {
335 : #ifdef UG_PARALLEL
336 : if(a.has_storage_type(PST_UNIQUE) && b.has_storage_type(PST_UNIQUE));
337 : else if(a.has_storage_type(PST_CONSISTENT) && b.has_storage_type(PST_CONSISTENT));
338 : else if (a.has_storage_type(PST_ADDITIVE) && b.has_storage_type(PST_ADDITIVE));
339 : else
340 : {
341 : a.change_storage_type(PST_ADDITIVE);
342 : b.change_storage_type(PST_ADDITIVE);
343 : }
344 : #endif
345 :
346 0 : for(size_t i = 0; i < a.size(); ++i)
347 : {
348 : // todo: move VecScaleAppend to ParallelVector
349 : VecScaleAdd(a[i], 1.0, a[i], s, b[i]);
350 : }
351 : return true;
352 : }
353 :
354 : /// computes the vector product
355 : number VecProd(vector_type& a, vector_type& b)
356 : {
357 : return a.dotprod(b);
358 : }
359 : };
360 :
361 : } // end namespace ug
362 :
363 : #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__GMRES__ */
|
