MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_ReducedHessianSecantUpdateLPBFGS_Strategy.cpp
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include "MoochoPack_ReducedHessianSecantUpdateLPBFGS_Strategy.hpp"
45 #include "MoochoPack_PBFGS_helpers.hpp"
46 #include "MoochoPack_NLPAlgo.hpp"
47 #include "MoochoPack_NLPAlgoState.hpp"
48 #include "ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp"
49 #include "ConstrainedOptPack/src/AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp"
50 #include "ConstrainedOptPack/src/AbstractLinAlgPack_BFGS_helpers.hpp"
51 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
52 #include "AbstractLinAlgPack_SpVectorOp.hpp"
53 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOpOut.hpp"
54 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
55 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
56 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymInitDiag.hpp"
57 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
58 #include "Midynamic_cast_verbose.h"
59 #include "MiWorkspacePack.h"
60 
61 namespace LinAlgOpPack {
63 }
64 
65 namespace MoochoPack {
66 
68  const proj_bfgs_updater_ptr_t& proj_bfgs_updater
69  ,size_type min_num_updates_proj_start
70  ,size_type max_num_updates_proj_start
71  ,size_type num_superbasics_switch_dense
72  ,size_type num_add_recent_updates
73  )
74  : proj_bfgs_updater_(proj_bfgs_updater)
75  , min_num_updates_proj_start_(min_num_updates_proj_start)
76  , max_num_updates_proj_start_(max_num_updates_proj_start)
77  , num_superbasics_switch_dense_(num_superbasics_switch_dense)
78  , num_add_recent_updates_(num_add_recent_updates)
79 {}
80 
82  DVectorSlice* s_bfgs, DVectorSlice* y_bfgs, bool first_update
83  ,std::ostream& out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s
84  ,MatrixOp *rHL_k
85  )
86 {
87  using std::setw;
88  using std::endl;
89  using std::right;
90  using Teuchos::dyn_cast;
91  namespace rcp = MemMngPack;
92  using Teuchos::RCP;
93  using LinAlgOpPack::V_MtV;
95  using AbstractLinAlgPack::norm_inf;
98  using Teuchos::Workspace;
100 
101  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
102  out << "\n*** (LPBFGS) Full limited memory BFGS to projected BFGS ...\n";
103  }
104 
105 #ifdef _WINDOWS
106  MHSB_t &rHL_super = dynamic_cast<MHSB_t&>(*rHL_k);
107 #else
108  MHSB_t &rHL_super = dyn_cast<MHSB_t>(*rHL_k);
109 #endif
110 
111  const size_type
112  n = algo->nlp().n(),
113  r = algo->nlp().r(),
114  n_pz = n-r;
115 
116  bool do_projected_rHL_RR = false;
117 
118  // See if we still have a limited memory BFGS update matrix
119  RCP<MatrixSymPosDefLBFGS> // We don't want this to be deleted until we are done with it
120  lbfgs_rHL_RR = Teuchos::rcp_const_cast<MatrixSymPosDefLBFGS>(
121  Teuchos::rcp_dynamic_cast<const MatrixSymPosDefLBFGS>(rHL_super.B_RR_ptr()) );
122 
123  if( lbfgs_rHL_RR.get() && rHL_super.Q_R().is_identity() ) {
124  //
125  // We have a limited memory BFGS matrix and have not started projected BFGS updating
126  // yet so lets determine if it is time to consider switching
127  //
128  // Check that the current update is sufficiently p.d. before we do anything
129  const value_type
130  sTy = dot(*s_bfgs,*y_bfgs),
131  yTy = dot(*y_bfgs,*y_bfgs);
132  if( !ConstrainedOptPack::BFGS_sTy_suff_p_d(
133  *s_bfgs,*y_bfgs,&sTy
134  , int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL )
135  && !proj_bfgs_updater().bfgs_update().use_dampening()
136  )
137  {
138  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
139  out << "\nWarning! use_damening == false so there is no way we can perform any"
140  " kind of BFGS update (projected or not) so we will skip it!\n";
141  }
142  quasi_newton_stats_(*s).set_k(0).set_updated_stats(
143  QuasiNewtonStats::INDEF_SKIPED );
144  return true;
145  }
146  // Consider if we can even look at the active set yet.
147  const bool consider_switch = lbfgs_rHL_RR->num_secant_updates() >= min_num_updates_proj_start();
148  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
149  out << "\nnum_previous_updates = " << lbfgs_rHL_RR->num_secant_updates()
150  << ( consider_switch ? " >= " : " < " )
151  << "min_num_updates_proj_start = " << min_num_updates_proj_start()
152  << ( consider_switch
153  ? "\nConsidering if we should switch to projected BFGS updating of superbasics ...\n"
154  : "\nNot time to consider switching to projected BFGS updating of superbasics yet!" );
155  }
156  if( consider_switch ) {
157  //
158  // Our job here is to determine if it is time to switch to projected projected
159  // BFGS updating.
160  //
161  if( act_set_stats_(*s).updated_k(-1) ) {
162  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
163  out << "\nDetermining if projected BFGS updating of superbasics should begin ...\n";
164  }
165  // Determine if the active set has calmed down enough
166  const SpVector
167  &nu_km1 = s->nu().get_k(-1);
168  const SpVectorSlice
169  nu_indep = nu_km1(s->var_indep());
170  const bool
171  act_set_calmed_down
172  = PBFGSPack::act_set_calmed_down(
173  act_set_stats_(*s).get_k(-1)
174  ,proj_bfgs_updater().act_set_frac_proj_start()
175  ,olevel,out
176  ),
177  max_num_updates_exceeded
178  = lbfgs_rHL_RR->m_bar() >= max_num_updates_proj_start();
179  do_projected_rHL_RR = act_set_calmed_down || max_num_updates_exceeded;
180  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
181  if( act_set_calmed_down ) {
182  out << "\nThe active set has calmed down enough so lets further consider switching to\n"
183  << "projected BFGS updating of superbasic variables ...\n";
184  }
185  else if( max_num_updates_exceeded ) {
186  out << "\nThe active set has not calmed down enough but num_previous_updates = "
187  << lbfgs_rHL_RR->m_bar() << " >= max_num_updates_proj_start = " << max_num_updates_proj_start()
188  << "\nso we will further consider switching to projected BFGS updating of superbasic variables ...\n";
189  }
190  else {
191  out << "\nIt is not time to switch to projected BFGS so just keep performing full limited memory BFGS for now ...\n";
192  }
193  }
194  if( do_projected_rHL_RR ) {
195  //
196  // Determine the set of initially fixed and free independent variables.
197  //
198  typedef Workspace<size_type> i_x_t;
199  typedef Workspace<ConstrainedOptPack::EBounds> bnd_t;
200  i_x_t i_x_free(wss,n_pz);
201  i_x_t i_x_fixed(wss,n_pz);
202  bnd_t bnd_fixed(wss,n_pz);
203  i_x_t l_x_fixed_sorted(wss,n_pz);
204  size_type n_pz_X = 0, n_pz_R = 0;
205  // rHL = rHL_scale * I
206  value_type
207  rHL_scale = sTy > 0.0 ? yTy/sTy : 1.0;
208  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
209  out << "\nScaling for diagonal intitial rHL = rHL_scale*I, rHL_scale = " << rHL_scale << std::endl;
210  }
211  value_type sRTBRRsR = 0.0, sRTyR = 0.0, sXTBXXsX = 0.0, sXTyX = 0.0;
212  // Get the elements in i_x_free[] for variables that are definitely free
213  // and initialize s_R'*B_RR*s_R and s_R'*y_R
214  PBFGSPack::init_i_x_free_sRTsR_sRTyR(
215  nu_indep, *s_bfgs, *y_bfgs
216  , &n_pz_R, &i_x_free[0], &sRTBRRsR, &sRTyR );
217  sRTBRRsR *= rHL_scale;
218  Workspace<value_type> rHL_XX_diag_ws(wss,nu_indep.nz());
219  DVectorSlice rHL_XX_diag(&rHL_XX_diag_ws[0],rHL_XX_diag_ws.size());
220  rHL_XX_diag = rHL_scale;
221  // Sort fixed variables according to |s_X(i)^2*B_XX(i,i)|/|sRTBRRsR| + |s_X(i)*y_X(i)|/|sRTyR|
222  // and initialize s_X'*B_XX*s_X and s_X*y_X
223  PBFGSPack::sort_fixed_max_cond_viol(
224  nu_indep,*s_bfgs,*y_bfgs,rHL_XX_diag,sRTBRRsR,sRTyR
225  ,&sXTBXXsX,&sXTyX,&l_x_fixed_sorted[0]);
226  // Pick initial set of i_x_free[] and i_x_fixed[] (sorted!)
227  PBFGSPack::choose_fixed_free(
228  proj_bfgs_updater().project_error_tol()
229  ,proj_bfgs_updater().super_basic_mult_drop_tol(),nu_indep
230  ,*s_bfgs,*y_bfgs,rHL_XX_diag,&l_x_fixed_sorted[0]
231  ,olevel,out,&sRTBRRsR,&sRTyR,&sXTBXXsX,&sXTyX
232  ,&n_pz_X,&n_pz_R,&i_x_free[0],&i_x_fixed[0],&bnd_fixed[0] );
233  if( n_pz_R < n_pz ) {
234  //
235  // We are ready to transition to projected BFGS updating!
236  //
237  // Determine if we are be using dense or limited memory BFGS?
238  const bool
239  low_num_super_basics = n_pz_R <= num_superbasics_switch_dense();
240  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ) {
241  out << "\nSwitching to projected BFGS updating ..."
242  << "\nn_pz_R = " << n_pz_R << ( low_num_super_basics ? " <= " : " > " )
243  << " num_superbasics_switch_dense = " << num_superbasics_switch_dense()
244  << ( low_num_super_basics
245  ? "\nThere are not too many superbasic variables so use dense projected BFGS ...\n"
246  : "\nThere are too many superbasic variables so use limited memory projected BFGS ...\n"
247  );
248  }
249  // Create new matrix to use for rHL_RR initialized to rHL_RR = rHL_scale*I
250  RCP<MatrixSymSecant>
251  rHL_RR = NULL;
252  if( low_num_super_basics ) {
253  rHL_RR = new MatrixSymPosDefCholFactor(
254  NULL // Let it allocate its own memory
255  ,NULL // ...
256  ,n_pz // maximum size
257  ,lbfgs_rHL_RR->maintain_original()
258  ,lbfgs_rHL_RR->maintain_inverse()
259  );
260  }
261  else {
262  rHL_RR = new MatrixSymPosDefLBFGS(
263  n_pz, lbfgs_rHL_RR->m()
264  ,lbfgs_rHL_RR->maintain_original()
265  ,lbfgs_rHL_RR->maintain_inverse()
266  ,lbfgs_rHL_RR->auto_rescaling()
267  );
268  }
269  rHL_RR->init_identity( n_pz_R, rHL_scale );
270  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
271  out << "\nrHL_RR after intialized to rHL_RR = rHL_scale*I but before performing current BFGS update:\nrHL_RR =\n"
272  << dynamic_cast<MatrixOp&>(*rHL_RR); // I know this is okay!
273  }
274  // Initialize rHL_XX = rHL_scale*I
275 #ifdef _WINDOWS
276  MatrixSymInitDiag
277  &rHL_XX = dynamic_cast<MatrixSymInitDiag&>(
278  const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr()));
279 #else
280  MatrixSymInitDiag
281  &rHL_XX = dyn_cast<MatrixSymInitDiag>(
282  const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr()));
283 #endif
284  rHL_XX.init_identity( n_pz_X, rHL_scale );
285  // Reinitialize rHL
286  rHL_super.initialize(
287  n_pz, n_pz_R, &i_x_free[0], &i_x_fixed[0], &bnd_fixed[0]
288  ,Teuchos::rcp_const_cast<const MatrixSymWithOpFactorized>(
289  Teuchos::rcp_dynamic_cast<MatrixSymWithOpFactorized>(rHL_RR))
290  ,NULL,BLAS_Cpp::no_trans,rHL_super.B_XX_ptr()
291  );
292  //
293  // Perform the current BFGS update first
294  //
295  MatrixSymOp
296  &rHL_RR_op = dynamic_cast<MatrixSymOp&>(*rHL_RR);
297  const GenPermMatrixSlice
298  &Q_R = rHL_super.Q_R(),
299  &Q_X = rHL_super.Q_X();
300  // Get projected BFGS update vectors y_bfgs_R, s_bfgs_R
301  Workspace<value_type>
302  y_bfgs_R_ws(wss,Q_R.cols()),
303  s_bfgs_R_ws(wss,Q_R.cols()),
304  y_bfgs_X_ws(wss,Q_X.cols()),
305  s_bfgs_X_ws(wss,Q_X.cols());
306  DVectorSlice y_bfgs_R(&y_bfgs_R_ws[0],y_bfgs_R_ws.size());
307  DVectorSlice s_bfgs_R(&s_bfgs_R_ws[0],s_bfgs_R_ws.size());
308  DVectorSlice y_bfgs_X(&y_bfgs_X_ws[0],y_bfgs_X_ws.size());
309  DVectorSlice s_bfgs_X(&s_bfgs_X_ws[0],s_bfgs_X_ws.size());
310  V_MtV( &y_bfgs_R, Q_R, BLAS_Cpp::trans, *y_bfgs ); // y_bfgs_R = Q_R'*y_bfgs
311  V_MtV( &s_bfgs_R, Q_R, BLAS_Cpp::trans, *s_bfgs ); // s_bfgs_R = Q_R'*s_bfgs
312  V_MtV( &y_bfgs_X, Q_X, BLAS_Cpp::trans, *y_bfgs ); // y_bfgs_X = Q_X'*y_bfgs
313  V_MtV( &s_bfgs_X, Q_X, BLAS_Cpp::trans, *s_bfgs ); // s_bfgs_X = Q_X'*s_bfgs
314  // Update rHL_RR
315  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
316  out << "\nPerform current BFGS update on " << n_pz_R << " x " << n_pz_R << " projected "
317  << "reduced Hessian for the superbasic variables where B = rHL_RR...\n";
318  }
319  QuasiNewtonStats quasi_newton_stats;
320  proj_bfgs_updater().bfgs_update().perform_update(
321  &s_bfgs_R(),&y_bfgs_R(),false,out,olevel,algo->algo_cntr().check_results()
322  ,&rHL_RR_op, &quasi_newton_stats );
323  // Perform previous updates if possible
324  if( lbfgs_rHL_RR->m_bar() ) {
325  const size_type num_add_updates = std::_MIN(num_add_recent_updates(),lbfgs_rHL_RR->m_bar());
326  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
327  out << "\nAdd the min(num_previous_updates,num_add_recent_updates) = min(" << lbfgs_rHL_RR->m_bar()
328  << "," << num_add_recent_updates() << ") = " << num_add_updates << " most recent BFGS updates if possible ...\n";
329  }
330  // Now add previous update vectors
331  const value_type
332  project_error_tol = proj_bfgs_updater().project_error_tol();
333  const DMatrixSlice
334  S = lbfgs_rHL_RR->S(),
335  Y = lbfgs_rHL_RR->Y();
336  size_type k = lbfgs_rHL_RR->k_bar(); // Location in S and Y of most recent update vectors
337  for( size_type l = 1; l <= num_add_updates; ++l, --k ) {
338  if(k == 0) k = lbfgs_rHL_RR->m_bar(); // see MatrixSymPosDefLBFGS
339  // Check to see if this update satisfies the required conditions.
340  // Start with the condition on s'*y since this are cheap to check.
341  V_MtV( &s_bfgs_X(), Q_X, BLAS_Cpp::trans, S.col(k) ); // s_bfgs_X = Q_X'*s_bfgs
342  V_MtV( &y_bfgs_X(), Q_X, BLAS_Cpp::trans, Y.col(k) ); // y_bfgs_X = Q_X'*y_bfgs
343  sRTyR = dot( s_bfgs_R, y_bfgs_R );
344  sXTyX = dot( s_bfgs_X, y_bfgs_X );
345  bool
346  sXTyX_cond = ::fabs(sXTyX/sRTyR) <= project_error_tol,
347  do_update = sXTyX_cond,
348  sXTBXXsX_cond = false;
349  if( sXTyX_cond ) {
350  // Check the second more expensive condition
351  V_MtV( &s_bfgs_R(), Q_R, BLAS_Cpp::trans, S.col(k) ); // s_bfgs_R = Q_R'*s_bfgs
352  V_MtV( &y_bfgs_R(), Q_R, BLAS_Cpp::trans, Y.col(k) ); // y_bfgs_R = Q_R'*y_bfgs
353  sRTBRRsR = transVtMtV( s_bfgs_R, rHL_RR_op, BLAS_Cpp::no_trans, s_bfgs_R );
354  sXTBXXsX = rHL_scale * dot( s_bfgs_X, s_bfgs_X );
355  sXTBXXsX_cond = sXTBXXsX/sRTBRRsR <= project_error_tol;
356  do_update = sXTBXXsX_cond && sXTyX_cond;
357  }
358  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
359  out << "\n---------------------------------------------------------------------"
360  << "\nprevious update " << l
361  << "\n\nChecking projection error:\n"
362  << "\n|s_X'*y_X|/|s_R'*y_R| = |" << sXTyX << "|/|" << sRTyR
363  << "| = " << ::fabs(sXTyX/sRTyR)
364  << ( sXTyX_cond ? " <= " : " > " ) << " project_error_tol = "
365  << project_error_tol;
366  if( sXTyX_cond ) {
367  out << "\n(s_X'*rHL_XX*s_X/s_R'*rHL_RR*s_R) = (" << sXTBXXsX
368  << ") = " << (sXTBXXsX/sRTBRRsR)
369  << ( sXTBXXsX_cond ? " <= " : " > " ) << " project_error_tol = "
370  << project_error_tol;
371  }
372  out << ( do_update
373  ? "\n\nAttemping to add this previous update where B = rHL_RR ...\n"
374  : "\n\nCan not add this previous update ...\n" );
375  }
376  if( do_update ) {
377  // ( rHL_RR, s_bfgs_R, y_bfgs_R ) -> rHL_RR (this should not throw an exception!)
378  try {
379  proj_bfgs_updater().bfgs_update().perform_update(
380  &s_bfgs_R(),&y_bfgs_R(),false,out,olevel,algo->algo_cntr().check_results()
381  ,&rHL_RR_op, &quasi_newton_stats );
382  }
383  catch( const MatrixSymSecant::UpdateSkippedException& excpt ) {
384  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
385  out << "\nOops! The " << l << "th most recent BFGS update was rejected?:\n"
386  << excpt.what() << std::endl;
387  }
388  }
389  }
390  }
391  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
392  out << "\n---------------------------------------------------------------------\n";
393  }
394  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
395  out << "\nrHL_RR after adding previous BFGS updates:\nrHL_BRR =\n"
396  << dynamic_cast<MatrixOp&>(*rHL_RR);
397  }
398  }
399  else {
400  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
401  out << "\nThere were no previous update vectors to add!\n";
402  }
403  }
404  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
405  out << "\nFull rHL after complete reinitialization:\nrHL =\n" << *rHL_k;
406  }
407  quasi_newton_stats_(*s).set_k(0).set_updated_stats(
408  QuasiNewtonStats::REINITIALIZED );
409  return true;
410  }
411  else {
412  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
413  out << "\nn_pz_R == n_pz = " << n_pz_R << ", No variables would be removed from "
414  << "the superbasis so just keep on performing limited memory BFGS for now ...\n";
415  }
416  do_projected_rHL_RR = false;
417  }
418  }
419  }
420  }
421  // If we have not switched to PBFGS then just update the full limited memory BFGS matrix
422  if(!do_projected_rHL_RR) {
423  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
424  out << "\nPerform current BFGS update on " << n_pz << " x " << n_pz << " full "
425  << "limited memory reduced Hessian B = rHL ...\n";
426  }
427  proj_bfgs_updater().bfgs_update().perform_update(
428  s_bfgs,y_bfgs,first_update,out,olevel,algo->algo_cntr().check_results()
429  ,lbfgs_rHL_RR.get()
430  ,&quasi_newton_stats_(*s).set_k(0)
431  );
432  return true;
433  }
434  }
435  else {
436  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
437  out << "\nWe have already switched to projected BFGS updating ...\n";
438  }
439  }
440  //
441  // If we get here then we must have switched to
442  // projected updating so lets just pass it on!
443  //
444  return proj_bfgs_updater().perform_update(
445  s_bfgs,y_bfgs,first_update,out,olevel,algo,s,rHL_k);
446 }
447 
448 void ReducedHessianSecantUpdateLPBFGS_Strategy::print_step( std::ostream& out, const std::string& L ) const
449 {
450  out
451  << L << "*** Perform limited memory LBFGS updating initially then switch to dense\n"
452  << L << "*** projected BFGS updating when appropriate.\n"
453  << L << "ToDo: Finish implementation!\n";
454 }
455 
456 } // end namespace MoochoPack
457 
458 #endif // 0
bool perform_update(DVectorSlice *s_bfgs, DVectorSlice *y_bfgs, bool first_update, std::ostream &out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s, MatrixOp *rHL_k)
T_To & dyn_cast(T_From &from)
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
ReducedHessianSecantUpdateLPBFGS_Strategy(const proj_bfgs_updater_ptr_t &proj_bfgs_updater=NULL, size_type min_num_updates_proj_start=0, size_type max_num_updates_proj_start=999999, size_type num_superbasics_switch_dense=500, size_type num_add_recent_updates=10)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_t size_type
void Vp_StMtV(VectorMutable *v_lhs, value_type alpha, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2, value_type beta=1.0)
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
void print_step(std::ostream &out, const std::string &leading_str) const
TEUCHOSCORE_LIB_DLL_EXPORT Teuchos::RCP< WorkspaceStore > get_default_workspace_store()