MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_DecompositionSystemHandlerVarReductPerm_Strategy.cpp
Go to the documentation of this file.
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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef MOOCHO_NO_BASIS_PERM_DIRECT_SOLVERS
00043 
00044 #include <ostream>
00045 #include <typeinfo>
00046 
00047 #include "MoochoPack_DecompositionSystemHandlerVarReductPerm_Strategy.hpp"
00048 #include "MoochoPack_Exceptions.hpp"
00049 #include "MoochoPack_moocho_algo_conversion.hpp"
00050 #include "IterationPack_print_algorithm_step.hpp"
00051 #include "ConstrainedOptPack_DecompositionSystemVarReductPerm.hpp"
00052 #include "NLPInterfacePack_NLPFirstOrder.hpp"
00053 #include "NLPInterfacePack_NLPVarReductPerm.hpp"
00054 #include "AbstractLinAlgPack_PermutationOut.hpp"
00055 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00056 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00057 #include "AbstractLinAlgPack_VectorMutable.hpp"
00058 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00059 #include "AbstractLinAlgPack_VectorOut.hpp"
00060 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00061 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00062 #include "Teuchos_dyn_cast.hpp"
00063 #include "Teuchos_Assert.hpp"
00064 
00065 namespace MoochoPack {
00066 
00067 // Constructors / initializers
00068 
00069 DecompositionSystemHandlerVarReductPerm_Strategy::DecompositionSystemHandlerVarReductPerm_Strategy()
00070   :select_new_decomposition_(false)
00071 {}
00072 
00073 // Overridden from DecompositionSystemHandler_Strategy
00074 
00075 bool DecompositionSystemHandlerVarReductPerm_Strategy::update_decomposition(
00076   NLPAlgo                                &algo
00077   ,NLPAlgoState                          &s
00078   ,NLPFirstOrder                         &nlp
00079   ,EDecompSysTesting                     decomp_sys_testing
00080   ,EDecompSysPrintLevel                  decomp_sys_testing_print_level
00081   ,bool                                  *new_decomp_selected
00082   )
00083 {
00084   using Teuchos::dyn_cast;
00085 
00086   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00087   std::ostream& out = algo.track().journal_out();
00088 
00089   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00090     out << "\nUpdating the decomposition ...\n";
00091   }
00092 
00093   DecompositionSystemVarReductPerm   &decomp_sys_perm = dyn_cast<DecompositionSystemVarReductPerm>(s.decomp_sys());
00094   NLPVarReductPerm                   &nlp_vrp         = dyn_cast<NLPVarReductPerm>(nlp);
00095 
00096   const size_type
00097     n  = nlp.n(),
00098     m  = nlp.m();
00099   size_type
00100     r  = s.decomp_sys().equ_decomp().size();
00101 
00102   bool decomp_updated     = false;
00103   bool get_new_basis      = false;
00104   bool new_basis_selected = false;
00105 
00106   if( select_new_decomposition_ ) {
00107     if( olevel >= PRINT_ALGORITHM_STEPS )
00108       out << "\nSome client called select_new_decomposition() so we will select a new basis ...\n";
00109     get_new_basis = true;
00110     select_new_decomposition_ = false;
00111   }
00112   else if( !decomp_sys_perm.has_basis() ) {
00113     if( olevel >= PRINT_ALGORITHM_STEPS )
00114       out << "\nDecompositionSystemVarReductPerm object currently does not have a basis so we must select one ...\n";
00115     get_new_basis = true;
00116   }
00117 
00118   // Get the iteration quantity container objects
00119   IterQuantityAccess<index_type>
00120     &num_basis_iq = s.num_basis();
00121   IterQuantityAccess<VectorMutable>
00122     &x_iq   = s.x(),
00123     &nu_iq  = s.nu();
00124   IterQuantityAccess<MatrixOp>
00125     *Gc_iq  = m  > 0                                   ? &s.Gc() : NULL,
00126     *Z_iq   = ( n > m && r > 0 )    || get_new_basis   ? &s.Z()  : NULL,
00127     *Y_iq   = ( r > 0 )             || get_new_basis   ? &s.Y()  : NULL,
00128     *Uz_iq  = ( m  > 0 && m  > r )  || get_new_basis   ? &s.Uz() : NULL,
00129     *Uy_iq  = ( m  > 0 && m  > r )  || get_new_basis   ? &s.Uy() : NULL;
00130   IterQuantityAccess<MatrixOpNonsing>
00131     *R_iq   = ( m > 0 )                                ? &s.R()  : NULL;
00132   
00133   //
00134   // Update (or select a new) range/null decomposition
00135   //
00136 
00137   // Determine if we will test the decomp_sys or not
00138   const DecompositionSystem::ERunTests
00139     ds_test_what = ( ( decomp_sys_testing == DST_TEST
00140                || ( decomp_sys_testing == DST_DEFAULT
00141                 && algo.algo_cntr().check_results() ) )
00142              ? DecompositionSystem::RUN_TESTS
00143              : DecompositionSystem::NO_TESTS );
00144     
00145   // Determine the output level for decomp_sys        
00146   DecompositionSystem::EOutputLevel ds_olevel;
00147   switch(olevel) {
00148     case PRINT_NOTHING:
00149     case PRINT_BASIC_ALGORITHM_INFO:
00150       ds_olevel = DecompositionSystem::PRINT_NONE;
00151       break;
00152     case PRINT_ALGORITHM_STEPS:
00153     case PRINT_ACTIVE_SET:
00154       ds_olevel = DecompositionSystem::PRINT_BASIC_INFO;
00155       break;
00156     case PRINT_VECTORS:
00157       ds_olevel = DecompositionSystem::PRINT_VECTORS;
00158       break;
00159     case PRINT_ITERATION_QUANTITIES:
00160       ds_olevel = DecompositionSystem::PRINT_EVERY_THING;
00161       break;
00162     default:
00163       TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
00164   };
00165 
00166   if( !get_new_basis ) {
00167     // Form the decomposition of Gc and Gh and update the decomposition system matrices
00168     if( olevel >= PRINT_ALGORITHM_STEPS ) {
00169       out << "\nUpdating the range/null decompostion matrices ...\n";
00170     }
00171     try {
00172       s.decomp_sys().update_decomp(
00173         static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
00174         ,ds_olevel                         // olevel
00175         ,ds_test_what                      // test_what
00176         ,Gc_iq->get_k(0)                   // Gc
00177         ,&Z_iq->set_k(0)                   // Z
00178         ,&Y_iq->set_k(0)                   // Y
00179         ,&R_iq->set_k(0)                   // R
00180         ,Uz_iq ? &Uz_iq->set_k(0) : NULL   // Uz
00181         ,Uy_iq ? &Uy_iq->set_k(0) : NULL   // Uy
00182         ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this!
00183         );
00184       s.equ_decomp(   s.decomp_sys().equ_decomp()   );
00185       s.equ_undecomp( s.decomp_sys().equ_undecomp() );
00186       decomp_updated = true;
00187     }
00188     catch( const DecompositionSystem::SingularDecomposition& except) {
00189       if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
00190         out
00191           << "\nOops!  This decomposition was singular; must select a new basis!\n"
00192           << except.what() << std::endl;
00193     }
00194   }
00195 
00196   if( !decomp_updated ) {
00197     if( s.get_P_var_current().get() == NULL ) {
00198       Teuchos::RCP<Permutation>
00199         P_var = nlp_vrp.factory_P_var()->create(),
00200         P_equ = nlp_vrp.factory_P_equ()->create();
00201       Range1D
00202         var_dep,
00203         equ_decomp;
00204       nlp_vrp.get_basis(
00205         P_var.get(), &var_dep, P_equ.get(), &equ_decomp );
00206       s.set_P_var_current( P_var );
00207       s.set_P_equ_current( P_equ );
00208     }
00209     Teuchos::RCP<Permutation>
00210       P_var = nlp_vrp.factory_P_var()->create(),
00211       P_equ = nlp_vrp.factory_P_equ()->create();
00212     Range1D
00213       var_dep,
00214       equ_decomp;
00215     bool nlp_selected_basis = false;
00216     bool do_permute_x = true;
00217     if( nlp_vrp.nlp_selects_basis() ) {
00218       // The nlp may select the new (or first) basis.
00219         
00220       // If this is the first basis, the NLPVarReductPerm interface specifies that it
00221       // will already be set for the nlp.  Check to see if this is the first basis
00222       // and if not, ask the nlp to give you the next basis.
00223       // I must form a loop here to deal with the
00224       // possibility that the basis the nlp selects will be singular.
00225       if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
00226         out
00227           << "\nThe NLP will attempt to select a basis "
00228           << "(k = " << s.k() << ")...\n";
00229       // If decomp_sys_per->has_basis() == false, the first execution of the while()
00230       // statement will not execute get_next_basis(...).    
00231       bool very_first_basis = !decomp_sys_perm.has_basis();
00232       bool a_basis_was_singular = false;
00233       if(very_first_basis)
00234         nlp_vrp.get_basis(
00235           P_var.get(), &var_dep, P_equ.get(), &equ_decomp );
00236       while( very_first_basis
00237            || nlp_vrp.get_next_basis(
00238              P_var.get(), &var_dep, P_equ.get(), &equ_decomp )
00239         )
00240       {
00241         try {
00242           very_first_basis = false;
00243           decomp_sys_perm.set_decomp(
00244             static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
00245             ,ds_olevel                         // olevel
00246             ,ds_test_what                      // test_what
00247             ,*P_var                            // P_var
00248             ,var_dep                           // var_dep
00249             ,P_equ.get()                       // P_equ
00250             ,&equ_decomp                       // equ_decomp
00251             ,Gc_iq->get_k(0)                   // Gc
00252             ,&Z_iq->set_k(0)                   // Z
00253             ,&Y_iq->set_k(0)                   // Y
00254             ,&R_iq->set_k(0)                   // R
00255             ,Uz_iq ? &Uz_iq->set_k(0) : NULL   // Uz
00256             ,Uy_iq ? &Uy_iq->set_k(0) : NULL   // Uy
00257             ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this to MATRICES_INDEP_IMPS
00258             );
00259           // If you get here the basis was not singular.
00260           nlp_selected_basis = true;
00261           break; // break out of the while(...) loop
00262         }
00263         // Catch the singularity exceptions and loop around
00264         catch( const DecompositionSystem::SingularDecomposition& except )
00265         {
00266           if( olevel >= PRINT_BASIC_ALGORITHM_INFO )
00267             out
00268               << "\nOops!  This decomposition was singular; ask the NLP for another basis!\n"
00269               << except.what() << std::endl;
00270           a_basis_was_singular = true;
00271         }
00272         // Any other exception gets thrown clean out of here.
00273       }
00274       do_permute_x =  !( very_first_basis && !a_basis_was_singular );
00275       if( olevel >= PRINT_BASIC_ALGORITHM_INFO && !nlp_selected_basis )
00276         out
00277           << "\nThe NLP was unable to provide a nonsigular basis "
00278           << "(k = " << s.k() << ")\n";
00279     }
00280     if(!nlp_selected_basis) {
00281       // If you get into here then the nlp could not select a nonsingular
00282       // basis so we will let the decomposition system select a basis.
00283       // and give it to the nlp.
00284       
00285       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ) {
00286         out
00287           << "\nThe decomposition system object is selecting the basis "
00288           << "(k = " << s.k() << ")...\n";
00289       }
00290       decomp_sys_perm.select_decomp(
00291         static_cast<int>(olevel) >= static_cast<int>(PRINT_BASIC_ALGORITHM_INFO) ? &out : NULL // out
00292         ,ds_olevel                                  // olevel
00293         ,ds_test_what                               // test_what
00294         ,nu_iq.updated_k(0)?&nu_iq.get_k(0):NULL    // nu
00295         ,&Gc_iq->get_k(0)                           // Gc
00296         ,P_var.get()                                // P_var
00297         ,&var_dep                                   // var_dep
00298         ,P_equ.get()                                // P_equ
00299         ,&equ_decomp                                // equ_decomp
00300         ,&Z_iq->set_k(0)                            // Z
00301         ,&Y_iq->set_k(0)                            // Y
00302         ,&R_iq->set_k(0)                            // R
00303         ,Uz_iq ? &Uz_iq->set_k(0) : NULL            // Uz
00304         ,Uy_iq ? &Uy_iq->set_k(0) : NULL            // Uy
00305         ,DecompositionSystem::MATRICES_ALLOW_DEP_IMPS // ToDo: Change this to MATRICES_INDEP_IMPS
00306         );
00307       nlp_vrp.set_basis(  *P_var, var_dep, P_equ.get(), &equ_decomp );
00308     }
00309     
00310     // If you get here then no unexpected exceptions where thrown and a new
00311     // basis has been selected.
00312         
00313     new_basis_selected = true;
00314     r                  = s.decomp_sys().equ_decomp().size();
00315 
00316     // Record this basis change
00317 
00318     const int
00319       last_updated_k = num_basis_iq.last_updated();
00320     const index_type
00321       num_basis = ( last_updated_k != IterQuantity::NONE_UPDATED ? num_basis_iq.get_k(last_updated_k) : 0 ) + 1;
00322     num_basis_iq.set_k(0) = num_basis;
00323 
00324     s.equ_decomp(   decomp_sys_perm.equ_decomp()   );
00325     s.equ_undecomp( decomp_sys_perm.equ_undecomp() );
00326           
00327     s.set_P_var_last( s.get_P_var_current() );
00328     s.set_P_equ_last( s.get_P_equ_current() );
00329 
00330     s.set_P_var_current( P_var );
00331     s.set_P_equ_current( P_equ );
00332 
00333     if( olevel >= PRINT_VECTORS ) {
00334       out << "\nNew permutations:"
00335         << "\nP_var_current() =\n" << s.P_var_current()
00336         << "\nP_equ_current() =\n" << s.P_equ_current();
00337     }
00338           
00339     if( do_permute_x ) {
00340       // Sort x according to this new basis.
00341       VectorMutable &x = x_iq.get_k(0);
00342       s.P_var_last().permute( BLAS_Cpp::trans, &x ); // Permute back to original order
00343       if( olevel >= PRINT_VECTORS ) {
00344         out << "\nx resorted to the original order\n" << x;
00345       }
00346       s.P_var_current().permute( BLAS_Cpp::no_trans, &x ); // Permute to new (current) order
00347       if( olevel >= PRINT_VECTORS ) {
00348         out << "\nx resorted to new basis \n" << x;
00349       }
00350     }
00351           
00352     // Set the new range and null spaces (these will update all of the set vectors!)
00353     s.set_space_range( decomp_sys_perm.space_range() );
00354     s.set_space_null(  decomp_sys_perm.space_null()  );
00355 
00356   }
00357 
00358   *new_decomp_selected = new_basis_selected;
00359 
00360   return true;
00361 }
00362 
00363 void DecompositionSystemHandlerVarReductPerm_Strategy::print_update_decomposition(
00364   const NLPAlgo                          &algo
00365   ,const NLPAlgoState                        &s
00366   ,std::ostream                           &out
00367   ,const std::string                      &L
00368   ) const
00369 {
00370   out
00371     << L << "*** Updating or selecting a new decomposition using a variable reduction\n"
00372     << L << "*** range/null decomposition object.\n"
00373     << L << "if decomp_sys does not support the DecompositionSystemVarReductPerm interface then throw exception!\n"
00374     << L << "if nlp does not support the NLPVarReductPerm interface then throw exception!\n"
00375     << L << "decomp_updated     = false\n"
00376     << L << "get_new_basis      = false\n"
00377     << L << "new_basis_selected = false\n"
00378     << L << "if( select_new_decomposition == true ) then\n"
00379     << L << "  get_new_basis            = true\n"
00380     << L << "  select_new_decomposition = false\n"
00381     << L << "end\n"
00382     << L << "if (decomp_sys does not have a basis) then\n"
00383     << L << "  get_new_basis = true\n"
00384     << L << "end\n"
00385     << L << "if (get_new_basis == true) then\n"
00386     << L << "  begin update decomposition\n"
00387     << L << "  (class = \'" << typeName(s.decomp_sys()) << "\')\n"
00388     ;
00389   s.decomp_sys().print_update_decomp( out, L + "    " );
00390   out
00391     << L << "  end update decomposition\n"
00392     << L << "if SingularDecomposition exception was not thrown then\n"
00393     << L << "  decomp_updated = true\n"
00394     << L << "end\n"
00395     << L << "if (decomp_updated == false) then\n"
00396     << L << "  nlp_selected_basis = false\n"
00397     << L << "  if (nlp.selects_basis() == true) then\n"
00398     << L << "    for each basis returned from nlp.get_basis(...) or nlp.get_next_basis()\n"
00399     << L << "      decomp_sys.set_decomp(Gc_k...) -> Z_k,Y_k,R_k,Uz_k,Uy_k \n"
00400     << L << "      if SingularDecompositon exception was not thrown then\n"
00401     << L << "        nlp_selected_basis = true\n"
00402     << L << "        exit loop\n"
00403     << L << "      end\n"
00404     << L << "    end\n"
00405     << L << "  end\n"
00406     << L << "  if (nlp_selected_basis == false) then\n"
00407     << L << "    decomp_sys.select_decomp(Gc_k...) -> P_var,P_equ,Z,Y,R,Uz,Uy\n"
00408     << L << "                                              and permute Gc\n"
00409     << L << "  end\n"
00410     << L << "  *** If you get here then no unexpected exceptions were thrown and a new\n"
00411     << L << "  *** basis has been selected\n"
00412     << L << "  num_basis_k = num_basis_k(last_updated) + 1\n"
00413     << L << "  P_var_last = P_var_current\n"
00414     << L << "  P_equ_last = P_equ_current\n"
00415     << L << "  P_var_current = P_var\n"
00416     << L << "  P_equ_current = P_equ\n"
00417     << L << "  Resort x_k according to P_var_current\n"
00418     << L << "end\n"
00419     ;
00420 }
00421 
00422 // Overridden from DecompositionSystemHandlerSelectNew_Strategy
00423 
00424 void DecompositionSystemHandlerVarReductPerm_Strategy::select_new_decomposition( bool select_new_decomposition )
00425 {
00426   select_new_decomposition_ = select_new_decomposition;
00427 }
00428 
00429 } // end namespace MoochoPack
00430 
00431 #endif // MOOCHO_NO_BASIS_PERM_DIRECT_SOLVERS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines