MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MoochoPack_DecompositionSystemHandlerVarReductPerm_Strategy.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 #ifndef MOOCHO_NO_BASIS_PERM_DIRECT_SOLVERS
43 
44 #include <ostream>
45 #include <typeinfo>
46 
62 #include "Teuchos_dyn_cast.hpp"
63 #include "Teuchos_Assert.hpp"
64 
65 namespace MoochoPack {
66 
67 // Constructors / initializers
68 
70  :select_new_decomposition_(false)
71 {}
72 
73 // Overridden from DecompositionSystemHandler_Strategy
74 
76  NLPAlgo &algo
77  ,NLPAlgoState &s
78  ,NLPFirstOrder &nlp
79  ,EDecompSysTesting decomp_sys_testing
80  ,EDecompSysPrintLevel decomp_sys_testing_print_level
81  ,bool *new_decomp_selected
82  )
83 {
84  using Teuchos::dyn_cast;
85 
86  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
87  std::ostream& out = algo.track().journal_out();
88 
89  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
90  out << "\nUpdating the decomposition ...\n";
91  }
92 
93  DecompositionSystemVarReductPerm &decomp_sys_perm = dyn_cast<DecompositionSystemVarReductPerm>(s.decomp_sys());
94  NLPVarReductPerm &nlp_vrp = dyn_cast<NLPVarReductPerm>(nlp);
95 
96  const size_type
97  n = nlp.n(),
98  m = nlp.m();
99  size_type
100  r = s.decomp_sys().equ_decomp().size();
101 
102  bool decomp_updated = false;
103  bool get_new_basis = false;
104  bool new_basis_selected = false;
105 
107  if( olevel >= PRINT_ALGORITHM_STEPS )
108  out << "\nSome client called select_new_decomposition() so we will select a new basis ...\n";
109  get_new_basis = true;
111  }
112  else if( !decomp_sys_perm.has_basis() ) {
113  if( olevel >= PRINT_ALGORITHM_STEPS )
114  out << "\nDecompositionSystemVarReductPerm object currently does not have a basis so we must select one ...\n";
115  get_new_basis = true;
116  }
117 
118  // Get the iteration quantity container objects
119  IterQuantityAccess<index_type>
120  &num_basis_iq = s.num_basis();
121  IterQuantityAccess<VectorMutable>
122  &x_iq = s.x(),
123  &nu_iq = s.nu();
124  IterQuantityAccess<MatrixOp>
125  *Gc_iq = m > 0 ? &s.Gc() : NULL,
126  *Z_iq = ( n > m && r > 0 ) || get_new_basis ? &s.Z() : NULL,
127  *Y_iq = ( r > 0 ) || get_new_basis ? &s.Y() : NULL,
128  *Uz_iq = ( m > 0 && m > r ) || get_new_basis ? &s.Uz() : NULL,
129  *Uy_iq = ( m > 0 && m > r ) || get_new_basis ? &s.Uy() : NULL;
130  IterQuantityAccess<MatrixOpNonsing>
131  *R_iq = ( m > 0 ) ? &s.R() : NULL;
132 
133  //
134  // Update (or select a new) range/null decomposition
135  //
136 
137  // Determine if we will test the decomp_sys or not
138  const DecompositionSystem::ERunTests
139  ds_test_what = ( ( decomp_sys_testing == DST_TEST
140  || ( decomp_sys_testing == DST_DEFAULT
141  && algo.algo_cntr().check_results() ) )
142  ? DecompositionSystem::RUN_TESTS
143  : DecompositionSystem::NO_TESTS );
144 
145  // Determine the output level for decomp_sys
146  DecompositionSystem::EOutputLevel ds_olevel;
147  switch(olevel) {
148  case PRINT_NOTHING:
150  ds_olevel = DecompositionSystem::PRINT_NONE;
151  break;
153  case PRINT_ACTIVE_SET:
154  ds_olevel = DecompositionSystem::PRINT_BASIC_INFO;
155  break;
156  case PRINT_VECTORS:
158  break;
160  ds_olevel = DecompositionSystem::PRINT_EVERY_THING;
161  break;
162  default:
163  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
164  };
165 
166  if( !get_new_basis ) {
167  // Form the decomposition of Gc and Gh and update the decomposition system matrices
168  if( olevel >= PRINT_ALGORITHM_STEPS ) {
169  out << "\nUpdating the range/null decompostion matrices ...\n";
170  }
171  try {
172  s.decomp_sys().update_decomp(
173  static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
174  ,ds_olevel // olevel
175  ,ds_test_what // test_what
176  ,Gc_iq->get_k(0) // Gc
177  ,&Z_iq->set_k(0) // Z
178  ,&Y_iq->set_k(0) // Y
179  ,&R_iq->set_k(0) // R
180  ,Uz_iq ? &Uz_iq->set_k(0) : NULL // Uz
181  ,Uy_iq ? &Uy_iq->set_k(0) : NULL // Uy
182  ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this!
183  );
184  s.equ_decomp( s.decomp_sys().equ_decomp() );
185  s.equ_undecomp( s.decomp_sys().equ_undecomp() );
186  decomp_updated = true;
187  }
188  catch( const DecompositionSystem::SingularDecomposition& except) {
189  if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
190  out
191  << "\nOops! This decomposition was singular; must select a new basis!\n"
192  << except.what() << std::endl;
193  }
194  }
195 
196  if( !decomp_updated ) {
197  if( s.get_P_var_current().get() == NULL ) {
199  P_var = nlp_vrp.factory_P_var()->create(),
200  P_equ = nlp_vrp.factory_P_equ()->create();
201  Range1D
202  var_dep,
203  equ_decomp;
204  nlp_vrp.get_basis(
205  P_var.get(), &var_dep, P_equ.get(), &equ_decomp );
206  s.set_P_var_current( P_var );
207  s.set_P_equ_current( P_equ );
208  }
210  P_var = nlp_vrp.factory_P_var()->create(),
211  P_equ = nlp_vrp.factory_P_equ()->create();
212  Range1D
213  var_dep,
214  equ_decomp;
215  bool nlp_selected_basis = false;
216  bool do_permute_x = true;
217  if( nlp_vrp.nlp_selects_basis() ) {
218  // The nlp may select the new (or first) basis.
219 
220  // If this is the first basis, the NLPVarReductPerm interface specifies that it
221  // will already be set for the nlp. Check to see if this is the first basis
222  // and if not, ask the nlp to give you the next basis.
223  // I must form a loop here to deal with the
224  // possibility that the basis the nlp selects will be singular.
225  if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
226  out
227  << "\nThe NLP will attempt to select a basis "
228  << "(k = " << s.k() << ")...\n";
229  // If decomp_sys_per->has_basis() == false, the first execution of the while()
230  // statement will not execute get_next_basis(...).
231  bool very_first_basis = !decomp_sys_perm.has_basis();
232  bool a_basis_was_singular = false;
233  if(very_first_basis)
234  nlp_vrp.get_basis(
235  P_var.get(), &var_dep, P_equ.get(), &equ_decomp );
236  while( very_first_basis
237  || nlp_vrp.get_next_basis(
238  P_var.get(), &var_dep, P_equ.get(), &equ_decomp )
239  )
240  {
241  try {
242  very_first_basis = false;
243  decomp_sys_perm.set_decomp(
244  static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
245  ,ds_olevel // olevel
246  ,ds_test_what // test_what
247  ,*P_var // P_var
248  ,var_dep // var_dep
249  ,P_equ.get() // P_equ
250  ,&equ_decomp // equ_decomp
251  ,Gc_iq->get_k(0) // Gc
252  ,&Z_iq->set_k(0) // Z
253  ,&Y_iq->set_k(0) // Y
254  ,&R_iq->set_k(0) // R
255  ,Uz_iq ? &Uz_iq->set_k(0) : NULL // Uz
256  ,Uy_iq ? &Uy_iq->set_k(0) : NULL // Uy
257  ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this to MATRICES_INDEP_IMPS
258  );
259  // If you get here the basis was not singular.
260  nlp_selected_basis = true;
261  break; // break out of the while(...) loop
262  }
263  // Catch the singularity exceptions and loop around
264  catch( const DecompositionSystem::SingularDecomposition& except )
265  {
266  if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
267  out
268  << "\nOops! This decomposition was singular; ask the NLP for another basis!\n"
269  << except.what() << std::endl;
270  a_basis_was_singular = true;
271  }
272  // Any other exception gets thrown clean out of here.
273  }
274  do_permute_x = !( very_first_basis && !a_basis_was_singular );
275  if( olevel >= PRINT_BASIC_ALGORITHM_INFO && !nlp_selected_basis )
276  out
277  << "\nThe NLP was unable to provide a nonsigular basis "
278  << "(k = " << s.k() << ")\n";
279  }
280  if(!nlp_selected_basis) {
281  // If you get into here then the nlp could not select a nonsingular
282  // basis so we will let the decomposition system select a basis.
283  // and give it to the nlp.
284 
285  if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ) {
286  out
287  << "\nThe decomposition system object is selecting the basis "
288  << "(k = " << s.k() << ")...\n";
289  }
290  decomp_sys_perm.select_decomp(
291  static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
292  ,ds_olevel // olevel
293  ,ds_test_what // test_what
294  ,nu_iq.updated_k(0)?&nu_iq.get_k(0):NULL // nu
295  ,&Gc_iq->get_k(0) // Gc
296  ,P_var.get() // P_var
297  ,&var_dep // var_dep
298  ,P_equ.get() // P_equ
299  ,&equ_decomp // equ_decomp
300  ,&Z_iq->set_k(0) // Z
301  ,&Y_iq->set_k(0) // Y
302  ,&R_iq->set_k(0) // R
303  ,Uz_iq ? &Uz_iq->set_k(0) : NULL // Uz
304  ,Uy_iq ? &Uy_iq->set_k(0) : NULL // Uy
305  ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this to MATRICES_INDEP_IMPS
306  );
307  nlp_vrp.set_basis( *P_var, var_dep, P_equ.get(), &equ_decomp );
308  }
309 
310  // If you get here then no unexpected exceptions where thrown and a new
311  // basis has been selected.
312 
313  new_basis_selected = true;
314  r = s.decomp_sys().equ_decomp().size();
315 
316  // Record this basis change
317 
318  const int
319  last_updated_k = num_basis_iq.last_updated();
320  const index_type
321  num_basis = ( last_updated_k != IterQuantity::NONE_UPDATED ? num_basis_iq.get_k(last_updated_k) : 0 ) + 1;
322  num_basis_iq.set_k(0) = num_basis;
323 
324  s.equ_decomp( decomp_sys_perm.equ_decomp() );
325  s.equ_undecomp( decomp_sys_perm.equ_undecomp() );
326 
327  s.set_P_var_last( s.get_P_var_current() );
328  s.set_P_equ_last( s.get_P_equ_current() );
329 
330  s.set_P_var_current( P_var );
331  s.set_P_equ_current( P_equ );
332 
333  if( olevel >= PRINT_VECTORS ) {
334  out << "\nNew permutations:"
335  << "\nP_var_current() =\n" << s.P_var_current()
336  << "\nP_equ_current() =\n" << s.P_equ_current();
337  }
338 
339  if( do_permute_x ) {
340  // Sort x according to this new basis.
341  VectorMutable &x = x_iq.get_k(0);
342  s.P_var_last().permute( BLAS_Cpp::trans, &x ); // Permute back to original order
343  if( olevel >= PRINT_VECTORS ) {
344  out << "\nx resorted to the original order\n" << x;
345  }
346  s.P_var_current().permute( BLAS_Cpp::no_trans, &x ); // Permute to new (current) order
347  if( olevel >= PRINT_VECTORS ) {
348  out << "\nx resorted to new basis \n" << x;
349  }
350  }
351 
352  // Set the new range and null spaces (these will update all of the set vectors!)
353  s.set_space_range( decomp_sys_perm.space_range() );
354  s.set_space_null( decomp_sys_perm.space_null() );
355 
356  }
357 
358  *new_decomp_selected = new_basis_selected;
359 
360  return true;
361 }
362 
364  const NLPAlgo &algo
365  ,const NLPAlgoState &s
366  ,std::ostream &out
367  ,const std::string &L
368  ) const
369 {
370  out
371  << L << "*** Updating or selecting a new decomposition using a variable reduction\n"
372  << L << "*** range/null decomposition object.\n"
373  << L << "if decomp_sys does not support the DecompositionSystemVarReductPerm interface then throw exception!\n"
374  << L << "if nlp does not support the NLPVarReductPerm interface then throw exception!\n"
375  << L << "decomp_updated = false\n"
376  << L << "get_new_basis = false\n"
377  << L << "new_basis_selected = false\n"
378  << L << "if( select_new_decomposition == true ) then\n"
379  << L << " get_new_basis = true\n"
380  << L << " select_new_decomposition = false\n"
381  << L << "end\n"
382  << L << "if (decomp_sys does not have a basis) then\n"
383  << L << " get_new_basis = true\n"
384  << L << "end\n"
385  << L << "if (get_new_basis == true) then\n"
386  << L << " begin update decomposition\n"
387  << L << " (class = \'" << typeName(s.decomp_sys()) << "\')\n"
388  ;
389  s.decomp_sys().print_update_decomp( out, L + " " );
390  out
391  << L << " end update decomposition\n"
392  << L << "if SingularDecomposition exception was not thrown then\n"
393  << L << " decomp_updated = true\n"
394  << L << "end\n"
395  << L << "if (decomp_updated == false) then\n"
396  << L << " nlp_selected_basis = false\n"
397  << L << " if (nlp.selects_basis() == true) then\n"
398  << L << " for each basis returned from nlp.get_basis(...) or nlp.get_next_basis()\n"
399  << L << " decomp_sys.set_decomp(Gc_k...) -> Z_k,Y_k,R_k,Uz_k,Uy_k \n"
400  << L << " if SingularDecompositon exception was not thrown then\n"
401  << L << " nlp_selected_basis = true\n"
402  << L << " exit loop\n"
403  << L << " end\n"
404  << L << " end\n"
405  << L << " end\n"
406  << L << " if (nlp_selected_basis == false) then\n"
407  << L << " decomp_sys.select_decomp(Gc_k...) -> P_var,P_equ,Z,Y,R,Uz,Uy\n"
408  << L << " and permute Gc\n"
409  << L << " end\n"
410  << L << " *** If you get here then no unexpected exceptions were thrown and a new\n"
411  << L << " *** basis has been selected\n"
412  << L << " num_basis_k = num_basis_k(last_updated) + 1\n"
413  << L << " P_var_last = P_var_current\n"
414  << L << " P_equ_last = P_equ_current\n"
415  << L << " P_var_current = P_var\n"
416  << L << " P_equ_current = P_equ\n"
417  << L << " Resort x_k according to P_var_current\n"
418  << L << "end\n"
419  ;
420 }
421 
422 // Overridden from DecompositionSystemHandlerSelectNew_Strategy
423 
425 {
427 }
428 
429 } // end namespace MoochoPack
430 
431 #endif // MOOCHO_NO_BASIS_PERM_DIRECT_SOLVERS
AbstractLinAlgPack::size_type size_type
std::string typeName(const T &t)
bool update_decomposition(NLPAlgo &algo, NLPAlgoState &s, NLPFirstOrder &nlp, EDecompSysTesting decomp_sys_testing, EDecompSysPrintLevel decomp_sys_testing_print_level, bool *new_decomp_selected)
rSQP Algorithm control class.
T * get() const
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.
T_To & dyn_cast(T_From &from)
Reduced space SQP state encapsulation interface.
std::ostream * out
AlgorithmTracker & track()
void set_space_null(const vec_space_ptr_t &space_null)
Set the VectorSpace of the null space (pz).
void print_update_decomposition(const NLPAlgo &algo, const NLPAlgoState &s, std::ostream &out, const std::string &leading_spaces) const
int n
void set_space_range(const vec_space_ptr_t &space_range)
Set the VectorSpace of the range space (py).
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)