C++ SmartPtr::MakeNewContainer方法代码示例

bool InexactSearchDirCalculator::ComputeSearchDirection()
{
   DBG_START_METH("InexactSearchDirCalculator::ComputeSearchDirection",
      dbg_verbosity);

   // First check if the iterates have converged to a locally
   // infeasible point
   Number curr_scaled_Ac_norm = InexCq().curr_scaled_Ac_norm();
   Jnlst().Printf(J_DETAILED, J_SOLVE_PD_SYSTEM, "curr_scaled_Ac_norm = %e\n", curr_scaled_Ac_norm);
   Number curr_inf = IpCq().curr_primal_infeasibility(NORM_2);
   // ToDo work on termination criteria
   if( curr_scaled_Ac_norm <= local_inf_Ac_tol_ && curr_inf > 1e-4 )
   {
      THROW_EXCEPTION(LOCALLY_INFEASIBLE, "The scaled norm of Ac is satisfying tolerance");
   }

   bool compute_normal = false;
   switch( decomposition_type_ )
   {
      case ALWAYS:
         compute_normal = true;
         break;
      case ADAPTIVE:
         compute_normal = InexData().next_compute_normal();
         break;
      case SWITCH_ONCE:
         compute_normal = InexData().next_compute_normal() || InexData().compute_normal();
         break;
   }

   SmartPtr<Vector> normal_x;
   SmartPtr<Vector> normal_s;
   bool retval;
   SmartPtr<IteratesVector> delta;
   SmartPtr<const IteratesVector> curr = IpData().curr();
   SmartPtr<IteratesVector> rhs;
   SmartPtr<Vector> tmp;

   // Now we set up the primal-dual system for computing the
   // tangential step and the search direction for the multipliers.
   // This is taken from IpPDSearchDirCal.cpp (rev 549).
   // We do not need entries for the variable bound multipliers

   // Upper part of right-hand-side vector is same for both systems
   rhs = curr->MakeNewContainer();
   tmp = curr->x()->MakeNew();
   tmp->AddOneVector(-1., *IpCq().curr_grad_lag_with_damping_x(), 0.);
   rhs->Set_x(*tmp);
   tmp = curr->s()->MakeNew();
   tmp->AddOneVector(-1., *IpCq().curr_grad_lag_with_damping_s(), 0.);
   rhs->Set_s(*tmp);
   tmp = curr->v_L()->MakeNew();
   tmp->AddOneVector(-1., *IpCq().curr_relaxed_compl_s_L(), 0.);
   rhs->Set_v_L(*tmp);
   tmp = curr->v_U()->MakeNew();
   tmp->AddOneVector(-1., *IpCq().curr_relaxed_compl_s_U(), 0.);
   rhs->Set_v_U(*tmp);

   // Loop through algorithms
   bool done = false;
   while( !done )
   {

      InexData().set_compute_normal(compute_normal);
      InexData().set_next_compute_normal(compute_normal);

      if( !compute_normal )
      {
         normal_x = NULL;
         normal_s = NULL;
      }
      else
      {
         retval = normal_step_calculator_->ComputeNormalStep(normal_x, normal_s);
         if( !retval )
         {
            return false;
         }
         // output
         if( Jnlst().ProduceOutput(J_VECTOR, J_SOLVE_PD_SYSTEM) )
         {
            Jnlst().Printf(J_VECTOR, J_SOLVE_PD_SYSTEM, "Normal step (without slack scaling):\n");
            normal_x->Print(Jnlst(), J_VECTOR, J_SOLVE_PD_SYSTEM, "normal_x");
            normal_s->Print(Jnlst(), J_VECTOR, J_SOLVE_PD_SYSTEM, "normal_s");
         }
      }

      // Lower part of right-hand-side vector is different for each system
      if( !compute_normal )
      {
         tmp = curr->y_c()->MakeNew();
         tmp->AddOneVector(-1., *IpCq().curr_c(), 0.);
         rhs->Set_y_c(*tmp);
         tmp = curr->y_d()->MakeNew();
         tmp->AddOneVector(-1., *IpCq().curr_d_minus_s(), 0.);
         rhs->Set_y_d(*tmp);
      }
      else
      {
         rhs->Set_y_c(*IpCq().curr_jac_c_times_vec(*normal_x));
//.........这里部分代码省略.........

  bool
  CGPenaltyLSAcceptor::TrySecondOrderCorrection(
    Number alpha_primal_test,
    Number& alpha_primal,
    SmartPtr<IteratesVector>& actual_delta)
  {
    DBG_START_METH("CGPenaltyLSAcceptor::TrySecondOrderCorrection",
                   dbg_verbosity);
    if (max_soc_==0) {
      return false;
    }
    bool accept = false;
    Index count_soc = 0;
    Number theta_soc_old = 0.;
    Number theta_soc_old2 = 0.;
    Number theta_trial = IpCq().trial_constraint_violation();
    Number theta_trial2 = IpCq().curr_primal_infeasibility(NORM_2);
    Number alpha_primal_soc = alpha_primal;
    // delta_y_c and delta_y_d are the steps used in the right hand
    // side for the SOC step
    SmartPtr<const Vector> delta_y_c = IpData().delta()->y_c();
    SmartPtr<const Vector> delta_y_d = IpData().delta()->y_d();
    SmartPtr<Vector> c_soc = IpCq().curr_c()->MakeNewCopy();
    SmartPtr<Vector> dms_soc = IpCq().curr_d_minus_s()->MakeNewCopy();
    while (count_soc<max_soc_ && !accept &&
           (count_soc==0 || (theta_trial<=kappa_soc_*theta_soc_old ||
                             theta_trial2<=kappa_soc_*theta_soc_old2)) ) {
      theta_soc_old = theta_trial;
      theta_soc_old2 = theta_trial2;
      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH,
                     "Trying second order correction number %d\n",
                     count_soc+1);
      // Compute SOC constraint violation
      /*
      Number c_over_r = 0.;
      if (IpData().BiggerJacPert()){
      c_over_r = IpCq().curr_cg_pert_fact();
      }*/
      c_soc->AddTwoVectors(1.0, *IpCq().trial_c(),
                           -CGPenData().CurrPenaltyPert(), *delta_y_c,
                           alpha_primal_soc);
      dms_soc->AddTwoVectors(1.0, *IpCq().trial_d_minus_s(),
                             -CGPenData().CurrPenaltyPert(), *delta_y_d,
                             alpha_primal_soc);
      // Compute the SOC search direction
      SmartPtr<IteratesVector> delta_soc =
        actual_delta->MakeNewIteratesVector(true);
      SmartPtr<IteratesVector> rhs = actual_delta->MakeNewContainer();
      rhs->Set_x(*IpCq().curr_grad_lag_with_damping_x());
      rhs->Set_s(*IpCq().curr_grad_lag_with_damping_s());
      rhs->Set_y_c(*c_soc);
      rhs->Set_y_d(*dms_soc);
      rhs->Set_z_L(*IpCq().curr_relaxed_compl_x_L());
      rhs->Set_z_U(*IpCq().curr_relaxed_compl_x_U());
      rhs->Set_v_L(*IpCq().curr_relaxed_compl_s_L());
      rhs->Set_v_U(*IpCq().curr_relaxed_compl_s_U());
      pd_solver_->Solve(-1.0, 0.0, *rhs, *delta_soc, true);
      // Update the delta_y_c and delta_y_d vectors in case we do
      // additional SOC steps
      delta_y_c = ConstPtr(delta_soc->y_c());
      delta_y_d = ConstPtr(delta_soc->y_d());
      // Compute step size
      alpha_primal_soc =
        IpCq().primal_frac_to_the_bound(IpData().curr_tau(),
                                        *delta_soc->x(),
                                        *delta_soc->s());
      // Check if trial point is acceptable
      try {
        // Compute the primal trial point
        IpData().SetTrialPrimalVariablesFromStep(alpha_primal_soc, *delta_soc->x(), *delta_soc->s());

        // in acceptance tests, use original step size!
        accept = CheckAcceptabilityOfTrialPoint(alpha_primal_test);
      }
      catch (IpoptNLP::Eval_Error& e) {
        e.ReportException(Jnlst(), J_DETAILED);
        Jnlst().Printf(J_WARNING, J_MAIN, "Warning: SOC step rejected due to evaluation error\n");
        IpData().Append_info_string("e");
        accept = false;
        // There is no point in continuing SOC procedure
        break;
      }
      if (accept) {
        Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "Second order correction step accepted with %d corrections.\n", count_soc+1);
        // Accept all SOC quantities
        alpha_primal = alpha_primal_soc;
        actual_delta = delta_soc;
      }
      else {
        count_soc++;
        theta_trial = IpCq().trial_constraint_violation();
        theta_trial2 = IpCq().trial_primal_infeasibility(NORM_2);
      }
    }
    if (accept) {
      ls_counter_ = 0;
    }
    return accept;
  }

//.........这里部分代码省略.........
    max_correction = correct_bound_multiplier(
                       *IpData().trial()->z_L(),
                       *IpCq().trial_slack_x_L(),
                       *IpCq().trial_compl_x_L(),
                       new_z_L);
    if (max_correction>0.) {
      Jnlst().Printf(J_DETAILED, J_MAIN,
                     "Some value in z_L becomes too large - maximal correction = %8.2e\n",
                     max_correction);
      corrected = true;
    }
    SmartPtr<const Vector> new_z_U;
    max_correction = correct_bound_multiplier(
                       *IpData().trial()->z_U(),
                       *IpCq().trial_slack_x_U(),
                       *IpCq().trial_compl_x_U(),
                       new_z_U);
    if (max_correction>0.) {
      Jnlst().Printf(J_DETAILED, J_MAIN,
                     "Some value in z_U becomes too large - maximal correction = %8.2e\n",
                     max_correction);
      corrected = true;
    }
    SmartPtr<const Vector> new_v_L;
    max_correction = correct_bound_multiplier(
                       *IpData().trial()->v_L(),
                       *IpCq().trial_slack_s_L(),
                       *IpCq().trial_compl_s_L(),
                       new_v_L);
    if (max_correction>0.) {
      Jnlst().Printf(J_DETAILED, J_MAIN,
                     "Some value in v_L becomes too large - maximal correction = %8.2e\n",
                     max_correction);
      corrected = true;
    }
    SmartPtr<const Vector> new_v_U;
    max_correction = correct_bound_multiplier(
                       *IpData().trial()->v_U(),
                       *IpCq().trial_slack_s_U(),
                       *IpCq().trial_compl_s_U(),
                       new_v_U);
    if (max_correction>0.) {
      Jnlst().Printf(J_DETAILED, J_MAIN,
                     "Some value in v_U becomes too large - maximal correction = %8.2e\n",
                     max_correction);
      corrected = true;
    }
    SmartPtr<IteratesVector> trial = IpData().trial()->MakeNewContainer();
    trial->Set_bound_mult(*new_z_L, *new_z_U, *new_v_L, *new_v_U);
    IpData().set_trial(trial);

    if (corrected) {
      IpData().Append_info_string("z");
    }

    // Accept the step
    IpData().AcceptTrialPoint();

    // If we want to recalculate the multipliers (e.g., as least
    // square estimates), call the calculator for that
    if (recalc_y_) {
      // There is no point in doing this if there are no constraints
      if (IpData().curr()->y_c()->Dim()+IpData().curr()->y_d()->Dim()==0) {
        recalc_y_ = false;
      }
    }
    if (recalc_y_ && IpCq().curr_constraint_violation()<recalc_y_feas_tol_) {
      if (Jnlst().ProduceOutput(J_MOREDETAILED, J_MAIN)) {
        Jnlst().Printf(J_MOREDETAILED, J_MAIN,
                       "dual infeasisibility before least square multiplier update = %e\n",
                       IpCq().curr_dual_infeasibility(NORM_MAX));
      }
      IpData().Append_info_string("y ");
      DBG_ASSERT(IsValid(eq_multiplier_calculator_));
      if (IpData().curr()->y_c()->Dim()+IpData().curr()->y_d()->Dim()>0) {
        SmartPtr<Vector> y_c = IpData().curr()->y_c()->MakeNew();
        SmartPtr<Vector> y_d = IpData().curr()->y_d()->MakeNew();
        bool retval =
          eq_multiplier_calculator_->CalculateMultipliers(*y_c, *y_d);
        if (retval) {
          SmartPtr<const IteratesVector> curr = IpData().curr();
          SmartPtr<IteratesVector> iterates = curr->MakeNewContainer();
          iterates->Set_x(*curr->x());
          iterates->Set_s(*curr->s());
          iterates->Set_z_L(*curr->z_L());
          iterates->Set_z_U(*curr->z_U());
          iterates->Set_v_L(*curr->v_L());
          iterates->Set_v_U(*curr->v_U());
          iterates->Set_y_c(*y_c);
          iterates->Set_y_d(*y_d);
          IpData().set_trial(iterates);
          IpData().AcceptTrialPoint();
        }
        else {
          Jnlst().Printf(J_DETAILED, J_MAIN,
                         "Recalculation of y multipliers skipped because eq_mult_calc returned false.\n");
        }
      }
    }
  }