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

  /* The functions SchurSolve do IFT step, if S_==NULL, and DenseGenSchurDriver otherwise. */
  bool DenseGenSchurDriver::SchurSolve(SmartPtr<IteratesVector> lhs, // new left hand side will be stored here
				       SmartPtr<const IteratesVector> rhs, // rhs r_s
				       SmartPtr<Vector> delta_u,  // should be (u_p - u_0) WATCH OUT FOR THE SIGN! I like it this way, so that u_0+delta_u = u_p, but victor always used it the other way round, so be careful. At the end, delta_nu is saved in here.
				       SmartPtr<IteratesVector> sol) // the vector K^(-1)*r_s which usually should have been computed before.

  {
    DBG_START_METH("DenseGenSchurDriver::SchurSolve", dbg_verbosity);
    DBG_ASSERT(IsValid(S_));
    bool retval;

    // set up rhs of equation (3.48a)
    SmartPtr<Vector> delta_rhs = delta_u->MakeNew();
    data_B()->Multiply(*sol, *delta_rhs);
    delta_rhs->Print(Jnlst(),J_VECTOR,J_USER1,"delta_rhs");
    delta_rhs->Scal(-1.0);
    delta_rhs->Axpy(1.0, *delta_u);
    delta_rhs->Print(Jnlst(),J_VECTOR,J_USER1,"rhs 3.48a");

    // solve equation (3.48a) for delta_nu
    SmartPtr<DenseVector> delta_nu = dynamic_cast<DenseVector*>(GetRawPtr(delta_rhs))->MakeNewDenseVector();
    delta_nu->Copy(*delta_rhs);
    S_->LUSolveVector(*delta_nu); // why is LUSolveVector not bool??
    delta_nu->Print(Jnlst(),J_VECTOR,J_USER1,"delta_nu");

    // solve equation (3.48b) for lhs (=delta_s)
    SmartPtr<IteratesVector> new_rhs = lhs->MakeNewIteratesVector();
    data_A()->TransMultiply(*delta_nu, *new_rhs);
    new_rhs->Axpy(-1.0, *rhs);
    new_rhs->Scal(-1.0);
    new_rhs->Print(Jnlst(),J_VECTOR,J_USER1,"new_rhs");
    retval = backsolver_->Solve(lhs, ConstPtr(new_rhs));

    return retval;
  }

  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;
  }

//.........这里部分代码省略.........

   // 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));
         tmp = normal_s->MakeNew();
         tmp->AddTwoVectors(1., *IpCq().curr_jac_d_times_vec(*normal_x), -1., *normal_s, 0.);
         rhs->Set_y_d(*tmp);

      }

      InexData().set_normal_x(normal_x);
      InexData().set_normal_s(normal_s);

      delta = rhs->MakeNewIteratesVector();
      retval = inexact_pd_solver_->Solve(*rhs, *delta);

      // Determine if acceptable step has been computed
      if( !compute_normal && (!retval || InexData().next_compute_normal()) )
      {
         // If normal step has not been computed and step is not satisfactory, try computing normal step
         InexData().set_compute_normal(true);
         compute_normal = true;
      }
      else
      {
         // If normal step has been computed, stop anyway
         done = true;
      }
   }

   if( retval )
   {
      // Store the search directions in the IpData object
      IpData().set_delta(delta);
      if( InexData().compute_normal() )
      {
         IpData().Append_info_string("NT ");
      }
      else
      {
         IpData().Append_info_string("PD ");
      }
   }

   return retval;
}

bool ProbingMuOracle::CalculateMu(
   Number  mu_min,
   Number  mu_max,
   Number& new_mu
   )
{
   DBG_START_METH("ProbingMuOracle::CalculateMu",
      dbg_verbosity);

   /////////////////////////////////////
   // Compute the affine scaling step //
   /////////////////////////////////////

   Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, "Solving the Primal Dual System for the affine step\n");
   // First get the right hand side
   SmartPtr<IteratesVector> rhs = IpData().curr()->MakeNewContainer();

   rhs->Set_x(*IpCq().curr_grad_lag_x());
   rhs->Set_s(*IpCq().curr_grad_lag_s());
   rhs->Set_y_c(*IpCq().curr_c());
   rhs->Set_y_d(*IpCq().curr_d_minus_s());
   rhs->Set_z_L(*IpCq().curr_compl_x_L());
   rhs->Set_z_U(*IpCq().curr_compl_x_U());
   rhs->Set_v_L(*IpCq().curr_compl_s_L());
   rhs->Set_v_U(*IpCq().curr_compl_s_U());

   // Get space for the affine scaling step
   SmartPtr<IteratesVector> step = rhs->MakeNewIteratesVector(true);

   // Now solve the primal-dual system to get the affine step.  We
   // allow a somewhat inexact solution here
   bool allow_inexact = true;
   bool retval = pd_solver_->Solve(-1.0, 0.0, *rhs, *step, allow_inexact);
   if( !retval )
   {
      Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, "The linear system could not be solved for the affine step!\n");
      return false;
   }

   DBG_PRINT_VECTOR(2, "step", *step);

   /////////////////////////////////////////////////////////////
   // Use Mehrotra's formula to compute the barrier parameter //
   /////////////////////////////////////////////////////////////

   // First compute the fraction-to-the-boundary step sizes
   Number alpha_primal_aff = IpCq().primal_frac_to_the_bound(1.0, *step->x(), *step->s());

   Number alpha_dual_aff = IpCq().dual_frac_to_the_bound(1.0, *step->z_L(), *step->z_U(), *step->v_L(), *step->v_U());

   Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, "  The affine maximal step sizes are\n"
      "   alpha_primal_aff = %23.16e\n"
      "   alpha_dual_aff = %23.16e\n", alpha_primal_aff, alpha_dual_aff);

   // now compute the average complementarity at the affine step
   // ToDo shoot for mu_min instead of 0?
   Number mu_aff = CalculateAffineMu(alpha_primal_aff, alpha_dual_aff, *step);
   Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, "  The average complementariy at the affine step is %23.16e\n", mu_aff);

   // get the current average complementarity
   Number mu_curr = IpCq().curr_avrg_compl();
   Jnlst().Printf(J_DETAILED, J_BARRIER_UPDATE, "  The average complementariy at the current point is %23.16e\n",
      mu_curr);
   DBG_ASSERT(mu_curr > 0.);

   // Apply Mehrotra's rule
   Number sigma = pow((mu_aff / mu_curr), 3);
   // Make sure, sigma is not too large
   sigma = Min(sigma, sigma_max_);

   Number mu = sigma * mu_curr;

   // Store the affine search direction (in case it is needed in the
   // line search for a corrector step)
   IpData().set_delta_aff(step);
   IpData().SetHaveAffineDeltas(true);

   char ssigma[40];
   sprintf(ssigma, " sigma=%8.2e", sigma);
   IpData().Append_info_string(ssigma);
   //sprintf(ssigma, " xi=%8.2e ", IpCq().curr_centrality_measure());
   //IpData().Append_info_string(ssigma);

   new_mu = Max(Min(mu, mu_max), mu_min);
   return true;
}