MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_EvalNewPointStd_Step.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 #include <ostream>
00043 #include <typeinfo>
00044 
00045 #include "MoochoPack_EvalNewPointStd_Step.hpp"
00046 #include "MoochoPack_Exceptions.hpp"
00047 #include "MoochoPack_moocho_algo_conversion.hpp"
00048 #include "IterationPack_print_algorithm_step.hpp"
00049 #include "NLPInterfacePack_NLPFirstOrder.hpp"
00050 #include "ConstrainedOptPack_DecompositionSystemVarReduct.hpp"
00051 #include "AbstractLinAlgPack_MatrixSymIdent.hpp"
00052 #include "AbstractLinAlgPack_PermutationOut.hpp"
00053 #include "AbstractLinAlgPack_MatrixOpNonsing.hpp"
00054 #include "AbstractLinAlgPack_MatrixOpOut.hpp"
00055 #include "AbstractLinAlgPack_VectorMutable.hpp"
00056 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00057 #include "AbstractLinAlgPack_VectorOut.hpp"
00058 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00059 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00060 #include "Teuchos_dyn_cast.hpp"
00061 #include "Teuchos_Assert.hpp"
00062 
00063 #ifdef TEUCHOS_DEBUG
00064 #include "DenseLinAlgPack_PermVecMat.hpp"
00065 #endif
00066 
00067 namespace MoochoPack {
00068 
00069 EvalNewPointStd_Step::EvalNewPointStd_Step(
00070   const decomp_sys_handler_ptr_t                              &decomp_sys_handler
00071   ,const deriv_tester_ptr_t                                   &deriv_tester
00072   ,const bounds_tester_ptr_t                                  &bounds_tester
00073   ,const decomp_sys_tester_ptr_t                              &decomp_sys_tester
00074   ,EFDDerivTesting                                            fd_deriv_testing
00075   ,DecompositionSystemHandler_Strategy::EDecompSysTesting     decomp_sys_testing
00076   ,DecompositionSystemHandler_Strategy::EDecompSysPrintLevel  decomp_sys_testing_print_level
00077   )
00078   :decomp_sys_handler_(decomp_sys_handler)
00079   ,deriv_tester_(deriv_tester)
00080   ,bounds_tester_(bounds_tester)
00081   ,decomp_sys_tester_(decomp_sys_tester)
00082   ,fd_deriv_testing_(fd_deriv_testing)
00083   ,decomp_sys_testing_(decomp_sys_testing)
00084   ,decomp_sys_testing_print_level_(decomp_sys_testing_print_level)
00085 {}
00086 
00087 bool EvalNewPointStd_Step::do_step(
00088   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00089   ,poss_type assoc_step_poss
00090   )
00091 {
00092   using Teuchos::rcp;
00093   using Teuchos::dyn_cast;
00094   using AbstractLinAlgPack::assert_print_nan_inf;
00095   using IterationPack::print_algorithm_step;
00096   using NLPInterfacePack::NLPFirstOrder;
00097 
00098   NLPAlgo         &algo   = rsqp_algo(_algo);
00099   NLPAlgoState    &s      = algo.rsqp_state();
00100   NLPFirstOrder   &nlp    = dyn_cast<NLPFirstOrder>(algo.nlp());
00101 
00102   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00103   EJournalOutputLevel ns_olevel = algo.algo_cntr().null_space_journal_output_level();
00104   std::ostream& out = algo.track().journal_out();
00105 
00106   // print step header.
00107   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00108     using IterationPack::print_algorithm_step;
00109     print_algorithm_step( algo, step_poss, type, assoc_step_poss, out );
00110   }
00111 
00112   if(!nlp.is_initialized())
00113     nlp.initialize(algo.algo_cntr().check_results());
00114 
00115   Teuchos::VerboseObjectTempState<NLP>
00116     nlpOutputTempState(rcp(&nlp,false),Teuchos::getFancyOStream(rcp(&out,false)),convertToVerbLevel(olevel));
00117 
00118   const size_type
00119     n  = nlp.n(),
00120     nb = nlp.num_bounded_x(),
00121     m  = nlp.m();
00122   size_type
00123     r  = 0;
00124 
00125   // Get the iteration quantity container objects
00126   IterQuantityAccess<value_type>
00127     &f_iq   = s.f();
00128   IterQuantityAccess<VectorMutable>
00129     &x_iq   = s.x(),
00130     *c_iq   = m > 0  ? &s.c() : NULL,
00131     &Gf_iq  = s.Gf();
00132   IterQuantityAccess<MatrixOp>
00133     *Gc_iq  = m  > 0 ? &s.Gc() : NULL,
00134     *Z_iq   = NULL,
00135     *Y_iq   = NULL,
00136     *Uz_iq  = NULL,
00137     *Uy_iq  = NULL;
00138   IterQuantityAccess<MatrixOpNonsing>
00139     *R_iq   = NULL;
00140 
00141   MatrixOp::EMatNormType mat_nrm_inf = MatrixOp::MAT_NORM_INF;
00142   const bool calc_matrix_norms = algo.algo_cntr().calc_matrix_norms();
00143   const bool calc_matrix_info_null_space_only = algo.algo_cntr().calc_matrix_info_null_space_only();
00144   
00145   if( x_iq.last_updated() == IterQuantity::NONE_UPDATED ) {
00146     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00147       out << "\nx is not updated for any k so set x_k = nlp.xinit() ...\n";
00148     }
00149     x_iq.set_k(0) = nlp.xinit();
00150   }
00151   
00152   // Validate x
00153   if( nb && algo.algo_cntr().check_results() ) {
00154     assert_print_nan_inf(
00155       x_iq.get_k(0), "x_k", true
00156       , int(olevel) >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
00157       );
00158     if( nlp.num_bounded_x() > 0 ) {
00159       if(!bounds_tester().check_in_bounds(
00160            int(olevel)  >= int(PRINT_ALGORITHM_STEPS) ? &out : NULL
00161            ,int(olevel) >= int(PRINT_VECTORS)                // print_all_warnings
00162            ,int(olevel) >= int(PRINT_ITERATION_QUANTITIES)  // print_vectors
00163            ,nlp.xl(),        "xl"
00164            ,nlp.xu(),        "xu"
00165            ,x_iq.get_k(0),   "x_k"
00166            ))
00167       {
00168         TEUCHOS_TEST_FOR_EXCEPTION(
00169           true, TestFailed
00170           ,"EvalNewPointStd_Step::do_step(...) : Error, "
00171           "the variables bounds xl <= x_k <= xu where violated!" );
00172       }
00173     }
00174   }
00175 
00176   Vector &x = x_iq.get_k(0);
00177 
00178   Range1D  var_dep(Range1D::INVALID), var_indep(Range1D::INVALID);
00179   if( s.get_decomp_sys().get() ) {
00180     const ConstrainedOptPack::DecompositionSystemVarReduct
00181       *decomp_sys_vr = dynamic_cast<ConstrainedOptPack::DecompositionSystemVarReduct*>(&s.decomp_sys());
00182     if(decomp_sys_vr) {
00183       var_dep   = decomp_sys_vr->var_dep();
00184       var_indep = decomp_sys_vr->var_indep();
00185     }
00186     s.var_dep(var_dep);
00187     s.var_indep(var_indep);
00188   }
00189 
00190   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00191     out << "\n||x_k||inf            = " << x.norm_inf();
00192     if( var_dep.size() )
00193       out << "\n||x(var_dep)_k||inf   = " << x.sub_view(var_dep)->norm_inf();
00194     if( var_indep.size() )
00195       out << "\n||x(var_indep)_k||inf = " << x.sub_view(var_indep)->norm_inf();
00196     out << std::endl;
00197   }
00198   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00199     out << "\nx_k = \n" << x;
00200     if( var_dep.size() )
00201       out << "\nx(var_dep)_k = \n" << *x.sub_view(var_dep);
00202   }
00203   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00204     if( var_indep.size() )
00205       out << "\nx(var_indep)_k = \n" << *x.sub_view(var_indep);
00206   }
00207 
00208   // Set the references to the current point's quantities to be updated
00209   const bool f_k_updated  = f_iq.updated_k(0);
00210   const bool Gf_k_updated = Gf_iq.updated_k(0);
00211   const bool c_k_updated  = m  > 0 ? c_iq->updated_k(0)  : false;
00212   const bool Gc_k_updated = m  > 0 ? Gc_iq->updated_k(0) : false;
00213   nlp.unset_quantities();
00214   if(!f_k_updated) nlp.set_f( &f_iq.set_k(0) );
00215   if(!Gf_k_updated) nlp.set_Gf( &Gf_iq.set_k(0) );
00216   if( m > 0 ) {
00217     if(!c_k_updated) nlp.set_c( &c_iq->set_k(0) );
00218     if(!Gc_k_updated) nlp.set_Gc( &Gc_iq->set_k(0) );
00219   }
00220 
00221   // Calculate Gc at x_k
00222   bool new_point = true;
00223   if(m > 0) {
00224     if(!Gc_k_updated) nlp.calc_Gc( x, new_point );
00225     new_point = false;
00226   }
00227 
00228   //
00229   // Update (or select a new) range/null decomposition
00230   //
00231   bool new_decomp_selected = false;
00232   if( m > 0 ) {
00233     
00234     // Update the range/null decomposition
00235     decomp_sys_handler().update_decomposition(
00236       algo, s, nlp, decomp_sys_testing(), decomp_sys_testing_print_level()
00237       ,&new_decomp_selected
00238       );
00239 
00240     r  = s.equ_decomp().size();
00241 
00242     Z_iq   = ( n > m && r > 0 )      ? &s.Z()  : NULL;
00243     Y_iq   = ( r > 0 )               ? &s.Y()  : NULL;
00244     Uz_iq  = ( m  > 0 && m  > r )    ? &s.Uz() : NULL;
00245     Uy_iq  = ( m  > 0 && m  > r )    ? &s.Uy() : NULL;
00246     R_iq   = ( m > 0 )               ? &s.R()  : NULL;
00247 
00248     // Determine if we will test the decomp_sys or not
00249     const DecompositionSystem::ERunTests
00250       ds_test_what = ( ( decomp_sys_testing() == DecompositionSystemHandler_Strategy::DST_TEST
00251                  || ( decomp_sys_testing() == DecompositionSystemHandler_Strategy::DST_DEFAULT
00252                   && algo.algo_cntr().check_results() ) )
00253                ? DecompositionSystem::RUN_TESTS
00254                : DecompositionSystem::NO_TESTS );
00255     
00256     // Determine the output level for decomp_sys        
00257     DecompositionSystem::EOutputLevel ds_olevel;
00258     switch(olevel) {
00259       case PRINT_NOTHING:
00260       case PRINT_BASIC_ALGORITHM_INFO:
00261         ds_olevel = DecompositionSystem::PRINT_NONE;
00262         break;
00263       case PRINT_ALGORITHM_STEPS:
00264       case PRINT_ACTIVE_SET:
00265         ds_olevel = DecompositionSystem::PRINT_BASIC_INFO;
00266         break;
00267       case PRINT_VECTORS:
00268         ds_olevel = DecompositionSystem::PRINT_VECTORS;
00269         break;
00270       case PRINT_ITERATION_QUANTITIES:
00271         ds_olevel = DecompositionSystem::PRINT_EVERY_THING;
00272         break;
00273       default:
00274         TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
00275     };
00276 
00277     // Test the decomposition system
00278     if( ds_test_what == DecompositionSystem::RUN_TESTS ) {
00279       // Set the output level
00280       if( decomp_sys_tester().print_tests() == DecompositionSystemTester::PRINT_NOT_SELECTED ) {
00281         DecompositionSystemTester::EPrintTestLevel  ds_olevel;
00282         switch(olevel) {
00283           case PRINT_NOTHING:
00284           case PRINT_BASIC_ALGORITHM_INFO:
00285             ds_olevel = DecompositionSystemTester::PRINT_NONE;
00286             break;
00287           case PRINT_ALGORITHM_STEPS:
00288           case PRINT_ACTIVE_SET:
00289             ds_olevel = DecompositionSystemTester::PRINT_BASIC;
00290             break;
00291           case PRINT_VECTORS:
00292             ds_olevel = DecompositionSystemTester::PRINT_MORE;
00293             break;
00294           case PRINT_ITERATION_QUANTITIES:
00295             ds_olevel = DecompositionSystemTester::PRINT_ALL;
00296             break;
00297           default:
00298             TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
00299         }
00300         decomp_sys_tester().print_tests(ds_olevel);
00301         decomp_sys_tester().dump_all( olevel == PRINT_ITERATION_QUANTITIES );
00302       }
00303       // Run the tests
00304       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00305         out << "\nTesting the range/null decompostion ...\n";
00306       }
00307       const bool
00308         decomp_sys_passed = decomp_sys_tester().test_decomp_system(
00309           s.decomp_sys()
00310           ,Gc_iq->get_k(0)                   // Gc
00311           ,Z_iq ? &Z_iq->get_k(0) : NULL     // Z
00312           ,&Y_iq->get_k(0)                   // Y
00313           ,&R_iq->get_k(0)                   // R
00314           ,m > r  ? &Uz_iq->get_k(0) : NULL  // Uz
00315           ,m > r  ? &Uy_iq->get_k(0) : NULL  // Uy
00316           ,( olevel >= PRINT_ALGORITHM_STEPS ) ? &out : NULL
00317           );
00318       TEUCHOS_TEST_FOR_EXCEPTION(
00319         !decomp_sys_passed, TestFailed
00320         ,"EvalNewPointStd_Step::do_step(...) : Error, "
00321         "the tests of the decomposition system failed!" );
00322     }
00323   }
00324   else {
00325     // Unconstrained problem
00326     Z_iq = &s.Z();
00327     dyn_cast<MatrixSymIdent>(Z_iq->set_k(0)).initialize( nlp.space_x() );
00328     s.equ_decomp(Range1D::Invalid);
00329     s.equ_undecomp(Range1D::Invalid);
00330   }
00331 
00332   // Calculate the rest of the quantities.  If decomp_sys is a variable
00333   // reduction decomposition system object, then nlp will be hip to the
00334   // basis selection and will permute these quantities to that basis.
00335   // Note that x will already be permuted to the current basis.
00336   if(!Gf_k_updated) { nlp.calc_Gf( x, new_point ); new_point = false; }
00337   if( m && (!c_k_updated || new_decomp_selected ) ) {
00338     if(c_k_updated) nlp.set_c( &c_iq->set_k(0) ); // This was not set earlier!
00339     nlp.calc_c( x, false);
00340   }
00341   if(!f_k_updated) {
00342     nlp.calc_f( x, false);
00343   }
00344   nlp.unset_quantities();
00345   
00346   // Check for NaN and Inf
00347   assert_print_nan_inf(f_iq.get_k(0),"f_k",true,&out); 
00348   if(m)
00349     assert_print_nan_inf(c_iq->get_k(0),"c_k",true,&out); 
00350   assert_print_nan_inf(Gf_iq.get_k(0),"Gf_k",true,&out); 
00351   
00352   // Print the iteration quantities before we test the derivatives for debugging
00353 
00354   // Update the selection of dependent and independent variables
00355   if( s.get_decomp_sys().get() ) {
00356     const ConstrainedOptPack::DecompositionSystemVarReduct
00357       *decomp_sys_vr = dynamic_cast<ConstrainedOptPack::DecompositionSystemVarReduct*>(&s.decomp_sys());
00358     if(decomp_sys_vr) {
00359       var_dep   = decomp_sys_vr->var_dep();
00360       var_indep = decomp_sys_vr->var_indep();
00361     }
00362   }
00363   
00364   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00365     out << "\nPrinting the updated iteration quantities ...\n";
00366   }
00367 
00368   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00369     out << "\nf_k                      = "     << f_iq.get_k(0);
00370     out << "\n||Gf_k||inf              = "     << Gf_iq.get_k(0).norm_inf();
00371     if( var_dep.size() )
00372       out << "\n||Gf_k(var_dep)_k||inf   = " << Gf_iq.get_k(0).sub_view(var_dep)->norm_inf();
00373     if( var_indep.size() )
00374       out << "\n||Gf_k(var_indep)_k||inf = " << Gf_iq.get_k(0).sub_view(var_indep)->norm_inf();
00375     if(m) {
00376       out << "\n||c_k||inf               = " << c_iq->get_k(0).norm_inf();
00377       if( calc_matrix_norms && !calc_matrix_info_null_space_only )
00378         out << "\n||Gc_k||inf              = " << Gc_iq->get_k(0).calc_norm(mat_nrm_inf).value;
00379       if( n > r && calc_matrix_norms && !calc_matrix_info_null_space_only )
00380         out << "\n||Z||inf                 = " << Z_iq->get_k(0).calc_norm(mat_nrm_inf).value;
00381       if( r && calc_matrix_norms && !calc_matrix_info_null_space_only )
00382         out << "\n||Y||inf                 = " << Y_iq->get_k(0).calc_norm(mat_nrm_inf).value;
00383       if( r && calc_matrix_norms && !calc_matrix_info_null_space_only  )
00384         out << "\n||R||inf                 = " << R_iq->get_k(0).calc_norm(mat_nrm_inf).value;
00385       if( algo.algo_cntr().calc_conditioning() && !calc_matrix_info_null_space_only ) {
00386         out << "\ncond_inf(R)              = " << R_iq->get_k(0).calc_cond_num(mat_nrm_inf).value;
00387       }
00388       if( m > r && calc_matrix_norms && !calc_matrix_info_null_space_only ) {
00389         out << "\n||Uz_k||inf              = " << Uz_iq->get_k(0).calc_norm(mat_nrm_inf).value;
00390         out << "\n||Uy_k||inf              = " << Uy_iq->get_k(0).calc_norm(mat_nrm_inf).value;
00391       }
00392     }
00393     out << std::endl;
00394   }
00395 
00396   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ITERATION_QUANTITIES) ) {
00397     if(m)
00398       out << "\nGc_k =\n" << Gc_iq->get_k(0);
00399     if( n > r )
00400       out << "\nZ_k =\n" << Z_iq->get_k(0);
00401     if(r) {
00402       out << "\nY_k =\n" << Y_iq->get_k(0);
00403       out << "\nR_k =\n" << R_iq->get_k(0);
00404     }
00405     if( m > r ) {
00406       out << "\nUz_k =\n" << Uz_iq->get_k(0);
00407       out << "\nUy_k =\n" << Uy_iq->get_k(0);
00408     }
00409   }
00410   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00411     out << "\nGf_k =\n" << Gf_iq.get_k(0);
00412     if( var_dep.size() )
00413       out << "\nGf(var_dep)_k =\n " << *Gf_iq.get_k(0).sub_view(var_dep);
00414   }
00415   if( static_cast<int>(ns_olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00416     if( var_indep.size() )
00417       out << "\nGf(var_indep)_k =\n" << *Gf_iq.get_k(0).sub_view(var_indep);
00418   }
00419   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) {
00420     if(m)
00421       out << "\nc_k = \n" << c_iq->get_k(0);
00422   }
00423   
00424   // Check the derivatives if we are checking the results
00425   if(   fd_deriv_testing() == FD_TEST
00426     || ( fd_deriv_testing() == FD_DEFAULT && algo.algo_cntr().check_results() )  )
00427   {
00428     
00429     if( olevel >= PRINT_ALGORITHM_STEPS ) {
00430       out << "\n*** Checking derivatives by finite differences\n";
00431     }
00432 
00433     const bool
00434       nlp_passed = deriv_tester().finite_diff_check(
00435         &nlp
00436         ,x
00437         ,nb ? &nlp.xl() : NULL
00438         ,nb ? &nlp.xu() : NULL
00439         ,m  ? &Gc_iq->get_k(0) : NULL
00440         ,&Gf_iq.get_k(0)
00441         ,olevel >= PRINT_VECTORS
00442         ,( olevel >= PRINT_ALGORITHM_STEPS ) ? &out : NULL
00443         );
00444     TEUCHOS_TEST_FOR_EXCEPTION(
00445       !nlp_passed, TestFailed
00446       ,"EvalNewPointStd_Step::do_step(...) : Error, "
00447       "the tests of the nlp derivatives failed!" );
00448   }
00449 
00450   return true;
00451 }
00452 
00453 void EvalNewPointStd_Step::print_step(
00454    const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00455   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00456   ) const
00457 {
00458   const NLPAlgo       &algo = rsqp_algo(_algo);
00459   const NLPAlgoState  &s    = algo.rsqp_state();
00460   const NLP           &nlp  = algo.nlp();
00461   const size_type
00462     m = nlp.m();
00463   out
00464     << L << "*** Evaluate the new point and update the range/null decomposition\n"
00465     << L << "if nlp is not initialized then initialize the nlp\n"
00466     << L << "if x is not updated for any k then set x_k = xinit\n";
00467   if(m) {
00468     out
00469       << L << "if Gc_k is not updated Gc_k = Gc(x_k) <: space_x|space_c\n"
00470       << L << "For Gc_k = [ Gc_k(:,equ_decomp), Gc_k(:,equ_undecomp) ] where:\n"
00471       << L << "  Gc_k(:,equ_decomp) <: space_x|space_c(equ_decomp) has full column rank r\n"
00472       << L << "Find:\n"
00473       << L << "  Z_k  <: space_x|space_null    s.t. Gc_k(:,equ_decomp)' * Z_k = 0\n"
00474       << L << "  Y_k  <: space_x|space_range   s.t. [Z_k Y_k] is nonsigular \n"
00475       << L << "  R_k  <: space_c(equ_decomp)|space_range\n"
00476       << L << "                                s.t. R_k = Gc_k(:,equ_decomp)' * Y_k\n"
00477       << L << "  if m > r : Uz_k <: space_c(equ_undecomp)|space_null\n"
00478       << L << "                                s.t. Uz_k = Gc_k(:,equ_undecomp)' * Z_k\n"
00479       << L << "  if m > r : Uy_k <: space_c(equ_undecomp)|space_range\n"
00480       << L << "                                s.t. Uy_k = Gc_k(:,equ_undecomp)' * Y_k\n"
00481       << L << "begin update decomposition (class \'" << typeName(decomp_sys_handler()) << "\')\n"
00482       ;
00483     decomp_sys_handler().print_update_decomposition( algo, s, out, L + "  " );
00484     out
00485       << L << "end update decomposition\n"
00486       << L << "if ( (decomp_sys_testing==DST_TEST)\n"
00487       << L << "  or (decomp_sys_testing==DST_DEFAULT and check_results==true)\n"
00488       << L << "  ) then\n"
00489       << L << "  check properties for Z_k, Y_k, R_k, Uz_k and Uy_k\n"
00490       << L << "end\n"
00491       ;
00492   }
00493   else {
00494     out 
00495       << L << "Z_k = eye(space_x)\n";
00496   }
00497   if(m) {
00498     out
00499       << L << "end\n";
00500   }
00501   out
00502     << L << "Gf_k = Gf(x_k) <: space_x\n"
00503     << L << "if m > 0 and c_k is not updated c_k = c(x_k) <: space_c\n"
00504     << L << "if f_k is not updated f_k = f(x_k) <: REAL\n"
00505     << L << "if ( (fd_deriv_testing==FD_TEST)\n"
00506     << L << "  or (fd_deriv_testing==FD_DEFAULT and check_results==true)\n"
00507     << L << "  ) then\n"
00508     << L << "  check Gc_k (if m > 0) and Gf_k by finite differences.\n"
00509     << L << "end\n"
00510     ;
00511 }
00512 
00513 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines