MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_PostEvalNewPointBarrier_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 // 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 #include <ostream>
00030 #include <typeinfo>
00031 #include <iostream>
00032 
00033 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00034 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00035 #include "AbstractLinAlgPack_assert_print_nan_inf.hpp"
00036 #include "AbstractLinAlgPack_VectorOut.hpp"
00037 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00038 #include "IterationPack_print_algorithm_step.hpp"
00039 #include "NLPInterfacePack_NLPFirstOrder.hpp"
00040 #include "MoochoPack_IpState.hpp"
00041 #include "MoochoPack_PostEvalNewPointBarrier_Step.hpp"
00042 #include "MoochoPack_moocho_algo_conversion.hpp"
00043 
00044 #include "OptionsFromStreamPack_StringToBool.hpp"
00045 
00046 #include "Teuchos_dyn_cast.hpp"
00047 
00048 namespace MoochoPack {
00049 
00050 bool PostEvalNewPointBarrier_Step::do_step(
00051   Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00052   ,poss_type assoc_step_poss
00053   )
00054   {
00055   using Teuchos::dyn_cast;
00056   using IterationPack::print_algorithm_step;
00057   using AbstractLinAlgPack::inv_of_difference;
00058   using AbstractLinAlgPack::correct_upper_bound_multipliers;
00059   using AbstractLinAlgPack::correct_lower_bound_multipliers;
00060   using LinAlgOpPack::Vp_StV;
00061 
00062   NLPAlgo            &algo   = dyn_cast<NLPAlgo>(_algo);
00063   IpState             &s      = dyn_cast<IpState>(_algo.state());
00064   NLP                 &nlp    = algo.nlp();
00065   
00066   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00067   std::ostream& out = algo.track().journal_out();
00068   
00069   if(!nlp.is_initialized())
00070     nlp.initialize(algo.algo_cntr().check_results());
00071 
00072   // print step header.
00073   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00074     {
00075     using IterationPack::print_algorithm_step;
00076     print_algorithm_step( _algo, step_poss, type, assoc_step_poss, out );
00077     }
00078 
00079   IterQuantityAccess<VectorMutable> &x_iq = s.x();
00080   IterQuantityAccess<MatrixSymDiagStd> &Vl_iq = s.Vl();
00081   IterQuantityAccess<MatrixSymDiagStd> &Vu_iq = s.Vu();
00082 
00084   // Calculate invXl = diag(1/(x-xl)) 
00085   //  and invXu = diag(1/(xu-x)) matrices
00087 
00088   // get references to x, invXl, and invXu
00089   VectorMutable& x = x_iq.get_k(0);
00090   MatrixSymDiagStd& invXu = s.invXu().set_k(0);
00091   MatrixSymDiagStd& invXl = s.invXl().set_k(0);
00092   
00093   //std::cout << "xu=\n";
00094   //nlp.xu().output(std::cout);
00095 
00096   inv_of_difference(1.0, nlp.xu(), x, &invXu.diag());
00097   inv_of_difference(1.0, x, nlp.xl(), &invXl.diag());
00098 
00099   //std::cout << "invXu=\v";
00100   //invXu.output(std::cout);
00101 
00102   //std::cout << "\ninvXl=\v";
00103   //invXl.output(std::cout);
00104   
00105   // Check for divide by zeros - 
00106     using AbstractLinAlgPack::assert_print_nan_inf;
00107     assert_print_nan_inf(invXu.diag(), "invXu", true, &out); 
00108     assert_print_nan_inf(invXl.diag(), "invXl", true, &out); 
00109   // These should never go negative either - could be a useful check
00110 
00111   // Initialize Vu and Vl if necessary
00112   if (  Vu_iq.last_updated() == IterQuantity::NONE_UPDATED )
00113     {
00114     VectorMutable& vu = Vu_iq.set_k(0).diag();    
00115     vu = 0;
00116     Vp_StV(&vu, s.barrier_parameter().get_k(-1), invXu.diag());
00117 
00118     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00119       {
00120       out << "\nInitialize Vu with barrier_parameter * invXu ...\n";
00121       }
00122     }
00123 
00124 if (  Vl_iq.last_updated() == IterQuantity::NONE_UPDATED  )
00125     {
00126     VectorMutable& vl = Vl_iq.set_k(0).diag();
00127     vl = 0;
00128     Vp_StV(&vl, s.barrier_parameter().get_k(-1), invXl.diag());
00129 
00130     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) 
00131       {
00132       out << "\nInitialize Vl with barrier_parameter * invXl ...\n";
00133       }
00134     }
00135 
00136   if (s.num_basis().updated_k(0))
00137     {
00138     // Basis changed, reorder Vl and Vu
00139     if (Vu_iq.updated_k(-1))
00140       { Vu_iq.set_k(0,-1); }
00141     if (Vl_iq.updated_k(-1))
00142       { Vl_iq.set_k(0,-1); }
00143       
00144     VectorMutable& vu = Vu_iq.set_k(0).diag();
00145     VectorMutable& vl = Vl_iq.set_k(0).diag();
00146 
00147     s.P_var_last().permute( BLAS_Cpp::trans, &vu ); // Permute back to original order
00148     s.P_var_last().permute( BLAS_Cpp::trans, &vl ); // Permute back to original order
00149 
00150     if( (int)olevel >= (int)PRINT_VECTORS ) 
00151       {
00152       out << "\nx resorted vl and vu to the original order\n" << x;
00153       }
00154 
00155     s.P_var_current().permute( BLAS_Cpp::no_trans, &vu ); // Permute to new (current) order
00156     s.P_var_current().permute( BLAS_Cpp::no_trans, &vl ); // Permute to new (current) order
00157 
00158     if( (int)olevel >= (int)PRINT_VECTORS ) 
00159       {
00160       out << "\nx resorted to new basis \n" << x;
00161       }
00162     }
00163 
00164   correct_upper_bound_multipliers(nlp.xu(), +NLP::infinite_bound(), &Vu_iq.get_k(0).diag());
00165   correct_lower_bound_multipliers(nlp.xl(), -NLP::infinite_bound(), &Vl_iq.get_k(0).diag());
00166   
00167   if( (int)olevel >= (int)PRINT_VECTORS ) 
00168     {
00169     out << "x=\n"  << s.x().get_k(0);
00170     out << "xl=\n" << nlp.xl();
00171     out << "vl=\n" << s.Vl().get_k(0).diag();
00172     out << "xu=\n" << nlp.xu();
00173     out << "vu=\n" << s.Vu().get_k(0).diag();
00174     }
00175 
00176   // Update general algorithm bound multipliers
00177   VectorMutable& v = s.nu().set_k(0);
00178   v = Vu_iq.get_k(0).diag();
00179   Vp_StV(&v,-1.0,Vl_iq.get_k(0).diag());
00180 
00181   // Print vector information
00182   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_VECTORS) ) 
00183     {
00184     out << "invXu_k.diag()=\n" << invXu.diag();
00185     out << "invXl_k.diag()=\n" << invXl.diag();
00186     out << "Vu_k.diag()=\n"    << Vu_iq.get_k(0).diag();
00187     out << "Vl_k.diag()=\n"    << Vl_iq.get_k(0).diag();
00188     out << "nu_k=\n"           << s.nu().get_k(0);
00189     }
00190 
00191   return true;
00192   }
00193 
00194 
00195 void PostEvalNewPointBarrier_Step::print_step(
00196   const Algorithm& _algo, poss_type step_poss, IterationPack::EDoStepType type
00197   ,poss_type assoc_step_poss, std::ostream& out, const std::string& L
00198   ) const
00199   {
00200   //const NLPAlgo   &algo = rsqp_algo(_algo);
00201   //const NLPAlgoState  &s    = algo.rsqp_state();
00202   out << L << "# Evaluate information specific to primal / dual barrier algorithms (Post EvalNewPoint)\n"
00203     << L << "invXl_k = diag(i, 1/(x(i)-xl))"
00204     << L << "invXu_k = diag(i, 1/(xu-x(i)))\n"
00205     << L << "if (Vu_k not updated) then\n"
00206     << L << "   Vu_k = mu_k*invXu_k\n"
00207     << L << "end\n"
00208     << L << "if (Vl_k not updated) then\n"
00209     << L << "   Vl_k = mu_k*invXl_k\n"
00210     << L << "end\n"
00211     << L << "nu_k_k = Vu_k.diag() - Vl_k.diag()\n";
00212   }
00213 
00214 } // end namespace MoochoPack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines