42 #ifndef MOOCHO_NO_BASIS_PERM_DIRECT_SOLVERS
65 namespace MoochoPack {
70 :select_new_decomposition_(false)
81 ,
bool *new_decomp_selected
90 out <<
"\nUpdating the decomposition ...\n";
93 DecompositionSystemVarReductPerm &decomp_sys_perm =
dyn_cast<DecompositionSystemVarReductPerm>(s.decomp_sys());
94 NLPVarReductPerm &nlp_vrp =
dyn_cast<NLPVarReductPerm>(nlp);
100 r = s.decomp_sys().equ_decomp().size();
102 bool decomp_updated =
false;
103 bool get_new_basis =
false;
104 bool new_basis_selected =
false;
108 out <<
"\nSome client called select_new_decomposition() so we will select a new basis ...\n";
109 get_new_basis =
true;
112 else if( !decomp_sys_perm.has_basis() ) {
114 out <<
"\nDecompositionSystemVarReductPerm object currently does not have a basis so we must select one ...\n";
115 get_new_basis =
true;
119 IterQuantityAccess<index_type>
120 &num_basis_iq = s.num_basis();
121 IterQuantityAccess<VectorMutable>
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;
138 const DecompositionSystem::ERunTests
139 ds_test_what = ( ( decomp_sys_testing ==
DST_TEST
141 && algo.algo_cntr().check_results() ) )
142 ? DecompositionSystem::RUN_TESTS
143 : DecompositionSystem::NO_TESTS );
146 DecompositionSystem::EOutputLevel ds_olevel;
150 ds_olevel = DecompositionSystem::PRINT_NONE;
154 ds_olevel = DecompositionSystem::PRINT_BASIC_INFO;
160 ds_olevel = DecompositionSystem::PRINT_EVERY_THING;
166 if( !get_new_basis ) {
169 out <<
"\nUpdating the range/null decompostion matrices ...\n";
172 s.decomp_sys().update_decomp(
180 ,Uz_iq ? &Uz_iq->set_k(0) : NULL
181 ,Uy_iq ? &Uy_iq->set_k(0) : NULL
182 ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS
184 s.equ_decomp( s.decomp_sys().equ_decomp() );
185 s.equ_undecomp( s.decomp_sys().equ_undecomp() );
186 decomp_updated =
true;
188 catch(
const DecompositionSystem::SingularDecomposition& except) {
191 <<
"\nOops! This decomposition was singular; must select a new basis!\n"
192 << except.what() << std::endl;
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();
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 );
210 P_var = nlp_vrp.factory_P_var()->create(),
211 P_equ = nlp_vrp.factory_P_equ()->create();
215 bool nlp_selected_basis =
false;
216 bool do_permute_x =
true;
217 if( nlp_vrp.nlp_selects_basis() ) {
227 <<
"\nThe NLP will attempt to select a basis "
228 <<
"(k = " << s.
k() <<
")...\n";
231 bool very_first_basis = !decomp_sys_perm.has_basis();
232 bool a_basis_was_singular =
false;
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 )
242 very_first_basis =
false;
243 decomp_sys_perm.set_decomp(
255 ,Uz_iq ? &Uz_iq->set_k(0) : NULL
256 ,Uy_iq ? &Uy_iq->set_k(0) : NULL
257 ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS
260 nlp_selected_basis =
true;
264 catch(
const DecompositionSystem::SingularDecomposition& except )
268 <<
"\nOops! This decomposition was singular; ask the NLP for another basis!\n"
269 << except.what() << std::endl;
270 a_basis_was_singular =
true;
274 do_permute_x = !( very_first_basis && !a_basis_was_singular );
277 <<
"\nThe NLP was unable to provide a nonsigular basis "
278 <<
"(k = " << s.
k() <<
")\n";
280 if(!nlp_selected_basis) {
287 <<
"\nThe decomposition system object is selecting the basis "
288 <<
"(k = " << s.
k() <<
")...\n";
290 decomp_sys_perm.select_decomp(
294 ,nu_iq.updated_k(0)?&nu_iq.get_k(0):NULL
303 ,Uz_iq ? &Uz_iq->set_k(0) : NULL
304 ,Uy_iq ? &Uy_iq->set_k(0) : NULL
305 ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS
307 nlp_vrp.set_basis( *P_var, var_dep, P_equ.
get(), &equ_decomp );
313 new_basis_selected =
true;
314 r = s.decomp_sys().equ_decomp().size();
319 last_updated_k = num_basis_iq.last_updated();
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;
324 s.equ_decomp( decomp_sys_perm.equ_decomp() );
325 s.equ_undecomp( decomp_sys_perm.equ_undecomp() );
327 s.set_P_var_last( s.get_P_var_current() );
328 s.set_P_equ_last( s.get_P_equ_current() );
330 s.set_P_var_current( P_var );
331 s.set_P_equ_current( P_equ );
334 out <<
"\nNew permutations:"
335 <<
"\nP_var_current() =\n" << s.P_var_current()
336 <<
"\nP_equ_current() =\n" << s.P_equ_current();
341 VectorMutable &x = x_iq.get_k(0);
344 out <<
"\nx resorted to the original order\n" << x;
348 out <<
"\nx resorted to new basis \n" << x;
358 *new_decomp_selected = new_basis_selected;
367 ,
const std::string &L
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"
382 << L <<
"if (decomp_sys does not have a basis) then\n"
383 << L <<
" get_new_basis = true\n"
385 << L <<
"if (get_new_basis == true) then\n"
386 << L <<
" begin update decomposition\n"
387 << L <<
" (class = \'" <<
typeName(s.decomp_sys()) <<
"\')\n"
389 s.decomp_sys().print_update_decomp( out, L +
" " );
391 << L <<
" end update decomposition\n"
392 << L <<
"if SingularDecomposition exception was not thrown then\n"
393 << L <<
" decomp_updated = true\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"
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"
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"
431 #endif // MOOCHO_NO_BASIS_PERM_DIRECT_SOLVERS
DecompositionSystemHandlerVarReductPerm_Strategy()
Constructor.
AbstractLinAlgPack::size_type size_type
std::string typeName(const T &t)
RTOp_index_type index_type
bool select_new_decomposition_
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.
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.
AlgorithmTracker & track()
void set_space_null(const vec_space_ptr_t &space_null)
Set the VectorSpace of the null space (pz).
void select_new_decomposition(bool select_new_decomposition)
void print_update_decomposition(const NLPAlgo &algo, const NLPAlgoState &s, std::ostream &out, const std::string &leading_spaces) const
RangePack::Range1D Range1D
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)