MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MoochoPack_ReducedHessianExactStd_Step.cpp
Go to the documentation of this file.
1 #if 0
2 
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
7 // Copyright (2003) Sandia Corporation
8 //
9 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10 // license for use of this work by or on behalf of the U.S. Government.
11 //
12 // Redistribution and use in source and binary forms, with or without
13 // modification, are permitted provided that the following conditions are
14 // met:
15 //
16 // 1. Redistributions of source code must retain the above copyright
17 // notice, this list of conditions and the following disclaimer.
18 //
19 // 2. Redistributions in binary form must reproduce the above copyright
20 // notice, this list of conditions and the following disclaimer in the
21 // documentation and/or other materials provided with the distribution.
22 //
23 // 3. Neither the name of the Corporation nor the names of the
24 // contributors may be used to endorse or promote products derived from
25 // this software without specific prior written permission.
26 //
27 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 //
39 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
40 //
41 // ***********************************************************************
42 // @HEADER
43 
44 #include <sstream>
45 #include <typeinfo>
46 #include <iomanip>
47 
52 #include "ConstrainedOptPack/src/VectorWithNorms.h"
54 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_MatrixSymOp.hpp"
61 #include "Midynamic_cast_verbose.h"
62 
63 namespace MoochoPack {
64 
66  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
67  , poss_type assoc_step_poss)
68 {
69  using Teuchos::dyn_cast;
72  typedef AbstractLinAlgPack::MatrixSymDenseInitialize MatrixSymDenseInitialize;
73  typedef AbstractLinAlgPack::MatrixSymOp MatrixSymOp;
74  using ConstrainedOptPack::NLPSecondOrder;
75 
76  NLPAlgo &algo = rsqp_algo(_algo);
77  NLPAlgoState &s = algo.rsqp_state();
78  NLPSecondOrder
79 #ifdef _WINDOWS
80  &nlp = dynamic_cast<NLPSecondOrder&>(algo.nlp());
81 #else
82  &nlp = dyn_cast<NLPSecondOrder>(algo.nlp());
83 #endif
84  MatrixSymOp
85  *HL_sym_op = dynamic_cast<MatrixSymOp*>(&s.HL().get_k(0));
86 
87  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
88  std::ostream& out = algo.track().journal_out();
89 
90  // print step header.
91  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
93  print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
94  }
95 
96  // problem size
97  size_type n = nlp.n(),
98  r = nlp.r(),
99  nind = n - r;
100 
101  // Compute HL first (You may want to move this into its own step later?)
102 
103  if( !s.lambda().updated_k(-1) ) {
104  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
105  out << "Initializing lambda_km1 = nlp.get_lambda_init ... \n";
106  }
107  nlp.get_init_lagrange_mult( &s.lambda().set_k(-1).v(), NULL );
108  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
109  out << "||lambda_km1||inf = " << s.lambda().get_k(-1).norm_inf() << std::endl;
110  }
111  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
112  out << "lambda_km1 = \n" << s.lambda().get_k(-1)();
113  }
114  }
115 
116  nlp.set_HL( HL_sym_op );
117  nlp.calc_HL( s.x().get_k(0)(), s.lambda().get_k(-1)(), false );
118 
119  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
120  s.HL().get_k(0).output( out << "\nHL_k = \n" );
121  }
122 
123  // If rHL has already been updated for this iteration then just leave it.
124  if( !s.rHL().updated_k(0) ) {
125 
126  if( !HL_sym_op ) {
127  std::ostringstream omsg;
128  omsg
129  << "ReducedHessianExactStd_Step::do_step(...) : Error, "
130  << "The matrix HL with the concrete type "
131  << typeName(s.HL().get_k(0)) << " does not support the "
132  << "MatrixSymOp iterface";
133  throw std::logic_error( omsg.str() );
134  }
135 
136  MatrixSymDenseInitialize
137  *rHL_sym_init = dynamic_cast<MatrixSymDenseInitialize*>(&s.rHL().set_k(0));
138  if( !rHL_sym_init ) {
139  std::ostringstream omsg;
140  omsg
141  << "ReducedHessianExactStd_Step::do_step(...) : Error, "
142  << "The matrix rHL with the concrete type "
143  << typeName(s.rHL().get_k(0)) << " does not support the "
144  << "MatrixSymDenseInitialize iterface";
145  throw std::logic_error( omsg.str() );
146  }
147 
148  // Compute the dense reduced Hessian
149  DMatrix rHL_sym_store(nind,nind);
150  DMatrixSliceSym rHL_sym(rHL_sym_store(),BLAS_Cpp::lower);
151  Mp_StMtMtM( &rHL_sym, 1.0, MatrixSymOp::DUMMY_ARG, *HL_sym_op
152  , s.Z().get_k(0), BLAS_Cpp::no_trans, 0.0 );
153 
154  if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
155  out << "\nLower triangular partion of dense reduced Hessian (ignore nonzeros above diagonal):\nrHL_dense = \n" << rHL_sym_store();
156  }
157 
158  // Set the reduced Hessain
159  rHL_sym_init->initialize( rHL_sym );
160 
161  if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES ) {
162  s.rHL().get_k(0).output( out << "\nrHL_k = \n" );
163  }
164  }
165 
166  return true;
167 }
168 
170  const Algorithm& algo, poss_type step_poss, IterationPack::EDoStepType type
171  , poss_type assoc_step_poss, std::ostream& out, const std::string& L ) const
172 {
173  out
174  << L << "*** Calculate the exact reduced Hessian of the Lagrangian\n"
175  << L << "if lambda_km1 is not updated then\n"
176  << L << " lambda_km1 = nlp.get_lambda_init\n"
177  << L << "end\n"
178  << L << "HL_k = HL(x_k,lambda_km1) <: R^(n+m) -> R^(n x n)\n"
179  << L << "if rHL_k is not updated then\n"
180  << L << " rHL_dense = Z_k' * HL_k * Z_k (MatrixSymOp interface for HL_k)\n"
181  << L << " rHL_k = rHL_dense (MatrixSymDenseInitialize interface for rHL_k)\n"
182  << L << "end\n";
183 }
184 
185 } // end namespace MoochoPack
186 
187 #endif // 0
AbstractLinAlgPack::size_type size_type
std::string typeName(const T &t)
Interface adding operations specific for a symmetric matrix {abstract}.
Not transposed.
DMatrixSliceSym nonconst_sym(DMatrixSlice gms, BLAS_Cpp::Uplo uplo)
Return a symmetric matrix.
EJournalOutputLevel
enum for journal output.
void Mp_StMtMtM(MatrixSymOp *sym_lhs, value_type alpha, MatrixSymOp::EMatRhsPlaceHolder dummy_place_holder, const MatrixSymOp &M, const MatrixOp &mwo_rhs, BLAS_Cpp::Transp mwo_rhs_trans, value_type beta=1.0)
sym_lhs = alpha * op(mwo_rhs') * M * op(mwo_rhs) + beta * sym_lhs
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 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
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
DenseLinAlgPack::DMatrixSliceSym DMatrixSliceSym
Mix-in Interface for initializing a matrix with a dense symmetric matrix.
int n
NLPAlgo & rsqp_algo(Algorithm &algo)
Convert from a Algorithm to a NLPAlgo.