Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
solverunconstrained.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2008-02-14
7 
8  Copyright (C) 2008 Universite Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef SOLVERUNCONSTRAINED_H
30 #define SOLVERUNCONSTRAINED_H
31 
32 #include <boost/numeric/ublas/banded.hpp>
33 #include <boost/numeric/ublas/symmetric.hpp>
34 #include <boost/numeric/ublas/io.hpp>
35 #include <boost/lambda/lambda.hpp>
36 #include <boost/lambda/if.hpp>
37 
38 #include <clapack.h>
39 
40 #include <feel/feelcore/feel.hpp>
41 #include <feel/feelopt/problem.hpp>
43 
44 namespace Feel
45 {
51 template<
52 typename Data,
53  template<class> class Problem = problem
54  >
56 {
57 public:
58 
61 
63  typedef Problem<Data> problem_type;
64 
65 
66 
67  enum
68  {
69  _E_n = problem_type::_E_n,
70  _E_g = 0,
71  _E_h = 0,
72  _E_nA = problem_type::_E_nA,
73  _E_nL = problem_type::_E_nL,
74  _E_nAL = problem_type::_E_nAL
75  };
76 
78  typedef typename problem_type::ad_0_type ad_0_type;
79 
81  typedef typename problem_type::ad_1_type ad_1_type;
82 
84  typedef typename problem_type::ad_2_type ad_2_type;
85 
87  typedef typename ad_0_type::value_type value_type;
88 
90  typedef typename problem_type::f_type f_type;
91 
93  //typedef typename problem_type::value_type value_type;
94 
96  typedef typename problem_type::vector_type vector_type;
97 
99  typedef typename problem_type::matrix_type matrix_type;
100  typedef typename problem_type::symmetric_matrix_type symmetric_matrix_type;
101  typedef ublas::banded_matrix<double> banded_matrix_type;
102 
103  typedef directionalScalingMatrix<value_type> dsm_type;
108  struct COptions
109  {
110  COptions()
111  :
112  Delta_init( 1.0 ),
113  zeta_min( 0.9 ),
114  CGtol( 1e-12 ),
115  deps( 1e-9 ),
116  max_TR_iter( 100 ),
117  TR_tol( 1e-7 ),
118  allow_Trust_radius_calculations( true ),
119  rho_decrease( 0.3 ),
120  rho_big( 0.9 ),
121  rho_increase_big( 1.5 ),
122  rho_small( 0.3 ),
123  rho_increase_small( 1.1 )
124  {
125  // do nothing here
126  }
127  value_type Delta_init;
128  value_type zeta_min;
129  value_type CGtol;
130  value_type deps;
131  value_type max_TR_iter;
132  value_type TR_tol;
133  bool allow_Trust_radius_calculations;
134 
135  value_type rho_decrease;
136  value_type rho_big;
137  value_type rho_increase_big;
138  value_type rho_small;
139  value_type rho_increase_small;
140  };
141 
144 
149  class Stats
150  {
151  public:
152 
157  Stats( solver_type* __s, bool __c = false )
158  :
159  M_s( __s ),
160  M_collect( __c ),
161  M_x ( _E_n ),
162  M_l ( _E_n ),
163  M_u ( _E_n )
164  {
165  // nothing to do here
166  }
167 
169  bool collectStats() const
170  {
171  return M_collect;
172  }
173 
175  void clear();
176 
177  void push( const vector_type & __x, vector_type const&, vector_type const& );
178 
179  void push( value_type norm_Tgrad_fx,
180  int n_CGiter, int n_restarts, int n_indef,
181  int n_crosses_def, int n_crosses_indef,
182  int n_truss_exit_def, int n_truss_exit_indef,
183  value_type Delta, value_type ared_til, value_type phi_til, value_type rho );
184  void show() const;
185 
186  private:
187 
188  solver_type *M_s;
189  bool M_collect;
190  mutable int iter;
191 
192  vector_type M_x;
193  vector_type M_l;
194  vector_type M_u;
195 
196  std::vector<value_type> norm_Tgrad_fx_hstr, norm_err;
197  std::vector<int> n_CGiter_hstr, n_restarts_hstr, n_indef_hstr,
198  n_crosses_def_hstr, n_crosses_indef_hstr,
199  n_truss_exit_def_hstr, n_truss_exit_indef_hstr;
200  std::vector<value_type> Delta_hstr, ared_til_hstr, phi_til_hstr, rho_hstr;
201 
202 
203 
204  };
205 
208 
210  SolverUnconstrained( value_type x_definitions[_E_n][3] );
211 
213  SolverUnconstrained( Data const& __data );
214 
217 
219  void redefine_problem( value_type x_definitions[_E_n][3] )
220  {
221  M_prob.define_problem( x_definitions );
222  }
223 
225  bool optimize( vector_type& __x );
226 
227  statistics_type const& stats() const
228  {
229  return M_solver_stats;
230  }
231 
232  problem_type& problem()
233  {
234  return M_prob;
235  }
236 
237 private:
238 
239  // make it private to avoid having it called
241 
242 
244  void read_options();
245 
247  void makeCauchyStep( vector_type & _x, value_type _Delta,
248  f_type&,
249  vector_type & _Tgrad_fx,
250  banded_matrix_type & _Hg,
251  symmetric_matrix_type & _Thess_fxT, symmetric_matrix_type & _Htil,
252  vector_type & _neg_grad_fx );
253 
255  void makeStep( vector_type & _x, vector_type & _s, value_type _Delta,
256  f_type&,
257  vector_type & _Tgrad_fx,
258  banded_matrix_type & _Hg,
259  symmetric_matrix_type & _Thess_fxT, symmetric_matrix_type & _Htil );
260 
262  value_type norm_Theta_x_grad_fx( vector_type & _x, value_type _Delta );
263 
264 
266  void CGstep( vector_type & _x, value_type _Delta, vector_type & _sCG, value_type &norm_s_til,
267  int &_CGiter,
268  int &_n_restarts, int &_n_indef,
269  int &_n_crosses_def, int &_n_crosses_indef,
270  int &_n_truss_exit_def, int &_n_truss_exit_indef,
271  value_type &_s_til_x_G_til_x_s_til, value_type &phi_til );
272 
281  value_type tau( vector_type const& _s, vector_type const& _d, value_type const& _Delta );
282 
293  value_type xi( vector_type const& __s, vector_type const& __d, value_type const& __xi1 );
294 
295 
297  void lambda_LS( vector_type & _x, vector_type & _lambda_l, vector_type & _lambda_u );
298 
299 private:
300 
301  // problem specification and data
302  problem_type M_prob;
303 
304  options_type M_options;
305 
306  statistics_type M_solver_stats;
307 
308  dsm_type M_theta;
309 
310 };
311 
312 // SolverUnconstrained
313 
314 //#define DEBUG_OUT_CG
315 //#define DEBUG_OUT_TR
316 
317 template<typename Data,template<class> class Problem>
319  :
320  M_prob( x_definitions ),
321  M_solver_stats ( this, true ),
322  M_theta( _E_n )
323 {
324  read_options();
325 }
326 template<typename Data,template<class> class Problem>
328  :
329  M_prob( __data ),
330  M_solver_stats ( this, true ),
331  M_theta( M_prob.lowerBounds(), M_prob.upperBounds() )
332 {
333  read_options();
334 }
335 
336 template<typename Data,template<class> class Problem>
338 {
339 }
340 
341 template<typename Data,template<class> class Problem>
342 bool
344 {
345  M_solver_stats.clear();
346 
347  // Controlling parameters
348  // trust region radius
349  value_type Delta = M_options.Delta_init;
350 
351  vector_type x_new ( x.size() );
352  vector_type stot ( x.size() );
353  value_type norm_Tgrad_fx = this->norm_Theta_x_grad_fx( x, Delta );
354  value_type norm_s_til = 0;
355 
356  M_prob.copy_x0_to_x( x );
357 
358 
359  // FIXME: if ( M_options.verbose() )
360  //M_prob.print_complete( x );
361 
362  if ( M_solver_stats.collectStats() )
363  M_solver_stats.push( norm_Tgrad_fx, 0, 0, 0, 0, 0, 0, 0, Delta, 0, 0, 0 );
364 
365 
366  try
367  {
368  int iter = 0;
369 
370  //
371  // Apply bound constrained Trust region algorithm
372  //
373  while ( iter == 0 ||
374  ( iter < M_options.max_TR_iter && norm_Tgrad_fx > M_options.TR_tol ) )
375  {
376  iter++;
377  int n_CGiter, n_restarts, n_indef,
378  n_crosses_def, n_crosses_indef, n_truss_exit_def, n_truss_exit_indef;
379  value_type _s_til_x_G_til_x_s_til, phi_til, rho, Delta_used = Delta;
380 
381  DVLOG(2) << "\n===================== iter = " << iter << " ===========================";
382  //DVLOG(2) << "\nx = " << x << "\n";
383  DVLOG(2) << "\n -> norm_Tgrad_fx = " << norm_Tgrad_fx;
384 
385 
386 
391  this->CGstep( x, Delta, stot, norm_s_til,
392  n_CGiter,
393  n_restarts, n_indef,
394  n_crosses_def, n_crosses_indef,
395  n_truss_exit_def, n_truss_exit_indef,
396  _s_til_x_G_til_x_s_til, phi_til );
397 
398  //
399  x_new = x + stot;
400 
401  f_type __fx_new;
402  M_prob.evaluate( x_new, __fx_new, diff_order<0>() );
403  f_type __fx;
404  M_prob.evaluate( x, __fx, diff_order<0>() );
405 
406  // compute actual merit function reduction
407  value_type ared_til = __fx_new.value( 0 ) - __fx.value( 0 ) + 0.5 * _s_til_x_G_til_x_s_til;
408 
409  rho = ared_til / phi_til;
410 
412  // Trust region radius calculation:
413  if ( M_options.allow_Trust_radius_calculations )
414  {
415  if ( rho <= 1e-8 )
416  Delta = M_options.rho_decrease * Delta;
417 
418  else
419  {
420  x = x + stot;
421 
422  if ( rho > M_options.rho_big )
423  Delta = std::max( M_options.rho_increase_big*norm_s_til, Delta );
424 
425  else if ( rho > M_options.rho_small )
426  Delta = std::max( M_options.rho_increase_small*norm_s_til, Delta );
427  }
428 
429  norm_Tgrad_fx = this->norm_Theta_x_grad_fx( x, Delta );
430  }
431 
432  else
433  {
434  x = x + stot;
435  norm_Tgrad_fx = this->norm_Theta_x_grad_fx( x, Delta );
436  }
437 
439 
440  if ( M_solver_stats.collectStats() )
441  {
442  M_solver_stats.push( norm_Tgrad_fx,
443  n_CGiter, n_restarts, n_indef,
444  n_crosses_def, n_crosses_indef,
445  n_truss_exit_def, n_truss_exit_indef,
446  Delta_used, ared_til, phi_til, rho );
447  }
448 
449  if ( norm_Tgrad_fx > 1e-5 )
450  M_prob.setAccuracy( std::min( 1e-1, norm_Tgrad_fx ) );
451 
452  DVLOG(2) << "norm_Tgrad_fx = " << norm_Tgrad_fx << "\n";
453  }
454  }
455 
456  catch ( std::exception const& __ex )
457  {
458  f_type __fx_new;
459  M_prob.evaluate ( x, __fx_new, diff_order<2>() );
460 
461  if ( norm_inf( __fx_new.gradient( 0 ) ) > 1e-10 )
462  throw __ex;
463  }
464 
465  vector_type __l( _E_n ), __u( _E_n );
466  lambda_LS( x, __l, __u );
467 
468  if ( M_solver_stats.collectStats() )
469  M_solver_stats.push( x, __l, __u );
470 
471  return true;
472 }
473 #if 0
474 template<typename Data,template<class> class Problem> void
476 {
477  M_solver_stats.show();
478 
479  vector_type lambda_l ( _E_n ), lambda_u ( _E_n );
480  this->lambda_LS( x, lambda_l, lambda_u );
481  M_prob.print_stationary_x( x, lambda_l, lambda_u );
482 }
483 #endif
484 template<typename Data,template<class> class Problem>
485 void
486 SolverUnconstrained<Data,Problem>::read_options()
487 {
488  value_type num;
489  char str[20];
490  std::ifstream fin( "t_options.inp" );
491 
492  if ( fin.fail() )
493  {
494  // no option file : give up
495  return;
496  }
497 
498  while ( !fin.eof() )
499  {
500  fin >> str >> num;
501 
502  if ( !strcmp( str, "Delta_init" ) ) M_options.Delta_init = num;
503 
504  if ( !strcmp( str, "zeta_min" ) ) M_options.zeta_min = num;
505 
506  if ( !strcmp( str, "CGtol" ) ) M_options.CGtol = num;
507 
508  if ( !strcmp( str, "deps" ) ) M_options.deps = num;
509 
510  if ( !strcmp( str, "max_TR_iter" ) ) M_options.max_TR_iter = num;
511 
512  if ( !strcmp( str, "TR_tol" ) ) M_options.TR_tol = num;
513 
514  if ( !strcmp( str, "allow_Trust_radius_calculations" ) ) M_options.allow_Trust_radius_calculations = num;
515 
516  if ( !strcmp( str, "rho_decrease" ) ) M_options.rho_decrease = num;
517 
518  if ( !strcmp( str, "rho_big" ) ) M_options.rho_big = num;
519 
520  if ( !strcmp( str, "rho_increase_big" ) ) M_options.rho_increase_big = num;
521 
522  if ( !strcmp( str, "rho_small" ) ) M_options.rho_small = num;
523 
524  if ( !strcmp( str, "rho_increase_small" ) ) M_options.rho_increase_small = num;
525  }
526 
527  fin.close();
528 
529 }
530 
531 
532 template<typename Data,template<class> class Problem>
533 void
534 SolverUnconstrained<Data,Problem>::makeCauchyStep( vector_type & _x, value_type _Delta,
535  f_type& __fx,
536  vector_type & _Tgrad_fx,
537  banded_matrix_type & _Hg,
538  symmetric_matrix_type & _Thess_fxT, symmetric_matrix_type & _Htil,
539  vector_type & _neg_grad_fx )
540 {
541  M_prob.evaluate ( _x, __fx, diff_order<2>() );
542 
543 
544  _neg_grad_fx = - __fx.gradient( 0 );
545 
546  makeStep( _x, _neg_grad_fx, _Delta, __fx, _Tgrad_fx,
547  _Hg, _Thess_fxT, _Htil );
548 }
549 
550 template<typename Data,template<class> class Problem>
551 void
552 SolverUnconstrained<Data,Problem>::makeStep( vector_type & _x, vector_type & _s, value_type _Delta,
553  f_type& __fx,
554  vector_type & _Tgrad_fx,
555  banded_matrix_type & _Hg,
556  symmetric_matrix_type & _Thess_fxT, symmetric_matrix_type & _Htil )
557 {
558  M_theta.update( _Delta, _x, _s );
559 
560  _Tgrad_fx = prod( M_theta(), __fx.gradient( 0 ) );
561 
562  banded_matrix_type __diag_grad_fx( _E_n, _E_n, 0, 0 );
563 
564  matrix<value_type> __m( outer_prod( _Tgrad_fx,
565  scalar_vector<value_type>( _Tgrad_fx.size(), 1 ) ) );
566  __diag_grad_fx = banded_adaptor<matrix<value_type> >( __m, 0, 0 );
567 
568  _Hg = prod( M_theta.jacobian(), __diag_grad_fx );
569 
570  _Thess_fxT = prod( M_theta(), prod( __fx.hessian( 0 ), M_theta() ) );
571  _Htil = _Thess_fxT + _Hg;
572 }
573 
574 template<typename Data,template<class> class Problem>
575 typename SolverUnconstrained<Data,Problem>::value_type
576 SolverUnconstrained<Data,Problem>::norm_Theta_x_grad_fx( vector_type & _x, value_type _Delta )
577 {
578  f_type __fx;
579  M_prob.evaluate ( _x, __fx, diff_order<1>() );
580 
581  vector_type __Tgrad_fx( _E_n );
582  vector_type __neg_grad_fx( _E_n );
583  __neg_grad_fx = - __fx.gradient( 0 );
584 
585  M_theta.update( _Delta, _x, __neg_grad_fx, dsm_type::NO_JACOBIAN );
586 
587  __Tgrad_fx = prod ( M_theta(), __fx.gradient( 0 ) );
588 
589  return norm_2( __Tgrad_fx );
590 }
591 
592 template<typename Data,template<class> class Problem>
593 typename SolverUnconstrained<Data,Problem>::value_type
594 SolverUnconstrained<Data,Problem>::tau( vector_type const& _s,
595  vector_type const& _d,
596  value_type const& _Delta )
597 {
598  value_type a = inner_prod( _d, _d );
599  GST_SMART_ASSERT( a != 0 ).error( "a is 0 and should not be. It will cause a divide by zero." );
600  value_type b = 2.0 * inner_prod( _s, _d );
601  value_type c = inner_prod( _s, _s ) - _Delta*_Delta;
602  return ( -b + sqrt( b*b - 4.0*a*c ) ) / ( 2.0*a );
603 }
604 
605 template<typename Data,template<class> class Problem>
606 typename SolverUnconstrained<Data,Problem>::value_type
607 SolverUnconstrained<Data,Problem>::xi( vector_type const& __s,
608  vector_type const& __d,
609  value_type const& _xi1 )
610 {
611 
612  vector_type __frac = - element_div ( __s, __d );
613 
614  namespace lambda = boost::lambda;
615 
616  value_type __xi = _xi1;
617  std::for_each( __frac.begin(), __frac.end(),
618  lambda::if_then( lambda::_1 > 0 && lambda::_1 < lambda::var( __xi ),
619  lambda::var( __xi ) = lambda::_1 ) );
620  return __xi;
621 }
622 
623 template<typename Data,template<class> class Problem>
624 void
625 SolverUnconstrained<Data,Problem>::CGstep( vector_type & _x, value_type _Delta, vector_type & _sCG,
626  value_type &norm_s_til,
627  int &_CGiter,
628  int &_n_restarts, int &_n_indef,
629  int &_n_crosses_def, int &_n_crosses_indef,
630  int &_n_truss_exit_def,
631  int &_n_truss_exit_indef,
632  value_type &_s_til_x_G_til_x_s_til,
633  value_type &phi_til )
634 {
635  bool _restart = false;
636  bool _done = false;
637  _n_restarts = 0;
638  _CGiter = 0;
639  _n_indef = 0;
640  _n_crosses_def = 0;
641  _n_crosses_indef = 0;
642  _n_truss_exit_def = 0;
643  _n_truss_exit_indef = 0;
644 
645  f_type __fx;
646 
647  vector_type _neg_grad_fx ( _x.size() );
648  vector_type _Tgrad_fx ( _x.size() );
649  vector_type _Htil_d ( _x.size() ); // Htil * dCG
650  vector_type _Htil_s_til ( _x.size() ); // Htil * s_til
651 
652  vector_type _s_til ( _x.size() );
653  vector_type _s_til_old ( _x.size() );
654  vector_type _s_til_eps ( _x.size() );
655  vector_type _rCG ( _x.size() );
656  vector_type _rCG_old ( _x.size() );
657  vector_type _dCG ( _x.size() );
658 
659  value_type _alpha = 0;
660  value_type _gamma; // d' Htil d
661  value_type _beta;
662  value_type _xi = 0;
663  value_type _tau;
664 
665  // diagonal matrices
666  banded_matrix_type _Hg ( _x.size(), _x.size(), 0, 0 ); // Gtil
667 
668  symmetric_matrix_type _Thess_fxT ( _E_nA, _E_nA ); // Btil
669  symmetric_matrix_type _Htil ( _E_nA, _E_nA ); // Htil
670 
671  DVLOG(2) << "\n\n[value_type SolverUnconstrained<Data,Problem>::CGstep]...\n";
672 
673  // INITIALIZE:
674  this->makeCauchyStep( _x, _Delta, __fx, _Tgrad_fx, _Hg, _Thess_fxT, _Htil, _neg_grad_fx );
675 
676  DVLOG(2) << "Trust region active (C) : " << M_theta.isTrustRegionActive() << "\n";
677 
678  value_type _norm_gtil = norm_2( _Tgrad_fx );
679  _s_til = zero_vector<value_type>( _s_til.size() );
680  _s_til_old = _s_til;
681  _rCG = -_Tgrad_fx;
682  _rCG_old = _rCG;
683  _dCG = _rCG;
684 
685  size_t _inner_iter = 0;
686 
687  _CGiter++;
688 
689  while ( !_done )
690  {
691  if ( _restart )
692  {
693  _n_restarts++;
694  _CGiter++;
695  // RE-INITIALIZE:
696  _inner_iter = 0;
697 
698  _s_til_eps = _s_til + M_options.deps * _dCG;
699 
700  this->makeStep( _x, _s_til_eps, _Delta,
701  __fx,
702  _Tgrad_fx, _Hg, _Thess_fxT, _Htil );
703 
704  _norm_gtil = norm_2( _Tgrad_fx );
705 
706  _Htil_s_til = prod( _Htil, _s_til );
707 
708  _s_til_old = _s_til;
709  _rCG = -_Tgrad_fx - _Htil_s_til;
710  _rCG_old = _rCG ;
711  _dCG = _rCG;
712 
713  _restart = false;
714  }
715 
716  //
717  // STEP 1:
718  //
719  _gamma = inner_prod( _dCG, prod( _Htil,_dCG ) );
720 
721  if ( _gamma <= 0 )
722  {
723  _n_indef++;
724  _tau = this->tau( _s_til, _dCG, _Delta );
725 
726  if ( ( !M_theta.isTrustRegionActive() ) || ( _CGiter == 1 ) )
727  {
728  _xi = _tau;
729  }
730 
731  else if ( _inner_iter == 0 )
732  {
733  _xi = this->xi( _s_til_eps, _dCG, _tau );
734  }
735 
736  else
737  {
738  _xi = this->xi( _s_til, _dCG, _tau );
739  }
740 
741  _s_til_old = _s_til;
742  _s_til += _xi * _dCG;
743 
744  if ( _xi < _tau )
745  {
746  _n_crosses_indef++;
747  _restart = true;
748  }
749 
750  else if ( M_theta.isTrustRegionActive() )
751  {
752  _sCG = _s_til;
753 
754  _n_truss_exit_indef++;
755  _done = true;
756  }
757  }
758 
759  //
760  // STEP 2:
761  //
762  if ( !_restart && !_done )
763  {
764  _alpha = inner_prod( _rCG, _rCG ) / _gamma;
765 
766  if ( ( !M_theta.isTrustRegionActive() ) || ( _CGiter == 1 ) )
767  _xi = _alpha;
768 
769  else if ( _inner_iter == 0 )
770  _xi = this->xi( _s_til_eps, _dCG, _alpha );
771 
772  else
773  _xi = this->xi( _s_til, _dCG, _alpha );
774 
775  _s_til_old = _s_til;
776  _s_til += _xi * _dCG;
777 
778  if ( norm_2( _s_til ) > _Delta )
779  {
780  _tau = this->tau( _s_til_old, _dCG, _Delta );
781 
782  _sCG = _s_til_old + _tau * _dCG;
783 
784  _n_truss_exit_def++;
785  _done = true;
786  }
787  }
788 
789  //
790  // STEP 3:
791  //
792  if ( !_restart && !_done )
793  {
794  if ( _xi >= _alpha )
795  {
796  _rCG_old = _rCG;
797  _rCG -= _alpha * prod( _Htil, _dCG );
798  }
799 
800  else if ( M_theta.isTrustRegionActive() )
801  {
802  _n_crosses_def++;
803  _restart = true;
804  }
805  }
806 
807  if ( norm_2( _rCG ) / _norm_gtil < M_options.CGtol )
808  {
809  _sCG = _s_til;
810 
811  DVLOG(2) << "\n\nNormal CG exit 1: ||rCG||/||g_til|| = " << norm_2( _rCG ) / _norm_gtil << "\n";
812 
813  _done = true;
814  }
815 
816  if ( _CGiter >= _x.size() )
817  {
818  _sCG = _s_til;
819 
820  DVLOG(2) << "\n\nNormal CG exit 2: _CGiter = " << _CGiter << "\n";
821 
822  _done = true;
823  }
824 
825  //
826  // STEP 4:
827  //
828  if ( !_restart && !_done )
829  {
830  _beta = inner_prod( _rCG, _rCG ) / inner_prod( _rCG_old, _rCG_old );
831 
832  DVLOG(2) << "\nbeta = " << _beta << "\n";
833 
834  _dCG = _rCG + _beta * _dCG;
835 
836  _CGiter++;
837  _inner_iter++;
838  }
839  }
840 
841  _s_til_x_G_til_x_s_til = inner_prod( _sCG, prod( _Hg, _sCG ) );
842 
843  phi_til = inner_prod( _Tgrad_fx, _sCG ) + 0.5 * inner_prod( _sCG, prod( _Htil,_sCG ) );
844 
845  norm_s_til = norm_2( _sCG );
846 
847  // restoring to original space
848  _sCG = prod( M_theta(), _sCG );
849 }
850 
851 
852 template<typename Data, template<class> class Problem>
853 void
854 SolverUnconstrained<Data,Problem>::lambda_LS( vector_type & _x, vector_type & _lambda_l, vector_type & _lambda_u )
855 {
856  //symmetric_matrix_type __A ( _E_nL, _E_nL );
857  matrix_type __A ( _E_nL, _E_nL );
858  vector_type __b ( _E_nL );
859 
860  // At = [ -I (l-x) 0; I 0 (x-u) ] --> AtA = [ I+(l-x)^2 -I; -I I+(x-u)^2 ]
861  __A = zero_matrix<double>( _E_nL, _E_nL );
862 
863  for ( int __i = 0; __i < _E_n; ++__i )
864  {
865  __A( __i, __i ) = 1.0 + ( M_prob.x_l( __i )-_x ( __i ) )*( M_prob.x_l( __i )-_x ( __i ) );
866 
867  __A( __i, _E_n+__i ) = -1.0;
868  __A( _E_n+__i,__i ) = -1.0;
869 
870  __A( ( _E_n+__i ), _E_n+__i ) = 1.0+( _x ( __i )-M_prob.x_u( __i ) )*( _x ( __i )-M_prob.x_u( __i ) );
871  }
872 
873  // b = [ \nabla f; 0; 0 ] --> Atb = [ -\nabla f; \nabla f ]
874  f_type __f_x;
875  M_prob.evaluate( _x, __f_x, diff_order<1>() );
876 
877  for ( int i = 0; i < _E_n; i++ )
878  {
879  value_type val = __f_x.gradient( 0, i );
880  __b ( i ) = val;
881  __b ( _E_n+i ) = -val;
882  }
883 
884  // after solve __b = [ _lambda_l; _lambda_] = A\b
885  //char __uplo = 'U';
886  int __info = 0;
887  int __nrhs = 1;
888  int __N = _E_nL;
889  //dppsv_ ( &__uplo, &__N, &__nrhs, __A.data(), __b.data(), &__ldb, &__info );
890  //dppsv_ ( &__uplo, &__N, &__nrhs, __A.data().begin(), __b.data().begin(), &__N, &__info );
891 
892  ublas::vector<int> __itype( __N );
893  //dspsv_( &__uplo, &__N, &__nrhs, __A.data().begin(), __itype.data().begin(), __b.data().begin(), &__N, &__info );
894  dgesv_( &__N,&__nrhs,__A.data().begin(),&__N,__itype.data().begin(),__b.data().begin(),&__N,&__info );
895 
896  if ( __info!= 0 )
897  {
898  std::ostringstream __err;
899 
900  __err << "[" << __PRETTY_FUNCTION__ << "] dppsv failed: " << __info << "\n"
901  << "A = " << __A << "\n"
902  << "B = " << __b << "\n";
903 
904  throw std::out_of_range( __err.str() );
905  }
906 
907  for ( int i = 0; i < _E_n; i++ )
908  {
909  _lambda_l ( i ) = __b ( i );
910  _lambda_u ( i ) = __b ( _E_n+i );
911  }
912 }
913 
914 
915 //
916 // Statistics class implementation
917 //
918 
919 template<typename Data, template<class> class Problem>
920 void
922 {
923  norm_Tgrad_fx_hstr.clear();
924  norm_err.clear();
925  n_CGiter_hstr.clear();
926  n_restarts_hstr.clear();
927  n_indef_hstr.clear();
928  n_crosses_def_hstr.clear();
929  n_crosses_indef_hstr.clear();
930  n_truss_exit_def_hstr.clear();
931  n_truss_exit_indef_hstr.clear();
932  Delta_hstr.clear();
933  ared_til_hstr.clear();
934  phi_til_hstr.clear();
935  rho_hstr.clear();
936 }
937 template<typename Data, template<class> class Problem>
938 void
940 {
941  M_x = __x;
942  M_l = __l;
943  M_u = __u;
944 }
945 
946 template<typename Data, template<class> class Problem>
947 void
949  int n_CGiter, int n_restarts, int n_indef,
950  int n_crosses_def, int n_crosses_indef,
951  int n_truss_exit_def, int n_truss_exit_indef,
952  value_type Delta,
953  value_type ared_til,
954  value_type phi_til,
955  value_type rho )
956 {
957  //push_x_in_hstr( _x );
958  norm_Tgrad_fx_hstr.push_back( norm_Tgrad_fx );
959  n_CGiter_hstr.push_back( n_CGiter );
960  n_restarts_hstr.push_back( n_restarts );
961  n_indef_hstr.push_back( n_indef );
962  n_crosses_def_hstr.push_back( n_crosses_def );
963  n_crosses_indef_hstr.push_back( n_crosses_indef );
964  n_truss_exit_def_hstr.push_back( n_truss_exit_def );
965  n_truss_exit_indef_hstr.push_back( n_truss_exit_indef );
966  Delta_hstr.push_back( Delta );
967  ared_til_hstr.push_back( ared_til );
968  phi_til_hstr.push_back( phi_til );
969  rho_hstr.push_back( rho );
970 }
971 
972 template<typename Data, template<class> class Problem>
973 void
974 SolverUnconstrained<Data,Problem>::Stats::show() const
975 {
976  M_s->problem().print_stationary_x ( M_x, M_l, M_u );
977 
978  iter = norm_Tgrad_fx_hstr.size();
979  std::cerr << "\n\nScaled Trust-Region Method Statistics:\n";
980  std::cerr << "\nNumber of outter iterations: " << iter-1 << "\n";
981 
982  std::cerr << "\n";
983  std::cerr << " k CGiters restarts indefs crosses trustEx ||g_til|| Delta^k ared_til^k phi_til^k rho^k \n";
984  std::cerr << "--------------------------------------------------------------------------------------------------------\n";
985 
986  for ( int k = 0; k < iter; k++ )
987  {
988  fprintf( stderr, "%3d ", k );
989 
990  if ( k == 0 )
991  std::cerr << " N/A N/A N/A N/A N/A";
992 
993  else
994  {
995  std::cerr << " " << n_CGiter_hstr[k]
996  << " " << n_restarts_hstr[k]
997  << " " << n_indef_hstr[k]
998  << " " << n_crosses_def_hstr[k] + n_crosses_indef_hstr[k]
999  << " " << n_truss_exit_def_hstr[k] + n_truss_exit_indef_hstr[k] << " ";
1000  }
1001 
1002  fprintf ( stderr, " %13.3e %9f", norm_Tgrad_fx_hstr[ k ], Delta_hstr[k] );
1003 
1004  if ( k == 0 )
1005  std::cerr << " N/A N/A N/A";
1006 
1007  else
1008  {
1009  fprintf ( stderr, " %11.3e", ared_til_hstr[ k ] );
1010  fprintf ( stderr, " %11.3e", phi_til_hstr[ k ] );
1011  fprintf ( stderr, " %11f", rho_hstr[ k ] );
1012  }
1013 
1014  std::cerr << std::endl;
1015  }
1016 
1017  //M_s->lambda_LS(
1018 }
1019 
1020 } // Feel
1021 #endif

Generated on Sun Oct 20 2013 08:25:05 for Feel++ by doxygen 1.8.4