Loading [MathJax]/extensions/tex2jax.js
ug4
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
line_search.h
Go to the documentation of this file.
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__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
50
51#endif
52
53
54
55namespace ug{
56
59// Line Search
62
69template <typename TVector>
71{
72 public:
73 typedef TVector vector_type;
74
75 public:
77 virtual void set_offset(std::string offset) = 0;
78
93 vector_type& d, number defect) = 0;
94
96
103 virtual std::string config_string() const = 0;
104
106 virtual ~ILineSearch() {}
107
108#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
109 virtual void setNewtonUpdater( SmartPtr<NewtonUpdaterGeneric<TVector> > nU ) = 0;
110#endif
111};
112
114template <typename TVector>
115class StandardLineSearch : public ILineSearch<TVector>
116{
117 public:
118 // type of
119 typedef TVector vector_type;
120
121 public:
124 : m_maxSteps(10), m_lambdaStart(1.0), m_lambdaReduce(0.5), m_alpha(0.25),
125 m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(false), m_bCheckAll(false), m_offset("")
126#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
127 ,
128 m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
129#endif
130 {};
131
133 StandardLineSearch(int maxSteps, number lambdaStart, number lambdaReduce, bool bAcceptBest)
134 : m_maxSteps(maxSteps), m_lambdaStart(lambdaStart), m_lambdaReduce(lambdaReduce), m_alpha(0.25),
135 m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(bAcceptBest), m_bCheckAll(false), m_offset("")
136#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
137 ,
138 m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
139#endif
140 {};
141
143 StandardLineSearch(int maxSteps, number lambdaStart, number lambdaReduce, bool bAcceptBest, bool bCheckAll)
144 : m_maxSteps(maxSteps), m_lambdaStart(lambdaStart), m_lambdaReduce(lambdaReduce), m_alpha(0.25),
145 m_maxDefect(1e+10), m_verbose(true), m_bAcceptBest(bAcceptBest), m_bCheckAll(bCheckAll), m_offset("")
146#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
147 ,
148 m_newtonUpdater(new NewtonUpdaterGeneric<vector_type>{})
149#endif
150 {};
151
153 virtual std::string config_string() const
154 {
155 std::stringstream ss;
156 ss << "StandardLineSearch( maxSteps = " << m_maxSteps << ", lambdaStart = " << m_lambdaStart << ", lambdaReduce = " << m_lambdaReduce << ", accept best = " <<
157 (m_bAcceptBest ? "true" : "false") << " check all = " << (m_bCheckAll ? "true" : "false");
158 return ss.str();
159
160 }
161
163 void set_maximum_steps(int steps) {m_maxSteps = steps;}
164
166 void set_lambda_start(number start) {m_lambdaStart = start;}
167
169 void set_reduce_factor(number factor) {m_lambdaReduce = factor;}
170
172 void set_suff_descent_factor(number factor) {m_alpha = factor;}
173
175 void set_accept_best(bool bAcceptBest) {m_bAcceptBest = bAcceptBest;}
176
178 void set_check_all(bool bCheckAll) {m_bCheckAll = bCheckAll;}
179
181 void set_maximum_defect(number maxDef) {m_maxDefect = maxDef;}
182
184 void set_verbose(bool level) {m_verbose = level;}
185
187 virtual void set_offset(std::string offset) {m_offset = offset;};
188
192 vector_type& d, number defect)
193 {
194 PROFILE_BEGIN_GROUP(StandardLineSearch_search, ""); // group?
195 // clone pattern for s
196 s.resize(u.size());
197
198 // start factor
199 number lambda = m_lambdaStart;
200
201 // some values
202 number norm, norm_old = defect;
203 bool converged = false;
204 std::vector<number> vRho;
205
206 // remember u
207 s = u;
208
209 // print heading line
210 if(m_verbose)
211 UG_LOG(m_offset << " ++++ Line Search:\n"
212 << " + Iter lambda Defect Rate \n");
213
214
215 // loop line search steps
216 for(int k = 1; k <= m_maxSteps; ++k)
217 {
218#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
219 // try on line u := u - lambda*p
220
221 //if( m_newtonUpdater != NULL )
222 //{
223 bool acceptedNewtonUpdate = m_newtonUpdater->updateSolution(u, 1.0, u, (-1)*lambda, p);
224
225// if( ! m_newtonUpdater->updateNewton(u, 1.0, u, (-1)*lambda, p) )
226 if( ! acceptedNewtonUpdate )
227// if( m_newtonUpdater->updateNewton(u,p) )
228 {
229 UG_LOG("Update in Line Search did not work.\n");
230 norm = std::numeric_limits<number>::max();
231 //VecScaleAdd(u, 1.0, u, (-1)*lambda, p);
232 //return false;
233 }
234 else
235 {
236 // compute new Defect
237 spOp->prepare(u);
238 spOp->apply(d, u);
239
240 // compute new Residuum
241 norm = d.norm();
242 }
243 //}
244// else
245// {
246// UG_LOG("Please set a Newton updater in Line search");
247// }
248#else
249 // try on line u := u - lambda*p
250 VecScaleAdd(u, 1.0, u, (-1)*lambda, p);
251
252 // compute new Defect
253 spOp->prepare(u);
254 spOp->apply(d, u);
255
256 // compute new Residuum
257 norm = d.norm();
258
259#endif
260
261
262 // compute reduction
263 vRho.push_back(norm/norm_old);
264
265 // print rate
266 if(m_verbose)
267 UG_LOG(m_offset << " + " << std::setw(4)
268 << k << ": " << std::setw(11)
269 << std::resetiosflags( ::std::ios::scientific )<<
270 lambda << " "
271 << std::scientific << norm << " " << vRho.back() <<"\n");
272
273 // check if reduction fits
274 if(vRho.back() <= 1 - m_alpha * std::fabs(lambda))
275 {
276 converged = true;
277 if(!m_bCheckAll) break;
278 }
279
280 lambda *= m_lambdaReduce;
281
282 // check if maximum number of steps reached
283 if(k == m_maxSteps)
284 {
285 // if not accept best, line search failed
286 if(!m_bAcceptBest)
287 {
288 UG_LOG(m_offset << " ++++ Line Search did not converge.\n");
289 return false;
290 }
291
292 // search minimum
293 size_t best = 0;
294 number rho_min = vRho.front();
295 for(size_t i = 1; i < vRho.size(); ++i)
296 {
297 if(rho_min > vRho[i])
298 {
299 rho_min = vRho[i];
300 best = i;
301 }
302 }
303
304 /* check if best is converging (i.e. rho < 1)
305 if(vRho[best] >= 1)
306 {
307 UG_LOG(m_offset << " ++++ Accept Best: No try with "
308 "Rate < 1, cannot accept any line search step.\n");
309 UG_LOG(m_offset << " ++++ Line Search did not converge.\n");
310 return false;
311 }*/
312
313 // accept best
314 if(m_verbose)
315 UG_LOG(m_offset << " ++++ Accept Best: Accepting step " <<
316 best+1 << ", Rate = "<< vRho[best] <<".\n");
317
318#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
319 // try on line u := u - lambda*p
320
321 //if( m_newtonUpdater != NULL )
322 //{
323 if( ! m_newtonUpdater->updateSolution(u, 1.0, s, (-1)*m_lambdaStart*std::pow(m_lambdaReduce, (number)best), p) )
324 {
325 UG_LOG("Update in Line Search kmax did not work.\n");
326
327 norm = std::numeric_limits<number>::max();
328 //return false;
329 }
330 else
331 {
332 spOp->prepare(u);
333 spOp->apply(d, u);
334
335 // compute new Residuum
336 norm = d.norm();
337 }
338 //}
339// else
340// {
341// UG_LOG("Please set a Newton updater in Line search");
342// }
343#else
344 // try on line u := u - lambda*p
345 VecScaleAdd(u, 1.0, s, (-1)*m_lambdaStart*std::pow(m_lambdaReduce, (number)best), p);
346 // compute new Defect
347 spOp->prepare(u);
348 spOp->apply(d, u);
349
350 // compute new Residuum
351 norm = d.norm();
352#endif
353
354
355 // check if defect-norm is smaller than maximum allowed defect value
356 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 // reset u and eventual local variables
370 m_newtonUpdater->resetSolution(u,s);
371#else
372 // reset u
373 u = s;
374#endif
375
376 }
377
378 // print end line
379 if(m_verbose)
380 {
381 //only for rate < 1, we call it "Line Search converged"
382 if(converged)
383 UG_LOG(m_offset << " ++++ Line Search converged.\n");
384 }
385
386
387#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
388 if( ! m_newtonUpdater->tellAndFixUpdateEvents(u) )
389 {
390 UG_LOG("unable to fix local Newton updates" << std::endl );
391 return false;
392 }
393#endif
394 // we're done
395 return true;
396 }
397
398
399#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
400
401 virtual void setNewtonUpdater( SmartPtr<NewtonUpdaterGeneric<TVector> > nU )
402 {
403 m_newtonUpdater = nU;
404 }
405
406#endif
407
408 protected:
411
412 protected:
415
418
421
424
427
430
433
436
438 std::string m_offset;
439
440#if ENABLE_NESTED_NEWTON_RESOLFUNC_UPDATE
441private:
443#endif
444};
445
446} // end namespace ug
447
448#endif /* __H__UG__LIB_DISC__OPERATOR__NON_LINEAR_OPERATOR__LINE_SEARCH__ */
parameterString p
Definition smart_pointer.h:108
Definition line_search.h:71
virtual std::string config_string() const =0
returns information about configuration parameters
virtual void set_offset(std::string offset)=0
set string to be printed before each output of line search
virtual ~ILineSearch()
virtual destructor
Definition line_search.h:106
TVector vector_type
Definition line_search.h:73
virtual bool search(SmartPtr< IOperator< vector_type > > spOp, vector_type &u, vector_type &p, vector_type &d, number defect)=0
describes a mapping X->Y
Definition operator.h:86
Definition newtonUpdaterGeneric.h:25
standard implementation of the line search based on the "sufficient descent"
Definition line_search.h:116
void set_reduce_factor(number factor)
sets factor by which line search factor is multiplied in each step
Definition line_search.h:169
std::string m_offset
number of spaces inserted before output
Definition line_search.h:438
number m_lambdaStart
initial step length scaling
Definition line_search.h:417
int m_maxSteps
maximum number of steps to be performed
Definition line_search.h:414
void set_check_all(bool bCheckAll)
sets iff all the max_steps line search steps must be tested even if the sufficient descent is achieve...
Definition line_search.h:178
void set_suff_descent_factor(number factor)
sets the factor controlling the sufficient descent
Definition line_search.h:172
vector_type s
solution in line direction
Definition line_search.h:410
number m_alpha
sufficient descent factor
Definition line_search.h:423
virtual std::string config_string() const
returns information about configuration parameters
Definition line_search.h:153
TVector vector_type
Definition line_search.h:119
number m_maxDefect
maximum allowed defect
Definition line_search.h:426
void set_maximum_defect(number maxDef)
sets maximum allowed norm of the defect (an exception is thrown if this value if exceeded)
Definition line_search.h:181
virtual void set_offset(std::string offset)
set string to be printed before each output of line search
Definition line_search.h:187
void set_maximum_steps(int steps)
sets maximum number of line search steps
Definition line_search.h:163
void set_verbose(bool level)
sets if info should be printed
Definition line_search.h:184
void set_lambda_start(number start)
sets start factor
Definition line_search.h:166
StandardLineSearch(int maxSteps, number lambdaStart, number lambdaReduce, bool bAcceptBest)
constructor
Definition line_search.h:133
StandardLineSearch()
default constructor (setting default values)
Definition line_search.h:123
void set_accept_best(bool bAcceptBest)
sets iff after max_steps the best try is used
Definition line_search.h:175
bool m_bAcceptBest
accept best
Definition line_search.h:432
virtual bool search(SmartPtr< IOperator< vector_type > > spOp, vector_type &u, vector_type &p, vector_type &d, number defect)
Definition line_search.h:190
bool m_bCheckAll
check all
Definition line_search.h:435
bool m_verbose
verbosity level
Definition line_search.h:429
StandardLineSearch(int maxSteps, number lambdaStart, number lambdaReduce, bool bAcceptBest, bool bCheckAll)
constructor
Definition line_search.h:143
number m_lambdaReduce
reduction factor for the step length
Definition line_search.h:420
#define UG_LOG(msg)
Definition log.h:367
double number
Definition types.h:124
void VecScaleAdd(vector_t &vOut, typename vector_t::value_type s1, const vector_t &v1, typename vector_t::value_type s2, const vector_t &v2)
Scales two Vectors, adds them and returns the sum in a third vector.
Definition math_vector_functions_common_impl.hpp:265
if(!(yy_init))
Definition lexer.cpp:997
the ug namespace
#define PROFILE_BEGIN_GROUP(name, groups)
Definition profiler.h:255
function ProblemDisc new(problemDesc, dom)