C++ AbstractLinAlgPack::num_bounded方法代码示例

QPSolverStats::ESolutionType
QPSolverRelaxedQPKWIK::imp_solve_qp(
  std::ostream* out, EOutputLevel olevel, ERunTests test_what
  ,const Vector& g, const MatrixSymOp& G
  ,value_type etaL
  ,const Vector* dL, const Vector* dU
  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
  ,const Vector* eL, const Vector* eU
  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
  ,value_type* obj_d
  ,value_type* eta, VectorMutable* d
  ,VectorMutable* nu
  ,VectorMutable* mu, VectorMutable* Ed
  ,VectorMutable* lambda, VectorMutable* Fd
  )
{
  using Teuchos::dyn_cast;
  using DenseLinAlgPack::nonconst_tri_ele;
  using LinAlgOpPack::dot;
  using LinAlgOpPack::V_StV;
  using LinAlgOpPack::assign;
  using LinAlgOpPack::V_StV;
  using LinAlgOpPack::V_MtV;
  using AbstractLinAlgPack::EtaVector;
  using AbstractLinAlgPack::transVtMtV;
  using AbstractLinAlgPack::num_bounded;
  using ConstrainedOptPack::MatrixExtractInvCholFactor;

  // /////////////////////////
  // Map to QPKWIK input

  // Validate that rHL is of the proper type.
  const MatrixExtractInvCholFactor &cG
    = dyn_cast<const MatrixExtractInvCholFactor>(G);

  // Determine the number of sparse bounds on variables and inequalities.
  // By default set for the dense case
  const value_type inf = this->infinite_bound();
  const size_type
    nd              = d->dim(),
    m_in            = E  ? b->dim()                  : 0,
    m_eq            = F  ? f->dim()                  : 0,
    nvarbounds      = dL ? num_bounded(*dL,*dU,inf)  : 0,
    ninequbounds    = E  ? num_bounded(*eL,*eU,inf)  : 0,
    nequalities     = F  ? f->dim()                  : 0;

  // Determine if this is a QP with a structure different from the
  // one just solved.
  
  const bool same_qp_struct = (  N_ == nd && M1_ == nvarbounds && M2_ == ninequbounds && M3_ == nequalities );

  /////////////////////////////////////////////////////////////////
  // Set the input parameters to be sent to QPKWIKNEW

  // N
  N_ = nd;

  // M1
  M1_ = nvarbounds;

  // M2
  M2_ = ninequbounds;

  // M3
  M3_ = nequalities;

  // GRAD
  GRAD_ = VectorDenseEncap(g)();

  // UINV_AUG
  //
  // UINV_AUG = [ sqrt(bigM)  0  ]
  //            [ 0           L' ]
  //
  UINV_AUG_.resize(N_+1,N_+1);
  cG.extract_inv_chol( &nonconst_tri_ele( UINV_AUG_(2,N_+1,2,N_+1), BLAS_Cpp::upper ) );
  UINV_AUG_(1,1) = 1.0 / ::sqrt( NUMPARAM_[2] );
  UINV_AUG_.col(1)(2,N_+1) = 0.0;
  UINV_AUG_.row(1)(2,N_+1) = 0.0;

  // LDUINV_AUG
  LDUINV_AUG_ = UINV_AUG_.rows();

  // IBND, BL , BU, A, LDA, YPY

  IBND_INV_.resize( nd + m_in);
  std::fill( IBND_INV_.begin(), IBND_INV_.end(), 0 ); // Initialize the zero
  IBND_.resize( my_max( 1, M1_ + M2_ ) );
  BL_.resize( my_max( 1, M1_ + M2_ ) );
  BU_.resize( my_max( 1, M1_ + M2_ + M3_ ) );
  LDA_ = my_max( 1, M2_ + M3_ );
  A_.resize( LDA_, (  M2_ + M3_ > 0 ? N_ : 1 ) );
  YPY_.resize( my_max( 1, M1_ + M2_ ) );
  if(M1_)
    YPY_(1,M1_) = 0.0; // Must be for this QP interface

  // Initialize variable bound constraints
  if( dL ) {
    VectorDenseEncap dL_de(*dL);
    VectorDenseEncap dU_de(*dU);
//.........这里部分代码省略.........

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
//.........这里部分代码省略.........