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
solverconstrained.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 CONSTRAINED_SOLVER_HPP
30 #define CONSTRAINED_SOLVER_HPP
31 
32 #include <feel/feelopt/problem.hpp>
33 
34 namespace Feel
35 {
41 template<
42 typename Data,
43  template<class> class Problem = problem
44  >
46 {
47 public:
48 
51 
53  typedef Problem<Data> problem_type;
54 
55  enum
56  {
57  _E_n = problem_type::_E_n,
58  _E_f = problem_type::_E_f,
59  _E_g = problem_type::_E_g,
60  _E_h = problem_type::_E_h,
61  _E_nA = problem_type::_E_nA,
62  _E_nL = problem_type::_E_nL,
63  _E_nAL = problem_type::_E_nAL
64  };
65 
67  typedef typename problem_type::ad_0_type ad_0_type;
68 
70  typedef typename problem_type::ad_1_type ad_1_type;
71 
73  typedef typename problem_type::ad_2_type ad_2_type;
74 
76  typedef typename ad_0_type::value_type value_type;
77 
78  typedef typename problem_type::f_type f_type;
79  typedef typename problem_type::g_type g_type;
80  typedef typename problem_type::h_type h_type;
81 
87  struct COptions
88  {
89  COptions()
90  :
91  eta ( 1e-8 ),
92  tau ( 0.995 ),
93  theta_N ( 0.2 ),
94  zeta_N ( 0.8 ),
95  eTOL ( 1e-6 ),
96  rho_N ( 0.3 ),
97 
98  Delta_init ( 1 ),
99  s_theta_init ( 10 ),
100  eMU_init ( 1e-1 ),
101  MU_init ( 1e-1 ),
102  nu0_init ( 10 ),
103 
104  rho_decrease ( 0.1 ),
105  rho_big ( 0.9 ),
106  rho_increase_big ( 5.0 ),
107  rho_small ( 0.3 ),
108  rho_increase_small ( 2.0 ),
109 
110  Primal_Dual ( true ),
111 
112  max_TR_iter ( 4 ),
113  max_hom_iter ( 8 )
114 
115  {
116  // nothing to do here
117  }
118  double eta;
119  double tau;
120  double theta_N;
121  double zeta_N;
122  double eTOL;
123  double rho_N;
124 
125  double Delta_init;
126  double s_theta_init;
127  double eMU_init;
128  double MU_init;
129  double nu0_init;
130 
131 
132 
133  double rho_decrease;
134  double rho_big;
135  double rho_increase_big;
136  double rho_small;
137  double rho_increase_small;
138 
139  bool Primal_Dual;
140 
141  int max_TR_iter;
142  int max_hom_iter;
143 
144  };
145 
148 
153  struct Stats
154  {
159  Stats( solver_type* __s, bool __c = false )
160  :
161  M_s( __s ),
162  M_collect( __c )
163  {
164  // nothing to do here
165  }
166 
168  bool collectStats() const
169  {
170  return M_collect;
171  }
172 
174  void clear();
175 
177  void push_x( const double *__x );
178 
180  void push_all( double mu, double E, double E0,
181  int n_CGiter,
182  bool truss_exit, bool SOC,
183  double Delta, double ared, double pred, double gamma );
184 
189  void show( bool verbose = false ) const;
190 
191  private:
193  Stats();
194 
195  private:
196 
197 
198  solver_type* M_s;
199 
200  bool M_collect;
201 
202  mutable int iter;
203  std::vector<double> x_hstr, mu_hstr, E_hstr, E0_hstr;
204  std::vector<int> n_CGiter_hstr;
205  std::vector<bool> truss_exit_hstr, SOC_hstr;
206  std::vector<double> Delta_hstr, ared_hstr, pred_hstr, gamma_hstr;
207  };
208 
211 
213  /*
214  \brief default constructor
215  takes the bounds for the control variables
216  \todo change this in the future by adding a class for control variables
217  */
218  SolverConstrained( double x_definitions[_E_n][3] );
219 
222 
224  void redefine_problem( double x_definitions[_E_n][3] )
225  {
226  M_prob.define_problem( x_definitions );
227  }
228 
233  bool optimize( double __x[_E_n] );
234 
237  {
238  return M_prob;
239  }
240 
242  statistics_type const& stats() const
243  {
244  return M_solver_stats;
245  }
246 
247 private:
248 
250  void read_options();
251 
252  void initialize_solver_with_x( double *_x, double *_s, double *_lambda_h, double *_lambda_g, double MU );
253 
254  void update_opt_objs( double *_x, double *_s, double *_lambda_h,
255  double *_lambda_g, double MU, bool update_Hess );
256 
257  void solve_AtA_system( double _r[ _E_n + _E_g+2*_E_n ], double _b[ _E_h + _E_g+2*_E_n ],
258  double _g[ _E_n + _E_g+2*_E_n ], double _l[ _E_h + _E_g+2*_E_n ] );
259 
260  void dogleg_solve( double *_vCP_til, double *_vN_til, double *_s, double Delta,
261  double *_v_til );
262 
263  double model_m( double *_s, double *_v_til );
264 
265  double get_vpred( double *_s, double *_v_til );
266 
268  double get_normal_steps( double *_s, double Delta,
269  double *_dx_normal, double *_ds_normal, double *_v_til );
270 
271  void get_Pr( double *_r, double *_Pr );
272 
273  bool slack_feasible( double *_w_til, double *_v_til );
274 
275  double get_q( double *_v_til, double *_w_til, double MU );
276 
278  double get_tangen_steps( double *_s, double Delta, double *_v_til, double MU,
279  double *_dx_tangen, double *_ds_tangen,
280  bool &truss_exit, int &n_CGiter );
281 
282  double get_steps( double *_s, double _Delta, double MU,
283  double *_dx, double *_ds, double &_vpred, double &_q,
284  bool &truss_exit, int &n_CGiter );
285 
286  double get_phi( double *_x, double *_s, double MU, double nu, double *rhs );
287  double get_phi( double *_x, double *_s, double MU, double nu );
288  double get_phi( double *_s, double MU, double nu );
289 
290  double get_SOC( double *_s, double *_rhs_xs_p_d, double *_xy, double *_sy );
291 
292  //
293  // Evaluations
294  //
295 
297  void evaluate_Ag( double *x );
299  void evaluate_Ah( double *x );
301  void evaluate_A_hat( double *_x, double *_s );
303  void evaluate_Hessian_L( double *_x, double *_lambda_h, double *_lambda_g );
305  void evaluate_G( double *_x, double *_s, double *_lambda_h, double *_lambda_g,
306  double MU, bool evaluate__Hessians_flag );
307  double evaluate_E_current( double *_s, double *_lambda_h, double *_lambda_g, double MU );
309  void evaluate_gradf( double *_x );
311  void evaluate_fxgxhx( double *_x );
312 
313  //
314  // Print
315  //
316 
318  void print_Ag();
320  void print_Ah();
322  void print_lambda_LS( double *_lambda_g, double *_lambda_h );
324  void print_G();
326  void print_A_hat();
328  void print_Hessian_L();
329 
332 
333 private:
334 
335  problem_type M_prob;
336 
338  options_type options;
339 
341  statistics_type M_solver_stats;
342 
343  double fx, gradf[ _E_n ];
344  double gx[ _E_g+2*_E_n ];
345  double hx[ _E_h ];
346  double Ag[ _E_n ][ _E_g+2*_E_n ];
347  double Ah[ _E_n ][ _E_h ];
348  double A_hat[ _E_n + _E_g+2*_E_n ][ _E_h + _E_g+2*_E_n ];
349 
350  double Hessian_L[ _E_n ][ _E_n ];
351  double G[ _E_n+_E_g+2*_E_n ][ _E_n+_E_g+2*_E_n ];
352 
353  static const double Max_Allowed_Delta = 1e8;
354 };
355 
356 
357 // SolverConstrained
358 
359 //#define DEBUG_OUT_CG
360 //#define DEBUG_OUT_TR
361 //#define DEBUG_OUT_HM
362 //#define DEBUG_x_OUT
363 
364 template<typename Data,template<class> class Problem>
366  :
367  M_prob( x_definitions ),
368  M_solver_stats( this, true ),
369  fx( 0 ),
370  gradf(),
371  gx(),
372  hx(),
373  Ag(),
374  Ah(),
375  A_hat(),
376  Hessian_L(),
377  G()
378 {
379  std::cerr << "SolverConstrained::SolverConstrained()\n";
380  read_options();
381 }
382 
383 template<typename Data,template<class> class Problem>
385 {
386 }
387 
388 template<typename Data,template<class> class Problem>
390 {
391  double num;
392  char str[20];
393  std::ifstream fin( "N_options.inp" );
394 
395  if ( fin.fail() )
396  return;
397 
398  options.Primal_Dual = true;
399 
400  while ( !fin.eof() )
401  {
402  fin >> str >> num;
403 
404  if ( !strcmp( str, "Delta_init" ) ) options.Delta_init = num;
405 
406  if ( !strcmp( str, "eta" ) ) options.eta = num;
407 
408  if ( !strcmp( str, "tau" ) ) options.tau = num;
409 
410  if ( !strcmp( str, "theta_N" ) ) options.theta_N = num;
411 
412  if ( !strcmp( str, "zeta_N" ) ) options.zeta_N = num;
413 
414  if ( !strcmp( str, "eTOL" ) ) options.eTOL = num;
415 
416  if ( !strcmp( str, "rho_N" ) ) options.rho_N = num;
417 
418  if ( !strcmp( str, "s_theta_init" ) ) options.s_theta_init = num;
419 
420  if ( !strcmp( str, "eMU_init" ) ) options.eMU_init = num;
421 
422  if ( !strcmp( str, "MU_init" ) ) options.MU_init = num;
423 
424  if ( !strcmp( str, "nu0_init" ) ) options.nu0_init = num;
425 
426  if ( !strcmp( str, "Primal_Dual" ) ) options.Primal_Dual = ( bool )num;
427 
428  if ( !strcmp( str, "max_TR_iter" ) ) options.max_TR_iter = ( int )num;
429 
430  if ( !strcmp( str, "max_hom_iter" ) ) options.max_hom_iter = ( int )num;
431 
432  if ( !strcmp( str, "rho_decrease" ) ) options.rho_decrease = num;
433 
434  if ( !strcmp( str, "rho_big" ) ) options.rho_big = num;
435 
436  if ( !strcmp( str, "rho_increase_big" ) ) options.rho_increase_big = num;
437 
438  if ( !strcmp( str, "rho_small" ) ) options.rho_small = num;
439 
440  if ( !strcmp( str, "rho_increase_small" ) ) options.rho_increase_small = num;
441  }
442 
443  fin.close();
444 
445 }
446 
447 template<typename Data,template<class> class Problem>
448 void SolverConstrained<Data,Problem>::print_Hessian_L()
449 {
450  int n = M_prob.n();
451 
452  for ( int i = 0; i < n; i++ )
453  {
454  printf( "\n" );
455 
456  for ( int j = 0; j < n; j++ )
457  fprintf( stderr, "%10.3f", Hessian_L[i][j] );
458  }
459 
460  std::cerr << "\n";
461 }
462 
463 template<typename Data,template<class> class Problem>
464 void SolverConstrained<Data,Problem>::evaluate_Hessian_L( double *_x, double *_lambda_h, double *_lambda_g )
465 {
466  f_type __fx;
467  g_type __gx;
468  h_type __hx;
469  M_prob.evaluate ( _x, __fx, diff_order<2>() );
470  M_prob.evaluate ( _x, __gx, diff_order<2>() );
471  M_prob.evaluate ( _x, __hx, diff_order<2>() );
472 
473  for ( int i = 0; i < _E_n; i++ )
474  {
475  // symetry
476  for ( int j = i; j < _E_n; j++ )
477  {
478  double val = 0;
479 
480  for ( int __k = 0; __k < _E_f; ++__k )
481  {
482  val += __fx.hessian( __k, i, j );
483  }
484 
485  for ( int m = 0; m < _E_h; m++ )
486  {
487  val += _lambda_h[ m ] * __hx.hessian( m,i,j );
488  }
489 
490  for ( int m = 0; m < _E_g; m++ )
491  {
492  val += _lambda_g[ m ] * __gx.hessian( m,i,j );
493  }
494 
495  Hessian_L[ i ][ j ] = val;
496  Hessian_L[ j ][ i ] = val;
497  }
498  }
499 }
500 
501 template<typename Data,template<class> class Problem>
502 void SolverConstrained<Data,Problem>::print_G()
503 {
504  int n = M_prob.n(), nv = n+_E_g+2*n;
505 
506  for ( int i = 0; i < nv; i++ )
507  {
508  printf( "\n" );
509 
510  for ( int j = 0; j < nv; j++ )
511  fprintf( stderr, "%12.5f", G[i][j] );
512  }
513 
514  std::cerr << "\n";
515 }
516 
517 template<typename Data,template<class> class Problem>
518 void SolverConstrained<Data,Problem>::evaluate_G( double *_x, double *_s, double *_lambda_h,
519  double *_lambda_g, double MU, bool evaluate_Hessians_flag )
520 {
521  int n = M_prob.n();
522  int nv = n+_E_g+2*n;
523 
524  if ( evaluate_Hessians_flag )
525  this->evaluate_Hessian_L( _x, _lambda_h, _lambda_g );
526 
527  for ( int i = 0; i < nv; i++ )
528  for ( int j = 0; j < nv; j++ )
529  G[ i ][ j ] = 0;
530 
531  for ( int i = 0; i < n; i++ )
532  for ( int j = 0; j < n; j++ )
533  G[ i ][ j ] = Hessian_L[ i ][ j ];
534 
535  for ( int i = 0; i < _E_g+2*n; i++ )
536  if ( options.Primal_Dual && _lambda_g[ i ] > 0 )
537  G[ n+i ][ n+i ] = _lambda_g[ i ] * _s[ i ];
538 
539  else
540  G[ n+i ][ n+i ] = MU;
541 
542  //this->print_G(); exit(0);
543 
544  /*
545  std::cerr << "\nHL = ";
546  this->print_Hessian_L();
547  std::cerr << "\nlambda_g = ";
548  __print_vec( _E_g+2*n, _lambda_g );
549  std::cerr << "\ns = ";
550  __print_vec( _E_g+2*n, _s );
551  std::cerr << "\nG = ";
552  this->print_G();
553  exit(0);
554  */
555 }
556 
557 
558 template<typename Data,template<class> class Problem>
559 void SolverConstrained<Data,Problem>::evaluate_Ag( double *_x )
560 {
561  g_type __gx;
562  M_prob.evaluate ( _x, __gx, diff_order<1>() );
563 
564  int n = M_prob.n();
565 
566  assert ( n == _E_n );
567 
568  // Ag = [ g_0 ... g_n -I I ]
569  for ( int i = 0; i < _E_n; i++ )
570  for ( int m = 0; m < _E_g+2*_E_n; m++ )
571  Ag[ i ][ m ] = 0;
572 
573  for ( int i = 0; i < _E_n; i++ )
574  {
575  for ( int m = 0; m < _E_g; m++ )
576  {
577  Ag[ i ][ m ] = __gx.gradient( m,i );
578  }
579 
580  int __index = _E_g+i;
581  Ag[ i ][ __index ] = -1.0;
582  Ag[ i ][ __index + n] = +1.0;
583  }
584 }
585 
586 template<typename Data,template<class> class Problem>
587 void SolverConstrained<Data,Problem>::evaluate_Ah( double *_x )
588 {
589  h_type __hx;
590  M_prob.evaluate ( _x, __hx, diff_order<1>() );
591 
592  int n = M_prob.n();
593 
594  // Ah = [ h_0 ... h_n ]
595  for ( int i = 0; i < n; i++ )
596  for ( int m = 0; m < _E_h; m++ )
597  Ah[ i ][ m ] = 0;
598 
599  for ( int i = 0; i < n; i++ )
600  for ( int m = 0; m < _E_h; m++ )
601  Ah[ i ][ m ] = __hx.gradient( m,i );
602 }
603 
604 
605 template<typename Data,template<class> class Problem>
606 void SolverConstrained<Data,Problem>::print_Ag()
607 {
608  int n = M_prob.n();
609 
610  for ( int i = 0; i < n; i++ )
611  {
612  printf( "\n" );
613 
614  for ( int j = 0; j < _E_g+2*n; j++ )
615  fprintf( stderr, "%10.1f", Ag[i][j] );
616  }
617 
618  std::cerr << "\n";
619 }
620 
621 template<typename Data,template<class> class Problem>
622 void SolverConstrained<Data,Problem>::print_Ah()
623 {
624  int n = M_prob.n();
625 
626  for ( int i = 0; i < n; i++ )
627  {
628  printf( "\n" );
629 
630  for ( int j = 0; j < _E_h; j++ )
631  fprintf( stderr, "%10.1f", Ah[i][j] );
632  }
633 
634  std::cerr << "\n";
635 }
636 
637 template<typename Data,template<class> class Problem>
638 void SolverConstrained<Data,Problem>::evaluate_fxgxhx( double *_x )
639 {
640  f_type __fx;
641  g_type __gx;
642  h_type __hx;
643 
644  // get the order 0 information only
645  M_prob.evaluate ( _x, __fx, diff_order<0>() );
646  M_prob.evaluate ( _x, __gx, diff_order<0>() );
647  M_prob.evaluate ( _x, __hx, diff_order<0>() );
648 
649  fx = __fx.value( 0 );
650 
651  for ( uint __i = 0; __i < __gx.value().size(); ++__i )
652  gx[__i] = __gx.value( __i );
653 
654  for ( uint __i = 0; __i < _E_h; ++__i )
655  hx[__i] = __hx.value( __i );
656 }
657 
658 
659 template<typename Data,template<class> class Problem>
660 void
661 SolverConstrained<Data,Problem>::initialize_solver_with_x( double *_x,
662  double *_s,
663  double *_lambda_h, double *_lambda_g,
664  double MU )
665 {
666  g_type __gx;
667  M_prob.evaluate ( _x, __gx, diff_order<0>() );
668 
669  for ( uint __i = 0; __i < __gx.value().size(); ++__i )
670  {
671  gx[__i] = __gx.value( __i );
672  _s[ __i ] = std::max( -gx[ __i ], options.s_theta_init );
673  }
674 
675  /*
676  for( int m = _E_g; m < ns; m++ )
677  _s[ m ] = options.s_theta_init;
678  */
679 
680  this->update_opt_objs( _x, _s, _lambda_h, _lambda_g, MU, true );
681 
682 }
683 
684 template<typename Data,template<class> class Problem>
685 void
686 SolverConstrained<Data,Problem>::evaluate_A_hat( double *_x, double *_s )
687 {
688  int n = M_prob.n();
689 
690  this->evaluate_Ag( _x );
691  this->evaluate_Ah( _x );
692 
693  // A_hat = [ Ah Ag ; 0 S ]
694  for ( int i = 0; i < n+_E_g+2*n; i++ )
695  for ( int m = 0; m < _E_h+_E_g+2*n; m++ )
696  A_hat[ i ][ m ] = 0;
697 
698  for ( int i = 0; i < n; i++ )
699  {
700  for ( int m = 0; m < _E_h; m++ )
701  A_hat[ i ][ m ] = Ah[ i ][ m ];
702 
703  for ( int m = 0; m < _E_g+2*n; m++ )
704  A_hat[ i ][ _E_h+m ] = Ag[ i ][ m ];
705  }
706 
707  for ( int m = 0; m < _E_g+2*n; m++ )
708  A_hat[ n+m ][ _E_h+m ] = _s[ m ];
709 
710 }
711 
712 template<typename Data,template<class> class Problem>
713 void SolverConstrained<Data,Problem>::print_A_hat()
714 {
715  int n = M_prob.n();
716 
717  for ( int i = 0; i < n+_E_g+2*n; i++ )
718  {
719  printf( "\n" );
720 
721  for ( int j = 0; j < _E_h+_E_g+2*n; j++ )
722  fprintf( stderr, "%14.9f", A_hat[i][j] );
723  }
724 
725  std::cerr << "\n";
726 }
727 
728 template<typename Data,template<class> class Problem>
729 void SolverConstrained<Data,Problem>::solve_AtA_system( double _r[ _E_n + _E_g+2*_E_n ], double _b[ _E_h + _E_g+2*_E_n ],
730  double _g[ _E_n + _E_g+2*_E_n], double _l[ _E_h + _E_g+2*_E_n ] )
731 {
732 
733  int n = M_prob.n();
734  int n_A_hat = n + _E_g+2*n;
735  int m_A_hat = _E_h + _E_g+2*n;
736  int n_big = n_A_hat + m_A_hat;
737  int big_N = _E_n + _E_g+2*_E_n + _E_h + _E_g+2*_E_n;
738 
739  double bigA[ big_N*big_N ], rhs[ big_N ];
740 
741  for ( int i = 0; i < n_big*n_big; i++ )
742  bigA[ i ] = 0;
743 
744  for ( int i = 0; i < n_A_hat; i++ )
745  {
746  bigA[ i + i*n_big ] = 1;
747  rhs[ i ] = _r[ i ];
748  }
749 
750  double val;
751 
752  for ( int i = 0; i < n_A_hat; i++ )
753  {
754  for ( int j = 0; j < m_A_hat; j++ )
755  {
756  val = A_hat[ i ][ j ];
757  bigA[ ( i ) + ( n_A_hat+j )*n_big ] = val;
758  bigA[ ( n_A_hat+j ) + ( i )*n_big ] = val;
759  }
760  }
761 
762  for ( int j = 0; j < m_A_hat; j++ )
763  rhs[ n_A_hat + j ] = -_b[ j ];
764 
765  /*
766  this->print_A_hat();
767  cerr << "\nA = [\n";
768  for( int i = 0; i < n_big; i++ )
769  {
770  for( int j = 0; j < n_big; j++ )
771  fprintf( stderr, "%10.1f", bigA[ i + j*n_big ]);
772  cerr << "\n";
773  }
774  cerr << "]\n\nb = [";
775  __print_vec( n_big, rhs );
776  cerr << "]\n";
777  */
778 
779  double sol[ big_N ];
780 
781 #warning disabled Lusolve
782  //__LUSolve( n_big, bigA, rhs, sol, STORAGE_GENERAL );
783 
784  for ( int i = 0; i < n_A_hat; i++ )
785  _g[ i ] = sol[ i ];
786 
787  for ( int i = 0; i < m_A_hat; i++ )
788  _l[ i ] = sol[ n_A_hat+i ];
789 
790 }
791 
792 template<typename Data,template<class> class Problem>
793 void SolverConstrained<Data,Problem>::evaluate_gradf( double *_x )
794 {
795  f_type __fx;
796  M_prob.evaluate ( _x, __fx, diff_order<1>() );
797 
798  for ( uint __i = 0; __i < _E_n; ++__i )
799  {
800  gradf[__i] = __fx.gradient( 0, __i );
801  }
802 }
803 
804 template<typename Data,template<class> class Problem>
805 void SolverConstrained<Data,Problem>::update_opt_objs( double *_x, double *_s,
806  double *_lambda_h, double *_lambda_g,
807  double MU, bool update_Hess )
808 {
809  int n = M_prob.n();
810  double _r[ _E_n + _E_g+2*_E_n ], _b[ _E_h + _E_g+2*_E_n ],
811  _dummy[ _E_n + _E_g+2*_E_n ];
812  double _lambda_LS[ _E_h + _E_g+2*_E_n ];
813 
814  this->evaluate_fxgxhx( _x );
815  this->evaluate_gradf( _x );
816 
817  for ( int i = 0; i < n; i++ )
818  _r[ i ] = -gradf[ i ];
819 
820  for ( int i = n; i < n+_E_g+2*n; i++ )
821  _r[ i ] = MU;
822 
823  for ( int i = 0; i < _E_h + _E_g+2*n; i++ )
824  _b[ i ] = 0;
825 
826  this->evaluate_A_hat( _x, _s );
827 
828  this->solve_AtA_system( _r, _b, _dummy, _lambda_LS );
829 
830  for ( int i = 0; i < _E_h; i++ )
831  _lambda_h[ i ] = _lambda_LS[ i ];
832 
833  for ( int i = 0; i < _E_g+2*n; i++ )
834  _lambda_g[ i ] = _lambda_LS[ _E_h+i ];
835 
836  this->evaluate_G( _x, _s, _lambda_h, _lambda_g, MU, update_Hess );
837 
838  /*
839  std::cerr << "\nlambda_h = ";
840  __print_vec( _E_h, _lambda_h );
841  std::cerr << "\nlambda_g = ";
842  __print_vec( _E_g+2*n, _lambda_g );
843  */
844 
845 }
846 
847 template<typename Data,template<class> class Problem>
848 void SolverConstrained<Data,Problem>::print_lambda_LS( double *_lambda_h, double *_lambda_g )
849 {
850  int n = M_prob.n();
851  std::cerr << "\nlambda_h = ";
852  __print_vec( _E_h, _lambda_h );
853  std::cerr << "\nlambda_g = ";
854  __print_vec( _E_g+2*n, _lambda_g );
855  std::cerr << "\n";
856 }
857 
858 template<typename Data,template<class> class Problem>
859 double SolverConstrained<Data,Problem>::evaluate_E_current( double *_s, double *_lambda_h, double *_lambda_g, double MU )
860 {
861  int n = M_prob.n();
862  double Ahlh[ _E_n ], Aglg[ _E_n ], dfAhlhAglg[ _E_n ], Slg_mue[ _E_g+2*_E_n], gps[ _E_g+2*_E_n];
863  double E = 0;
864 
865  for ( int i = 0; i < n; i++ )
866  {
867  Ahlh[ i ] = 0;
868 
869  for ( int m = 0; m < _E_h; m++ )
870  Ahlh[ i ] += Ah[ i ][ m ] * _lambda_h[ m ];
871 
872  Aglg[ i ] = 0;
873 
874  for ( int m = 0; m < _E_g+2*n; m++ )
875  Aglg[ i ] += Ag[ i ][ m ] * _lambda_g[ m ];
876 
877  dfAhlhAglg[ i ] = gradf[ i ] + Ahlh[ i ] + Aglg[ i ];
878  }
879 
880  E = std::max( E, __Norm_infty( n, dfAhlhAglg ) );
881 
882  for ( int m = 0; m < _E_g+2*n; m++ )
883  Slg_mue[ m ] = _s[ m ] * _lambda_g[ m ] - MU;
884 
885  E = std::max( E, __Norm_infty( _E_g+2*n, Slg_mue ) );
886 
887  E = std::max( E, __Norm_infty( _E_h, hx ) );
888 
889  for ( int m = 0; m < _E_g+2*n; m++ )
890  gps[ m ] = gx[ m ] + _s[ m ];
891 
892  E = std::max( E, __Norm_infty( _E_g+2*n, gps ) );
893 
894  /*
895  this->print_Ah();
896  this->print_Ag();
897  this->print_A_hat();
898  this->print_G();
899  __print_vec( _E_g+2*n, _s );
900  __print_vec( _E_g+2*n, _lambda_g );
901  std::cerr << "\n|| df + Ahlh + Aglg ||_infty = " << __Norm_infty( n, dfAhlhAglg );
902  std::cerr << "\nmu = " << MU;
903  std::cerr << "\n|| Slg - mue ||_infty = " << __Norm_infty( _E_g+2*n, Slg_mue );
904  std::cerr << "\n|| hx ||_infty = " << __Norm_infty( _E_h, hx );
905  std::cerr << "\n|| gx + s ||_infty = " << __Norm_infty( _E_g+2*n, gps );
906  std::cerr << "\n";
907  */
908 
909 
910  return E;
911 }
912 
913 template<typename Data,template<class> class Problem>
914 double SolverConstrained<Data,Problem>::model_m( double *_s, double *_v_til )
915 {
916  int n = M_prob.n(), nv = n+_E_g+2*n;
917  double A_hatTxv_til[ _E_h + _E_g+2*_E_n ], rhs[ _E_h + _E_g+2*_E_n ];
918 
919  for ( int i = 0; i < _E_h+_E_g+2*n; i++ )
920  {
921  A_hatTxv_til[ i ] = 0;
922 
923  for ( int j = 0; j < nv; j++ )
924  A_hatTxv_til[ i ] += A_hat[ j ][ i ] * _v_til[ j ];
925  }
926 
927  for ( int m = 0; m < _E_h; m++ )
928  rhs[ m ] = hx[ m ];
929 
930  for ( int m = 0; m < _E_g+2*n; m++ )
931  rhs[ _E_h + m ] = gx[ m ] + _s[ m ];
932 
933  double model = 2. * __Dot( _E_h+_E_g+2*n, rhs, A_hatTxv_til )
934  + __Dot( _E_h+_E_g+2*n, A_hatTxv_til, A_hatTxv_til );
935 
936  /*
937  this->print_A_hat();
938  __print_vec( nv, _v_til );
939  __print_vec( _E_h+_E_g+2*n, A_hatTxv_til );
940  std::cerr << "\n m = " << model;
941  exit(0);
942  */
943 
944 
945  return model;
946 }
947 
948 template<typename Data,template<class> class Problem>
949 void SolverConstrained<Data,Problem>::dogleg_solve( double *_vCP_til, double *_vN_til, double *_s, double _Delta,
950  double *_v_til )
951 {
952  int n = M_prob.n(), nv = n+_E_g+2*n;
953  double v_diff[ _E_n+_E_g+2*_E_n ];
954  double zeta_N = options.zeta_N, tau = options.tau;
955  double theta1 = 1, theta2 = 1, theta3 = 1;
956 
957  theta1 = std::min( theta1, zeta_N * _Delta / __Norm_2( nv, _vN_til ) );
958 
959  for ( int i = n; i < nv; i++ )
960  if ( _vN_til[ i ] < 0 )
961  theta1 = std::min( theta1, -tau/2./_vN_til[ i ] );
962 
963  assert( theta1 > 0 );
964 
965  if ( theta1 == 1 )
966  for ( int i = 0; i < nv; i++ )
967  _v_til[ i ] = theta1*_vN_til[ i ];
968 
969  else
970  {
971  for ( int i = 0; i < nv; i++ )
972  v_diff[ i ] = _vN_til[ i ] - _vCP_til[ i ];
973 
974  double a = __Dot( nv, v_diff, v_diff );
975  double b = 2. * __Dot( nv, v_diff, _vCP_til );
976  double c = __Dot( nv, _vCP_til, _vCP_til ) - ( zeta_N*_Delta )*( zeta_N*_Delta );
977  double root = b*b - 4.*a*c;
978  bool exist = false;
979 
980  if ( root >= 0 )
981  {
982  double r2 = ( ( -b + std::sqrt( root ) ) / ( 2.*a ) );
983 
984  if ( r2 > 1 )
985  {
986  exist = true;
987  theta2 = 1;
988  }
989 
990  else if ( r2 > 0 )
991  {
992  exist = true;
993  theta2 = r2;
994  }
995 
996  if ( exist )
997  for ( int i = n; i < nv; i++ )
998  if ( v_diff[ i ] < 0 )
999  theta2 = std::min( theta2, ( -tau/2. - _vCP_til[ i ] ) / ( v_diff[ i ] ) );
1000 
1001  if ( theta2 < 0 )
1002  exist = false;
1003  }
1004 
1005  if ( !exist )
1006  {
1007  theta3 = std::min( theta3, zeta_N * _Delta / __Norm_2( nv, _vCP_til ) );
1008 
1009  for ( int i = n; i < nv; i++ )
1010  if ( _vCP_til[ i ] < 0 )
1011  theta3 = std::min( theta3, -tau/2./_vCP_til[ i ] );
1012 
1013  assert( theta3 > 0 );
1014 
1015  for ( int i = 0; i < nv; i++ )
1016  _v_til[ i ] = theta3 * _vCP_til[ i ];
1017  }
1018 
1019  else
1020  for ( int i = 0; i < nv; i++ )
1021  _v_til[ i ] = ( 1-theta2 ) * _vCP_til[ i ] + theta2 * _vN_til[ i ];
1022 
1023  double theta1xvN_til[ _E_n+_E_g+2*_E_n ];
1024 
1025  for ( int i = 0; i < nv; i++ )
1026  theta1xvN_til[ i ] = theta1 * _vN_til[ i ];
1027 
1028  if ( this->model_m( _s, _v_til ) >= this->model_m( _s, theta1xvN_til ) )
1029  for ( int i = 0; i < nv; i++ )
1030  _v_til[ i ] = theta1xvN_til[ i ];
1031  }
1032 
1033  /*
1034  std::cerr << "\ntheta1 = " << theta1;
1035  std::cerr << "\ntheta2 = " << theta2;
1036  std::cerr << "\ntheta3 = " << theta3;
1037  std::cerr << "\nNorm(_v_til) = " << __Norm_2( nv, _v_til );
1038  std::cerr << "\n_v_til = ";
1039  __print_vec( nv, _v_til );
1040  std::cerr << "\ndx_n = ";
1041  exit(0);
1042  */
1043 
1044 }
1045 
1046 template<typename Data,template<class> class Problem>
1047 double SolverConstrained<Data,Problem>::get_vpred( double *_rhs, double *_v_til )
1048 {
1049  int n = M_prob.n(), nv = n + _E_g+2*n;
1050  double A_hatTxv_til[ _E_h+_E_g+2*_E_n ], rhspAv[ _E_h + _E_g+2*_E_n ];
1051 
1052  for ( int i = 0; i < _E_h+_E_g+2*n; i++ )
1053  {
1054  A_hatTxv_til[ i ] = 0;
1055 
1056  for ( int j = 0; j < nv; j++ )
1057  A_hatTxv_til[ i ] += A_hat[ j ][ i ] * _v_til[ j ];
1058 
1059  rhspAv[ i ] = _rhs[ i ] + A_hatTxv_til[ i ];
1060  }
1061 
1062  return __Norm_2( _E_h+_E_g+2*n, _rhs ) - __Norm_2( _E_h+_E_g+2*n, rhspAv );
1063 }
1064 
1065 template<typename Data,template<class> class Problem>
1066 double SolverConstrained<Data,Problem>::get_normal_steps( double *_s, double _Delta,
1067  double *_dx_normal, double *_ds_normal,
1068  double *_v_til )
1069 {
1070  int n = M_prob.n(), nv = n + _E_g+2*n;
1071  double vCP_til[ _E_n + _E_g+2*_E_n ], vN_til[ _E_n + _E_g+2*_E_n ],
1072  rhs[ _E_h + _E_g+2*_E_n ], A_hatxrhs[ _E_n + _E_g+2*_E_n ],
1073  A_hatTxA_hat[ _E_h + _E_g+2*_E_n ][ _E_h + _E_g+2*_E_n ],
1074  A_hatTxA_hat2[ _E_h + _E_g+2*_E_n ][ _E_h + _E_g+2*_E_n ],
1075  A_hat2xrhs [ _E_h+_E_g+2*_E_n ];
1076  double zeros[ _E_n + _E_g+2*_E_n ], dummy[ _E_n + _E_g+2*_E_n ], x_til[ _E_h + _E_g+2*_E_n ];
1077 
1078  for ( int m = 0; m < _E_h; m++ )
1079  rhs[ m ] = hx[ m ];
1080 
1081  for ( int m = 0; m < _E_g+2*n; m++ )
1082  rhs[ _E_h + m ] = gx[ m ] + _s[ m ];
1083 
1084  // A^T * A
1085  for ( int mi = 0; mi < _E_h+_E_g+2*n; mi++ )
1086  {
1087  for ( int mj = 0; mj < _E_h+_E_g+2*n; mj++ )
1088  {
1089  A_hatTxA_hat[ mi ][ mj ] = 0;
1090 
1091  for ( int i = 0; i < nv; i++ )
1092  A_hatTxA_hat[ mi ][ mj ] += A_hat[ i ][ mi ] * A_hat[ i ][ mj ];
1093  }
1094  }
1095 
1096  // (A^T * A )^2
1097  for ( int mi = 0; mi < _E_h+_E_g+2*n; mi++ )
1098  {
1099  for ( int mj = 0; mj < _E_h+_E_g+2*n; mj++ )
1100  {
1101  A_hatTxA_hat2[ mi ][ mj ] = 0;
1102 
1103  for ( int i = 0; i < _E_h+_E_g+2*n; i++ )
1104  A_hatTxA_hat2[ mi ][ mj ] += A_hatTxA_hat[ i ][ mi ] * A_hatTxA_hat[ i ][ mj ];
1105  }
1106  }
1107 
1108  for ( int i = 0; i < nv; i++ )
1109  {
1110  zeros[ i ] = 0;
1111  A_hatxrhs[ i ] = 0;
1112 
1113  for ( int m = 0; m < _E_h+_E_g+2*n; m++ )
1114  A_hatxrhs[ i ] += A_hat[ i ][ m ] * rhs[ m ];
1115  }
1116 
1117  for ( int i = 0; i < _E_h+_E_g+2*n; i++ )
1118  {
1119  A_hat2xrhs[ i ] = 0;
1120 
1121  for ( int m = 0; m < _E_h+_E_g+2*n; m++ )
1122  A_hat2xrhs[ i ] += A_hatTxA_hat2[ i ][ m ] * rhs[ m ];
1123  }
1124 
1125  double norm = __Norm_2( nv, A_hatxrhs );
1126  double alpha = norm*norm / __Dot( _E_h+_E_g+2*n, rhs, A_hat2xrhs );
1127 
1128  for ( int i = 0; i < nv; i++ )
1129  vCP_til[ i ] = -alpha * A_hatxrhs[ i ];
1130 
1131  this->solve_AtA_system( zeros, rhs, dummy, x_til );
1132 
1133  for ( int i = 0; i < n+_E_g+2*n; i++ )
1134  {
1135  vN_til[ i ] = 0;
1136 
1137  for ( int m = 0; m < _E_h+_E_g+2*n; m++ )
1138  vN_til[ i ] += -A_hat[ i ][ m ] * x_til[ m ];
1139  }
1140 
1141  this->dogleg_solve( vCP_til, vN_til, _s, _Delta, _v_til );
1142 
1143  for ( int i = 0; i < n; i++ )
1144  _dx_normal[ i ] = _v_til[ i ];
1145 
1146  for ( int i = 0; i < _E_g+2*n; i++ )
1147  _ds_normal[ i ] = _s[ i ] * _v_til[ n+i ];
1148 
1149  /*
1150  std::cerr << "\nalpha = " << alpha;
1151  std::cerr << "\nATA = \n";
1152  for( int i = 0; i < _E_h+_E_g+2*n; i++ )
1153  {
1154  printf("\n");
1155  for( int j = 0; j < _E_h+_E_g+2*n; j++ )
1156  fprintf( stderr, "%10.1f", A_hatTxA_hat2[i][j]);
1157  }
1158  std::cerr << "\n";
1159  std::cerr << "\nhx = [";
1160  __print_vec( _E_h, hx );
1161  std::cerr << "]\n";
1162  std::cerr << "\ngx = [";
1163  __print_vec( _E_g+2*n, gx );
1164  std::cerr << "]\n";
1165  std::cerr << "\ns = [";
1166  __print_vec( _E_g+2*n, _s );
1167  std::cerr << "]\n";
1168  std::cerr << "\nA_hat = [";
1169  this->print_A_hat();
1170  std::cerr << "]\n";
1171  std::cerr << "\nrhs = [";
1172  __print_vec( _E_h+_E_g+2*n, rhs );
1173  std::cerr << "]\n";
1174  std::cerr << "\nA_hatxrhs = [";
1175  __print_vec( nv, A_hatxrhs );
1176  std::cerr << "]\n";
1177  __print_vec( n, _dx_normal );
1178  std::cerr << "\nds_n = ";
1179  __print_vec( _E_g+2*n, _ds_normal );
1180  exit(0);
1181  */
1182 
1183  return this->get_vpred( rhs, _v_til );
1184 
1185 }
1186 
1187 template<typename Data,template<class> class Problem>
1188 double SolverConstrained<Data,Problem>::get_q( double *_v_til, double *_w_til, double MU )
1189 {
1190  int n = M_prob.n(), nd = n + _E_g+2*n;
1191  double d_til[ _E_n+_E_g+2*_E_n ], gr[ _E_n+_E_g+2*_E_n ], Gd_til[ _E_n+_E_g+2*_E_n ];
1192 
1193  __SumVecs( nd, _v_til, _w_til, d_til );
1194 
1195  for ( int i = 0; i < n; i++ )
1196  gr[ i ] = gradf[ i ];
1197 
1198  for ( int i = n; i < nd; i++ )
1199  gr[ i ] = -MU;
1200 
1201  for ( int i = 0; i < nd; i++ )
1202  {
1203  Gd_til[ i ] = 0;
1204 
1205  for ( int j = 0; j < nd; j++ )
1206  Gd_til[ i ] += G[ i ][ j ] * d_til[ j ];
1207  }
1208 
1209  return __Dot( nd, gr, d_til ) + 0.5 * __Dot( nd, d_til, Gd_til );
1210 }
1211 
1212 template<typename Data,template<class> class Problem>
1213 void SolverConstrained<Data,Problem>::get_Pr( double *_r, double *_Pr )
1214 {
1215  int n = M_prob.n(), nw = n + _E_g+2*n;
1216  double zeros[ _E_n + _E_g+2*_E_n ], dummy[ _E_n + _E_g+2*_E_n ],
1217  Atr[ _E_h + _E_g+2*_E_n ], x_Atr[ _E_h + _E_g+2*_E_n ];
1218 
1219  for ( int i = 0; i < nw; i++ )
1220  {
1221  zeros[ i ] = 0;
1222  }
1223 
1224  for ( int i = 0; i < _E_h+_E_g+2*n; i++ )
1225  {
1226  Atr[ i ] = 0;
1227 
1228  for ( int j = 0; j < nw; j++ )
1229  Atr[ i ] += A_hat[ j ][ i ] * _r[ j ];
1230  }
1231 
1232  this->solve_AtA_system( zeros, Atr, dummy, x_Atr );
1233 
1234  for ( int i = 0; i < nw; i++ )
1235  {
1236  _Pr[ i ] = _r[i];
1237 
1238  for ( int m = 0; m < _E_h+_E_g+2*n; m++ )
1239  _Pr[ i ] += -A_hat[ i ][ m ] * x_Atr[ m ];
1240  }
1241 }
1242 
1243 template<typename Data,template<class> class Problem>
1244 bool SolverConstrained<Data,Problem>::slack_feasible( double *_w_til, double *_v_til )
1245 {
1246  int n = M_prob.n();
1247  double tau = options.tau;
1248  bool feasible = true;
1249 
1250  for ( int i = n; i < n+_E_g+2*n; i++ )
1251  {
1252  if ( _w_til[ i ] + _v_til[ i ] + tau < 1e-13 )
1253  feasible = false;
1254 
1255  //std::cerr << "\n" << i << ": " << _w_til[ i ] + _v_til[ i ] + tau;
1256  }
1257 
1258  return feasible;
1259 }
1260 
1261 template<typename Data,template<class> class Problem>
1262 double SolverConstrained<Data,Problem>::get_tangen_steps( double *_s, double _Delta, double *_v_til, double MU,
1263  double *_dx_tangen, double *_ds_tangen,
1264  bool &truss_exit, int &n_CGiter )
1265 {
1266  int n = M_prob.n();
1267  int nw = n + _E_g+2*n;
1268  double last_slack_feasible_rg;
1269  double w_til[ _E_n+_E_g+2*_E_n ], wP_til[ _E_n+_E_g+2*_E_n ],
1270  last_slack_feasible_w_til[ _E_n+_E_g+2*_E_n ],
1271  r[ _E_n+_E_g+2*_E_n ], rP[ _E_n+_E_g+2*_E_n ], Gp[ _E_n+_E_g+2*_E_n ],
1272  g[ _E_n+_E_g+2*_E_n ], gP[ _E_n+_E_g+2*_E_n ],
1273  p[ _E_n+_E_g+2*_E_n ], pP[ _E_n+_E_g+2*_E_n ],
1274  last_slack_feasible_p[ _E_n+_E_g+2*_E_n ];
1275 
1276  truss_exit = false;
1277  n_CGiter = 0;
1278 
1279  for ( int i = 0; i < nw; i++ )
1280  {
1281  Gp[ i ] = 0;
1282 
1283  for ( int j = 0; j < nw; j++ )
1284  Gp[ i ] += G[ i ][ j ] * _v_til[ j ];
1285  }
1286 
1287  for ( int i = 0; i < n; i++ )
1288  r[ i ] = gradf[ i ] + Gp[ i ];
1289 
1290  for ( int i = n; i < nw; i++ )
1291  r[ i ] = -MU + Gp[ i ];
1292 
1293  get_Pr( r, g );
1294 
1295  /*
1296  this->print_A_hat();
1297  __print_vec( nw, r );
1298  __print_vec( nw, g );
1299  this->print_G();
1300  */
1301 
1302  for ( int i = 0; i < nw; i++ )
1303  {
1304  w_til[ i ] = 0;
1305  wP_til[ i ] = 0;
1306  last_slack_feasible_w_til[ i ] = 0;
1307  p[ i ] = -g[ i ];
1308  pP[ i ] = -g[ i ];
1309  last_slack_feasible_p[ i ] = -g[ i ];
1310  }
1311 
1312  last_slack_feasible_rg = __Dot( nw, g, r );
1313 
1314  double tol = 0.01 * std::sqrt( __Dot( nw, g, r ) );
1315  int iter = 0;
1316  bool done = false;
1317  double gamma, alpha, beta;
1318 
1319  while ( iter < 2*( n-_E_h ) && !done )
1320  {
1321 #ifdef DEBUG_OUT_CG
1322  std::cerr << "\nCGiter = " << iter << ",\t r'g = " << __Dot( nw, g, r );
1323 #endif
1324 
1325  for ( int i = 0; i < nw; i++ )
1326  {
1327  Gp[ i ] = 0;
1328 
1329  for ( int j = 0; j < nw; j++ )
1330  Gp[ i ] += G[ i ][ j ] * p[ j ];
1331  }
1332 
1333  gamma = __Dot( nw, p, Gp );
1334 
1335  if ( gamma <= 0 )
1336  {
1337  double a = __Dot( nw, p, p );
1338  double b = 2. * __Dot( nw, w_til, p );
1339  double c = __Dot( nw, w_til, w_til ) + __Dot( nw, _v_til, _v_til ) - _Delta*_Delta;
1340  double thetaa = ( -b + std::sqrt( b*b - 4.*a*c ) ) / ( 2.*a );
1341  assert( thetaa > 0 );
1342 
1343  for ( int i = 0; i < nw; i++ )
1344  wP_til[ i ] = w_til[ i ] + thetaa * p[ i ];
1345 
1346  done = true;
1347  truss_exit = true;
1348 #ifdef DEBUG_OUT_CG
1349  std::cerr << "\nCG indefinitness EXIT: ||wP_til||^2 + ||v_til||^2 = "
1350  << __Dot( nw, wP_til, wP_til ) + __Dot( nw, _v_til, _v_til )
1351  << " ,\t Delta^2 = "
1352  << _Delta*_Delta << "\n";
1353 #endif
1354  }
1355 
1356  if ( !done )
1357  {
1358  alpha = __Dot( nw, r, g ) / gamma;
1359 
1360  for ( int i = 0; i < nw; i++ )
1361  wP_til[ i ] = w_til[ i ] + alpha * p[ i ];
1362 
1363  if ( __Dot( nw, wP_til, wP_til ) + __Dot( nw, _v_til, _v_til ) > _Delta*_Delta )
1364  {
1365  double a = __Dot( nw, p, p );
1366  double b = 2. * __Dot( nw, w_til, p );
1367  double c = __Dot( nw, w_til, w_til ) + __Dot( nw, _v_til, _v_til ) - _Delta*_Delta;
1368  double thetab = ( -b + std::sqrt( b*b - 4.*a*c ) ) / ( 2.*a );
1369  assert( thetab > 0 );
1370 
1371  for ( int i = 0; i < nw; i++ )
1372  wP_til[ i ] = w_til[ i ] + thetab * p[ i ];
1373 
1374  done = true;
1375  truss_exit = true;
1376 #ifdef DEBUG_OUT_CG
1377  std::cerr << "\nCG truss EXIT: ||wP_til||^2 + ||v_til||^2 = "
1378  << __Dot( nw, wP_til, wP_til ) + __Dot( nw, _v_til, _v_til )
1379  << " ,\t Delta^2 = "
1380  << _Delta*_Delta << "\n";
1381 #endif
1382  }
1383  }
1384 
1385  if ( !done )
1386  {
1387  for ( int i = 0; i < nw; i++ )
1388  rP[ i ] = r[ i ] + alpha * Gp[ i ];
1389 
1390  get_Pr( rP, gP );
1391 
1392  /*
1393  this->print_A_hat();
1394  __print_vec( nw, rP );
1395  __print_vec( nw, gP );
1396  */
1397  if ( __Dot( nw, gP, rP ) < tol )
1398  {
1399  done = true;
1400 #ifdef DEBUG_OUT_CG
1401  std::cerr << "\nCG tol EXIT: \trP'gP = "
1402  << __Dot( nw, gP, rP )
1403  << " \t ( tol = "
1404  << tol << " )\n";
1405 #endif
1406  }
1407  }
1408 
1409  if ( !done )
1410  {
1411  beta = __Dot( nw, gP, rP ) / __Dot( nw, g, r );
1412 
1413  for ( int i = 0; i < nw; i++ )
1414  {
1415  pP[ i ] = -gP[ i ] + beta * p[ i ];
1416  w_til[ i ] = wP_til[ i ];
1417  r[ i ] = rP[ i ];
1418  g[ i ] = gP[ i ];
1419  p[ i ] = pP[ i ];
1420  }
1421 
1422  if ( this->slack_feasible( wP_til, _v_til ) )
1423  {
1424  for ( int i = 0; i < nw; i++ )
1425  {
1426  last_slack_feasible_w_til[ i ] = wP_til[ i ];
1427  last_slack_feasible_p[ i ] = p[ i ];
1428  }
1429 
1430  last_slack_feasible_rg = __Dot( nw, g, r );
1431  }
1432  }
1433 
1434  iter++;
1435  }
1436 
1437  if ( !( this->slack_feasible( wP_til, _v_til ) ) )
1438  {
1439  double a = __Dot( nw, last_slack_feasible_p, last_slack_feasible_p );
1440  double b = 2. * __Dot( nw, last_slack_feasible_w_til, last_slack_feasible_p );
1441  double c = __Dot( nw, last_slack_feasible_w_til, last_slack_feasible_w_til )
1442  + __Dot( nw, _v_til, _v_til ) - _Delta*_Delta;
1443  double thetac = ( -b + std::sqrt( b*b - 4.*a*c ) ) / ( 2.*a );
1444  assert( thetac > 0 );
1445  double tau = options.tau;
1446  truss_exit = true;
1447 
1448  for ( int i = n; i < nw; i++ )
1449  if ( last_slack_feasible_p[ i ] < 0 )
1450  {
1451  thetac = std::min( thetac,
1452  ( -tau - _v_til[ i ] - last_slack_feasible_w_til[ i ] )
1453  / last_slack_feasible_p[ i ] );
1454  }
1455 
1456  assert( thetac > 0 );
1457 
1458  for ( int i = 0; i < nw; i++ )
1459  wP_til[ i ] = last_slack_feasible_w_til[ i ] + thetac * last_slack_feasible_p[ i ];
1460 
1461 #ifdef DEBUG_OUT_CG
1462  std::cerr << "\nCorrected slack feasibility violation of wP_til: tau = " << tau << ", wP_til_s + v_til_s = " ;
1463 
1464  for ( int i = n; i < nw; i++ )
1465  std::cerr << "\n " << wP_til[ i ] + _v_til[ i ];
1466 
1467  std::cerr<< "\nr'g = " << last_slack_feasible_rg << "\n";
1468 #endif
1469  }
1470 
1471  n_CGiter = iter;
1472 
1473  for ( int i = 0; i < n; i++ )
1474  _dx_tangen[ i ] = wP_til[ i ];
1475 
1476  for ( int i = 0; i < _E_g+2*n; i++ )
1477  _ds_tangen[ i ] = _s[ i ] * wP_til[ n+i ];
1478 
1479  /*
1480  __print_vec( nw, _v_til );
1481  std::cerr << "\nwPtil = ";
1482  __print_vec( nw, wP_til );
1483 
1484  this->print_G();
1485  __print_vec( nw, wP_til );
1486  for( int i = 0; i < nw; i++ )
1487  {
1488  Gp[ i ] = 0;
1489  for( int j = 0; j < nw; j++ )
1490  Gp[ i ] += G[ i ][ j ] * wP_til[ j ];
1491  }
1492  __print_vec( nw, Gp );
1493 
1494  __print_vec( nw, wP_til );
1495  __print_vec( nw, _v_til );
1496  std::cerr << "\ntol = "<< tol;
1497  this->print_G();
1498  __print_vec( nw, Gp );
1499  __print_vec( n, gradf );
1500  this->print_A_hat();
1501  __print_vec( nw, r );
1502  __print_vec( nw, g );
1503  //exit(0);
1504  */
1505 
1506 
1507  return this->get_q( _v_til, wP_til, MU );
1508 
1509 }
1510 
1511 template<typename Data,template<class> class Problem>
1512 double SolverConstrained<Data,Problem>::get_steps( double *_s, double _Delta, double MU,
1513  double *_dx, double *_ds, double &_vpred, double &_q,
1514  bool &truss_exit, int &n_CGiter )
1515 {
1516  int n = M_prob.n();
1517  double dx_normal[ _E_n ], ds_normal[ _E_g+2*_E_n ], v_til[ _E_n+_E_g+2*_E_n ];
1518  double dx_tangen[ _E_n ], ds_tangen[ _E_g+2*_E_n ], d[ _E_n+_E_g+2*_E_n ];
1519 
1520  _vpred = this->get_normal_steps( _s, _Delta, dx_normal, ds_normal, v_til );
1521  _q = this->get_tangen_steps( _s, _Delta, v_til, MU, dx_tangen, ds_tangen, truss_exit, n_CGiter );
1522 
1523  __SumVecs( n, dx_normal, dx_tangen, _dx );
1524  __SumVecs( _E_g+2*n, ds_normal, ds_tangen, _ds );
1525 
1526  /*
1527  std::cerr << "\n dx_normal = ";
1528  __print_vec( n, dx_normal );
1529  std::cerr << "\n ds_normal = ";
1530  __print_vec( _E_g+2*n, ds_normal );
1531  std::cerr << "\n dx_tangen = ";
1532  __print_vec( n, dx_tangen );
1533  std::cerr << "\n ds_tangen = ";
1534  __print_vec( _E_g+2*n, ds_tangen );
1535  exit(0);
1536  */
1537 
1538  for ( int i = 0; i < n; i++ )
1539  d[ i ] = _dx[ i ];
1540 
1541  for ( int i = 0; i < _E_g+2*n; i++ )
1542  d[ n+i ] = _ds[ i ];
1543 
1544  return __Norm_2( n+_E_g+2*n, d );
1545 }
1546 
1547 template<typename Data,template<class> class Problem>
1548 double SolverConstrained<Data,Problem>::get_phi( double *_x, double *_s, double MU,
1549  double nu, double *rhs )
1550 {
1551  double phi_val;
1552 
1553 #if 0
1554 
1555  int n = M_prob.n();
1556  double f_at_x, g_at_x[ _E_g+2*_E_n ], h_at_x[ _E_h ];
1557 
1558 
1559  M_prob.fn_double( _x, f_at_x );
1560  M_prob.gn_double( _x, g_at_x );
1561 
1562  for ( int i = 0; i < n; i++ )
1563  {
1564  g_at_x[ _E_g+i ] = M_prob.x_l( i ) - _x[i];
1565  g_at_x[ _E_g+n+i ] = _x[i] - M_prob.x_u( i );
1566  }
1567 
1568  M_prob.hn_double( _x, h_at_x );
1569 
1570  for ( int m = 0; m < _E_h; m++ )
1571  rhs[ m ] = h_at_x[ m ];
1572 
1573  for ( int m = 0; m < _E_g+2*n; m++ )
1574  rhs[ _E_h + m ] = g_at_x[ m ] + _s[ m ];
1575 
1576  phi_val = f_at_x;
1577 
1578  for ( int i = 0; i < _E_g+2*n; i++ )
1579  phi_val += -MU * std::log( _s[ i ] );
1580 
1581  phi_val += nu * __Norm_2( _E_g+2*n, rhs );
1582 #else
1583  f_type __fx;
1584  g_type __gx;
1585  h_type __hx;
1586 
1587  // get the order 0 information only
1588  M_prob.evaluate ( _x, __fx, diff_order<0>() );
1589  M_prob.evaluate ( _x, __gx, diff_order<0>() );
1590  M_prob.evaluate ( _x, __hx, diff_order<0>() );
1591 
1592  double __f_at_x = __fx.value( 0 );
1593  double __g_at_x[ _E_g+2*_E_n ];
1594 
1595  for ( uint __i = 0; __i < __gx.value().size(); ++__i )
1596  __g_at_x[__i] = __gx.value( __i );
1597 
1598  double __h_at_x[ _E_h ];
1599 
1600  for ( uint __i = 0; __i < _E_h; ++__i )
1601  __h_at_x[__i] = __hx.value( __i );
1602 
1603  for ( int m = 0; m < _E_h; m++ )
1604  rhs[ m ] = __h_at_x[ m ];
1605 
1606  for ( int m = 0; m < _E_g+2*_E_n; m++ )
1607  rhs[ _E_h + m ] = __g_at_x[ m ] + _s[ m ];
1608 
1609  phi_val = __f_at_x;
1610 
1611  for ( int i = 0; i < _E_g+2*_E_n; i++ )
1612  phi_val += -MU * std::log( _s[ i ] );
1613 
1614  phi_val += nu * __Norm_2( _E_g+2*_E_n, rhs );
1615 #endif
1616  return phi_val;
1617 }
1618 
1619 template<typename Data,template<class> class Problem>
1620 double SolverConstrained<Data,Problem>::get_phi( double *_x, double *_s, double MU, double nu )
1621 {
1622  double rhs[ _E_h + _E_g+2*_E_n ];
1623  double phi_val = this->get_phi( _x, _s, MU, nu, rhs );
1624 
1625  return phi_val;
1626 }
1627 
1628 template<typename Data,template<class> class Problem>
1629 double SolverConstrained<Data,Problem>::get_SOC( double *_s, double *_rhs_xs_pd,
1630  double *_xy, double *_sy )
1631 {
1632  int n = M_prob.n();
1633 
1634  // Second order corrections y:
1635  double zeros[ _E_n + _E_g+2*_E_n ], dummy[ _E_n + _E_g+2*_E_n ],
1636  x_til[ _E_h + _E_g+2*_E_n ], y[ _E_n + _E_g+2*_E_n ];
1637 
1638  for ( int i = 0; i < n+_E_g+2*_E_n; i++ )
1639  zeros[ i ] = 0;
1640 
1641  this->solve_AtA_system( zeros, _rhs_xs_pd, dummy, x_til );
1642 
1643  for ( int i = 0; i < n+_E_g+2*n; i++ )
1644  {
1645  y[ i ] = 0;
1646 
1647  for ( int m = 0; m < _E_h+_E_g+2*n; m++ )
1648  y[ i ] += A_hat[ i ][ m ] * x_til[ m ];
1649  }
1650 
1651  for ( int i = 0; i < n; i++ )
1652  _xy[ i ] = y[ i ];
1653 
1654  for ( int i = 0; i < _E_g+2*n; i++ )
1655  _sy[ i ] = y[ n+i ];
1656 
1657  return __Norm_2( n+_E_g+2*n, y );
1658 }
1659 
1660 template<typename Data,template<class> class Problem>
1661 double SolverConstrained<Data,Problem>::get_phi( double *_s, double MU, double nu )
1662 {
1663  int n = M_prob.n();
1664  double rhs[ _E_h + _E_g+2*_E_n ];
1665  double phi_val;
1666 
1667  for ( int m = 0; m < _E_h; m++ )
1668  rhs[ m ] = hx[ m ];
1669 
1670  for ( int m = 0; m < _E_g+2*n; m++ )
1671  rhs[ _E_h + m ] = gx[ m ] + _s[ m ];
1672 
1673  phi_val = fx;
1674 
1675  for ( int i = 0; i < _E_g+2*n; i++ )
1676  phi_val += -MU * std::log( _s[ i ] );
1677 
1678  phi_val += nu * __Norm_2( _E_g+2*n, rhs );
1679 
1680  /*
1681  std::cerr << "\nfx = " << fx << "\n";
1682  std::cerr << "\ns = \n";
1683  __print_vec( _E_g+2*n, _s );
1684  std::cerr << "\ngx = \n";
1685  __print_vec( _E_g+2*n, gx );
1686  std::cerr << "\nrhs( x, s ) = \n";
1687  __print_vec( _E_h+_E_g+2*n, rhs );
1688  std::cerr << "\nhn = \n";
1689  __print_vec( _E_h, hx );
1690  */
1691 
1692  return phi_val;
1693 }
1694 
1695 
1696 template<typename Data,template<class> class Problem>
1697 bool
1699 {
1700  int n = M_prob.n();
1701 
1702  M_solver_stats.clear();
1703 
1704 
1705  std::cout << "\n[SolverConstrained<Data,Problem>::optimize()]: Optimizing...\n";
1706 
1707  // Controlling parameters
1708  double Delta = std::min( options.Delta_init, ( double )Max_Allowed_Delta );
1709  double eMU = options.eMU_init;
1710  double MU = options.MU_init;
1711  double MU_last = MU;
1712  double nu = options.nu0_init;
1713  double rho_N = options.rho_N;
1714 
1715  double rhs_xs_p_d[ _E_h + _E_g+2*_E_n ];
1716 
1717  double dx[ _E_n ];
1718  double x_p_dx[ _E_n ];
1719  double s[ _E_g+2*_E_n ];
1720  double ds[ _E_g+2*_E_n ];
1721  double s_p_ds[ _E_g+2*_E_n ];
1722  double lambda_h[ _E_h ], lambda_g[ _E_g+2*_E_n ];
1723  double vpred, q, pred, ared, gamma, norm_d;
1724 
1725  M_prob.copy_x0_to_x( __x );
1726  this->initialize_solver_with_x( __x, s, lambda_h, lambda_g, MU );
1727 
1728  double E = this->evaluate_E_current( s, lambda_h, lambda_g, MU );
1729 
1730  int err_iter = 0, hom_iter = 0;
1731 
1732  while ( err_iter < 500 && hom_iter < options.max_hom_iter && E > options.eTOL )
1733  {
1734 
1735  double max_TR_iter = options.max_TR_iter;
1736 
1737  if ( hom_iter+1 >= options.max_hom_iter )
1738  max_TR_iter = 2 * max_TR_iter;
1739 
1740  err_iter++;
1741 
1742 #ifdef DEBUG_OUT_HM
1743  std::cerr << "\n================================================================================"
1744  << "\n homotopy iter = " << hom_iter
1745  << ", MU = " << MU
1746  << ", E(x,y,MU) = " << this->evaluate_E_current( s, lambda_h, lambda_g, MU )
1747  << "\n================================================================================\n";
1748 #endif
1749 
1750  int iter = 0;
1751 
1752  while ( iter < max_TR_iter && E > eMU )
1753  {
1754 
1755  E = this->evaluate_E_current( s, lambda_h, lambda_g, MU );
1756 
1757  double Delta_used = Delta;
1758  bool truss_exit = false, SOC = false;
1759  int n_CGiter = 0;
1760 
1761  bool evaluate_Hess = false;
1762 
1763  if ( iter%2 == 0 )
1764  evaluate_Hess = true;
1765 
1766  norm_d = this->get_steps( s, Delta, MU, dx, ds, vpred, q, truss_exit, n_CGiter );
1767 
1768  if ( q > 0 )
1769  nu = std::max( 1.5*nu, q / ( ( 1 - rho_N ) * vpred ) );
1770 
1771  pred = -q + nu * vpred;
1772  __SumVecs( n, __x, dx, x_p_dx );
1773  __SumVecs( _E_g+2*n, s, ds, s_p_ds );
1774  ared = get_phi( s, MU, nu ) - get_phi( x_p_dx, s_p_ds, MU, nu, rhs_xs_p_d );
1775  gamma = ared / pred;
1776 
1777 #ifdef DEBUG_OUT_TR
1778  std::cerr << "\n-------------------------\n"
1779  << "\nTRiter = " << iter
1780  << "\n E(x,s,MU) = " << E
1781  << "\n gamma = " << gamma
1782 
1783  << "\n vpred = " << vpred
1784  << "\n pred = " << pred
1785  << "\n nu = " << nu
1786  << "\n ared = " << ared
1787  << "\n norm_d = " << norm_d
1788  << "\n Delta = " << Delta
1789 
1790  ;
1791  __print_vec( n, dx );
1792  __print_vec( _E_g+2*n, ds );
1793 #endif
1794 
1795  if ( gamma >= options.eta )
1796  {
1797  __CopyVecs( n, x_p_dx, __x );
1798  __CopyVecs( _E_g+2*n, s_p_ds, s );
1799 
1800  if ( gamma >= options.rho_big )
1801  Delta = std::max( options.rho_increase_big * norm_d , Delta );
1802 
1803  else if ( gamma >= options.rho_small )
1804  Delta = std::max( options.rho_increase_small * norm_d , Delta );
1805 
1806  Delta = std::min( Delta, ( double )Max_Allowed_Delta );
1807  }
1808 
1809  else
1810  {
1811  SOC = false;
1812  double xy[ _E_n ], sy[ _E_g+2*_E_n ];
1813  double norm_xypsy = this->get_SOC( s_p_ds, rhs_xs_p_d, xy, sy );
1814 #ifdef DEBUG_OUT_TR
1815  std::cerr << "\n SECOND ORDER CORRECTION:"
1816  << "\n ||xy||+||sy|| = " << norm_xypsy;
1817 #endif
1818 
1819  if ( norm_xypsy > 0 )
1820  {
1821  bool sy_feasible = true;
1822  double tau = options.tau;
1823 
1824  for ( int i = 0; i < _E_g+2*n; i++ )
1825  if ( ds[ i ] + sy[ i ] + tau*s[ i ] < 0 )
1826  sy_feasible = false;
1827 
1828  if ( sy_feasible )
1829  {
1830  __SumVecs( n, __x, dx, x_p_dx );
1831  __SumVecs( _E_g+2*n, s, ds, s_p_ds );
1832 
1833  __SumVecs( n, x_p_dx, xy, x_p_dx );
1834  __SumVecs( _E_g+2*n, s_p_ds, sy, s_p_ds );
1835 
1836  ared = get_phi( s, MU, nu ) - get_phi( x_p_dx, s_p_ds, MU, nu );
1837  gamma = ared / pred;
1838 #ifdef DEBUG_OUT_TR
1839  std::cerr << "\n gamma = " << gamma;
1840 #endif
1841 
1842  if ( gamma >= options.eta )
1843  {
1844  SOC = true;
1845  __CopyVecs( n, x_p_dx, __x );
1846  __CopyVecs( _E_g+2*n, s_p_ds, s );
1847  // (Delta = Delta)
1848  }
1849  }
1850  }
1851 
1852  if ( !SOC )
1853  Delta = options.rho_decrease * Delta;
1854  }
1855 
1856  if ( gamma >= options.eta ) // (successful step)
1857  {
1858  if ( !options.Primal_Dual )
1859  MU_last = MU;
1860 
1861  this->update_opt_objs( __x, s, lambda_h, lambda_g, MU_last, evaluate_Hess );
1862 
1863  M_solver_stats.push_all( MU, E, this->evaluate_E_current( s, lambda_h, lambda_g, 0 ),
1864  n_CGiter, truss_exit, SOC,
1865  Delta_used, ared, pred, gamma );
1866  iter++;
1867  MU_last = MU;
1868  }
1869 
1870 #ifdef DEBUG_OUT_TR
1871  std::cerr << "\n E(x+Dx,s+ds,MU) = " << E << " ( MU = " << MU << ")\n";
1872 #endif
1873  }
1874 
1875  MU = MU * options.theta_N;
1876 
1877  // be careful with reducing eMU because we reach Floating point exception
1878  // so cut out to 0 when reaches 1e-10
1879  if ( eMU < 1e-10 )
1880  eMU = 0;
1881 
1882  else
1883  eMU = eMU * options.theta_N * eMU;
1884 
1885  // reinitialize Delta
1886  Delta = options.Delta_init;
1887 
1888  // reinitialize Delta
1889  nu = options.nu0_init;
1890 
1891  // reevaluate E to start over the homotopy loop
1892  E = this->evaluate_E_current( s, lambda_h, lambda_g, 0 );
1893 
1894 #ifdef DEBUG_OUT_HM
1895  std::cerr << "\n E(x+Dx,s+ds,0) = " << E << "\n";
1896 #endif
1897  // go to next homotopy loop
1898  hom_iter++;
1899  }
1900 
1901  M_solver_stats.show();
1902  //std::cerr << "Number of homotopy levels: " << hom_iter;
1903 
1904  /*
1905  std::cerr << "\n\nx = \n";
1906  __print_vec( n, __x );
1907  std::cerr << "\n\ns = \n";
1908  __print_vec( _E_g+2*n, s );
1909  std::cerr << "\n";
1910  */
1911 
1912  //M_prob.print_stationaryN_x( __x, s, lambda_h, lambda_g );
1913 
1914  //std::cerr << "\n";
1915 
1916  return true;
1917 }
1918 
1919 //
1920 // Stats
1921 //
1922 template<typename Data, template<class> class Problem>
1923 void
1925 {
1926  x_hstr.clear();
1927  mu_hstr.clear();
1928  E_hstr.clear();
1929  E0_hstr.clear();
1930 
1931  n_CGiter_hstr.clear();
1932  truss_exit_hstr.clear();
1933  SOC_hstr.clear();
1934 
1935  Delta_hstr.clear();
1936  ared_hstr.clear();
1937  pred_hstr.clear();
1938  gamma_hstr.clear();
1939 }
1940 
1941 template<typename Data, template<class> class Problem>
1942 void
1944 {
1945  for ( int i = 0; i < solver_type::_E_n; i++ )
1946  x_hstr.push_back( __x[i] );
1947 }
1948 
1949 template<typename Data, template<class> class Problem>
1950 void
1952  int n_CGiter,
1953  bool truss_exit, bool SOC,
1954  double Delta, double ared, double pred, double gamma )
1955 {
1956  //push_x_in_hstr( _x );
1957  mu_hstr.push_back( mu );
1958  E_hstr.push_back( E );
1959  E0_hstr.push_back( E0 );
1960  n_CGiter_hstr.push_back( n_CGiter );
1961  truss_exit_hstr.push_back( truss_exit );
1962  SOC_hstr.push_back( SOC );
1963  Delta_hstr.push_back( Delta );
1964  ared_hstr.push_back( ared );
1965  pred_hstr.push_back( pred );
1966  gamma_hstr.push_back( gamma );
1967 }
1968 
1969 
1970 template<typename Data, template<class> class Problem>
1971 void
1973 {
1974  iter = mu_hstr.size();
1975  std::cerr << "\nScaled Trust-Region Method Statistics:\n\n";
1976 
1977  std::cerr << " k mu E(x,s;0) E(x,s;mu) CGiters trustEx SOC Delta^k ared^k pred^k gamma^k \n";
1978  std::cerr << "------------------------------------------------------------------------------------------------\n";
1979 
1980  double mu_last = 0;
1981 
1982  for ( int k = 0; k < iter; k++ )
1983  {
1984  if ( mu_hstr[ k ] != mu_last )
1985  {
1986  fprintf( stderr, "%3d ", k+1 );
1987  fprintf ( stderr, " %3.1e", mu_hstr[ k ] );
1988  fprintf ( stderr, " %10.3e", E0_hstr[ k ] );
1989  fprintf ( stderr, " %10.3e", E_hstr[ k ] );
1990  fprintf ( stderr, " %5d", n_CGiter_hstr[ k ] );
1991 
1992  if ( truss_exit_hstr[ k ] == true )
1993  fprintf ( stderr, " -Y- " );
1994 
1995  else
1996  fprintf ( stderr, " N " );
1997 
1998  if ( SOC_hstr[ k ] == true )
1999  fprintf ( stderr, "-Y- " );
2000 
2001  else
2002  fprintf ( stderr, " N " );
2003 
2004  fprintf ( stderr, " %3.1e", Delta_hstr[ k ] );
2005  fprintf ( stderr, " %9.2e", ared_hstr[ k ] );
2006  fprintf ( stderr, " %9.2e", pred_hstr[ k ] );
2007  fprintf ( stderr, " %7.1f", gamma_hstr[ k ] );
2008  std::cerr << std::endl;
2009  mu_last = mu_hstr[ k ];
2010  }
2011 
2012  else if ( !verbose )
2013  {
2014  fprintf( stderr, "%3d ", k+1 );
2015  fprintf ( stderr, " " );
2016  fprintf ( stderr, " %10.3e", E0_hstr[ k ] );
2017  fprintf ( stderr, " %10.3e", E_hstr[ k ] );
2018  fprintf ( stderr, " %5d", n_CGiter_hstr[ k ] );
2019 
2020  if ( truss_exit_hstr[ k ] == true )
2021  fprintf ( stderr, " -Y- " );
2022 
2023  else
2024  fprintf ( stderr, " N " );
2025 
2026  if ( SOC_hstr[ k ] == true )
2027  fprintf ( stderr, "-Y- " );
2028 
2029  else
2030  fprintf ( stderr, " N " );
2031 
2032  fprintf ( stderr, " %3.1e", Delta_hstr[ k ] );
2033  fprintf ( stderr, " %9.2e", ared_hstr[ k ] );
2034  fprintf ( stderr, " %9.2e", pred_hstr[ k ] );
2035  fprintf ( stderr, " %7.1f", gamma_hstr[ k ] );
2036  std::cerr << std::endl;
2037  }
2038  }
2039 
2040  std::cerr << "\nNumber of accepted steps: " << iter << "\n";
2041 }
2042 } // Feel
2043 #endif
holds the statistics for SolverConstrained
Definition: solverconstrained.hpp:153
Stats(solver_type *__s, bool __c=false)
default constructor
Definition: solverconstrained.hpp:159
void clear()
clear all statistics
Definition: solverconstrained.hpp:1924
Polynomial< Poly, Type > dx(Polynomial< Poly, Type > const &p)
compute of a polynomial p
Definition: operations.hpp:63
Problem< Data > problem_type
problem data type
Definition: solverconstrained.hpp:53
bool optimize(double __x[_E_n])
optimizes the problem defined by problem_type
Definition: solverconstrained.hpp:1698
Expr< Val< typename mpl::if_< boost::is_arithmetic< ExprT1 >, mpl::identity< Cst< ExprT1 > >, mpl::identity< ExprT1 > >::type::type > > val(ExprT1 const &__e1)
precompute expression tensor
Definition: val.hpp:304
statistics_type const & stats() const
Definition: solverconstrained.hpp:242
number of multipliers
Definition: solverconstrained.hpp:62
size of the matrix
Definition: solverconstrained.hpp:61
problem_type::ad_0_type ad_0_type
automatic differentiation type of order one
Definition: solverconstrained.hpp:67
Definition: solverconstrained.hpp:45
void push_all(double mu, double E, double E0, int n_CGiter, bool truss_exit, bool SOC, double Delta, double ared, double pred, double gamma)
insert all statistics
Definition: solverconstrained.hpp:1951
problem_type & problem()
Definition: solverconstrained.hpp:236
ad_0_type::value_type value_type
numerical type
Definition: solverconstrained.hpp:76
bool collectStats() const
tells if we collect statistics regarding the solver convergence
Definition: solverconstrained.hpp:168
void show(bool verbose=false) const
show statistics
Definition: solverconstrained.hpp:1972
void redefine_problem(double x_definitions[_E_n][3])
redefined the control variables bounds and if they are avtive
Definition: solverconstrained.hpp:224
SolverConstrained< Data, Problem > solver_type
this solver type
Definition: solverconstrained.hpp:50
~SolverConstrained()
destructor
Definition: solverconstrained.hpp:384
number of objective functions
Definition: solverconstrained.hpp:58
number of equality constraints
Definition: solverconstrained.hpp:60
problem_type::ad_1_type ad_1_type
automatic differentiation type of order one
Definition: solverconstrained.hpp:70
COptions options_type
option data type
Definition: solverconstrained.hpp:147
holds the options for the SolverConstrained
Definition: solverconstrained.hpp:87
problem_type::ad_2_type ad_2_type
automatic differentiation type of order two
Definition: solverconstrained.hpp:73
size of the multipliers matrix
Definition: solverconstrained.hpp:63
number of inequality constraints
Definition: solverconstrained.hpp:59
Stats statistics_type
statistics data type
Definition: solverconstrained.hpp:210
number of control variables
Definition: solverconstrained.hpp:57
void push_x(const double *__x)
insert x in statistics
Definition: solverconstrained.hpp:1943

Generated on Sun Dec 22 2013 13:11:13 for Feel++ by doxygen 1.8.5