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