MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_ReducedHessianSecantUpdateLPBFGS_Strategy.cpp
Go to the documentation of this file.
00001 #if 0
00002 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00007 //                  Copyright (2003) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include "MoochoPack_ReducedHessianSecantUpdateLPBFGS_Strategy.hpp"
00045 #include "MoochoPack_PBFGS_helpers.hpp"
00046 #include "MoochoPack_NLPAlgo.hpp"
00047 #include "MoochoPack_NLPAlgoState.hpp"
00048 #include "ConstrainedOptPack_MatrixSymPosDefLBFGS.hpp"
00049 #include "ConstrainedOptPack/src/AbstractLinAlgPack_MatrixSymPosDefCholFactor.hpp"
00050 #include "ConstrainedOptPack/src/AbstractLinAlgPack_BFGS_helpers.hpp"
00051 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00052 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00053 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixOpOut.hpp"
00054 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_GenPermMatrixSlice.hpp"
00055 #include "AbstractLinAlgPack_GenPermMatrixSliceOp.hpp"
00056 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymInitDiag.hpp"
00057 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00058 #include "Midynamic_cast_verbose.h"
00059 #include "MiWorkspacePack.h"
00060 
00061 namespace LinAlgOpPack {
00062   using AbstractLinAlgPack::Vp_StMtV;
00063 }
00064 
00065 namespace MoochoPack {
00066 
00067 ReducedHessianSecantUpdateLPBFGS_Strategy::ReducedHessianSecantUpdateLPBFGS_Strategy(
00068   const proj_bfgs_updater_ptr_t&  proj_bfgs_updater
00069   ,size_type                      min_num_updates_proj_start
00070   ,size_type                      max_num_updates_proj_start
00071   ,size_type                      num_superbasics_switch_dense
00072   ,size_type                      num_add_recent_updates
00073   )
00074   : proj_bfgs_updater_(proj_bfgs_updater)
00075   , min_num_updates_proj_start_(min_num_updates_proj_start)
00076   , max_num_updates_proj_start_(max_num_updates_proj_start)
00077   , num_superbasics_switch_dense_(num_superbasics_switch_dense)
00078   , num_add_recent_updates_(num_add_recent_updates)
00079 {}
00080 
00081 bool ReducedHessianSecantUpdateLPBFGS_Strategy::perform_update(
00082   DVectorSlice* s_bfgs, DVectorSlice* y_bfgs, bool first_update
00083   ,std::ostream& out, EJournalOutputLevel olevel, NLPAlgo *algo, NLPAlgoState *s
00084   ,MatrixOp *rHL_k
00085   )
00086 {
00087   using std::setw;
00088   using std::endl;
00089   using std::right;
00090   using Teuchos::dyn_cast;
00091   namespace rcp = MemMngPack;
00092   using Teuchos::RCP;
00093   using LinAlgOpPack::V_MtV;
00094   using DenseLinAlgPack::dot;
00095   using AbstractLinAlgPack::norm_inf;
00096   using AbstractLinAlgPack::transVtMtV;
00097   typedef ConstrainedOptPack::MatrixHessianSuperBasic MHSB_t;
00098   using Teuchos::Workspace;
00099   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00100 
00101   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00102     out << "\n*** (LPBFGS) Full limited memory BFGS to projected BFGS ...\n";
00103   }
00104 
00105 #ifdef _WINDOWS
00106   MHSB_t &rHL_super = dynamic_cast<MHSB_t&>(*rHL_k);
00107 #else
00108   MHSB_t &rHL_super = dyn_cast<MHSB_t>(*rHL_k);
00109 #endif
00110 
00111   const size_type
00112     n    = algo->nlp().n(),
00113     r    = algo->nlp().r(),
00114     n_pz = n-r;
00115 
00116   bool do_projected_rHL_RR = false;
00117 
00118   // See if we still have a limited memory BFGS update matrix
00119   RCP<MatrixSymPosDefLBFGS> // We don't want this to be deleted until we are done with it
00120     lbfgs_rHL_RR = Teuchos::rcp_const_cast<MatrixSymPosDefLBFGS>(
00121       Teuchos::rcp_dynamic_cast<const MatrixSymPosDefLBFGS>(rHL_super.B_RR_ptr()) );
00122 
00123   if( lbfgs_rHL_RR.get() && rHL_super.Q_R().is_identity()  ) {
00124     //
00125     // We have a limited memory BFGS matrix and have not started projected BFGS updating
00126     // yet so lets determine if it is time to consider switching
00127     //
00128     // Check that the current update is sufficiently p.d. before we do anything
00129     const value_type
00130       sTy = dot(*s_bfgs,*y_bfgs),
00131       yTy = dot(*y_bfgs,*y_bfgs);
00132     if( !ConstrainedOptPack::BFGS_sTy_suff_p_d(
00133       *s_bfgs,*y_bfgs,&sTy
00134       ,  int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL )
00135       && !proj_bfgs_updater().bfgs_update().use_dampening()
00136       )
00137     {
00138       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00139         out << "\nWarning!  use_damening == false so there is no way we can perform any"
00140             " kind of BFGS update (projected or not) so we will skip it!\n";
00141       }
00142       quasi_newton_stats_(*s).set_k(0).set_updated_stats(
00143         QuasiNewtonStats::INDEF_SKIPED );
00144       return true;
00145     }
00146     // Consider if we can even look at the active set yet.
00147     const bool consider_switch  = lbfgs_rHL_RR->num_secant_updates() >= min_num_updates_proj_start();
00148     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00149       out << "\nnum_previous_updates = " << lbfgs_rHL_RR->num_secant_updates()
00150         << ( consider_switch ? " >= " : " < " )
00151         << "min_num_updates_proj_start = " << min_num_updates_proj_start()
00152         << ( consider_switch
00153            ? "\nConsidering if we should switch to projected BFGS updating of superbasics ...\n"
00154            : "\nNot time to consider switching to projected BFGS updating of superbasics yet!" );
00155     }
00156     if( consider_switch ) {
00157       // 
00158       // Our job here is to determine if it is time to switch to projected projected
00159       // BFGS updating.
00160       //
00161       if( act_set_stats_(*s).updated_k(-1) ) {
00162         if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00163           out << "\nDetermining if projected BFGS updating of superbasics should begin ...\n";
00164         }
00165         // Determine if the active set has calmed down enough
00166         const SpVector
00167           &nu_km1 = s->nu().get_k(-1);
00168         const SpVectorSlice
00169           nu_indep = nu_km1(s->var_indep());
00170         const bool 
00171           act_set_calmed_down
00172           = PBFGSPack::act_set_calmed_down(
00173             act_set_stats_(*s).get_k(-1)
00174             ,proj_bfgs_updater().act_set_frac_proj_start()
00175             ,olevel,out
00176             ),
00177           max_num_updates_exceeded
00178           = lbfgs_rHL_RR->m_bar() >= max_num_updates_proj_start();
00179         do_projected_rHL_RR = act_set_calmed_down || max_num_updates_exceeded;
00180         if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00181           if( act_set_calmed_down ) {
00182             out << "\nThe active set has calmed down enough so lets further consider switching to\n"
00183               << "projected BFGS updating of superbasic variables ...\n";
00184           }
00185           else if( max_num_updates_exceeded ) {
00186             out << "\nThe active set has not calmed down enough but num_previous_updates = "
00187               << lbfgs_rHL_RR->m_bar() << " >= max_num_updates_proj_start = " << max_num_updates_proj_start()
00188               << "\nso we will further consider switching to projected BFGS updating of superbasic variables ...\n";
00189           }
00190           else {
00191             out << "\nIt is not time to switch to projected BFGS so just keep performing full limited memory BFGS for now ...\n";
00192           }
00193         }
00194         if( do_projected_rHL_RR ) {
00195           //
00196           // Determine the set of initially fixed and free independent variables.
00197           //
00198           typedef Workspace<size_type>                              i_x_t;
00199           typedef Workspace<ConstrainedOptPack::EBounds>   bnd_t;
00200           i_x_t   i_x_free(wss,n_pz);
00201           i_x_t   i_x_fixed(wss,n_pz);
00202           bnd_t   bnd_fixed(wss,n_pz);
00203           i_x_t   l_x_fixed_sorted(wss,n_pz);
00204           size_type n_pz_X = 0, n_pz_R = 0;
00205           // rHL = rHL_scale * I
00206           value_type
00207             rHL_scale = sTy > 0.0 ? yTy/sTy : 1.0;
00208           if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00209             out << "\nScaling for diagonal intitial rHL = rHL_scale*I, rHL_scale = " << rHL_scale << std::endl;
00210           }
00211           value_type sRTBRRsR = 0.0, sRTyR = 0.0, sXTBXXsX = 0.0, sXTyX = 0.0;
00212           // Get the elements in i_x_free[] for variables that are definitely free
00213           // and initialize s_R'*B_RR*s_R and s_R'*y_R
00214           PBFGSPack::init_i_x_free_sRTsR_sRTyR(
00215             nu_indep, *s_bfgs, *y_bfgs
00216             , &n_pz_R, &i_x_free[0], &sRTBRRsR, &sRTyR );
00217           sRTBRRsR *= rHL_scale;
00218           Workspace<value_type> rHL_XX_diag_ws(wss,nu_indep.nz());
00219           DVectorSlice rHL_XX_diag(&rHL_XX_diag_ws[0],rHL_XX_diag_ws.size());
00220           rHL_XX_diag = rHL_scale;
00221           // Sort fixed variables according to |s_X(i)^2*B_XX(i,i)|/|sRTBRRsR| + |s_X(i)*y_X(i)|/|sRTyR|
00222           // and initialize s_X'*B_XX*s_X and s_X*y_X
00223           PBFGSPack::sort_fixed_max_cond_viol(
00224             nu_indep,*s_bfgs,*y_bfgs,rHL_XX_diag,sRTBRRsR,sRTyR
00225             ,&sXTBXXsX,&sXTyX,&l_x_fixed_sorted[0]);
00226           // Pick initial set of i_x_free[] and i_x_fixed[] (sorted!)
00227           PBFGSPack::choose_fixed_free(
00228             proj_bfgs_updater().project_error_tol()
00229             ,proj_bfgs_updater().super_basic_mult_drop_tol(),nu_indep
00230             ,*s_bfgs,*y_bfgs,rHL_XX_diag,&l_x_fixed_sorted[0]
00231             ,olevel,out,&sRTBRRsR,&sRTyR,&sXTBXXsX,&sXTyX
00232             ,&n_pz_X,&n_pz_R,&i_x_free[0],&i_x_fixed[0],&bnd_fixed[0] );
00233           if( n_pz_R < n_pz ) {
00234             //
00235             // We are ready to transition to projected BFGS updating!
00236             //
00237             // Determine if we are be using dense or limited memory BFGS?
00238             const bool
00239               low_num_super_basics = n_pz_R <= num_superbasics_switch_dense();
00240             if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ) {
00241               out << "\nSwitching to projected BFGS updating ..."
00242                 << "\nn_pz_R = " << n_pz_R << ( low_num_super_basics ? " <= " : " > " )
00243                 << " num_superbasics_switch_dense = " << num_superbasics_switch_dense()
00244                 << ( low_num_super_basics
00245                    ? "\nThere are not too many superbasic variables so use dense projected BFGS ...\n"
00246                    : "\nThere are too many superbasic variables so use limited memory projected BFGS ...\n"
00247                   );
00248             }
00249             // Create new matrix to use for rHL_RR initialized to rHL_RR = rHL_scale*I
00250             RCP<MatrixSymSecant>
00251               rHL_RR = NULL;
00252             if( low_num_super_basics ) {
00253               rHL_RR = new MatrixSymPosDefCholFactor(
00254                 NULL    // Let it allocate its own memory
00255                 ,NULL   // ...
00256                 ,n_pz   // maximum size
00257                 ,lbfgs_rHL_RR->maintain_original()
00258                 ,lbfgs_rHL_RR->maintain_inverse()
00259                 );
00260             }
00261             else {
00262               rHL_RR = new MatrixSymPosDefLBFGS(
00263                 n_pz, lbfgs_rHL_RR->m()
00264                 ,lbfgs_rHL_RR->maintain_original()
00265                 ,lbfgs_rHL_RR->maintain_inverse()
00266                 ,lbfgs_rHL_RR->auto_rescaling()
00267                 );
00268             }
00269             rHL_RR->init_identity( n_pz_R, rHL_scale );
00270             if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00271               out << "\nrHL_RR after intialized to rHL_RR = rHL_scale*I but before performing current BFGS update:\nrHL_RR =\n"
00272                 << dynamic_cast<MatrixOp&>(*rHL_RR); // I know this is okay!
00273             }
00274             // Initialize rHL_XX = rHL_scale*I
00275 #ifdef _WINDOWS
00276             MatrixSymInitDiag
00277               &rHL_XX = dynamic_cast<MatrixSymInitDiag&>(
00278                 const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr()));
00279 #else
00280             MatrixSymInitDiag
00281               &rHL_XX = dyn_cast<MatrixSymInitDiag>(
00282                 const_cast<MatrixSymOp&>(*rHL_super.B_XX_ptr()));
00283 #endif
00284             rHL_XX.init_identity( n_pz_X, rHL_scale );
00285             // Reinitialize rHL
00286             rHL_super.initialize(
00287               n_pz, n_pz_R, &i_x_free[0], &i_x_fixed[0], &bnd_fixed[0]
00288               ,Teuchos::rcp_const_cast<const MatrixSymWithOpFactorized>(
00289                 Teuchos::rcp_dynamic_cast<MatrixSymWithOpFactorized>(rHL_RR))
00290               ,NULL,BLAS_Cpp::no_trans,rHL_super.B_XX_ptr()
00291               );
00292             //
00293             // Perform the current BFGS update first
00294             //
00295             MatrixSymOp
00296               &rHL_RR_op = dynamic_cast<MatrixSymOp&>(*rHL_RR);
00297             const GenPermMatrixSlice
00298               &Q_R = rHL_super.Q_R(),
00299               &Q_X = rHL_super.Q_X();
00300             // Get projected BFGS update vectors y_bfgs_R, s_bfgs_R
00301             Workspace<value_type>
00302               y_bfgs_R_ws(wss,Q_R.cols()),
00303               s_bfgs_R_ws(wss,Q_R.cols()),
00304               y_bfgs_X_ws(wss,Q_X.cols()),
00305               s_bfgs_X_ws(wss,Q_X.cols());
00306             DVectorSlice y_bfgs_R(&y_bfgs_R_ws[0],y_bfgs_R_ws.size());
00307             DVectorSlice s_bfgs_R(&s_bfgs_R_ws[0],s_bfgs_R_ws.size());
00308             DVectorSlice y_bfgs_X(&y_bfgs_X_ws[0],y_bfgs_X_ws.size());
00309             DVectorSlice s_bfgs_X(&s_bfgs_X_ws[0],s_bfgs_X_ws.size());
00310             V_MtV( &y_bfgs_R, Q_R, BLAS_Cpp::trans, *y_bfgs );  // y_bfgs_R = Q_R'*y_bfgs
00311             V_MtV( &s_bfgs_R, Q_R, BLAS_Cpp::trans, *s_bfgs );  // s_bfgs_R = Q_R'*s_bfgs
00312             V_MtV( &y_bfgs_X, Q_X, BLAS_Cpp::trans, *y_bfgs );  // y_bfgs_X = Q_X'*y_bfgs
00313             V_MtV( &s_bfgs_X, Q_X, BLAS_Cpp::trans, *s_bfgs );  // s_bfgs_X = Q_X'*s_bfgs
00314             // Update rHL_RR
00315             if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00316               out << "\nPerform current BFGS update on " << n_pz_R << " x " << n_pz_R << " projected "
00317                 << "reduced Hessian for the superbasic variables where B = rHL_RR...\n";
00318             }
00319             QuasiNewtonStats quasi_newton_stats;
00320             proj_bfgs_updater().bfgs_update().perform_update(
00321               &s_bfgs_R(),&y_bfgs_R(),false,out,olevel,algo->algo_cntr().check_results()
00322               ,&rHL_RR_op, &quasi_newton_stats );
00323             // Perform previous updates if possible
00324             if( lbfgs_rHL_RR->m_bar() ) {
00325               const size_type num_add_updates = std::_MIN(num_add_recent_updates(),lbfgs_rHL_RR->m_bar());
00326               if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00327                 out << "\nAdd the min(num_previous_updates,num_add_recent_updates) = min(" << lbfgs_rHL_RR->m_bar()
00328                   << "," << num_add_recent_updates() << ") = " << num_add_updates << " most recent BFGS updates if possible ...\n";
00329               }
00330               // Now add previous update vectors
00331               const value_type
00332                 project_error_tol = proj_bfgs_updater().project_error_tol();
00333               const DMatrixSlice
00334                 S = lbfgs_rHL_RR->S(),
00335                 Y = lbfgs_rHL_RR->Y();
00336               size_type k = lbfgs_rHL_RR->k_bar();  // Location in S and Y of most recent update vectors
00337               for( size_type l = 1; l <= num_add_updates; ++l, --k ) {
00338                 if(k == 0) k = lbfgs_rHL_RR->m_bar();  // see MatrixSymPosDefLBFGS
00339                 // Check to see if this update satisfies the required conditions.
00340                 // Start with the condition on s'*y since this are cheap to check.
00341                 V_MtV( &s_bfgs_X(), Q_X, BLAS_Cpp::trans, S.col(k) ); // s_bfgs_X = Q_X'*s_bfgs
00342                 V_MtV( &y_bfgs_X(), Q_X, BLAS_Cpp::trans, Y.col(k) ); // y_bfgs_X = Q_X'*y_bfgs
00343                 sRTyR    = dot( s_bfgs_R, y_bfgs_R );
00344                 sXTyX    = dot( s_bfgs_X, y_bfgs_X );
00345                 bool
00346                   sXTyX_cond    = ::fabs(sXTyX/sRTyR) <= project_error_tol,
00347                   do_update     = sXTyX_cond,
00348                   sXTBXXsX_cond = false;
00349                 if( sXTyX_cond ) {
00350                   // Check the second more expensive condition
00351                   V_MtV( &s_bfgs_R(), Q_R, BLAS_Cpp::trans, S.col(k) ); // s_bfgs_R = Q_R'*s_bfgs
00352                   V_MtV( &y_bfgs_R(), Q_R, BLAS_Cpp::trans, Y.col(k) ); // y_bfgs_R = Q_R'*y_bfgs
00353                   sRTBRRsR = transVtMtV( s_bfgs_R, rHL_RR_op, BLAS_Cpp::no_trans, s_bfgs_R );
00354                   sXTBXXsX = rHL_scale * dot( s_bfgs_X, s_bfgs_X );
00355                   sXTBXXsX_cond = sXTBXXsX/sRTBRRsR <= project_error_tol;
00356                   do_update     = sXTBXXsX_cond && sXTyX_cond;
00357                 }
00358                 if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00359                   out << "\n---------------------------------------------------------------------"
00360                     << "\nprevious update " << l
00361                     << "\n\nChecking projection error:\n"
00362                     << "\n|s_X'*y_X|/|s_R'*y_R| = |" << sXTyX << "|/|" << sRTyR
00363                     << "| = " << ::fabs(sXTyX/sRTyR)
00364                     << ( sXTyX_cond ? " <= " : " > " ) << " project_error_tol = "
00365                     << project_error_tol;
00366                   if( sXTyX_cond ) {
00367                     out << "\n(s_X'*rHL_XX*s_X/s_R'*rHL_RR*s_R) = (" << sXTBXXsX
00368                         << ") = " << (sXTBXXsX/sRTBRRsR)
00369                         << ( sXTBXXsX_cond ? " <= " : " > " ) << " project_error_tol = "
00370                         << project_error_tol;
00371                   }
00372                   out << ( do_update
00373                     ? "\n\nAttemping to add this previous update where B = rHL_RR ...\n"
00374                     : "\n\nCan not add this previous update ...\n" );
00375                 }
00376                 if( do_update ) {
00377                   // ( rHL_RR, s_bfgs_R, y_bfgs_R ) -> rHL_RR (this should not throw an exception!)
00378                   try {
00379                     proj_bfgs_updater().bfgs_update().perform_update(
00380                       &s_bfgs_R(),&y_bfgs_R(),false,out,olevel,algo->algo_cntr().check_results()
00381                       ,&rHL_RR_op, &quasi_newton_stats );
00382                   }
00383                   catch( const MatrixSymSecant::UpdateSkippedException& excpt ) {
00384                     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00385                       out << "\nOops!  The " << l << "th most recent BFGS update was rejected?:\n"
00386                         << excpt.what() << std::endl;
00387                     }
00388                   }
00389                 }
00390               }
00391               if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00392                 out << "\n---------------------------------------------------------------------\n";
00393                 }
00394               if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00395                 out << "\nrHL_RR after adding previous BFGS updates:\nrHL_BRR =\n"
00396                   << dynamic_cast<MatrixOp&>(*rHL_RR);
00397               }
00398             }
00399             else {
00400               if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00401                 out << "\nThere were no previous update vectors to add!\n";
00402               }
00403             }
00404             if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00405               out << "\nFull rHL after complete reinitialization:\nrHL =\n" << *rHL_k;
00406             }
00407             quasi_newton_stats_(*s).set_k(0).set_updated_stats(
00408               QuasiNewtonStats::REINITIALIZED );
00409             return true;
00410           }
00411           else {
00412             if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00413               out << "\nn_pz_R == n_pz = " << n_pz_R << ", No variables would be removed from "
00414                 << "the superbasis so just keep on performing limited memory BFGS for now ...\n";
00415             }
00416             do_projected_rHL_RR = false;
00417           }
00418         }
00419       }
00420     }
00421     // If we have not switched to PBFGS then just update the full limited memory BFGS matrix
00422     if(!do_projected_rHL_RR) {
00423       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00424         out << "\nPerform current BFGS update on " << n_pz << " x " << n_pz << " full "
00425           << "limited memory reduced Hessian B = rHL ...\n";
00426       }
00427       proj_bfgs_updater().bfgs_update().perform_update(
00428         s_bfgs,y_bfgs,first_update,out,olevel,algo->algo_cntr().check_results()
00429         ,lbfgs_rHL_RR.get()
00430         ,&quasi_newton_stats_(*s).set_k(0)
00431         );
00432       return true;
00433     }
00434   }
00435   else {
00436     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00437       out << "\nWe have already switched to projected BFGS updating ...\n";
00438     }
00439   }
00440   //
00441   // If we get here then we must have switched to
00442   // projected updating so lets just pass it on!
00443   //
00444   return proj_bfgs_updater().perform_update(
00445     s_bfgs,y_bfgs,first_update,out,olevel,algo,s,rHL_k);
00446 }
00447 
00448 void ReducedHessianSecantUpdateLPBFGS_Strategy::print_step( std::ostream& out, const std::string& L ) const
00449 {
00450   out
00451     << L << "*** Perform limited memory LBFGS updating initially then switch to dense\n"
00452     << L << "*** projected BFGS updating when appropriate.\n"
00453     << L << "ToDo: Finish implementation!\n";
00454 }
00455 
00456 }  // end namespace MoochoPack
00457 
00458 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines