Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Martin Rupp
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__AUTO_LINEAR_OPERATOR__LINEAR_SOLVER__
34 : #define __H__UG__LIB_DISC__OPERATOR__AUTO_LINEAR_OPERATOR__LINEAR_SOLVER__
35 : #include <iostream>
36 : #include <string>
37 :
38 : #include "lib_algebra/operator/interface/preconditioned_linear_operator_inverse.h"
39 : #include "lib_algebra/operator/interface/linear_solver_profiling.h"
40 : #ifdef UG_PARALLEL
41 : #include "lib_algebra/parallelization/parallelization.h"
42 : #endif
43 : #include "common/util/ostream_util.h"
44 : namespace ug{
45 :
46 : // this is for matrix operators which recalculate their "real" operator
47 : // if necessary
48 : class UpdateableMatrixOperator
49 : {
50 : public:
51 : virtual ~UpdateableMatrixOperator() {}
52 : virtual void calculate_matrix() = 0;
53 : };
54 :
55 : /// linear solver using abstract preconditioner interface
56 : /**
57 : * This class is a linear iterating scheme, that uses any implementation
58 : * of the ILinearIterator interface to precondition the iteration.
59 : *
60 : * \tparam TAlgebra algebra type
61 : */
62 : template <typename TVector>
63 : class AutoLinearSolver
64 : : public IPreconditionedLinearOperatorInverse<TVector>
65 : {
66 : public:
67 : /// Vector type
68 : typedef TVector vector_type;
69 :
70 : /// Base type
71 : typedef IPreconditionedLinearOperatorInverse<vector_type> base_type;
72 :
73 : protected:
74 : using base_type::convergence_check;
75 : using base_type::linear_operator;
76 : using base_type::preconditioner;
77 : using base_type::write_debug;
78 :
79 : public:
80 0 : AutoLinearSolver //(double reductionAlwaysAccept, double worseThenAverage, double desiredDefect)
81 : (double desiredDefect, double desiredReduction)
82 0 : {
83 0 : m_bInited=-1;
84 : // m_N = 0;
85 : // m_reductionAlwaysAccept = reductionAlwaysAccept;
86 : // m_worseThenAverage = worseThenAverage;
87 0 : m_desiredDefect = desiredDefect;
88 0 : m_desiredReduction = desiredReduction;
89 0 : m_initCalled = m_initsDone = 0;
90 0 : m_savedTime = 0.0;
91 0 : }
92 0 : AutoLinearSolver()
93 0 : {
94 0 : m_bInited=-1;
95 : // m_N = 0;
96 0 : m_reductionAlwaysAccept = 0.001;
97 0 : m_worseThenAverage = 2.0;
98 0 : m_initCalled = m_initsDone = 0;
99 0 : m_savedTime = 0.0;
100 0 : }
101 :
102 : /// returns if parallel solving is supported
103 0 : virtual bool supports_parallel() const
104 : {
105 0 : if(preconditioner().valid())
106 0 : return preconditioner()->supports_parallel();
107 : else return true;
108 : }
109 :
110 : private:
111 : double m_reductionAlwaysAccept;
112 : double m_worseThenAverage;
113 : int m_bInited;
114 : // size_t m_N;
115 : SmartPtr<ILinearOperator<vector_type,vector_type> > pJ;
116 : vector_type m_u;
117 : double m_avgReduction;
118 : size_t m_initsDone;
119 : size_t m_initCalled;
120 :
121 : double m_lastInitTime;
122 : double m_reductionPerTime, m_lastCallTime, m_lastCallReduction;
123 : double m_desiredReduction, m_desiredDefect;
124 : double m_savedTime;
125 :
126 : public:
127 :
128 0 : void set_reduction_always_accept(double d)
129 : {
130 0 : m_reductionAlwaysAccept = d;
131 0 : }
132 :
133 0 : void set_reinit_when_worse_then_average(double d)
134 : {
135 0 : m_worseThenAverage = d;
136 0 : }
137 :
138 : /// returns the name of the solver
139 0 : virtual const char* name() const {return "Auto Iterative Linear Solver";}
140 :
141 :
142 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type,vector_type> > J, const vector_type& u)
143 : {
144 0 : m_u = u;
145 0 : return init_op(J);
146 : }
147 0 : virtual bool init(SmartPtr<ILinearOperator<vector_type,vector_type> > J)
148 : {
149 : m_u.resize(0);
150 0 : return init_op(J);
151 : }
152 :
153 0 : bool init_op(SmartPtr<ILinearOperator<vector_type,vector_type> > J)
154 : {
155 0 : ILinearOperatorInverse<vector_type, vector_type>::init(J);
156 : //UG_LOG("ALS:init\n");
157 0 : m_initCalled++;
158 0 : if(m_bInited == -1)
159 : {
160 0 : m_bInited=false;
161 0 : pJ = J;
162 0 : reinit();
163 : }
164 : else
165 : {
166 0 : m_bInited=false;
167 0 : pJ = J;
168 : // if(u.size() != m_N)
169 : // reinit(u);
170 : // else
171 : // m_u = u;
172 : }
173 :
174 0 : return true;
175 : }
176 :
177 0 : bool reinit()
178 : {
179 0 : if(m_bInited) return true;
180 0 : m_bInited = true;
181 : //UG_LOG("ALS:reinit\n");
182 :
183 : double tStart = get_clock_s();
184 : // this does not work with stuff.
185 0 : SmartPtr<UpdateableMatrixOperator> uo = pJ.template cast_dynamic<UpdateableMatrixOperator> ();
186 0 : if(uo.valid()) uo->calculate_matrix();
187 0 : if(m_u.size() != 0.0)
188 0 : { if(!base_type::init(pJ, m_u)) return false; }
189 0 : else if(!base_type::init(pJ)) return false;
190 0 : m_lastInitTime = get_clock_s()-tStart;
191 :
192 0 : m_initsDone++;
193 :
194 : // m_N = u.size();
195 :
196 0 : return true;
197 : }
198 :
199 0 : bool compute_correction(vector_type &c, vector_type &d)
200 : {
201 : // Compute a correction c := B*d using one iterative step
202 : // Internally the defect is updated d := d - A*c = b - A*(x+c)
203 0 : if(preconditioner().valid()) {
204 : LS_PROFILE_BEGIN(LS_ApplyPrecond);
205 :
206 0 : if(!preconditioner()->apply(c, d))
207 : {
208 : UG_LOG("ERROR in 'LinearSolver::apply': Iterator "
209 : "Operator applied incorrectly. Aborting.\n");
210 0 : return false;
211 : }
212 0 : linear_operator()->apply_sub(d, c);
213 : LS_PROFILE_END(LS_ApplyPrecond);
214 : }
215 :
216 : return true;
217 : }
218 :
219 0 : virtual bool apply(vector_type& x, const vector_type& b)
220 : {
221 : //UG_LOG("ALS:apply\n");
222 0 : SmartPtr<vector_type> spB = b.clone_without_values();
223 0 : *spB = b;
224 0 : return apply_return_defect(x, *spB);
225 : }
226 : /// solves the system and returns the last defect
227 0 : virtual bool apply_return_defect(vector_type& x, vector_type& b)
228 : {
229 : //UG_LOG("ALS:return_defect\n");
230 : LS_PROFILE_BEGIN(LS_ApplyReturnDefect);
231 :
232 : #ifdef UG_PARALLEL
233 : if(!b.has_storage_type(PST_ADDITIVE) || !x.has_storage_type(PST_CONSISTENT))
234 : UG_THROW("LinearSolver::apply: Inadequate storage format of Vectors.");
235 : #endif
236 :
237 : // rename b as d (for convenience)
238 : vector_type& d = b;
239 :
240 : // build defect: d := b - J(u)*x
241 : LS_PROFILE_BEGIN(LS_BuildDefect);
242 0 : linear_operator()->apply_sub(d, x);
243 : LS_PROFILE_END(LS_BuildDefect);
244 :
245 : // create correction
246 : LS_PROFILE_BEGIN(LS_CreateCorrection);
247 0 : SmartPtr<vector_type> spC = x.clone_without_values();
248 : vector_type& c = *spC;
249 : LS_PROFILE_END(LS_CreateCorrection);
250 :
251 : LS_PROFILE_BEGIN(LS_ComputeStartDefect);
252 0 : prepare_conv_check();
253 0 : convergence_check()->start(d);
254 : LS_PROFILE_END(LS_ComputeStartDefect);
255 :
256 : int loopCnt = 0;
257 : char ext[20]; snprintf(ext, 20, "_iter%03d", loopCnt);
258 0 : std::string name("LS_Defect_"); name.append(ext).append(".vec");
259 0 : write_debug(d, name.c_str());
260 0 : name = std::string("LS_Solution_"); name.append(ext).append(".vec");
261 0 : write_debug(x, name.c_str());
262 :
263 : // Iteration loop
264 :
265 0 : if(m_bInited == false)
266 : {
267 0 : vector_type x2 = x;
268 0 : vector_type d2 = d;
269 : try{
270 : double tStartIterationTime = get_clock_s();
271 0 : while(!convergence_check()->iteration_ended())
272 : {
273 : //double tComputeTime = get_clock_s();
274 0 : if( !compute_correction(c, d) ) return false;
275 0 : x += c;
276 0 : convergence_check()->update(d);
277 :
278 :
279 0 : double r = convergence_check()->rate();
280 0 : double reduction = convergence_check()->reduction();
281 0 : double defect = convergence_check()->defect();
282 0 : double steps = convergence_check()->step();
283 : // if(r > m_reductionAlwaysAccept
284 : // && (r >= 1 || r * m_worseThenAverage > m_avgReduction))
285 0 : double spentTime = get_clock_s() - tStartIterationTime;
286 :
287 0 : double tComputeTime = spentTime/steps; //get_clock_s()-tComputeTime;
288 :
289 : // double timeForOriginalAlgorithm = reduction/m_reductionPerTime + m_lastInitTime;
290 0 : double approxSteps =
291 0 : std::min(log(m_desiredDefect/defect)/log(r),
292 0 : log(m_desiredReduction/reduction)/log(r) );
293 0 : double approxRemainingTimeForSolution = approxSteps*tComputeTime;
294 : /*UG_LOG("\n");
295 : UG_LOG("m_lastCallReduction = " << reset_floats << m_lastCallReduction << "\n");
296 : UG_LOG("reduction = " << reset_floats << reduction << "\n");
297 : UG_LOG("reduction remain = " << reset_floats << m_lastCallReduction/reduction << "\n");
298 : UG_LOG("approxSteps = " << reset_floats << approxSteps << "\n");
299 : UG_LOG("approxTimeForSolution = " << approxTimeForSolution << "\n");
300 : UG_LOG("spentTime = " << spentTime << "\n");
301 : UG_LOG("tComputeTime = " << tComputeTime << "\n");
302 : UG_LOG("rate = " << r << "\n");
303 : UG_LOG("reduction = " << reduction << "\n");*/
304 0 : if(r > 1)
305 : {
306 0 : m_savedTime -= spentTime;
307 : UG_LOG("AutoLinearSolver: REINIT because reduction rate >= 1.")
308 0 : UG_LOG(" [ Inits called: " << m_initCalled << ", inits done: " << reset_floats << m_initsDone+1
309 : << " (" << (100.0*(m_initsDone+1))/m_initCalled << " %), Total Time saved: " << m_savedTime << " s]\n");
310 :
311 :
312 0 : reinit();
313 : break;
314 : }
315 0 : double approxResolveTime = m_lastCallTime + m_lastInitTime;
316 0 : if(approxResolveTime < approxRemainingTimeForSolution)
317 : {
318 0 : m_savedTime -= spentTime;
319 : UG_LOG("AutoLinearSolver: REINIT because\n" << reset_floats )
320 : UG_LOG(" approximated remaining Time for Solution with old preconditioner = " << approxRemainingTimeForSolution << " s\n");
321 : UG_LOG(" > approximated Time for Reinit preconditioner and solve = " << approxResolveTime << " s\n");
322 : UG_LOG(" with old preconditioner:\n");
323 0 : UG_LOG(" reduction rate = " << r << " avg reduction = " << convergence_check()->avg_rate() << "\n")
324 : UG_LOG(" steps = " << steps << ", reduction = " << reduction << "\n")
325 : UG_LOG(" spentTime = " << spentTime << ", time per Step = " << tComputeTime << "\n");
326 : UG_LOG(" approximation of solution with old:\n")
327 0 : UG_LOG(" reduction remain = " << reset_floats << m_lastCallReduction/reduction << "\n");
328 : UG_LOG(" approxSteps with last reduction rate = " << approxSteps << "\n");
329 0 : UG_LOG(" [ Inits called: " << m_initCalled << ", inits done: " << reset_floats << m_initsDone+1
330 : << " (" << (100.0*(m_initsDone+1))/m_initCalled << " %), Total Time saved: " << m_savedTime << " s]\n");
331 0 : reinit();
332 :
333 : break;
334 : }
335 : else
336 : {
337 0 : x2 = x;
338 0 : d2 = d;
339 : }
340 : }
341 0 : if(!m_bInited)
342 : {
343 0 : double spentTime = get_clock_s() - tStartIterationTime;
344 0 : m_savedTime += m_lastCallTime+m_lastInitTime-spentTime;
345 : UG_LOG("AutoLinearSolver solved with old preconditioner [")
346 0 : UG_LOG(" Inits called: " << m_initCalled << ", inits done: " << reset_floats << m_initsDone+1
347 : << " (" << (100.0*(m_initsDone+1))/m_initCalled << " %), Total Time saved: " << m_savedTime << " s]\n");
348 : UG_LOG(" time spent with old = " << reset_floats << spentTime << " s, ");
349 0 : UG_LOG(" last time init + solve = " << m_lastCallTime+m_lastInitTime << " s, ");
350 :
351 0 : if(m_lastCallTime+m_lastInitTime > spentTime)
352 0 : { UG_LOG("saved " << m_lastCallTime+m_lastInitTime - spentTime << " s!.\n"); }
353 : else
354 0 : { UG_LOG("spent too much time with old! " << spentTime - m_lastCallTime+m_lastInitTime<< " s!\n"); }
355 :
356 : }
357 :
358 : }
359 0 : catch(...)
360 : {
361 :
362 : }
363 0 : x = x2;
364 0 : d = d2;
365 : }
366 :
367 0 : if(!convergence_check()->iteration_ended())
368 : {
369 : double T = get_clock_s();
370 0 : convergence_check()->start(d);
371 0 : while(!convergence_check()->iteration_ended())
372 : {
373 0 : if( !compute_correction(c, d) ) return false;
374 0 : x += c;
375 0 : convergence_check()->update(d);
376 : }
377 0 : m_avgReduction = convergence_check()->avg_rate();
378 0 : T = get_clock_s() - T;
379 0 : m_reductionPerTime = convergence_check()->reduction()/T;
380 0 : m_lastCallTime = T;
381 0 : m_lastCallReduction = convergence_check()->reduction();
382 : }
383 :
384 :
385 : // write some information when ending the iteration
386 0 : if(!convergence_check()->post())
387 : {
388 : UG_LOG("ERROR in 'LinearSolver::apply': post-convergence-check "
389 : "signaled failure. Aborting.\n");
390 : return false;
391 : }
392 :
393 : // end profiling of whole function
394 : LS_PROFILE_END(LS_ApplyReturnDefect);
395 :
396 : // we're done
397 : return true;
398 : }
399 :
400 0 : void print_information()
401 : {
402 : UG_LOG("AutoLinearSolver:\n");
403 0 : UG_LOG(" avg reduction is " << m_avgReduction << "\n");
404 0 : UG_LOG(" Inits called: " << m_initCalled << ", inits done: " << reset_floats << m_initsDone
405 : << " (" << (100.0*(m_initsDone))/m_initCalled << " %)\n");
406 0 : UG_LOG(" m_reductionPerTime = " << m_reductionPerTime << " s\n");
407 0 : UG_LOG(" reductionTime for 0.1 reduction = " << 0.1/m_reductionPerTime << " s\n");
408 0 : UG_LOG(" m_lastInitTime = " << m_lastInitTime << " s\n");
409 0 : UG_LOG(" SAVED TIME: " << m_savedTime << " s\n");
410 0 : }
411 :
412 : protected:
413 : /// prepares the convergence check output
414 0 : void prepare_conv_check()
415 : {
416 0 : convergence_check()->set_name(name());
417 0 : convergence_check()->set_symbol('%');
418 0 : if(preconditioner().valid())
419 : {
420 : std::string s;
421 0 : if(preconditioner().valid())
422 0 : s = std::string(" (Precond: ") + preconditioner()->name() + ")";
423 : else
424 : s = " (No Preconditioner) ";
425 0 : convergence_check()->set_info(s);
426 : }
427 0 : }
428 : };
429 :
430 : } // end namespace ug
431 :
432 : #endif /* __H__UG__LIB_DISC__OPERATOR__LINEAR_OPERATOR__LINEAR_SOLVER__ */
|
