Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel, modifications nested Newton: Markus Knodel
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__NON_LINEAR_OPERATOR__LINE_SEARCH__
34 : #define __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__LINE_SEARCH__
35 :
36 : #include <ostream>
37 : #include <string>
38 : #include <vector>
39 : #include <cmath>
40 : #include <sstream>
41 : #include <limits>
42 :
43 : #include "common/common.h"
44 :
45 : //#include "lib_disc/operator/non_linear_operator/newton_solver/nestedNewtonRFSwitch.h"
46 :
47 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
48 :
49 : #include "lib_disc/operator/non_linear_operator/newton_solver/newtonUpdaterGeneric.h"
50 :
51 : //#endif
52 :
53 :
54 :
55 : namespace ug{
56 :
57 : ///////////////////////////////////////////////////////////
58 : ///////////////////////////////////////////////////////////
59 : // Line Search
60 : ///////////////////////////////////////////////////////////
61 : ///////////////////////////////////////////////////////////
62 :
63 : /** LineSearch
64 : *
65 : * This is the base class for a line search object. An instance is
66 : * passed to a newton solver to control the line search.
67 : *
68 : */
69 : template <typename TVector>
70 : class ILineSearch
71 : {
72 : public:
73 : typedef TVector vector_type;
74 :
75 : public:
76 : /// set string to be printed before each output of line search
77 : virtual void set_offset(std::string offset) = 0;
78 :
79 : /**
80 : * Performs a line search to a given direction.
81 : *
82 : * \param[in] Op Non-linear operator
83 : * \param[in] u current solution
84 : * \param[in] p search direction
85 : * \param[in,out] d defect
86 : * \param[in] defect norm of current defect
87 : *
88 : * \return true if line search successful
89 : * false if line search failed
90 : */
91 : virtual bool search(SmartPtr<IOperator<vector_type> > spOp,
92 : vector_type& u, vector_type& p,
93 : vector_type& d, number defect) = 0;
94 :
95 : /// returns information about configuration parameters
96 : /**
97 : * this should return necessary information about parameters and possibly
98 : * calling config_string of subcomponents.
99 : *
100 : * \returns std::string necessary information about configuration parameters
101 : */
102 :
103 : virtual std::string config_string() const = 0;
104 :
105 : /// virtual destructor
106 : virtual ~ILineSearch() {}
107 :
108 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
109 : virtual void setNewtonUpdater( SmartPtr<NewtonUpdaterGeneric<TVector> > nU ) = 0;
110 : //#endif
111 : virtual bool createNewtonUpdater() = 0;
112 :
113 : };
114 :
115 : /// standard implementation of the line search based on the "sufficient descent"
116 : template <typename TVector>
117 : class StandardLineSearch : public ILineSearch<TVector>
118 : {
119 : public:
120 : // type of
121 : typedef TVector vector_type;
122 :
123 : public:
124 : /// default constructor (setting default values)
125 0 : StandardLineSearch()
126 0 : : m_maxSteps(10), m_lambdaStart(1.0), m_lambdaReduce(0.5), m_alpha(0.25),
127 0 : m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(false), m_bCheckAll(false), m_offset(""),
128 0 : m_newtonUpdater(SPNULL)
129 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
130 : // ,
131 : // m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
132 : //#endif
133 0 : {};
134 :
135 : /// constructor
136 0 : StandardLineSearch(int maxSteps, number lambdaStart, number lambdaReduce, bool bAcceptBest)
137 0 : : m_maxSteps(maxSteps), m_lambdaStart(lambdaStart), m_lambdaReduce(lambdaReduce), m_alpha(0.25),
138 0 : m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(bAcceptBest), m_bCheckAll(false), m_offset(""),
139 0 : m_newtonUpdater(SPNULL)
140 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
141 : // ,
142 : // m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
143 : //#endif
144 0 : {};
145 :
146 : /// constructor
147 0 : StandardLineSearch(int maxSteps, number lambdaStart, number lambdaReduce, bool bAcceptBest, bool bCheckAll)
148 0 : : m_maxSteps(maxSteps), m_lambdaStart(lambdaStart), m_lambdaReduce(lambdaReduce), m_alpha(0.25),
149 0 : m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(bAcceptBest), m_bCheckAll(bCheckAll), m_offset(""),
150 0 : m_newtonUpdater(SPNULL)
151 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
152 : // ,
153 : // m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
154 : //#endif
155 0 : {};
156 :
157 : /// returns information about configuration parameters
158 0 : virtual std::string config_string() const
159 : {
160 0 : std::stringstream ss;
161 0 : ss << "StandardLineSearch( maxSteps = " << m_maxSteps << ", lambdaStart = " << m_lambdaStart << ", lambdaReduce = " << m_lambdaReduce << ", accept best = " <<
162 0 : (m_bAcceptBest ? "true" : "false") << " check all = " << (m_bCheckAll ? "true" : "false");
163 0 : return ss.str();
164 :
165 0 : }
166 :
167 : /// sets maximum number of line search steps
168 0 : void set_maximum_steps(int steps) {m_maxSteps = steps;}
169 :
170 : /// sets start factor
171 0 : void set_lambda_start(number start) {m_lambdaStart = start;}
172 :
173 : /// sets factor by which line search factor is multiplied in each step
174 0 : void set_reduce_factor(number factor) {m_lambdaReduce = factor;}
175 :
176 : /// sets the factor controlling the sufficient descent
177 0 : void set_suff_descent_factor(number factor) {m_alpha = factor;}
178 :
179 : /// sets iff after max_steps the best try is used
180 0 : void set_accept_best(bool bAcceptBest) {m_bAcceptBest = bAcceptBest;}
181 :
182 : /// sets iff all the max_steps line search steps must be tested even if the sufficient descent is achieved
183 0 : void set_check_all(bool bCheckAll) {m_bCheckAll = bCheckAll;}
184 :
185 : /// sets maximum allowed norm of the defect (an exception is thrown if this value if exceeded)
186 0 : void set_maximum_defect(number maxDef) {m_maxDefect = maxDef;}
187 :
188 : /// sets if info should be printed
189 0 : void set_verbose(bool level) {m_verbose = level;}
190 :
191 : /// \copydoc ILineSearch::set_offset
192 0 : virtual void set_offset(std::string offset) {m_offset = offset;};
193 :
194 : /// \copydoc ILineSearch::search
195 0 : virtual bool search(SmartPtr<IOperator<vector_type> > spOp,
196 : vector_type& u, vector_type& p,
197 : vector_type& d, number defect)
198 : {
199 : PROFILE_BEGIN_GROUP(StandardLineSearch_search, ""); // group?
200 : // clone pattern for s
201 0 : s.resize(u.size());
202 :
203 : // start factor
204 0 : number lambda = m_lambdaStart;
205 :
206 : // some values
207 : number norm, norm_old = defect;
208 : bool converged = false;
209 : std::vector<number> vRho;
210 :
211 : // remember u
212 0 : s = u;
213 :
214 : // print heading line
215 0 : if(m_verbose)
216 : UG_LOG(m_offset << " ++++ Line Search:\n"
217 : << " + Iter lambda Defect Rate \n");
218 :
219 :
220 : // loop line search steps
221 0 : for(int k = 1; k <= m_maxSteps; ++k)
222 : {
223 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
224 0 : if( m_newtonUpdater != SPNULL )
225 : {
226 : // try on line u := u - lambda*p
227 :
228 0 : bool acceptedNewtonUpdate = m_newtonUpdater->updateSolution(u, 1.0, u, (-1)*lambda, p);
229 :
230 0 : if( ! acceptedNewtonUpdate )
231 : {
232 : UG_LOG("Update in Line Search did not work.\n");
233 : norm = std::numeric_limits<number>::max();
234 : //VecScaleAdd(u, 1.0, u, (-1)*lambda, p);
235 : //return false;
236 : }
237 : else
238 : {
239 : // compute new Defect
240 0 : spOp->prepare(u);
241 0 : spOp->apply(d, u);
242 :
243 : // compute new Residuum
244 0 : norm = d.norm();
245 : }
246 : }
247 : else
248 : {
249 : //#else
250 : // try on line u := u - lambda*p
251 0 : VecScaleAdd(u, 1.0, u, (-1)*lambda, p);
252 :
253 : // compute new Defect
254 0 : spOp->prepare(u);
255 0 : spOp->apply(d, u);
256 :
257 : // compute new Residuum
258 0 : norm = d.norm();
259 : }
260 : //#endif
261 :
262 :
263 : // compute reduction
264 0 : vRho.push_back(norm/norm_old);
265 :
266 : // print rate
267 0 : if(m_verbose)
268 0 : UG_LOG(m_offset << " + " << std::setw(4)
269 : << k << ": " << std::setw(11)
270 : << std::resetiosflags( ::std::ios::scientific )<<
271 : lambda << " "
272 : << std::scientific << norm << " " << vRho.back() <<"\n");
273 :
274 : // check if reduction fits
275 0 : if(vRho.back() <= 1 - m_alpha * std::fabs(lambda))
276 : {
277 : converged = true;
278 0 : if(!m_bCheckAll) break;
279 : }
280 :
281 0 : lambda *= m_lambdaReduce;
282 :
283 : // check if maximum number of steps reached
284 0 : if(k == m_maxSteps)
285 : {
286 : // if not accept best, line search failed
287 0 : if(!m_bAcceptBest)
288 : {
289 : UG_LOG(m_offset << " ++++ Line Search did not converge.\n");
290 : return false;
291 : }
292 :
293 : // search minimum
294 : size_t best = 0;
295 0 : number rho_min = vRho.front();
296 0 : for(size_t i = 1; i < vRho.size(); ++i)
297 : {
298 0 : if(rho_min > vRho[i])
299 : {
300 : rho_min = vRho[i];
301 : best = i;
302 : }
303 : }
304 :
305 : /* check if best is converging (i.e. rho < 1)
306 : if(vRho[best] >= 1)
307 : {
308 : UG_LOG(m_offset << " ++++ Accept Best: No try with "
309 : "Rate < 1, cannot accept any line search step.\n");
310 : UG_LOG(m_offset << " ++++ Line Search did not converge.\n");
311 : return false;
312 : }*/
313 :
314 : // accept best
315 0 : if(m_verbose)
316 0 : UG_LOG(m_offset << " ++++ Accept Best: Accepting step " <<
317 : best+1 << ", Rate = "<< vRho[best] <<".\n");
318 :
319 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
320 0 : if( m_newtonUpdater != SPNULL )
321 : {
322 : // try on line u := u - lambda*p
323 :
324 0 : if( ! m_newtonUpdater->updateSolution(u, 1.0, s, (-1)*m_lambdaStart*std::pow(m_lambdaReduce, (number)best), p) )
325 : {
326 : UG_LOG("Update in Line Search kmax did not work.\n");
327 :
328 : norm = std::numeric_limits<number>::max();
329 : //return false;
330 : }
331 : else
332 : {
333 0 : spOp->prepare(u);
334 0 : spOp->apply(d, u);
335 :
336 : // compute new Residuum
337 0 : norm = d.norm();
338 : }
339 : }
340 : else
341 : {
342 : //#else
343 : // try on line u := u - lambda*p
344 0 : VecScaleAdd(u, 1.0, s, (-1)*m_lambdaStart*std::pow(m_lambdaReduce, (number)best), p);
345 : // compute new Defect
346 0 : spOp->prepare(u);
347 0 : spOp->apply(d, u);
348 :
349 : // compute new Residuum
350 0 : norm = d.norm();
351 : }
352 : //#endif
353 :
354 :
355 : // check if defect-norm is smaller than maximum allowed defect value
356 0 : if (norm > m_maxDefect)
357 : {
358 : UG_LOG("ERROR in 'StandardLineSearch::search':"
359 : " maximum defect-limit is reached.\n");
360 : return false;
361 : }
362 :
363 : // break to finish
364 : break;
365 : }
366 :
367 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
368 :
369 0 : if( m_newtonUpdater != SPNULL )
370 : {
371 : // reset u and eventual local variables
372 0 : m_newtonUpdater->resetSolution(u,s);
373 : }
374 : else
375 : {
376 : //#else
377 : // reset u
378 0 : u = s;
379 : }
380 : //#endif
381 :
382 : }
383 :
384 : // print end line
385 0 : if(m_verbose)
386 : {
387 : //only for rate < 1, we call it "Line Search converged"
388 0 : if(converged)
389 : UG_LOG(m_offset << " ++++ Line Search converged.\n");
390 : }
391 :
392 :
393 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
394 0 : if( m_newtonUpdater != SPNULL )
395 : {
396 0 : if( ! m_newtonUpdater->tellAndFixUpdateEvents(u) )
397 : {
398 : UG_LOG("unable to fix local Newton updates" << std::endl );
399 : return false;
400 : }
401 : }
402 : //#endif
403 : // we're done
404 : return true;
405 0 : }
406 :
407 :
408 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
409 :
410 0 : virtual void setNewtonUpdater( SmartPtr<NewtonUpdaterGeneric<TVector> > nU )
411 : {
412 0 : m_newtonUpdater = nU;
413 0 : }
414 :
415 0 : virtual bool createNewtonUpdater()
416 : {
417 0 : if( m_newtonUpdater != SPNULL )
418 : {
419 0 : m_newtonUpdater = SmartPtr<NewtonUpdaterGeneric<TVector> >
420 0 : (new NewtonUpdaterGeneric<TVector>{});
421 :
422 0 : return true;
423 : }
424 :
425 : return false;
426 :
427 : }
428 :
429 :
430 : //#endif
431 :
432 : protected:
433 : /// solution in line direction
434 : vector_type s;
435 :
436 : protected:
437 : /// maximum number of steps to be performed
438 : int m_maxSteps;
439 :
440 : /// initial step length scaling
441 : number m_lambdaStart;
442 :
443 : /// reduction factor for the step length
444 : number m_lambdaReduce;
445 :
446 : /// sufficient descent factor
447 : number m_alpha;
448 :
449 : /// maximum allowed defect
450 : number m_maxDefect;
451 :
452 : /// verbosity level
453 : bool m_verbose;
454 :
455 : /// accept best
456 : bool m_bAcceptBest;
457 :
458 : /// check all
459 : bool m_bCheckAll;
460 :
461 : /// number of spaces inserted before output
462 : std::string m_offset;
463 :
464 : //#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
465 : private:
466 : SmartPtr<NewtonUpdaterGeneric<TVector> > m_newtonUpdater;
467 : //#endif
468 : };
469 :
470 : } // end namespace ug
471 :
472 : #endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__LINE_SEARCH__ */
|
