MoochoPack_CheckConvergenceIP_Strategy.cpp

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

Generated on Tue Jul 13 09:29:31 2010 for MoochoPack : Framework for Large-Scale Optimization Algorithms by  doxygen 1.4.7