MoochoPack_DecompositionSystemHandlerVarReductPerm_Strategy.cpp

00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef MOOCHO_NO_BASIS_PERM_DIRECT_SOLVERS
00030 
00031 #include <ostream>
00032 #include <typeinfo>
00033 
00034 #include "MoochoPack_DecompositionSystemHandlerVarReductPerm_Strategy.hpp"
00035 #include "MoochoPack_Exceptions.hpp"
00036 #include "MoochoPack_moocho_algo_conversion.hpp"
00037 #include "IterationPack_print_algorithm_step.hpp"
00038 #include "ConstrainedOptPack_DecompositionSystemVarReductPerm.hpp"
00039 #include "NLPInterfacePack_NLPFirstOrder.hpp"
00040 #include "NLPInterfacePack_NLPVarReductPerm.hpp"
00041 #include "AbstractLinAlgPack_PermutationOut.hpp"
00042 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00043 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00044 #include "AbstractLinAlgPack_VectorMutable.hpp"
00045 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00046 #include "AbstractLinAlgPack_VectorOut.hpp"
00047 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00048 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00049 #include "Teuchos_dyn_cast.hpp"
00050 #include "Teuchos_TestForException.hpp"
00051 
00052 namespace MoochoPack {
00053 
00054 // Constructors / initializers
00055 
00056 DecompositionSystemHandlerVarReductPerm_Strategy::DecompositionSystemHandlerVarReductPerm_Strategy()
00057   :select_new_decomposition_(false)
00058 {}
00059 
00060 // Overridden from DecompositionSystemHandler_Strategy
00061 
00062 bool DecompositionSystemHandlerVarReductPerm_Strategy::update_decomposition(
00063   NLPAlgo                                &algo
00064   ,NLPAlgoState                          &s
00065   ,NLPFirstOrder                         &nlp
00066   ,EDecompSysTesting                     decomp_sys_testing
00067   ,EDecompSysPrintLevel                  decomp_sys_testing_print_level
00068   ,bool                                  *new_decomp_selected
00069   )
00070 {
00071   using Teuchos::dyn_cast;
00072 
00073   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00074   std::ostream& out = algo.track().journal_out();
00075 
00076   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00077     out << "\nUpdating the decomposition ...\n";
00078   }
00079 
00080   DecompositionSystemVarReductPerm   &decomp_sys_perm = dyn_cast<DecompositionSystemVarReductPerm>(s.decomp_sys());
00081   NLPVarReductPerm                   &nlp_vrp         = dyn_cast<NLPVarReductPerm>(nlp);
00082 
00083   const size_type
00084     n  = nlp.n(),
00085     m  = nlp.m();
00086   size_type
00087     r  = s.decomp_sys().equ_decomp().size();
00088 
00089   bool decomp_updated     = false;
00090   bool get_new_basis      = false;
00091   bool new_basis_selected = false;
00092 
00093   if( select_new_decomposition_ ) {
00094     if( olevel >= PRINT_ALGORITHM_STEPS )
00095       out << "\nSome client called select_new_decomposition() so we will select a new basis ...\n";
00096     get_new_basis = true;
00097     select_new_decomposition_ = false;
00098   }
00099   else if( !decomp_sys_perm.has_basis() ) {
00100     if( olevel >= PRINT_ALGORITHM_STEPS )
00101       out << "\nDecompositionSystemVarReductPerm object currently does not have a basis so we must select one ...\n";
00102     get_new_basis = true;
00103   }
00104 
00105   // Get the iteration quantity container objects
00106   IterQuantityAccess<index_type>
00107     &num_basis_iq = s.num_basis();
00108   IterQuantityAccess<VectorMutable>
00109     &x_iq   = s.x(),
00110     &nu_iq  = s.nu();
00111   IterQuantityAccess<MatrixOp>
00112     *Gc_iq  = m  > 0                                   ? &s.Gc() : NULL,
00113     *Z_iq   = ( n > m && r > 0 )    || get_new_basis   ? &s.Z()  : NULL,
00114     *Y_iq   = ( r > 0 )             || get_new_basis   ? &s.Y()  : NULL,
00115     *Uz_iq  = ( m  > 0 && m  > r )  || get_new_basis   ? &s.Uz() : NULL,
00116     *Uy_iq  = ( m  > 0 && m  > r )  || get_new_basis   ? &s.Uy() : NULL;
00117   IterQuantityAccess<MatrixOpNonsing>
00118     *R_iq   = ( m > 0 )                                ? &s.R()  : NULL;
00119   
00120   //
00121   // Update (or select a new) range/null decomposition
00122   //
00123 
00124   // Determine if we will test the decomp_sys or not
00125   const DecompositionSystem::ERunTests
00126     ds_test_what = ( ( decomp_sys_testing == DST_TEST
00127                || ( decomp_sys_testing == DST_DEFAULT
00128                 && algo.algo_cntr().check_results() ) )
00129              ? DecompositionSystem::RUN_TESTS
00130              : DecompositionSystem::NO_TESTS );
00131     
00132   // Determine the output level for decomp_sys        
00133   DecompositionSystem::EOutputLevel ds_olevel;
00134   switch(olevel) {
00135     case PRINT_NOTHING:
00136     case PRINT_BASIC_ALGORITHM_INFO:
00137       ds_olevel = DecompositionSystem::PRINT_NONE;
00138       break;
00139     case PRINT_ALGORITHM_STEPS:
00140     case PRINT_ACTIVE_SET:
00141       ds_olevel = DecompositionSystem::PRINT_BASIC_INFO;
00142       break;
00143     case PRINT_VECTORS:
00144       ds_olevel = DecompositionSystem::PRINT_VECTORS;
00145       break;
00146     case PRINT_ITERATION_QUANTITIES:
00147       ds_olevel = DecompositionSystem::PRINT_EVERY_THING;
00148       break;
00149     default:
00150       TEST_FOR_EXCEPT(true); // Should not get here!
00151   };
00152 
00153   if( !get_new_basis ) {
00154     // Form the decomposition of Gc and Gh and update the decomposition system matrices
00155     if( olevel >= PRINT_ALGORITHM_STEPS ) {
00156       out << "\nUpdating the range/null decompostion matrices ...\n";
00157     }
00158     try {
00159       s.decomp_sys().update_decomp(
00160         static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
00161         ,ds_olevel                         // olevel
00162         ,ds_test_what                      // test_what
00163         ,Gc_iq->get_k(0)                   // Gc
00164         ,&Z_iq->set_k(0)                   // Z
00165         ,&Y_iq->set_k(0)                   // Y
00166         ,&R_iq->set_k(0)                   // R
00167         ,Uz_iq ? &Uz_iq->set_k(0) : NULL   // Uz
00168         ,Uy_iq ? &Uy_iq->set_k(0) : NULL   // Uy
00169         ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this!
00170         );
00171       s.equ_decomp(   s.decomp_sys().equ_decomp()   );
00172       s.equ_undecomp( s.decomp_sys().equ_undecomp() );
00173       decomp_updated = true;
00174     }
00175     catch( const DecompositionSystem::SingularDecomposition& except) {
00176       if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
00177         out
00178           << "\nOops!  This decomposition was singular; must select a new basis!\n"
00179           << except.what() << std::endl;
00180     }
00181   }
00182 
00183   if( !decomp_updated ) {
00184     if( s.get_P_var_current().get() == NULL ) {
00185       Teuchos::RCP<Permutation>
00186         P_var = nlp_vrp.factory_P_var()->create(),
00187         P_equ = nlp_vrp.factory_P_equ()->create();
00188       Range1D
00189         var_dep,
00190         equ_decomp;
00191       nlp_vrp.get_basis(
00192         P_var.get(), &var_dep, P_equ.get(), &equ_decomp );
00193       s.set_P_var_current( P_var );
00194       s.set_P_equ_current( P_equ );
00195     }
00196     Teuchos::RCP<Permutation>
00197       P_var = nlp_vrp.factory_P_var()->create(),
00198       P_equ = nlp_vrp.factory_P_equ()->create();
00199     Range1D
00200       var_dep,
00201       equ_decomp;
00202     bool nlp_selected_basis = false;
00203     bool do_permute_x = true;
00204     if( nlp_vrp.nlp_selects_basis() ) {
00205       // The nlp may select the new (or first) basis.
00206         
00207       // If this is the first basis, the NLPVarReductPerm interface specifies that it
00208       // will already be set for the nlp.  Check to see if this is the first basis
00209       // and if not, ask the nlp to give you the next basis.
00210       // I must form a loop here to deal with the
00211       // possibility that the basis the nlp selects will be singular.
00212       if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
00213         out
00214           << "\nThe NLP will attempt to select a basis "
00215           << "(k = " << s.k() << ")...\n";
00216       // If decomp_sys_per->has_basis() == false, the first execution of the while()
00217       // statement will not execute get_next_basis(...).    
00218       bool very_first_basis = !decomp_sys_perm.has_basis();
00219       bool a_basis_was_singular = false;
00220       if(very_first_basis)
00221         nlp_vrp.get_basis(
00222           P_var.get(), &var_dep, P_equ.get(), &equ_decomp );
00223       while( very_first_basis
00224            || nlp_vrp.get_next_basis(
00225              P_var.get(), &var_dep, P_equ.get(), &equ_decomp )
00226         )
00227       {
00228         try {
00229           very_first_basis = false;
00230           decomp_sys_perm.set_decomp(
00231             static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
00232             ,ds_olevel                         // olevel
00233             ,ds_test_what                      // test_what
00234             ,*P_var                            // P_var
00235             ,var_dep                           // var_dep
00236             ,P_equ.get()                       // P_equ
00237             ,&equ_decomp                       // equ_decomp
00238             ,Gc_iq->get_k(0)                   // Gc
00239             ,&Z_iq->set_k(0)                   // Z
00240             ,&Y_iq->set_k(0)                   // Y
00241             ,&R_iq->set_k(0)                   // R
00242             ,Uz_iq ? &Uz_iq->set_k(0) : NULL   // Uz
00243             ,Uy_iq ? &Uy_iq->set_k(0) : NULL   // Uy
00244             ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this to MATRICES_INDEP_IMPS
00245             );
00246           // If you get here the basis was not singular.
00247           nlp_selected_basis = true;
00248           break; // break out of the while(...) loop
00249         }
00250         // Catch the singularity exceptions and loop around
00251         catch( const DecompositionSystem::SingularDecomposition& except )
00252         {
00253           if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
00254             out
00255               << "\nOops!  This decomposition was singular; ask the NLP for another basis!\n"
00256               << except.what() << std::endl;
00257           a_basis_was_singular = true;
00258         }
00259         // Any other exception gets thrown clean out of here.
00260       }
00261       do_permute_x =  !( very_first_basis && !a_basis_was_singular );
00262       if( olevel >= PRINT_BASIC_ALGORITHM_INFO && !nlp_selected_basis )
00263         out
00264           << "\nThe NLP was unable to provide a nonsigular basis "
00265           << "(k = " << s.k() << ")\n";
00266     }
00267     if(!nlp_selected_basis) {
00268       // If you get into here then the nlp could not select a nonsingular
00269       // basis so we will let the decomposition system select a basis.
00270       // and give it to the nlp.
00271       
00272       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ) {
00273         out
00274           << "\nThe decomposition system object is selecting the basis "
00275           << "(k = " << s.k() << ")...\n";
00276       }
00277       decomp_sys_perm.select_decomp(
00278         static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
00279         ,ds_olevel                                  // olevel
00280         ,ds_test_what                               // test_what
00281         ,nu_iq.updated_k(0)?&nu_iq.get_k(0):NULL    // nu
00282         ,&Gc_iq->get_k(0)                           // Gc
00283         ,P_var.get()                                // P_var
00284         ,&var_dep                                   // var_dep
00285         ,P_equ.get()                                // P_equ
00286         ,&equ_decomp                                // equ_decomp
00287         ,&Z_iq->set_k(0)                            // Z
00288         ,&Y_iq->set_k(0)                            // Y
00289         ,&R_iq->set_k(0)                            // R
00290         ,Uz_iq ? &Uz_iq->set_k(0) : NULL            // Uz
00291         ,Uy_iq ? &Uy_iq->set_k(0) : NULL            // Uy
00292         ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this to MATRICES_INDEP_IMPS
00293         );
00294       nlp_vrp.set_basis(  *P_var, var_dep, P_equ.get(), &equ_decomp );
00295     }
00296     
00297     // If you get here then no unexpected exceptions where thrown and a new
00298     // basis has been selected.
00299         
00300     new_basis_selected = true;
00301     r                  = s.decomp_sys().equ_decomp().size();
00302 
00303     // Record this basis change
00304 
00305     const int
00306       last_updated_k = num_basis_iq.last_updated();
00307     const index_type
00308       num_basis = ( last_updated_k != IterQuantity::NONE_UPDATED ? num_basis_iq.get_k(last_updated_k) : 0 ) + 1;
00309     num_basis_iq.set_k(0) = num_basis;
00310 
00311     s.equ_decomp(   decomp_sys_perm.equ_decomp()   );
00312     s.equ_undecomp( decomp_sys_perm.equ_undecomp() );
00313           
00314     s.set_P_var_last( s.get_P_var_current() );
00315     s.set_P_equ_last( s.get_P_equ_current() );
00316 
00317     s.set_P_var_current( P_var );
00318     s.set_P_equ_current( P_equ );
00319 
00320     if( olevel >= PRINT_VECTORS ) {
00321       out << "\nNew permutations:"
00322         << "\nP_var_current() =\n" << s.P_var_current()
00323         << "\nP_equ_current() =\n" << s.P_equ_current();
00324     }
00325           
00326     if( do_permute_x ) {
00327       // Sort x according to this new basis.
00328       VectorMutable &x = x_iq.get_k(0);
00329       s.P_var_last().permute( BLAS_Cpp::trans, &x ); // Permute back to original order
00330       if( olevel >= PRINT_VECTORS ) {
00331         out << "\nx resorted to the original order\n" << x;
00332       }
00333       s.P_var_current().permute( BLAS_Cpp::no_trans, &x ); // Permute to new (current) order
00334       if( olevel >= PRINT_VECTORS ) {
00335         out << "\nx resorted to new basis \n" << x;
00336       }
00337     }
00338           
00339     // Set the new range and null spaces (these will update all of the set vectors!)
00340     s.set_space_range( decomp_sys_perm.space_range() );
00341     s.set_space_null(  decomp_sys_perm.space_null()  );
00342 
00343   }
00344 
00345   *new_decomp_selected = new_basis_selected;
00346 
00347   return true;
00348 }
00349 
00350 void DecompositionSystemHandlerVarReductPerm_Strategy::print_update_decomposition(
00351   const NLPAlgo                          &algo
00352   ,const NLPAlgoState                        &s
00353   ,std::ostream                           &out
00354   ,const std::string                      &L
00355   ) const
00356 {
00357   out
00358     << L << "*** Updating or selecting a new decomposition using a variable reduction\n"
00359     << L << "*** range/null decomposition object.\n"
00360     << L << "if decomp_sys does not support the DecompositionSystemVarReductPerm interface then throw exception!\n"
00361     << L << "if nlp does not support the NLPVarReductPerm interface then throw exception!\n"
00362     << L << "decomp_updated     = false\n"
00363     << L << "get_new_basis      = false\n"
00364     << L << "new_basis_selected = false\n"
00365     << L << "if( select_new_decomposition == true ) then\n"
00366     << L << "  get_new_basis            = true\n"
00367     << L << "  select_new_decomposition = false\n"
00368     << L << "end\n"
00369     << L << "if (decomp_sys does not have a basis) then\n"
00370     << L << "  get_new_basis = true\n"
00371     << L << "end\n"
00372     << L << "if (get_new_basis == true) then\n"
00373     << L << "  begin update decomposition\n"
00374     << L << "  (class = \'" << typeName(s.decomp_sys()) << "\')\n"
00375     ;
00376   s.decomp_sys().print_update_decomp( out, L + "    " );
00377   out
00378     << L << "  end update decomposition\n"
00379     << L << "if SingularDecomposition exception was not thrown then\n"
00380     << L << "  decomp_updated = true\n"
00381     << L << "end\n"
00382     << L << "if (decomp_updated == false) then\n"
00383     << L << "  nlp_selected_basis = false\n"
00384     << L << "  if (nlp.selects_basis() == true) then\n"
00385     << L << "    for each basis returned from nlp.get_basis(...) or nlp.get_next_basis()\n"
00386     << L << "      decomp_sys.set_decomp(Gc_k...) -> Z_k,Y_k,R_k,Uz_k,Uy_k \n"
00387     << L << "      if SingularDecompositon exception was not thrown then\n"
00388     << L << "        nlp_selected_basis = true\n"
00389     << L << "        exit loop\n"
00390     << L << "      end\n"
00391     << L << "    end\n"
00392     << L << "  end\n"
00393     << L << "  if (nlp_selected_basis == false) then\n"
00394     << L << "    decomp_sys.select_decomp(Gc_k...) -> P_var,P_equ,Z,Y,R,Uz,Uy\n"
00395     << L << "                                              and permute Gc\n"
00396     << L << "  end\n"
00397     << L << "  *** If you get here then no unexpected exceptions were thrown and a new\n"
00398     << L << "  *** basis has been selected\n"
00399     << L << "  num_basis_k = num_basis_k(last_updated) + 1\n"
00400     << L << "  P_var_last = P_var_current\n"
00401     << L << "  P_equ_last = P_equ_current\n"
00402     << L << "  P_var_current = P_var\n"
00403     << L << "  P_equ_current = P_equ\n"
00404     << L << "  Resort x_k according to P_var_current\n"
00405     << L << "end\n"
00406     ;
00407 }
00408 
00409 // Overridden from DecompositionSystemHandlerSelectNew_Strategy
00410 
00411 void DecompositionSystemHandlerVarReductPerm_Strategy::select_new_decomposition( bool select_new_decomposition )
00412 {
00413   select_new_decomposition_ = select_new_decomposition;
00414 }
00415 
00416 } // end namespace MoochoPack
00417 
00418 #endif // MOOCHO_NO_BASIS_PERM_DIRECT_SOLVERS

Generated on Wed May 12 21:32:12 2010 for MoochoPack : Framework for Large-Scale Optimization Algorithms by  doxygen 1.4.7