MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_CheckConvergenceIP_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 #include <assert.h>
00043 
00044 #include <ostream>
00045 //#include <limits>
00046 //#include <sstream>
00047 
00048 #include "MoochoPack_CheckConvergenceIP_Strategy.hpp"
00049 #include "MoochoPack_IpState.hpp"
00050 #include "MoochoPack_moocho_algo_conversion.hpp"
00051 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00052 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00053 #include "AbstractLinAlgPack_VectorOut.hpp"
00054 #include "IterationPack_print_algorithm_step.hpp"
00055 #include "Teuchos_dyn_cast.hpp"
00056 
00057 namespace MoochoPack {
00058 
00059 CheckConvergenceIP_Strategy::CheckConvergenceIP_Strategy(
00060   EOptErrorCheck         opt_error_check
00061   ,EScaleKKTErrorBy      scale_opt_error_by
00062   ,EScaleKKTErrorBy      scale_feas_error_by
00063   ,EScaleKKTErrorBy      scale_comp_error_by
00064   ,bool                  scale_opt_error_by_Gf
00065   )
00066   :
00067   CheckConvergenceStd_Strategy(
00068     opt_error_check,
00069     scale_opt_error_by,
00070     scale_feas_error_by,
00071     scale_comp_error_by,
00072     scale_opt_error_by_Gf
00073     )
00074   {}
00075 
00076 bool CheckConvergenceIP_Strategy::Converged(
00077   Algorithm& _algo
00078   )
00079   {
00080   using Teuchos::dyn_cast;
00081   using AbstractLinAlgPack::num_bounded;
00082   using AbstractLinAlgPack::IP_comp_err_with_mu;
00083 
00084   // Calculate kkt errors and check for overall convergence
00085   //bool found_solution = CheckConvergenceStd_Strategy::Converged(_algo);
00086   bool found_solution = false;
00087 
00088   // Recalculate the complementarity error including mu
00089   
00090   // Get the iteration quantities
00091   IpState &s = dyn_cast<IpState>(*_algo.get_state());
00092   NLPAlgo& algo = rsqp_algo(_algo);
00093   NLP& nlp = algo.nlp();
00094   
00095   EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
00096   std::ostream& out = algo.track().journal_out();
00097 
00098   // Get necessary iteration quantities
00099   const value_type &mu_km1 = s.barrier_parameter().get_k(-1);
00100   const Vector& x_k = s.x().get_k(0);
00101   const VectorMutable& Gf_k = s.Gf().get_k(0);
00102   const Vector& rGL_k = s.rGL().get_k(0);
00103   const Vector& c_k = s.c().get_k(0);
00104   const Vector& vl_k = s.Vl().get_k(0).diag();
00105   const Vector& vu_k = s.Vu().get_k(0).diag();
00106   
00107   // Calculate the errors with Andreas' scaling
00108   value_type& opt_err = s.opt_kkt_err().set_k(0);
00109   value_type& feas_err = s.feas_kkt_err().set_k(0);
00110   value_type& comp_err = s.comp_kkt_err().set_k(0);
00111   value_type& comp_err_mu = s.comp_err_mu().set_k(0);
00112 
00113   // scaling
00114   value_type scale_1 = 1 + x_k.norm_1()/x_k.dim();
00115 
00116   Teuchos::RCP<VectorMutable> temp = Gf_k.clone();
00117   temp->axpy(-1.0, vl_k);
00118   temp->axpy(1.0, vu_k);
00119   value_type scale_2 = temp->norm_1();
00120   scale_2 += vl_k.norm_1() + vu_k.norm_1();
00121 
00122   *temp = nlp.infinite_bound();
00123   const size_type nlb = num_bounded(nlp.xl(), *temp, nlp.infinite_bound());
00124   *temp = -nlp.infinite_bound();
00125   const size_type nub = num_bounded(*temp, nlp.xu(), nlp.infinite_bound());
00126   scale_2 = 1 + scale_2/(1+nlp.m()+nlb+nub);
00127 
00128   // Calculate the opt_err
00129   opt_err = rGL_k.norm_inf() / scale_2;
00130 
00131   // Calculate the feas_err
00132   feas_err = c_k.norm_inf() / scale_1;
00133   
00134   // Calculate the comp_err
00135   if( (int)olevel >= (int)PRINT_VECTORS )
00136     {
00137     out << "\nx =\n"    << s.x().get_k(0);
00138     out << "\nxl =\n"   << nlp.xl();
00139     out << "\nvl =\n"   << s.Vl().get_k(0).diag();
00140     out << "\nxu =\n"   << nlp.xu();
00141     out << "\nvu =\n"   << s.Vu().get_k(0).diag();
00142     }
00143 
00144   comp_err = IP_comp_err_with_mu(
00145     0.0, nlp.infinite_bound(), s.x().get_k(0), nlp.xl(), nlp.xu()
00146     ,s.Vl().get_k(0).diag(), s.Vu().get_k(0).diag());
00147 
00148   comp_err_mu = IP_comp_err_with_mu(
00149     mu_km1, nlp.infinite_bound(), s.x().get_k(0), nlp.xl(), nlp.xu()
00150     ,s.Vl().get_k(0).diag(), s.Vu().get_k(0).diag());
00151 
00152   comp_err = comp_err / scale_2;
00153   comp_err_mu = comp_err_mu / scale_2;
00154 
00155   // check for convergence
00156   
00157   const value_type opt_tol = algo.algo_cntr().opt_tol();
00158   const value_type feas_tol = algo.algo_cntr().feas_tol();
00159   const value_type comp_tol = algo.algo_cntr().comp_tol();
00160 
00161   if (opt_err < opt_tol && feas_err < feas_tol && comp_err < comp_tol)
00162     {
00163     found_solution = true;
00164     }
00165 
00166   if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) || (int(olevel) >= int(PRINT_BASIC_ALGORITHM_INFO) && found_solution) )
00167     {
00168     out 
00169       << "\nopt_kkt_err_k   = " << opt_err << ( opt_err < opt_tol ? " < " : " > " )
00170       << "opt_tol = " << opt_tol
00171       << "\nfeas_kkt_err_k   = " << feas_err << ( feas_err < feas_tol ? " < " : " > " )
00172       << "feas_tol = " << feas_tol
00173       << "\ncomp_kkt_err_k   = " << comp_err << ( comp_err < comp_tol ? " < " : " > " )
00174       << "comp_tol = " << comp_tol
00175       << "\ncomp_err_mu      = " << comp_err_mu
00176       << "\nbarrier_parameter_k (mu_km1) = " << mu_km1
00177       << "comp_tol = " << comp_tol
00178       << std::endl;
00179     }
00180     
00181   return found_solution;
00182   }
00183 
00184 void CheckConvergenceIP_Strategy::print_step( const Algorithm& _algo, std::ostream& out, const std::string& L ) const
00185   {
00186   out 
00187     << L << "CheckConvergenceIP_Strategy\n"
00188     << L << " IP_found_solution = CheckConvergedStd_Strategy::Converged(_algo, reportFinalSolution)\n";
00189   
00190   CheckConvergenceStd_Strategy::print_step(_algo, out, L+"   ");
00191   
00192   out 
00193     << L << "*** recalculate comp_err\n"
00194     << L << "comp_err_k = 0.0"
00195     << L << "for all i = 1 to n\n"
00196     << L << "   comp_err_k = max( comp_err_k, vl_k(i)*(x_k(i)-xl_k(i))-mu_km1, vu_k(i)*(xu_k(i)-x_k(i))-mu_k )\n"
00197     << L << "next i\n"
00198     << L << "if IP_found_solution then\n"
00199     << L << "   IP_found_solution = false\n"
00200     << L << "   if (comp_err_k < comp_tol && mu_k < comp_tol) then\n"
00201     << L << "      IP_found_solution = true\n"
00202     << L << "   end\n"
00203     << L << "end\n"
00204     << L << "return IP_found_solution\n";
00205   }
00206 
00207 } // end namespace MoochoPack
00208 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines