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_SOLVER__CG__
34 : #define __H__UG__LIB_DISC__OPERATOR__LINEAR_SOLVER__CG__
35 :
36 : #include <iostream>
37 : #include <string>
38 :
39 : #include "lib_algebra/operator/interface/operator.h"
40 : #include "lib_algebra/operator/interface/preconditioned_linear_operator_inverse.h"
41 : #include "common/profiler/profiler.h"
42 : #include "lib_algebra/operator/interface/pprocess.h"
43 : #ifdef UG_PARALLEL
44 : #include "lib_algebra/parallelization/parallelization.h"
45 : #endif
46 :
47 : namespace ug{
48 :
49 : /// the CG method as a solver for linear operators
50 : /**
51 : * This class implements the CG - method for the solution of linear operator
52 : * problems like A*x = b, where the solution x = A^{-1} b is computed.
53 : *
54 : * For detailed description of the algorithm, please refer to:
55 : *
56 : * - Barrett, Berry, Chan, Demmel, Donatom Dongarra, Eijkhout, Pozo, Romine,
57 : * Van der Vorst, "Templates for the Solution of Linear Systems: Building
58 : * Blocks for Iterative Methods", p.13, Fig, 2.5
59 : *
60 : * - Saad, "Iterative Methods For Sparse Linear Systems", p277, Alg. 9.1
61 : *
62 : * \tparam TVector vector type
63 : */
64 : template <typename TVector>
65 : class CG
66 : : public IPreconditionedLinearOperatorInverse<TVector>
67 : {
68 : public:
69 : /// Vector type
70 : typedef TVector vector_type;
71 :
72 : /// Base type
73 : typedef IPreconditionedLinearOperatorInverse<vector_type> base_type;
74 :
75 : protected:
76 : using base_type::convergence_check;
77 : using base_type::linear_operator;
78 : using base_type::preconditioner;
79 : using base_type::write_debug;
80 :
81 : public:
82 : /// constructors
83 0 : CG() : base_type() {}
84 :
85 0 : CG(SmartPtr<ILinearIterator<vector_type,vector_type> > spPrecond)
86 0 : : base_type ( spPrecond ) {}
87 :
88 0 : CG(SmartPtr<ILinearIterator<vector_type,vector_type> > spPrecond, SmartPtr<IConvergenceCheck<vector_type> > spConvCheck)
89 0 : : base_type ( spPrecond, spConvCheck) {}
90 :
91 : /// name of solver
92 0 : virtual const char* name() const {return "CG";}
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 : PROFILE_BEGIN_GROUP(CG_apply_return_defect, "CG algebra");
106 : // check parallel storage types
107 : #ifdef UG_PARALLEL
108 : if(!b.has_storage_type(PST_ADDITIVE) || !x.has_storage_type(PST_CONSISTENT))
109 : UG_THROW("CG::apply_return_defect:"
110 : "Inadequate storage format of Vectors.");
111 : #endif
112 :
113 : // rename r as b (for convenience)
114 : vector_type& r = b;
115 :
116 : // Build defect: r := b - J(u)*x
117 0 : linear_operator()->apply_sub(r, x);
118 :
119 : // create help vector (h will be consistent r)
120 0 : SmartPtr<vector_type> spQ = r.clone_without_values(); vector_type& q = *spQ;
121 0 : SmartPtr<vector_type> spZ = x.clone_without_values(); vector_type& z = *spZ;
122 0 : SmartPtr<vector_type> spP = x.clone_without_values(); vector_type& p = *spP;
123 :
124 0 : write_debugXR(x, r, convergence_check()->step());
125 :
126 : // Preconditioning
127 0 : if(preconditioner().valid())
128 : {
129 0 : enter_precond_debug_section(convergence_check()->step());
130 : // apply z = M^-1 * s
131 0 : if(!preconditioner()->apply(z, r))
132 : {
133 : UG_LOG("ERROR in 'CG::apply_return_defect': "
134 : "Cannot apply preconditioner. Aborting.\n");
135 0 : this->leave_vector_debug_writer_section();
136 0 : return false;
137 : }
138 0 : this->leave_vector_debug_writer_section();
139 : }
140 0 : else z = r;
141 :
142 : // make z consistent
143 : #ifdef UG_PARALLEL
144 : if(!z.change_storage_type(PST_CONSISTENT))
145 : UG_THROW("CG::apply_return_defect: "
146 : "Cannot convert z to consistent vector.");
147 : #endif
148 :
149 : // post-process the correction
150 0 : m_corr_post_process.apply (z);
151 :
152 : // compute start defect
153 0 : prepare_conv_check();
154 0 : convergence_check()->start(r);
155 :
156 : // start search direction
157 0 : p = z;
158 :
159 : // start rho
160 : number rhoOld = VecProd(z, r), rho;
161 :
162 : // Iteration loop
163 0 : while(!convergence_check()->iteration_ended())
164 : {
165 : // Build q = A*p (q is additive afterwards)
166 0 : linear_operator()->apply(q, p);
167 :
168 : // lambda = (q,p)
169 : number lambda = VecProd(q, p);
170 :
171 : // check lambda
172 0 : if(lambda == 0.0)
173 : {
174 0 : if (p.size())
175 : {
176 : UG_LOG("ERROR in 'CG::apply_return_defect': lambda=" <<
177 : lambda<< " is not admitted. Aborting solver.\n");
178 : return false;
179 : }
180 : // in cases where a proc has no geometry, we do not want to fail here
181 : // so we set lambda = 1, this is not harmful
182 : else
183 : lambda = 1.0;
184 : }
185 :
186 : // alpha = rho / (q,p)
187 0 : const number alpha = rhoOld/lambda;
188 :
189 : // Update x := x + alpha*p
190 0 : VecScaleAdd(x, 1.0, x, alpha, p);
191 :
192 : // Update r := r - alpha*t
193 0 : VecScaleAdd(r, 1.0, r, -alpha, q);
194 :
195 0 : write_debugXR(x, r, convergence_check()->step());
196 :
197 : // Check convergence
198 0 : convergence_check()->update(r);
199 0 : if(convergence_check()->iteration_ended()) break;
200 :
201 : // Preconditioning
202 0 : if(preconditioner().valid())
203 : {
204 0 : enter_precond_debug_section(convergence_check()->step());
205 : // apply z = M^-1 * r
206 0 : if(!preconditioner()->apply(z, r))
207 : {
208 : UG_LOG("ERROR in 'CG::apply_return_defect': "
209 : "Cannot apply preconditioner. Aborting.\n");
210 0 : this->leave_vector_debug_writer_section();
211 0 : return false;
212 : }
213 0 : this->leave_vector_debug_writer_section();
214 : }
215 0 : else z = r;
216 :
217 : #ifdef UG_PARALLEL
218 : // make z consistent
219 : if(!z.change_storage_type(PST_CONSISTENT))
220 : UG_THROW("CG::apply_return_defect': "
221 : "Cannot convert z to consistent vector.");
222 : #endif
223 :
224 : // post-process the correction
225 0 : m_corr_post_process.apply (z);
226 :
227 : // new rho = (z,r)
228 : rho = VecProd(z, r);
229 :
230 : // new beta = rho / rhoOld
231 0 : const number beta = rho/rhoOld;
232 :
233 : // new direction p := beta * p + z
234 0 : VecScaleAdd(p, beta, p, 1.0, z);
235 :
236 : // remember old rho
237 : rhoOld = rho;
238 : }
239 :
240 : // post output
241 0 : return convergence_check()->post();
242 : }
243 :
244 : /// adds a post-process for the iterates
245 0 : void add_postprocess_corr (SmartPtr<IPProcessVector<vector_type> > p)
246 : {
247 0 : m_corr_post_process.add (p);
248 0 : }
249 :
250 : /// removes a post-process for the iterates
251 0 : void remove_postprocess_corr (SmartPtr<IPProcessVector<vector_type> > p)
252 : {
253 0 : m_corr_post_process.remove (p);
254 0 : }
255 :
256 : protected:
257 : /// adjust output of convergence check
258 0 : void prepare_conv_check()
259 : {
260 : // set iteration symbol and name
261 0 : convergence_check()->set_name(name());
262 0 : convergence_check()->set_symbol('%');
263 :
264 : // set preconditioner string
265 : std::string s;
266 0 : if(preconditioner().valid())
267 0 : s = std::string(" (Precond: ") + preconditioner()->name() + ")";
268 : else
269 : s = " (No Preconditioner) ";
270 0 : convergence_check()->set_info(s);
271 0 : }
272 :
273 : /// debugger output: solution and residual
274 0 : void write_debugXR(vector_type &x, vector_type &r, int loopCnt)
275 : {
276 0 : if(!this->vector_debug_writer_valid()) return;
277 : char ext[20]; snprintf(ext, 20, "_iter%03d", loopCnt);
278 0 : write_debug(r, std::string("CG_Residual") + ext + ".vec");
279 0 : write_debug(x, std::string("CG_Solution") + ext + ".vec");
280 : }
281 :
282 : /// debugger section for the preconditioner
283 0 : void enter_precond_debug_section(int loopCnt)
284 : {
285 0 : if(!this->vector_debug_writer_valid()) return;
286 : char ext[20]; snprintf(ext, 20, "_iter%03d", loopCnt);
287 0 : this->enter_vector_debug_writer_section(std::string("CG_Precond_") + ext);
288 : }
289 :
290 : protected:
291 : number VecProd(vector_type& a, vector_type& b)
292 : {
293 : return a.dotprod(b);
294 : }
295 :
296 : protected:
297 : /// postprocessor for the correction in the iterations
298 : /**
299 : * These postprocess operations are applied to the preconditioned
300 : * defect before the orthogonalization. The goal is to prevent the
301 : * useless kernel parts to prevail in the (floating point) arithmetics.
302 : */
303 : PProcessChain<vector_type> m_corr_post_process;
304 : };
305 :
306 : } // end namespace ug
307 :
308 : #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_SOLVER__CG__ */
|
