MoochoPack : Framework for Large-Scale Optimization Algorithms  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
MoochoPack_UpdateReducedSigma_Step.cpp
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 <typeinfo>
44 #include <iostream>
45 
46 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
47 #include "AbstractLinAlgPack_MatrixSymOpNonsing.hpp"
48 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
49 #include "AbstractLinAlgPack_MultiVectorMutable.hpp"
50 #include "AbstractLinAlgPack_VectorStdOps.hpp"
51 #include "AbstractLinAlgPack_VectorOut.hpp"
52 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
53 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
54 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
55 #include "ConstrainedOptPack_MatrixIdentConcat.hpp"
56 #include "IterationPack_print_algorithm_step.hpp"
57 #include "MoochoPack_IpState.hpp"
58 #include "MoochoPack_UpdateReducedSigma_Step.hpp"
59 #include "MoochoPack_moocho_algo_conversion.hpp"
60 
61 #include "OptionsFromStreamPack_StringToIntMap.hpp"
62 
63 #include "Teuchos_dyn_cast.hpp"
64 
65 namespace MoochoPack {
66 
68  const e_update_methods update_method
69  )
70  :
71  update_method_(update_method)
72  {}
73 
75  Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
76  ,poss_type assoc_step_poss
77  )
78  {
79  using Teuchos::dyn_cast;
80  using IterationPack::print_algorithm_step;
81 
82  NLPAlgo &algo = dyn_cast<NLPAlgo>(_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) )
90  {
91  using IterationPack::print_algorithm_step;
92  print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
93  }
94 
95  switch (update_method_)
96  {
97  case ALWAYS_EXPLICIT:
98  {
99  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
100  {
101  out << "\nupdate_method is always_explicit, forming Reduced Sigma Explicitly ...\n";
102  }
103  FormReducedSigmaExplicitly(algo,s,olevel,out);
104  }
105  break;
106  case BFGS_PRIMAL:
107  case BFGS_DUAL_NO_CORRECTION:
108  case BFGS_DUAL_EXPLICIT_CORRECTION:
109  case BFGS_DUAL_SCALING_CORRECTION:
110  {
111  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Option BFGS_primal not handled yet in UpdateReducedSigma\n");
112  }
113  break;
114  default:
115  TEUCHOS_TEST_FOR_EXCEPT(true); // local error ?
116  };
117 
118  if( (int)olevel >= (int)PRINT_ITERATION_QUANTITIES )
119  {
120  out << "\nrHB_k =\n" << s.rHB().get_k(0);
121  }
122 
123  return true;
124  }
125 
126 
127 void UpdateReducedSigma_Step::print_step(
128  const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
129  ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
130  ) const
131  {
132  out << L << "*** Update Z'*Sigma*Z\n"
133  << L << "if (update_method == always_explicit) then\n"
134  << L << " Sigma_k = invXl*Vl-invXu*Vu\n"
135  << L << " Sigma_I = Sigma_k.sub_view(Z.I_rng)\n"
136  << L << " Sigma_D_sqrt = (Sigma_k.sub_view(Z.D_rng))^1/2\n"
137  << L << " J = Sigma_D_sqrt*Z\n"
138  << L << " rHB_k = J'*J + Sigma_I\n"
139  << L << "elsif ( update_method == BFGS_??? ) then\n"
140  << L << " NOT IMPLEMENTED YET!\n"
141  << L << "end\n";
142  }
143 
144 void UpdateReducedSigma_Step::FormReducedSigmaExplicitly(
145  NLPAlgo& algo, IpState& s, EJournalOutputLevel olevel, std::ostream& out
146  )
147  {
148  namespace mmp = MemMngPack;
149  using Teuchos::dyn_cast;
152  using LinAlgOpPack::Mp_M;
153  using LinAlgOpPack::Mp_MtM;
154  using LinAlgOpPack::M_MtM;
155  using LinAlgOpPack::M_StM;
156  using LinAlgOpPack::V_MtV;
157  using LinAlgOpPack::assign;
158 
159  // Calculate Reduced Sigma directly from
160  // Sigma = invXl*Vl + invXu*Vu
161  // Z_kT*Sigma*Z_k
162 
163  // Get the iteration quantities
164  const MatrixIdentConcat &Z = dyn_cast<MatrixIdentConcat>(s.Z().get_k(0));
165  const MatrixSymDiagStd &invXl = s.invXl().get_k(0);
166  const MatrixSymDiagStd &invXu = s.invXu().get_k(0);
167  const MatrixSymDiagStd &Vl = s.Vl().get_k(0);
168  const MatrixSymDiagStd &Vu = s.Vu().get_k(0);
169 
170  MatrixSymDiagStd &Sigma = s.Sigma().set_k(0);
171 
172  MatrixSymOpNonsing& rHB = dyn_cast<MatrixSymOpNonsing>(s.rHB().set_k(0));
173  if (rHB.cols() != Z.cols())
174  {
175  // Initialize space in rHB
176  dyn_cast<MatrixSymInitDiag>(rHB).init_identity(Z.space_rows(), 0.0);
177  }
178 
179  // Calculate Sigma = invXl*Vl + invXu*Vu
180 
181  ele_wise_prod(1.0, invXl.diag(), Vl.diag(), &(Sigma.diag() = 0.0));
182 
183  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
184  out << "\n||Sigma_l||inf = " << Sigma.diag().norm_inf() << std::endl;
185  if( (int)olevel >= (int)PRINT_VECTORS )
186  out << "\nSigma_l =\n" << Sigma.diag();
187 
188  ele_wise_prod(1.0, invXu.diag(), Vu.diag(), &Sigma.diag() );
189 
190  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
191  out << "\n||Sigma_k||inf = ||Sigma_l + Sigma_u||inf = " << Sigma.diag().norm_inf() << std::endl;
192  if( (int)olevel >= (int)PRINT_VECTORS )
193  out << "\nSigma_k = Sigma_l + Sigma_u =\n" << Sigma.diag();
194 
195  // Calculate the cross term (Z'*Sigma*Ypy) first
196  VectorSpace::vec_mut_ptr_t temp = Z.space_cols().create_member(0.0);
197  ele_wise_prod(1.0, s.Ypy().get_k(0), Sigma.diag(), temp.get());
198  VectorMutable& w_sigma = s.w_sigma().set_k(0);
199  V_MtV(&w_sigma, Z, BLAS_Cpp::trans, *temp);
200 
201  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
202  out << "\n||w_sigma_k||inf = " << w_sigma.norm_inf() << std::endl;
203  if( (int)olevel >= (int)PRINT_VECTORS )
204  out << "\nw_sigma_k = \n" << w_sigma;
205 
206  // Calculate Reduced Sigma
207  // Try sigma^1/2 making use of dependent and independent variables
208 
210  Sigma_D_diag = Sigma.diag().sub_view(Z.D_rng()),
211  Sigma_I_diag = Sigma.diag().sub_view(Z.I_rng());
212  const size_type
213  Sigma_D_nz = Sigma_D_diag->nz(),
214  Sigma_I_nz = Sigma_I_diag->nz();
215 
216  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
217  {
218  out << "\nSigma_D.diag().nz() = " << Sigma_D_nz;
219  out << "\nSigma_I.diag().nz() = " << Sigma_I_nz << std::endl;
220  }
221 
222  if( Sigma_D_nz || Sigma_I_nz )
223  {
224  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
225  {
226  out << "\nForming explicit, nonzero rHB_k = Z_k'*Sigma_k*Z_k ...\n";
227  }
228  if( Sigma_D_nz )
229  {
230 
231  MatrixSymDiagStd Sigma_D_sqrt(Sigma_D_diag->clone());
232 
233  ele_wise_sqrt(&Sigma_D_sqrt.diag());
234 
235  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) )
236  {
237  out << "\nSigma_D_sqrt =\n" << Sigma_D_sqrt;
238  }
239 
241  J = Sigma_D_sqrt.space_cols().create_members(Z.cols());
242  M_MtM(
243  static_cast<MatrixOp*>(J.get())
244  ,Sigma_D_sqrt, BLAS_Cpp::no_trans, Z.D(), BLAS_Cpp::no_trans);
245 
246  LinAlgOpPack::syrk( *J, BLAS_Cpp::trans, 1.0, 0.0, &rHB );
247 
248  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) )
249  {
250  out << "\nJ =\n" << *J;
251  out << "\nJ'*J =\n" << rHB;
252  }
253 
254  }
255 
256  if( Sigma_I_nz )
257  {
258 
259  const MatrixSymDiagStd Sigma_I(
260  Teuchos::rcp_const_cast<VectorMutable>(Sigma_I_diag)
261  );
262 
263  if(Sigma_D_nz)
264  {
265  Mp_M( &rHB, Sigma_I, BLAS_Cpp::no_trans );
266  }
267  else
268  {
269  assign( &rHB, Sigma_I, BLAS_Cpp::no_trans );
270  }
271 
272  }
273  }
274  else
275  {
276  if( (int)olevel >= (int)PRINT_ALGORITHM_STEPS )
277  {
278  out << "\nSigma_k is zero so setting rHB_k = 0.0 ...\n";
279  }
280  rHB.zero_out();
281  }
282 
283  /*
284  // Try the full unspecialised calculation, this is expensive, but will
285  // serve as a debug for the more efficient calculations.
286 
287  VectorSpace::multi_vec_mut_ptr_t Sigma_Z = Z.space_cols().create_members(Z.cols());
288  M_MtM((MatrixOp*)Sigma_Z.get(), Sigma, BLAS_Cpp::no_trans, Z, BLAS_Cpp::no_trans);
289 
290  //std::cout << "Sigma_Z\n";
291  //Sigma_Z->output(std::cout);
292 
293  VectorSpace::multi_vec_mut_ptr_t ZT_Sigma_Z = Z.space_rows().create_members(Z.cols());
294  M_MtM((MatrixOp*)ZT_Sigma_Z.get(), (MatrixOp&)Z, BLAS_Cpp::trans, (MatrixOp&)*Sigma_Z, BLAS_Cpp::no_trans);
295 
296  std::cout << "ZT_Sigma_Z=\n";
297  ZT_Sigma_Z->output(std::cout);
298  */
299  }
300 
301 
302 namespace {
303 
304 const int local_num_options = 1;
305 
306 enum local_EOptions
307  {
308  UPDATE_METHOD
309  };
310 
311 const char* local_SOptions[local_num_options] =
312  {
313  "update_method",
314  };
315 
316 const int num_update_methods = 5;
317 
318 const char* s_update_methods[num_update_methods] =
319  {
320  "always_explicit",
321  "BFGS_primal",
322  "BFGS_dual_no_correction",
323  "BFGS_dual_explicit_correction",
324  "BFGS_dual_scaling_correction"
325  };
326 
327 }
328 
329 
330 UpdateReducedSigma_StepSetOptions::UpdateReducedSigma_StepSetOptions(
331  UpdateReducedSigma_Step* target
332  , const char opt_grp_name[] )
333  :
334  OptionsFromStreamPack::SetOptionsFromStreamNode(
335  opt_grp_name, local_num_options, local_SOptions ),
336  OptionsFromStreamPack::SetOptionsToTargetBase< UpdateReducedSigma_Step >( target )
337  {
338  }
339 
340 void UpdateReducedSigma_StepSetOptions::setOption(
341  int option_num, const std::string& option_value )
342  {
344 
345  typedef UpdateReducedSigma_Step target_t;
346  switch( (local_EOptions)option_num )
347  {
348  case UPDATE_METHOD:
349  {
350  StringToIntMap config_map( UpdateReducedSigma_opt_grp_name, num_update_methods, s_update_methods );
351  target().update_method( (UpdateReducedSigma_Step::e_update_methods) config_map( option_value.c_str() ) );
352  }
353  break;
354  default:
355  TEUCHOS_TEST_FOR_EXCEPT(true); // Local error only?
356  }
357  }
358 
359 } // end namespace MoochoPack
UpdateReducedSigma_Step(const e_update_methods update_method=ALWAYS_EXPLICIT)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
T_To & dyn_cast(T_From &from)
rSQP Algorithm control class.
T * get() const
void assign(MatrixOp *M_lhs, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
virtual std::ostream & journal_out() const
void M_StM(MatrixOp *M_lhs, value_type alpha, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
size_t size_type
void syrk(const MatrixOp &mwo_rhs, BLAS_Cpp::Transp M_trans, value_type alpha, value_type beta, MatrixSymOp *sym_lhs)
void ele_wise_prod(const value_type &alpha, const Vector &v_rhs1, const Vector &v_rhs2, VectorMutable *v_lhs)
AlgorithmTracker & track()
void Mp_M(MatrixOp *M_lhs, const MatrixOp &M_rhs, BLAS_Cpp::Transp trans_rhs)
AlgorithmState & state()
void ele_wise_sqrt(VectorMutable *z)
bool do_step(Algorithm &algo, poss_type step_poss, IterationPack::EDoStepType type, poss_type assoc_step_poss)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)