MoochoPack_CheckConvergenceIP_Strategy.cpp

// @HEADER
// ***********************************************************************
// 
// Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
//                  Copyright (2003) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#include <assert.h>

#include <ostream>
//#include <limits>
//#include <sstream>

#include "MoochoPack_CheckConvergenceIP_Strategy.hpp"
#include "MoochoPack_IpState.hpp"
#include "MoochoPack_moocho_algo_conversion.hpp"
#include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
#include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
#include "AbstractLinAlgPack_VectorOut.hpp"
#include "IterationPack_print_algorithm_step.hpp"
#include "Teuchos_dyn_cast.hpp"

namespace MoochoPack {

CheckConvergenceIP_Strategy::CheckConvergenceIP_Strategy(
  EOptErrorCheck         opt_error_check
  ,EScaleKKTErrorBy      scale_opt_error_by
  ,EScaleKKTErrorBy      scale_feas_error_by
  ,EScaleKKTErrorBy      scale_comp_error_by
  ,bool                  scale_opt_error_by_Gf
  )
  :
  CheckConvergenceStd_Strategy(
    opt_error_check,
    scale_opt_error_by,
    scale_feas_error_by,
    scale_comp_error_by,
    scale_opt_error_by_Gf
    )
  {}

bool CheckConvergenceIP_Strategy::Converged(
  Algorithm& _algo
  )
  {
  using Teuchos::dyn_cast;
  using AbstractLinAlgPack::num_bounded;
  using AbstractLinAlgPack::IP_comp_err_with_mu;

  // Calculate kkt errors and check for overall convergence
  //bool found_solution = CheckConvergenceStd_Strategy::Converged(_algo);
  bool found_solution = false;

  // Recalculate the complementarity error including mu
  
  // Get the iteration quantities
  IpState &s = dyn_cast<IpState>(*_algo.get_state());
  NLPAlgo& algo = rsqp_algo(_algo);
  NLP& nlp = algo.nlp();
  
  EJournalOutputLevel olevel = algo.algo_cntr().journal_output_level();
  std::ostream& out = algo.track().journal_out();

  // Get necessary iteration quantities
  const value_type &mu_km1 = s.barrier_parameter().get_k(-1);
  const Vector& x_k = s.x().get_k(0);
  const VectorMutable& Gf_k = s.Gf().get_k(0);
  const Vector& rGL_k = s.rGL().get_k(0);
  const Vector& c_k = s.c().get_k(0);
  const Vector& vl_k = s.Vl().get_k(0).diag();
  const Vector& vu_k = s.Vu().get_k(0).diag();
  
  // Calculate the errors with Andreas' scaling
  value_type& opt_err = s.opt_kkt_err().set_k(0);
  value_type& feas_err = s.feas_kkt_err().set_k(0);
  value_type& comp_err = s.comp_kkt_err().set_k(0);
  value_type& comp_err_mu = s.comp_err_mu().set_k(0);

  // scaling
  value_type scale_1 = 1 + x_k.norm_1()/x_k.dim();

  Teuchos::RCP<VectorMutable> temp = Gf_k.clone();
  temp->axpy(-1.0, vl_k);
  temp->axpy(1.0, vu_k);
  value_type scale_2 = temp->norm_1();
  scale_2 += vl_k.norm_1() + vu_k.norm_1();

  *temp = nlp.infinite_bound();
  const size_type nlb = num_bounded(nlp.xl(), *temp, nlp.infinite_bound());
  *temp = -nlp.infinite_bound();
  const size_type nub = num_bounded(*temp, nlp.xu(), nlp.infinite_bound());
  scale_2 = 1 + scale_2/(1+nlp.m()+nlb+nub);

  // Calculate the opt_err
  opt_err = rGL_k.norm_inf() / scale_2;

  // Calculate the feas_err
  feas_err = c_k.norm_inf() / scale_1;
  
  // Calculate the comp_err
  if( (int)olevel >= (int)PRINT_VECTORS )
    {
    out << "\nx =\n"    << s.x().get_k(0);
    out << "\nxl =\n"   << nlp.xl();
    out << "\nvl =\n"   << s.Vl().get_k(0).diag();
    out << "\nxu =\n"   << nlp.xu();
    out << "\nvu =\n"   << s.Vu().get_k(0).diag();
    }

  comp_err = IP_comp_err_with_mu(
    0.0, nlp.infinite_bound(), s.x().get_k(0), nlp.xl(), nlp.xu()
    ,s.Vl().get_k(0).diag(), s.Vu().get_k(0).diag());

  comp_err_mu = IP_comp_err_with_mu(
    mu_km1, nlp.infinite_bound(), s.x().get_k(0), nlp.xl(), nlp.xu()
    ,s.Vl().get_k(0).diag(), s.Vu().get_k(0).diag());

  comp_err = comp_err / scale_2;
  comp_err_mu = comp_err_mu / scale_2;

  // check for convergence
  
  const value_type opt_tol = algo.algo_cntr().opt_tol();
  const value_type feas_tol = algo.algo_cntr().feas_tol();
  const value_type comp_tol = algo.algo_cntr().comp_tol();

  if (opt_err < opt_tol && feas_err < feas_tol && comp_err < comp_tol)
    {
    found_solution = true;
    }

  if( int(olevel) >= int(PRINT_ALGORITHM_STEPS) || (int(olevel) >= int(PRINT_BASIC_ALGORITHM_INFO) && found_solution) )
    {
    out 
      << "\nopt_kkt_err_k   = " << opt_err << ( opt_err < opt_tol ? " < " : " > " )
      << "opt_tol = " << opt_tol
      << "\nfeas_kkt_err_k   = " << feas_err << ( feas_err < feas_tol ? " < " : " > " )
      << "feas_tol = " << feas_tol
      << "\ncomp_kkt_err_k   = " << comp_err << ( comp_err < comp_tol ? " < " : " > " )
      << "comp_tol = " << comp_tol
      << "\ncomp_err_mu      = " << comp_err_mu
      << "\nbarrier_parameter_k (mu_km1) = " << mu_km1
      << "comp_tol = " << comp_tol
      << std::endl;
    }
    
  return found_solution;
  }

void CheckConvergenceIP_Strategy::print_step( const Algorithm& _algo, std::ostream& out, const std::string& L ) const
  {
  out 
    << L << "CheckConvergenceIP_Strategy\n"
    << L << " IP_found_solution = CheckConvergedStd_Strategy::Converged(_algo, reportFinalSolution)\n";
  
  CheckConvergenceStd_Strategy::print_step(_algo, out, L+"   ");
  
  out 
    << L << "*** recalculate comp_err\n"
    << L << "comp_err_k = 0.0"
    << L << "for all i = 1 to n\n"
    << 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"
    << L << "next i\n"
    << L << "if IP_found_solution then\n"
    << L << "   IP_found_solution = false\n"
    << L << "   if (comp_err_k < comp_tol && mu_k < comp_tol) then\n"
    << L << "      IP_found_solution = true\n"
    << L << "   end\n"
    << L << "end\n"
    << L << "return IP_found_solution\n";
  }

} // end namespace MoochoPack