24 #ifndef _REDUCEDNEWTON__ALGORITHM_H_
25 #define _REDUCEDNEWTON__ALGORITHM_H_
49 template <
typename PROBLEM,
typename VECTOR>
145 {
return gradient*gradient_transposed;}
147 unsigned int _nonlinear_maxiter, _linear_maxiter, _line_maxiter;
148 double _nonlinear_tol, _nonlinear_global_tol, _linear_tol, _linear_global_tol, _linesearch_rho, _linesearch_c;
149 bool _compute_functionals_in_every_step;
150 std::string _postindex;
156 using namespace dealii;
160 template <
typename PROBLEM,
typename VECTOR>
163 param_reader.
SetSubsection(
"reducednewtonalgorithm parameters");
164 param_reader.
declare_entry(
"nonlinear_maxiter",
"10",Patterns::Integer(0));
165 param_reader.
declare_entry(
"nonlinear_tol",
"1.e-7",Patterns::Double(0));
166 param_reader.
declare_entry(
"nonlinear_global_tol",
"1.e-11",Patterns::Double(0));
168 param_reader.
declare_entry(
"linear_maxiter",
"40",Patterns::Integer(0));
169 param_reader.
declare_entry(
"linear_tol",
"1.e-10",Patterns::Double(0));
170 param_reader.
declare_entry(
"linear_global_tol",
"1.e-12",Patterns::Double(0));
172 param_reader.
declare_entry(
"line_maxiter",
"4",Patterns::Integer(0));
173 param_reader.
declare_entry(
"linesearch_rho",
"0.9",Patterns::Double(0));
174 param_reader.
declare_entry(
"linesearch_c",
"0.1",Patterns::Double(0));
176 param_reader.
declare_entry(
"compute_functionals_in_every_step",
"false",Patterns::Bool());
182 template <
typename PROBLEM,
typename VECTOR>
189 :
ReducedAlgorithm<PROBLEM, VECTOR>(OP,S,param_reader,Except,Output,base_priority)
192 param_reader.
SetSubsection(
"reducednewtonalgorithm parameters");
193 _nonlinear_maxiter = param_reader.
get_integer (
"nonlinear_maxiter");
194 _nonlinear_tol = param_reader.
get_double (
"nonlinear_tol");
195 _nonlinear_global_tol = param_reader.
get_double (
"nonlinear_global_tol");
197 _linear_maxiter = param_reader.
get_integer (
"linear_maxiter");
198 _linear_tol = param_reader.
get_double (
"linear_tol");
199 _linear_global_tol = param_reader.
get_double (
"linear_global_tol");
201 _line_maxiter = param_reader.
get_integer (
"line_maxiter");
202 _linesearch_rho = param_reader.
get_double (
"linesearch_rho");
203 _linesearch_c = param_reader.
get_double (
"linesearch_c");
205 _compute_functionals_in_every_step = param_reader.
get_bool (
"compute_functionals_in_every_step");
207 _postindex =
"_"+this->
GetProblem()->GetName();
212 template <
typename PROBLEM,
typename VECTOR>
220 template <
typename PROBLEM,
typename VECTOR>
228 this->GetReducedProblem()->ComputeReducedCostFunctional(q);
232 this->GetExceptionHandler()->HandleCriticalException(e,
"ReducedNewtonAlgorithm::NewtonResidual");
237 this->GetReducedProblem()->ComputeReducedGradient(q,gradient,gradient_transposed);
241 this->GetExceptionHandler()->HandleCriticalException(e,
"ReducedNewtonAlgorithm::NewtonResidual");
244 return sqrt(Residual(gradient,gradient_transposed));
249 template <
typename PROBLEM,
typename VECTOR>
259 std::stringstream out;
260 this->GetOutputHandler()->InitNewtonOut(out);
262 out <<
"**************************************************\n";
263 out <<
"* Starting Reduced Newton Algorithm *\n";
264 out <<
"* Solving : "<<this->GetProblem()->GetName()<<
"\t*\n";
268 this->GetReducedProblem()->StateSizeInfo(out);
269 out <<
"**************************************************";
270 this->GetOutputHandler()->Write(out,1+this->GetBasePriority(),1,1);
272 this->GetOutputHandler()->SetIterationNumber(iter,
"OptNewton"+_postindex);
274 this->GetOutputHandler()->Write(q,
"Control"+_postindex,
"control");
278 cost = this->GetReducedProblem()->ComputeReducedCostFunctional(q);
282 this->GetExceptionHandler()->HandleCriticalException(e,
"ReducedNewtonAlgorithm::Solve");
285 out<<
"CostFunctional: " << cost;
286 this->GetOutputHandler()->Write(out,2+this->GetBasePriority());
288 if (_compute_functionals_in_every_step ==
true)
292 this->GetReducedProblem()->ComputeReducedFunctionals(q);
296 this->GetExceptionHandler()->HandleCriticalException(e);
302 this->GetReducedProblem()->ComputeReducedGradient(q,gradient,gradient_transposed);
306 this->GetExceptionHandler()->HandleCriticalException(e,
"ReducedNewtonAlgorithm::Solve");
309 double res = Residual(gradient,gradient_transposed);
310 double firstres = res;
314 this->GetOutputHandler()->Write(gradient,
"NewtonResidual"+_postindex,
"control");
315 out<<
"\t Newton step: " <<iter<<
"\t Residual (abs.): "<<sqrt(res)<<
"\n";
316 out<<
"\t Newton step: " <<iter<<
"\t Residual (rel.): "<<std::scientific<<sqrt(res)/sqrt(res)<<
"\n";
317 this->GetOutputHandler()->Write(out,3+this->GetBasePriority());
320 unsigned int miniter = 0;
324 global_tol = std::max(_nonlinear_global_tol,global_tol);
325 while(( (res >= global_tol*global_tol) && (res >= _nonlinear_tol*_nonlinear_tol*firstres) ) || iter < miniter )
328 this->GetOutputHandler()->SetIterationNumber(iter,
"OptNewton"+_postindex);
330 if(iter > _nonlinear_maxiter)
338 liniter = SolveReducedLinearSystem(q,gradient,gradient_transposed, dq);
344 this->GetExceptionHandler()->HandleException(e,
"ReducedNewtonAlgorithm::Solve");
349 this->GetExceptionHandler()->HandleException(e,
"ReducedNewtonAlgorithm::Solve");
354 this->GetExceptionHandler()->HandleCriticalException(e,
"ReducedNewtonAlgorithm::Solve");
359 lineiter = ReducedNewtonLineSearch(dq,gradient,cost,q);
365 this->GetExceptionHandler()->HandleException(e,
"ReducedNewtonAlgorithm::Solve");
373 out<<
"CostFunctional: " << cost;
374 this->GetOutputHandler()->Write(out,3+this->GetBasePriority());
376 if (_compute_functionals_in_every_step ==
true)
380 this->GetReducedProblem()->ComputeReducedFunctionals(q);
384 this->GetExceptionHandler()->HandleCriticalException(e);
392 this->GetReducedProblem()->ComputeReducedGradient(q,gradient,gradient_transposed);
396 this->GetExceptionHandler()->HandleCriticalException(e,
"ReducedNewtonAlgorithm::Solve");
399 this->GetOutputHandler()->Write(q,
"Control"+_postindex,
"control");
400 this->GetOutputHandler()->Write(gradient,
"NewtonResidual"+_postindex,
"control");
402 res = Residual(gradient,gradient_transposed);
404 out<<
"\t Newton step: " <<iter<<
"\t Residual (rel.): "<<this->GetOutputHandler()->ZeroTolerance(sqrt(res)/sqrt(firstres),1.0)<<
"\t LinearIters ["<<liniter<<
"]\t LineSearch {"<<lineiter<<
"} ";
405 this->GetOutputHandler()->Write(out,3+this->GetBasePriority());
409 out<<
"CostFunctional: " << cost;
410 this->GetOutputHandler()->Write(out,2+this->GetBasePriority());
413 this->GetReducedProblem()->ComputeReducedFunctionals(q);
417 this->GetExceptionHandler()->HandleCriticalException(e,
"ReducedNewtonAlgorithm::Solve");
420 out <<
"**************************************************\n";
421 out <<
"* Stopping Reduced Newton Algorithm *\n";
422 out <<
"* after "<<std::setw(6)<<iter<<
" Iterations *\n";
424 out <<
"* with rel. Residual "<<std::scientific << std::setw(11) << this->GetOutputHandler()->ZeroTolerance(sqrt(res)/sqrt(firstres),1.0)<<
" *\n";
426 out <<
"**************************************************";
427 this->GetOutputHandler()->Write(out,1+this->GetBasePriority(),1,1);
433 template <
typename PROBLEM,
typename VECTOR>
439 std::stringstream out;
444 r_transposed = gradient_transposed;
445 d = gradient_transposed;
447 double res = Residual(r,r_transposed);
448 double firstres = res;
452 out <<
"Starting Reduced Linear Solver with Residual: "<<sqrt(res);
453 this->GetOutputHandler()->Write(out,4+this->GetBasePriority());
455 unsigned int iter = 0;
456 double cgalpha, cgbeta, oldres;
458 this->GetOutputHandler()->SetIterationNumber(iter,
"OptNewtonCg"+_postindex);
462 while(res>= std::min(0.25,sqrt(firstres))*firstres && res>=_linear_global_tol*_linear_global_tol)
465 this->GetOutputHandler()->SetIterationNumber(iter,
"OptNewtonCg"+_postindex);
466 if(iter > _linear_maxiter)
468 throw DOpEIterationException(
"Iteration count exceeded bounds!",
"ReducedNewtonAlgorithm::SolveReducedLinearSystem");
473 this->GetReducedProblem()->ComputeReducedHessianVector(q,d,Hd,Hd_transposed);
477 this->GetExceptionHandler()->HandleCriticalException(e,
"ReducedNewtonAlgorithm::SolveReducedLinearSystem");
480 cgalpha = res / (Hd*d);
493 r_transposed.add(cgalpha,Hd_transposed);
496 res = Residual(r,r_transposed);
503 this->GetReducedProblem()->ComputeReducedHessianVector(q,dq,Hd,Hd_transposed);
507 this->GetExceptionHandler()->HandleCriticalException(e);
510 r_transposed = gradient_transposed;
512 r_transposed.add(1.,Hd_transposed);
513 res = Residual(r,r_transposed);
516 out<<
"\t Cg step: " <<iter<<
"\t Residual: "<<sqrt(res);
517 this->GetOutputHandler()->Write(out,4+this->GetBasePriority());
519 cgbeta = res / oldres;
521 d.equ(-1,r_transposed);
529 template <
typename PROBLEM,
typename VECTOR>
535 double rho = _linesearch_rho;
536 double c = _linesearch_c;
539 bool force_linesearch =
false;
544 costnew = this->GetReducedProblem()->ComputeReducedCostFunctional(q);
549 force_linesearch =
true;
550 this->GetOutputHandler()->Write(
"Computing Cost Failed",4+this->GetBasePriority());
554 unsigned int iter =0;
556 double reduction = gradient*dq;
559 this->GetOutputHandler()->WriteError(
"Waring: computed direction doesn't seem to be a descend direction!");
563 if(_line_maxiter > 0)
565 if(fabs(reduction) < 1.e-10*cost)
567 if(std::isinf(costnew) || std::isnan(costnew) || (costnew >= cost + c*alpha*reduction) || force_linesearch)
569 this->GetOutputHandler()->Write(
"\t linesearch ",4+this->GetBasePriority());
570 while(std::isinf(costnew) || std::isnan(costnew) || (costnew >= cost + c*alpha*reduction) || force_linesearch)
573 if(iter > _line_maxiter)
577 throw DOpEException(
"Iteration count exceeded bounds while unable to compute the CostFunctional!",
"ReducedNewtonAlgorithm::ReducedNewtonLineSearch");
582 throw DOpEIterationException(
"Iteration count exceeded bounds!",
"ReducedNewtonAlgorithm::ReducedNewtonLineSearch");
585 force_linesearch =
false;
586 q.
add(alpha*(rho-1.),dq);
591 costnew = this->GetReducedProblem()->ComputeReducedCostFunctional(q);
596 force_linesearch =
true;
597 this->GetOutputHandler()->Write(
"Computing Cost Failed",4+this->GetBasePriority());
Definition: dopeexception.h:76
void PrintInfos(std::stringstream &out)
Definition: controlvector.cc:542
ReducedNewtonAlgorithm(PROBLEM *OP, ReducedProblemInterface< PROBLEM, VECTOR > *S, ParameterReader ¶m_reader, DOpEExceptionHandler< VECTOR > *Except=NULL, DOpEOutputHandler< VECTOR > *Output=NULL, int base_priority=0)
Definition: reducednewtonalgorithm.h:183
double get_double(const std::string &entry_name)
Definition: parameterreader.h:115
Definition: dopeexception.h:90
bool get_bool(const std::string &entry_name)
Definition: parameterreader.h:148
Definition: parameterreader.h:36
Definition: optproblemcontainer.h:70
~ReducedNewtonAlgorithm()
Definition: reducednewtonalgorithm.h:213
virtual double Residual(const ControlVector< VECTOR > &gradient, const ControlVector< VECTOR > &gradient_transposed)
Definition: reducednewtonalgorithm.h:143
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation=std::string())
Definition: parameterreader.h:98
Definition: controlvector.h:48
Definition: reducednewtonalgorithm.h:50
PROBLEM * GetProblem()
Definition: reducedalgorithm.h:187
virtual int Solve(ControlVector< VECTOR > &q, double global_tol=-1.)
Definition: reducednewtonalgorithm.h:250
virtual int SolveReducedLinearSystem(const ControlVector< VECTOR > &q, const ControlVector< VECTOR > &gradient, const ControlVector< VECTOR > &gradient_transposed, ControlVector< VECTOR > &dq)
Definition: reducednewtonalgorithm.h:434
Definition: reducedalgorithm.h:54
void add(double s, const ControlVector &dq)
Definition: controlvector.cc:363
void ReInit()
Definition: controlvector.cc:60
static void declare_params(ParameterReader ¶m_reader)
Definition: reducedalgorithm.h:241
virtual int ReducedNewtonLineSearch(const ControlVector< VECTOR > &dq, const ControlVector< VECTOR > &gradient, double &cost, ControlVector< VECTOR > &q)
Definition: reducednewtonalgorithm.h:530
int get_integer(const std::string &entry_name)
Definition: parameterreader.h:126
static void declare_params(ParameterReader ¶m_reader)
Definition: reducednewtonalgorithm.h:161
void SetSubsection(const std::string subsection)
Definition: parameterreader.h:93
Definition: reducedprobleminterface.h:338
Definition: dopeexception.h:35
double NewtonResidual(const ControlVector< VECTOR > &q)
Definition: reducednewtonalgorithm.h:221
Definition: optproblemcontainer.h:72