Line data Source code
1 : /*
2 : * Copyright (c) 2010-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__BICGSTAB__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__BICGSTAB__
35 :
36 : #include <iostream>
37 : #include <string>
38 : #include <sstream>
39 :
40 : #include "lib_algebra/operator/interface/operator.h"
41 : #include "lib_algebra/operator/interface/preconditioned_linear_operator_inverse.h"
42 : #include "lib_algebra/operator/interface/linear_solver_profiling.h"
43 : #include "lib_algebra/operator/interface/pprocess.h"
44 : #ifdef UG_PARALLEL
45 : #include "lib_algebra/parallelization/parallelization.h"
46 : #endif
47 : #include "common/util/string_util.h"
48 :
49 : namespace ug{
50 :
51 : /// the BiCGStab method as a solver for linear operators
52 : /**
53 : * This class implements the BiCGStab - method for the solution of linear
54 : * operator problems like A*x = b, where the solution x = A^{-1} b is computed.
55 : *
56 : * For detailed description of the algorithm, please refer to:
57 : *
58 : * - Barrett, Berry, Chan, Demmel, Donatom Dongarra, Eijkhout, Pozo, Romine,
59 : * Van der Vorst, "Templates for the Solution of Linear Systems: Building
60 : * Blocks for Iterative Methods", p.24, Fig, 2.10
61 : *
62 : * - Saad, "Iterative Methods For Sparse Linear Systems", p246, Alg. 7.7
63 : *
64 : * \tparam TVector vector type
65 : */
66 : template <typename TVector>
67 : class BiCGStab
68 : : public IPreconditionedLinearOperatorInverse<TVector>
69 : {
70 : public:
71 : /// Vector type
72 : typedef TVector vector_type;
73 :
74 : /// Base type
75 : typedef IPreconditionedLinearOperatorInverse<vector_type> base_type;
76 :
77 : protected:
78 : using base_type::convergence_check;
79 : using base_type::linear_operator;
80 : using base_type::preconditioner;
81 : using base_type::write_debug;
82 :
83 : public:
84 : /// constructors
85 0 : BiCGStab() :
86 : m_numRestarts(0), m_minOrtho(0.0)
87 : {};
88 :
89 0 : BiCGStab(SmartPtr<ILinearIterator<vector_type,vector_type> > spPrecond)
90 : : base_type ( spPrecond ),
91 0 : m_numRestarts(0), m_minOrtho(0.0)
92 0 : {}
93 :
94 0 : BiCGStab( SmartPtr<ILinearIterator<vector_type> > spPrecond,
95 : SmartPtr<IConvergenceCheck<vector_type> > spConvCheck)
96 : : base_type(spPrecond, spConvCheck),
97 0 : m_numRestarts(0), m_minOrtho(0.0)
98 0 : {};
99 :
100 : /// name of solver
101 0 : virtual const char* name() const {return "BiCGStab";}
102 :
103 : /// returns if parallel solving is supported
104 0 : virtual bool supports_parallel() const
105 : {
106 0 : if(preconditioner().valid())
107 0 : return preconditioner()->supports_parallel();
108 : return true;
109 : }
110 :
111 : // Solve J(u)*x = b, such that x = J(u)^{-1} b
112 0 : virtual bool apply_return_defect(vector_type& x, vector_type& b)
113 : {
114 : LS_PROFILE_BEGIN(LS_ApplyReturnDefect);
115 :
116 : // check correct storage type in parallel
117 : #ifdef UG_PARALLEL
118 : if(!b.has_storage_type(PST_ADDITIVE) || !x.has_storage_type(PST_CONSISTENT)){
119 : UG_LOG("BICGStab:b_storage_type"<<b.has_storage_type(PST_ADDITIVE)<<"\n");
120 : UG_LOG("BICGStab:x_storage_type"<<x.has_storage_type(PST_CONSISTENT)<<"\n");
121 : UG_THROW("BiCGStab: Inadequate storage format of Vectors.");
122 : }
123 : #endif
124 :
125 : // build defect: r := b - A*x
126 0 : linear_operator()->apply_sub(b, x);
127 : vector_type& r = b;
128 :
129 : // create vectors
130 0 : SmartPtr<vector_type> spR = r.clone_without_values(); vector_type& r0 = *spR;
131 0 : SmartPtr<vector_type> spP = r.clone_without_values(); vector_type& p = *spP;
132 0 : SmartPtr<vector_type> spV = r.clone_without_values(); vector_type& v = *spV;
133 0 : SmartPtr<vector_type> spT = r.clone_without_values(); vector_type& t = *spT;
134 0 : SmartPtr<vector_type> spS = r.clone_without_values(); vector_type& s = *spS;
135 0 : SmartPtr<vector_type> spQ = x.clone_without_values(); vector_type& q = *spQ;
136 :
137 : // prepare convergence check
138 : prepare_conv_check();
139 :
140 : // compute start defect norm
141 0 : convergence_check()->start(r);
142 :
143 : // convert b to unique (should already be unique due to norm calculation)
144 : #ifdef UG_PARALLEL
145 : if(!r.change_storage_type(PST_UNIQUE))
146 : UG_THROW("BiCGStab: Cannot convert b to a vector with the 'unique' parallel storage type.");
147 : #endif
148 :
149 : // needed variables
150 : number rho = 1, alpha = 1, omega = 1, norm_r0 = 0.0;
151 :
152 : // restart flag (set to true at first run)
153 : bool bRestart = true;
154 :
155 0 : write_debugXR(x, r, convergence_check()->step(), 'i');
156 :
157 : // Iteration loop
158 0 : while(!convergence_check()->iteration_ended())
159 : {
160 : // check for restart based on fixed step number restart
161 0 : if(m_numRestarts > 0 &&
162 0 : (convergence_check()->step() % m_numRestarts == 0))
163 : {
164 0 : std::stringstream ss; ss <<
165 0 : "Restarting: at every "<<m_numRestarts<<" Iterations";
166 0 : convergence_check()->print_line(ss.str());
167 : bRestart = true;
168 0 : }
169 :
170 : // check if start values have to be set
171 0 : if(bRestart)
172 : {
173 : // reset arbitrary vectors
174 0 : r0 = r;
175 :
176 : // make r additive unique
177 : #ifdef UG_PARALLEL
178 : if(!r0.change_storage_type(PST_UNIQUE))
179 : UG_THROW("BiCGStab: Cannot convert r to unique vector.");
180 : #endif
181 :
182 : // set start vectors and alpha, omega:
183 : // This will lead to p = r, since beta = 0.0
184 : p = 0.0; alpha = 0.0;
185 : v = 0.0; omega = 1.0;
186 :
187 : // set rhoOld to 1 (rhoOld = rho, see below)
188 : rho = 1.0;
189 :
190 : // remember start norm
191 0 : norm_r0 = convergence_check()->defect();
192 :
193 : // remove restart flag
194 : bRestart = false;
195 : }
196 :
197 : // remember current rho
198 : const number rhoOld = rho;
199 :
200 : // Compute rho new
201 0 : if (!r.size())
202 : rho = 1.0;
203 : else
204 : rho = VecProd(r0, r);
205 :
206 : // check for restart compare (r, r0) > m_minOrtho * ||r|| ||r0||
207 0 : const number norm_r = convergence_check()->defect();
208 0 : if(fabs(rho)/(norm_r * norm_r0) <= m_minOrtho)
209 : {
210 0 : std::stringstream ss; ss <<
211 0 : "Restarting: Min Orthogonality "<<m_minOrtho<<" missed: "
212 : <<"(r,r0)="<<fabs(rho)<<", ||r||="<<norm_r<<", ||r0||= "
213 : <<norm_r0;
214 0 : convergence_check()->print_line(ss.str());
215 : bRestart = true;
216 0 : }
217 :
218 : // check that rhoOld valid
219 0 : if(rhoOld == 0.0)
220 : {
221 : UG_LOG("BiCGStab: Method breakdown with rhoOld = "<<rhoOld<<
222 : ". Aborting iteration.\n");
223 : return false;
224 : }
225 :
226 : // Compute new beta
227 0 : const number beta = (rho/rhoOld) * (alpha/omega);
228 :
229 : // update p = r + beta * p - beta * omega * v
230 0 : VecScaleAdd(p, 1.0, r, beta, p, -beta*omega, v);
231 :
232 : // apply q = M^-1 * p ...
233 0 : if(preconditioner().valid())
234 : {
235 0 : enter_precond_debug_section(convergence_check()->step(), 'a');
236 0 : if(!preconditioner()->apply(q, p))
237 : {
238 : UG_LOG("BiCGStab: Cannot apply preconditioner. Aborting.\n");
239 0 : this->leave_vector_debug_writer_section();
240 0 : return false;
241 : }
242 0 : this->leave_vector_debug_writer_section();
243 : }
244 : // ... or copy q = p
245 : else
246 : {
247 0 : q = p;
248 :
249 : // make q consistent
250 : #ifdef UG_PARALLEL
251 : if(!q.change_storage_type(PST_CONSISTENT))
252 : UG_THROW("BiCGStab: Cannot convert q to consistent vector.");
253 : #endif
254 : }
255 :
256 : // post-process the correction
257 0 : m_corr_post_process.apply (q);
258 :
259 : // compute v := A*q
260 0 : linear_operator()->apply(v, q);
261 :
262 : // make v unique
263 : #ifdef UG_PARALLEL
264 : if(!v.change_storage_type(PST_UNIQUE))
265 : UG_THROW("BiCGStab: Cannot convert v to unique vector.");
266 : #endif
267 :
268 : // alpha = (v,r)
269 0 : if (!v.size())
270 : alpha = 1.0;
271 : else
272 : alpha = VecProd(v, r0);
273 :
274 : // check validity of alpha
275 0 : if(alpha == 0.0){
276 : UG_LOG("BiCGStab: Method breakdown: alpha = "<<alpha<<
277 : " is an invalid value. Aborting iteration.\n");
278 : return false;
279 : }
280 :
281 : // alpha = rho/(v,r)
282 0 : alpha = rho/alpha;
283 :
284 : // add: x := x + alpha * q
285 0 : VecScaleAdd(x, 1.0, x, alpha, q);
286 :
287 : // compute s = r - alpha*v
288 0 : VecScaleAdd(s, 1.0, r, -alpha, v);
289 :
290 : // check convergence
291 0 : convergence_check()->update(s);
292 :
293 0 : write_debugXR(x, s, convergence_check()->step(), 'a');
294 :
295 : // if finished: set output to last defect and exist loop
296 0 : if(convergence_check()->iteration_ended())
297 : {
298 0 : r = s; break;
299 : }
300 :
301 : // apply q = M^-1 * s ...
302 0 : if(preconditioner().valid())
303 : {
304 0 : enter_precond_debug_section(convergence_check()->step(), 'b');
305 0 : if(!preconditioner()->apply(q, s))
306 : {
307 : UG_LOG("BiCGStab: Cannot apply preconditioner. Aborting.\n");
308 0 : this->leave_vector_debug_writer_section();
309 0 : return false;
310 : }
311 0 : this->leave_vector_debug_writer_section();
312 : }
313 : // ... or set q:=s
314 : else
315 : {
316 0 : q = s;
317 :
318 : // make q consistent
319 : #ifdef UG_PARALLEL
320 : if(!q.change_storage_type(PST_CONSISTENT))
321 : UG_THROW("BiCGStab: Cannot convert q to consistent vector.");
322 : #endif
323 : }
324 :
325 : // post-process the correction
326 0 : m_corr_post_process.apply (q);
327 :
328 : // compute t := A*q
329 0 : linear_operator()->apply(t, q);
330 :
331 : // make t unique
332 : #ifdef UG_PARALLEL
333 : if(!t.change_storage_type(PST_UNIQUE))
334 : UG_THROW("BiCGStab: Cannot convert t to unique vector.");
335 : #endif
336 :
337 : // tt = (t,t)
338 : number tt;
339 0 : if (!t.size())
340 : tt = 1.0;
341 : else
342 : tt = VecProd(t, t);
343 :
344 : // omega = (s,t)
345 0 : if (!s.size())
346 : omega = 1.0;
347 : else
348 : omega = VecProd(s, t);
349 :
350 : // check tt
351 0 : if(tt == 0.0)
352 : {
353 : UG_LOG("BiCGStab: Method breakdown tt = "<<tt<<" is an "
354 : "invalid value. Aborting iteration.\n");
355 : return false;
356 : }
357 :
358 : // omega = (s,t)/(t,t)
359 0 : omega = omega/tt;
360 :
361 : // add: x := x + omega * q
362 0 : VecScaleAdd(x, 1.0, x, omega, q);
363 :
364 : // compute r = s - omega*t
365 0 : VecScaleAdd(r, 1.0, s, -omega, t);
366 :
367 : // check convergence
368 0 : convergence_check()->update(r);
369 :
370 0 : write_debugXR(x, r, convergence_check()->step(), 'b');
371 :
372 : // check values
373 0 : if(omega == 0.0)
374 : {
375 : UG_LOG("BiCGStab: Method breakdown with omega = "<<omega<<
376 : ". Aborting iteration.\n");
377 : return false;
378 : }
379 : }
380 :
381 : // print ending output
382 0 : return convergence_check()->post();
383 : }
384 :
385 :
386 : /// sets to restart at given number of iteration steps
387 0 : void set_restart(int numRestarts) {m_numRestarts = numRestarts;}
388 :
389 : /// sets to restart if given orthogonality missed
390 0 : void set_min_orthogonality(number minOrtho) {m_minOrtho = minOrtho;}
391 :
392 : /// adds a post-process for the iterates
393 0 : void add_postprocess_corr (SmartPtr<IPProcessVector<vector_type> > p)
394 : {
395 0 : m_corr_post_process.add (p);
396 0 : }
397 :
398 : /// removes a post-process for the iterates
399 0 : void remove_postprocess_corr (SmartPtr<IPProcessVector<vector_type> > p)
400 : {
401 0 : m_corr_post_process.remove (p);
402 0 : }
403 :
404 : protected:
405 : /// prepares the output of the convergence check
406 0 : void prepare_conv_check()
407 : {
408 : // set iteration symbol and name
409 0 : convergence_check()->set_name(name());
410 0 : convergence_check()->set_symbol('%');
411 :
412 : // set preconditioner string
413 : std::string s;
414 0 : if(preconditioner().valid())
415 0 : s = std::string(" (Precond: ") + preconditioner()->name() + ")";
416 : else
417 : s = " (No Preconditioner) ";
418 0 : convergence_check()->set_info(s);
419 0 : }
420 :
421 : /// debugger output: solution and residual
422 0 : void write_debugXR(vector_type &x, vector_type &r, int loopCnt, char phase)
423 : {
424 0 : if(!this->vector_debug_writer_valid()) return;
425 0 : std::string ext = GetStringPrintf("-%c_iter%03d", phase, loopCnt);
426 0 : write_debug(r, std::string("BiCGStab_Residual") + ext + ".vec");
427 0 : write_debug(x, std::string("BiCGStab_Solution") + ext + ".vec");
428 : }
429 :
430 : /// debugger section for the preconditioner
431 0 : void enter_precond_debug_section(int loopCnt, char phase)
432 : {
433 0 : if(!this->vector_debug_writer_valid()) return;
434 0 : std::string ext = GetStringPrintf("-%c_iter%03d", phase, loopCnt);
435 0 : this->enter_vector_debug_writer_section(std::string("BiCGStab_Precond") + ext);
436 : }
437 :
438 : public:
439 0 : virtual std::string config_string() const
440 : {
441 0 : std::stringstream ss;
442 0 : ss << "BiCGStab( restart = " << m_numRestarts << ", min_orthogonality = " << m_minOrtho << ")\n";
443 0 : ss << base_type::config_string_preconditioner_convergence_check();
444 0 : return ss.str();
445 0 : }
446 : protected:
447 : /// restarts at every numRestarts steps (numRestarts <= 0 --> never)
448 : int m_numRestarts;
449 :
450 : /// minimal value in (0,1) accepted for Orthoginality before restart
451 : number m_minOrtho;
452 :
453 : /// postprocessor for the correction in the iterations
454 : /**
455 : * These postprocess operations are applied to the preconditioned
456 : * defect before the orthogonalization. The goal is to prevent the
457 : * useless kernel parts to prevail in the (floating point) arithmetics.
458 : */
459 : PProcessChain<vector_type> m_corr_post_process;
460 : };
461 :
462 : } // end namespace ug
463 :
464 : #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__BICGSTAB__ */
|
