MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MoochoPack_TangentialStepIP_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 #include <iostream>
44 
48 #include "MoochoPack_IpState.hpp"
60 #include "Teuchos_dyn_cast.hpp"
61 #include "Teuchos_Assert.hpp"
62 
63 namespace MoochoPack {
64 
66  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
67  ,poss_type assoc_step_poss
68  )
69  {
70  using BLAS_Cpp::no_trans;
71  using Teuchos::dyn_cast;
73  using LinAlgOpPack::Vt_S;
75  using LinAlgOpPack::V_StV;
76  using LinAlgOpPack::V_MtV;
78  using LinAlgOpPack::M_StM;
81 
82  NLPAlgo &algo = rsqp_algo(_algo);
83  IpState &s = dyn_cast<IpState>(_algo.state());
84 
85  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
86  std::ostream& out = algo.track().journal_out();
87 
88  // print step header.
89  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
91  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
92  }
93 
94  // Compute qp_grad which is an approximation to rGf + Z'*(mu*(invXu*e-invXl*e) + no_cross_term
95  // minimize round off error by calc'ing Z'*(Gf + mu*(invXu*e-invXl*e))
96 
97  // qp_grad_k = Z'*(Gf + mu*(invXu*e-invXl*e))
98  const MatrixSymDiagStd &invXu = s.invXu().get_k(0);
99  const MatrixSymDiagStd &invXl = s.invXl().get_k(0);
100  const value_type &mu = s.barrier_parameter().get_k(0);
101  const MatrixOp &Z_k = s.Z().get_k(0);
102 
103  Teuchos::RCP<VectorMutable> rhs = s.Gf().get_k(0).clone();
104  Vp_StV( rhs.get(), mu, invXu.diag() );
105  Vp_StV( rhs.get(), -1.0*mu, invXl.diag() );
106 
107  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
108  {
109  out << "\n||Gf_k + mu_k*(invXu_k-invXl_k)||inf = " << rhs->norm_inf() << std::endl;
110  }
111  if( (int)olevel >= (int)PRINT_VECTORS)
112  {
113  out << "\nGf_k + mu_k*(invXu_k-invXl_k) =\n" << *rhs;
114  }
115 
116  VectorMutable &qp_grad_k = s.qp_grad().set_k(0);
117  V_MtV(&qp_grad_k, Z_k, BLAS_Cpp::trans, *rhs);
118 
119  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
120  {
121  out << "\n||Z_k'*(Gf_k + mu_k*(invXu_k-invXl_k))||inf = " << qp_grad_k.norm_inf() << std::endl;
122  }
123  if( (int)olevel >= (int)PRINT_VECTORS )
124  {
125  out << "\nZ_k'*(Gf_k + mu_k*(invXu_k-invXl_k)) =\n" << qp_grad_k;
126  }
127 
128  // error check for cross term
129  value_type &zeta = s.zeta().set_k(0);
130  const Vector &w_sigma = s.w_sigma().get_k(0);
131 
132  // need code to calculate damping parameter
133  zeta = 1.0;
134 
135  Vp_StV(&qp_grad_k, zeta, w_sigma);
136 
137  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
138  {
139  out << "\n||qp_grad_k||inf = " << qp_grad_k.norm_inf() << std::endl;
140  }
141  if( (int)olevel >= (int)PRINT_VECTORS )
142  {
143  out << "\nqp_grad_k =\n" << qp_grad_k;
144  }
145 
146  // build the "Hessian" term B = rHL + rHB
147  // should this be MatrixSymOpNonsing
148  const MatrixSymOp &rHL_k = s.rHL().get_k(0);
149  const MatrixSymOp &rHB_k = s.rHB().get_k(0);
150  MatrixSymOpNonsing &B_k = dyn_cast<MatrixSymOpNonsing>(s.B().set_k(0));
151  if (B_k.cols() != Z_k.cols())
152  {
153  // Initialize space in rHB
154  dyn_cast<MatrixSymInitDiag>(B_k).init_identity(Z_k.space_rows(), 0.0);
155  }
156 
157  // M_StM(&B_k, 1.0, rHL_k, no_trans);
158  assign(&B_k, rHL_k, BLAS_Cpp::no_trans);
159  if( (int)olevel >= (int)PRINT_VECTORS )
160  {
161  out << "\nB_k = rHL_k =\n" << B_k;
162  }
163  Mp_StM(&B_k, 1.0, rHB_k, BLAS_Cpp::no_trans);
164  if( (int)olevel >= (int)PRINT_VECTORS )
165  {
166  out << "\nB_k = rHL_k + rHB_k =\n" << B_k;
167  }
168 
169  // Solve the system pz = - inv(rHL) * qp_grad
170  VectorMutable &pz_k = s.pz().set_k(0);
171  V_InvMtV( &pz_k, B_k, no_trans, qp_grad_k );
172  Vt_S( &pz_k, -1.0 );
173 
174  // Zpz = Z * pz
175  V_MtV( &s.Zpz().set_k(0), s.Z().get_k(0), no_trans, pz_k );
176 
177  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
178  {
179  out << "\n||pz||inf = " << s.pz().get_k(0).norm_inf()
180  << "\nsum(Zpz) = " << AbstractLinAlgPack::sum(s.Zpz().get_k(0)) << std::endl;
181  }
182 
183  if( (int)olevel >= (int)PRINT_VECTORS )
184  {
185  out << "\npz_k = \n" << s.pz().get_k(0);
186  out << "\nnu_k = \n" << s.nu().get_k(0);
187  out << "\nZpz_k = \n" << s.Zpz().get_k(0);
188  out << std::endl;
189  }
190 
191  if(algo.algo_cntr().check_results())
192  {
193  assert_print_nan_inf(s.pz().get_k(0), "pz_k",true,&out);
194  assert_print_nan_inf(s.Zpz().get_k(0), "Zpz_k",true,&out);
195  }
196 
197  return true;
198  }
199 
201  , poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss
202  , std::ostream& out, const std::string& L ) const
203  {
204  out
205  << L << "*** Calculate the null space step by solving an unconstrainted QP\n"
206  << L << "zeta_k = 1.0\n"
207  << L << "qp_grad_k = Z_k'*(Gf_k + mu_k*(invXu_k-invXl_k)) + zeta_k*w_sigma_k\n"
208  << L << "B_k = rHL_k + rHB_k\n"
209  << L << "pz_k = -inv(B_k)*qp_grad_k\n"
210  << L << "Zpz_k = Z_k*pz_k\n"
211  ;
212  }
213 
214 } // end namespace MoochoPack
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
void assign(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = V_rhs.
rSQP Algorithm control class.
T * get() const
void V_StV(VectorMutable *v_lhs, value_type alpha, const V &V_rhs)
v_lhs = alpha * V_rhs.
void Mp_StM(DMatrixSlice *vs_lhs, value_type alpha, const MatrixOp &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1)
m_lhs += alpha * op(mwo_rhs1).
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
Transposed.
void V_InvMtV(DVectorSlice *vs_lhs, const MatrixOpNonsing &mwo_rhs1, BLAS_Cpp::Transp trans_rhs1, const DVectorSlice &vs_rhs2)
vs_lhs = inv(op(mwo_rhs1)) * vs_rhs2.
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.
void M_StM(MatrixOp *M_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
gm_lhs = alpha * M_rhs.
T_To & dyn_cast(T_From &from)
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.
void V_InvMtV(VectorMutable *v_lhs, const MatrixNonsing &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const Vector &v_rhs2)
v_lhs = inv(op(M_rhs1)) * v_rhs2
AlgorithmTracker & track()
void Mp_StM(MatrixOp *mwo_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
void print_step(const IterationPack::Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss, std::ostream &out, const std::string &leading_str) const
Called by Algorithm::print_algorithm() to print out what this step does in Matlab like format...
bool assert_print_nan_inf(const value_type &val, const char name[], bool throw_excpt, std::ostream *out)
This function asserts if a value_type scalare is a NaN or Inf and optionally prints out these entires...
const f_dbl_prec const f_int const f_int const f_int f_dbl_prec rhs[]
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
Acts as the central hub for an iterative algorithm.
value_type sum(const Vector &v_rhs)
result = sum( v_rhs(i), i = 1,,,dim )
NLPAlgo & rsqp_algo(Algorithm &algo)
Convert from a Algorithm to a NLPAlgo.