MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MoochoPack_CheckSkipBFGSUpdateStd_Step.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <math.h>
43 
44 #include <ostream>
45 
51 
52 namespace MoochoPack {
53 
55  value_type skip_bfgs_prop_const
56  )
57  :skip_bfgs_prop_const_(skip_bfgs_prop_const)
58 {}
59 
61  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
62  , poss_type assoc_step_poss
63  )
64 {
65  NLPAlgo &algo = rsqp_algo(_algo);
66  NLPAlgoState &s = algo.rsqp_state();
67 
68  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
69  std::ostream& out = algo.track().journal_out();
70 
71  // print step header.
72  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
74  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
75  }
76 
77  IterQuantityAccess<MatrixSymOp>
78  &rHL = s.rHL();
79  if( rHL.updated_k(-1) ) {
80  bool skip_update = true;
81  if( !s.Ypy().updated_k(-1) || !s.Zpz().updated_k(-1)
82  || !s.rGL().updated_k(-1) || !s.c().updated_k(-1) )
83  {
84  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
85  out
86  << "\n*** Warning, rHL_km1 is updated but one or more of the quantities Ypy_km1, Zpz_km1\n"
87  ",rGL_km1 and c_km1 are not updated. Therefore, there is not sufficient information\n"
88  "to determine if to skip the BFGS update or not. Check storage requirements for\n"
89  " the above quantities\n"
90  "The BFGS wupdate will be skipped in this case ...\n";
91  }
92  skip_update = true;
93  }
94  else {
95  // The information exists for the update so determine
96  // if we are in the region to perform the BFGS update.
97  //
98  // Check if we are to skip the update for this iteration
99  const value_type
100  nrm_rGL_km1 = s.rGL().get_k(-1).norm_2(),
101  nrm_c_km1 = s.c().get_k(-1).norm_2(),
102  nrm_Zpz_km1 = s.Zpz().get_k(-1).norm_2(),
103  nrm_Ypy_km1 = s.Ypy().get_k(-1).norm_2();
104  // ratio = (skip_bfgs_prop_const / sqrt(||rGL_km1|| + ||c_km1||)) * ( ||Zpz_km1|| / ||Ypy_km1|| )
105  value_type
106  ratio = ( skip_bfgs_prop_const() / ::sqrt( nrm_rGL_km1 + nrm_c_km1 ) )
107  * ( nrm_Zpz_km1 / nrm_Ypy_km1 );
108  // If ratio < 1.0 then skip the update
109  skip_update = ratio < 1.0;
110  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
111  out
112  << std::endl
113  << "ratio = (skip_bfgs_prop_const/sqrt(||rGL_km1||+||c_km1||))*(||Zpz_km1||/||Ypy_km1||)\n"
114  << " = (" << skip_bfgs_prop_const() << "/sqrt("<<nrm_rGL_km1<<"+"<<nrm_c_km1<<"))\n"
115  << " * ("<<nrm_Zpz_km1<<"/"<<nrm_Ypy_km1<<")\n"
116  << " = " << ratio << std::endl
117  << "ratio " << (skip_update ? '<' : '>' ) << " 1\n"
118  << (skip_update
119  ? "Skipping BFGS update ...\n"
120  : "Perform BFGS update if you can ...\n" );
121  }
122  }
123  if( skip_update ) {
124  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS ) {
125  out
126  << "\nrHL_k = rHL_km1\n";
127  }
128  const MatrixSymOp &rHL_km1 = rHL.get_k(-1);
129  rHL.set_k(0) = rHL_km1;
130  quasi_newton_stats_(s).set_k(0).set_updated_stats(
132  if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
133  out << "\nrHL_k =\n" << rHL.get_k(0);
134  }
135  }
136  }
137 
138  return true;
139 }
140 
142  , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
143  , std::ostream& out, const std::string& L ) const
144 {
145  out
146  << L << "*** Check if we should do the BFGS update\n"
147  << L << "if rHL_km1 is update then\n"
148  << L << " If Ypy_km1, Zpz_km1, rGL_km1, or c_km1 is not updated then\n"
149  << L << " *** Warning, insufficient information to determine if we should\n"
150  << L << " *** perform the update. Check for sufficient backward storage.\n"
151  << L << " rHL_k = rHL_km1\n"
152  << L << " else\n"
153  << L << " *** Check if we are in the proper region\n"
154  << L << " ratio = (skip_bfgs_prop_const/sqrt(norm(rGL_km1,2)+norm(c_km1,2)))\n"
155  << L << " * (norm(Zpz_km1,2)/norm(Ypy_km1,2) )\n"
156  << L << " if ratio < 1 then \n"
157  << L << " rHL_k = rHL_km1\n"
158  << L << " end\n"
159  << L << " end\n"
160  << L << "end\n";
161 }
162 
163 } // end namespace MoochoPack
rSQP Algorithm control class.
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
virtual std::ostream & journal_out() const
Return a reference to a std::ostream to be used to output debug information and the like...
EJournalOutputLevel
enum for journal output.
Reduced space SQP state encapsulation interface.
std::ostream * out
void print_algorithm_step(const Algorithm &algo, Algorithm::poss_type step_poss, EDoStepType type, Algorithm::poss_type assoc_step_poss, std::ostream &out)
Prints to 'out' the algorithm step.
CheckSkipBFGSUpdateStd_Step(value_type skip_bfgs_prop_const=10.0)
AlgorithmTracker & track()
void sqrt(DVectorSlice *vs_lhs, const DVectorSlice &vs_rhs)
vs_lhs = sqrt(vs_rhs)
void print_step(const Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
AbstractLinAlgPack::value_type value_type
Acts as the central hub for an iterative algorithm.
NLPAlgoState & rsqp_state()
<<std aggr>="">> members for algo_cntr
NLPAlgo & rsqp_algo(Algorithm &algo)
Convert from a Algorithm to a NLPAlgo.