MoochoPack : Framework for Large-Scale Optimization Algorithms Version of the Day
MoochoPack_CheckSkipBFGSUpdateStd_Step.cpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <math.h>
00043 
00044 #include <ostream>
00045 
00046 #include "MoochoPack_CheckSkipBFGSUpdateStd_Step.hpp"
00047 #include "MoochoPack_moocho_algo_conversion.hpp"
00048 #include "IterationPack_print_algorithm_step.hpp"
00049 #include "AbstractLinAlgPack_MatrixSymOp.hpp"
00050 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00051 
00052 namespace MoochoPack {
00053 
00054 CheckSkipBFGSUpdateStd_Step::CheckSkipBFGSUpdateStd_Step(
00055   value_type  skip_bfgs_prop_const
00056   )
00057   :skip_bfgs_prop_const_(skip_bfgs_prop_const)
00058 {}
00059 
00060 bool CheckSkipBFGSUpdateStd_Step::do_step(
00061   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00062   , poss_type assoc_step_poss
00063   )
00064 {
00065   NLPAlgo &algo = rsqp_algo(_algo);
00066   NLPAlgoState  &s    = algo.rsqp_state();
00067 
00068   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00069   std::ostream& out = algo.track().journal_out();
00070 
00071   // print step header.
00072   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00073     using IterationPack::print_algorithm_step;
00074     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00075   }
00076 
00077   IterQuantityAccess<MatrixSymOp>
00078     &rHL = s.rHL();
00079   if( rHL.updated_k(-1) ) {
00080     bool skip_update = true;
00081     if( !s.Ypy().updated_k(-1) || !s.Zpz().updated_k(-1)
00082       || !s.rGL().updated_k(-1) || !s.c().updated_k(-1) )
00083     {
00084       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00085         out
00086           << "\n*** Warning, rHL_km1 is updated but one or more of the quantities Ypy_km1, Zpz_km1\n"
00087             ",rGL_km1 and c_km1 are not updated.  Therefore, there is not sufficient information\n"
00088             "to determine if to skip the BFGS update or not.  Check storage requirements for\n"
00089             " the above quantities\n"
00090             "The BFGS wupdate will be skipped in this case ...\n";
00091       } 
00092       skip_update = true;
00093     }
00094     else {
00095       // The information exists for the update so determine
00096       // if we are in the region to perform the BFGS update.
00097       //
00098       // Check if we are to skip the update for this iteration
00099       const value_type
00100         nrm_rGL_km1 = s.rGL().get_k(-1).norm_2(),
00101         nrm_c_km1 = s.c().get_k(-1).norm_2(),
00102         nrm_Zpz_km1 = s.Zpz().get_k(-1).norm_2(),
00103         nrm_Ypy_km1 = s.Ypy().get_k(-1).norm_2();
00104       // ratio = (skip_bfgs_prop_const / sqrt(||rGL_km1|| + ||c_km1||)) * ( ||Zpz_km1|| / ||Ypy_km1|| )
00105       value_type
00106         ratio = ( skip_bfgs_prop_const() / ::sqrt( nrm_rGL_km1 + nrm_c_km1 ) )
00107           * ( nrm_Zpz_km1 / nrm_Ypy_km1 );
00108       // If ratio < 1.0 then skip the update
00109       skip_update = ratio < 1.0;
00110       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00111         out
00112           << std::endl
00113           << "ratio = (skip_bfgs_prop_const/sqrt(||rGL_km1||+||c_km1||))*(||Zpz_km1||/||Ypy_km1||)\n"
00114           << "      = (" << skip_bfgs_prop_const() << "/sqrt("<<nrm_rGL_km1<<"+"<<nrm_c_km1<<"))\n"
00115           << "        * ("<<nrm_Zpz_km1<<"/"<<nrm_Ypy_km1<<")\n"
00116           << "      = " << ratio << std::endl
00117           << "ratio " << (skip_update ? '<' : '>' ) << " 1\n"
00118           << (skip_update
00119             ? "Skipping BFGS update ...\n"
00120             : "Perform BFGS update if you can ...\n"  );
00121       }
00122     }
00123     if( skip_update ) {
00124       if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
00125         out
00126           << "\nrHL_k = rHL_km1\n";
00127       }
00128       const MatrixSymOp &rHL_km1 = rHL.get_k(-1);
00129       rHL.set_k(0) = rHL_km1;
00130       quasi_newton_stats_(s).set_k(0).set_updated_stats(
00131         QuasiNewtonStats::SKIPED );
00132       if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
00133         out << "\nrHL_k =\n" << rHL.get_k(0);
00134       }
00135     }
00136   }
00137   
00138   return true;
00139 }
00140 
00141 void CheckSkipBFGSUpdateStd_Step::print_step( const Algorithm& algo
00142   , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
00143   , std::ostream& out, const std::string& L ) const
00144 {
00145   out
00146     << L << "*** Check if we should do the BFGS update\n"
00147     << L << "if rHL_km1 is update then\n"
00148     << L << "    If Ypy_km1, Zpz_km1, rGL_km1, or c_km1 is not updated then\n"
00149     << L << "        *** Warning, insufficient information to determine if we should\n"
00150     << L << "        *** perform the update.  Check for sufficient backward storage.\n"
00151     << L << "        rHL_k = rHL_km1\n"
00152     << L << "    else\n"
00153     << L << "        *** Check if we are in the proper region\n"
00154     << L << "        ratio = (skip_bfgs_prop_const/sqrt(norm(rGL_km1,2)+norm(c_km1,2)))\n"
00155     << L << "                 * (norm(Zpz_km1,2)/norm(Ypy_km1,2) )\n"
00156     << L << "        if ratio < 1 then \n"
00157     << L << "            rHL_k = rHL_km1\n"
00158     << L << "        end\n"
00159     << L << "    end\n"
00160     << L << "end\n";
00161 }
00162 
00163 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends