MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MoochoPack_CheckDescentQuasiNormalStep_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 <ostream>
43 
54 #include "Teuchos_Assert.hpp"
55 
56 namespace MoochoPack {
57 
59  const calc_fd_prod_ptr_t& calc_fd_prod
60  )
61  :calc_fd_prod_(calc_fd_prod)
62 {}
63 
65  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
66  ,poss_type assoc_step_poss
67  )
68 {
69  using BLAS_Cpp::no_trans;
71  using LinAlgOpPack::V_MtV;
72 
73  NLPAlgo &algo = rsqp_algo(_algo);
74  NLPAlgoState &s = algo.rsqp_state();
75  NLP &nlp = algo.nlp();
76  const Range1D equ_decomp = s.equ_decomp();
77 
78  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
79  std::ostream& out = algo.track().journal_out();
80 
81  // print step header.
82  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
84  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
85  }
86 
87  const size_type
88  nb = nlp.num_bounded_x();
89 
90  // Get iteration quantities
91  IterQuantityAccess<VectorMutable>
92  &c_iq = s.c(),
93  &Ypy_iq = s.Ypy();
94  const Vector::vec_ptr_t
95  cd_k = c_iq.get_k(0).sub_view(equ_decomp);
96  const Vector
97  &Ypy_k = Ypy_iq.get_k(0);
98 
99  value_type descent_c = -1.0;
100  if( s.get_iter_quant_id( Gc_name ) != AlgorithmState::DOES_NOT_EXIST ) {
101  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
102  out << "\nGc_k exists; compute descent_c = c_k(equ_decomp)'*Gc_k(:,equ_decomp)'*Ypy_k ...\n";
103  }
104  const MatrixOp::mat_ptr_t
105  Gcd_k = s.Gc().get_k(0).sub_view(Range1D(),equ_decomp);
106  VectorSpace::vec_mut_ptr_t
107  t = cd_k->space().create_member();
108  V_MtV( t.get(), *Gcd_k, BLAS_Cpp::trans, Ypy_k );
109  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
110  out << "\nGc_k(:,equ_decomp)'*Ypy_k =\n" << *t;
111  }
112  descent_c = dot( *cd_k, *t );
113  }
114  else {
115  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
116  out << "\nGc_k does not exist; compute descent_c = c_k(equ_decomp)'*FDGc_k(:,equ_decomp)'*Ypy_k "
117  << "using finite differences ...\n";
118  }
119  VectorSpace::vec_mut_ptr_t
120  t = nlp.space_c()->create_member();
121  calc_fd_prod().calc_deriv_product(
122  s.x().get_k(0),nb?&nlp.xl():NULL,nb?&nlp.xu():NULL
123  ,Ypy_k,NULL,&c_iq.get_k(0),true,&nlp
124  ,NULL,t.get()
125  ,static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ? &out : NULL
126  );
127  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
128  out << "\nFDGc_k(:,equ_decomp)'*Ypy_k =\n" << *t->sub_view(equ_decomp);
129  }
130  descent_c = dot( *cd_k, *t->sub_view(equ_decomp) );
131  }
132 
133  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
134  out << "\ndescent_c = " << descent_c << std::endl;
135  }
136 
137  if( descent_c > 0.0 ) { // ToDo: add some allowance for > 0.0 for finite difference errors!
138  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
139  out << "\nError, descent_c > 0.0; this is not a descent direction\n"
140  << "Throw TestFailed and terminate the algorithm ...\n";
141  }
143  true, TestFailed
144  ,"CheckDescentQuasiNormalStep_Step::do_step(...) : Error, descent for the decomposed constraints "
145  "with respect to the quasi-normal step c_k(equ_decomp)'*FDGc_k(:,equ_decomp)'*Ypy_k = "
146  << descent_c << " > 0.0; This is not a descent direction!\n" );
147  }
148 
149  return true;
150 }
151 
153  const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
154  ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
155  ) const
156 {
157  out
158  << L << "*** Check for descent in the decomposed equality constraints for the quasi-normal step\n"
159  << L << "if Gc_k exists then\n"
160  << L << " descent_c = c_k(equ_decomp)'*Gc_k(:,equ_decomp)'*Ypy_k\n"
161  << L << "else\n"
162  << L << " descent_c = c_k(equ_decomp)'*FDGc(:,equ_decomp)'*Ypy_k (finite diff.)\n"
163  << L << "end\n"
164  << L << "if descent > 0.0 then\n"
165  << L << " throw TestFailed exception!\n"
166  << L << "end\n"
167  ;
168 }
169 
170 } // end namespace MoochoPack
AbstractLinAlgPack::size_type size_type
CheckDescentQuasiNormalStep_Step(const calc_fd_prod_ptr_t &calc_fd_prod)
Constructor.
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
Thrown if a runtime test failed.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
rSQP Algorithm control class.
const std::string Gc_name
Transposed.
Not transposed.
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.
AlgorithmTracker & track()
virtual iq_id_type get_iter_quant_id(const std::string &iq_name) const
Return the iteration quantity id (iq_id) for the iteration quantity.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
result = v_rhs1' * v_rhs2
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
AbstractLinAlgPack::value_type value_type
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
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.