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

  SmartPtr<const Vector>
  AugRestoSystemSolver::Neg_Omega_c_plus_D_c(
    const SmartPtr<const Vector>& sigma_tilde_n_c_inv,
    const SmartPtr<const Vector>& sigma_tilde_p_c_inv,
    const Vector* D_c,
    const Vector& any_vec_in_c)
  {
    DBG_START_METH("AugRestoSystemSolver::Neg_Omega_c_plus_D_c",dbg_verbosity);
    SmartPtr<Vector> retVec;
    if (IsValid(sigma_tilde_n_c_inv) || IsValid(sigma_tilde_p_c_inv) || D_c) {
      if (!neg_omega_c_plus_D_c_cache_.
          GetCachedResult3Dep(retVec, GetRawPtr(sigma_tilde_n_c_inv), GetRawPtr(sigma_tilde_p_c_inv), D_c)) {
        DBG_PRINT((1,"Not found in cache\n"));
        retVec = any_vec_in_c.MakeNew();

        Number fact1, fact2;
        SmartPtr<const Vector> v1;
        SmartPtr<const Vector> v2;

        if (IsValid(sigma_tilde_n_c_inv)) {
          v1 = sigma_tilde_n_c_inv;
          fact1 = -1.;
        }
        else {
          v1 = &any_vec_in_c;
          fact1 = 0.;
        }
        if (IsValid(sigma_tilde_p_c_inv)) {
          v2 = sigma_tilde_p_c_inv;
          fact2 = -1.;
        }
        else {
          v2 = &any_vec_in_c;
          fact2 = 0.;
        }
        retVec->AddTwoVectors(fact1, *v1, fact2, *v2, 0.);

        if (D_c) {
          retVec->Axpy(1.0, *D_c);
        }

        neg_omega_c_plus_D_c_cache_.
        AddCachedResult3Dep(retVec, GetRawPtr(sigma_tilde_n_c_inv), GetRawPtr(sigma_tilde_p_c_inv), D_c);
      }
    }
    return ConstPtr(retVec);
  }

  char CGPenaltyLSAcceptor::UpdatePenaltyParameter()
  {
    DBG_START_METH("CGPenaltyLSAcceptor::UpdatePenaltyParameter",
                   dbg_verbosity);
    char info_alpha_primal_char = 'n';
    // We use the new infeasibility here...
    Number trial_inf = IpCq().trial_primal_infeasibility(NORM_2);
    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
                   "trial infeasibility = %8.2\n", trial_inf);
    if (curr_eta_<0.) {
      // We need to initialize the eta tolerance
      curr_eta_ = Max(eta_min_, Min(gamma_tilde_,
                                    gamma_hat_*IpCq().curr_nlp_error()));
    }
    // Check if the penalty parameter is to be increased
    Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
                   "Starting tests for penalty parameter update:\n");
    bool increase = (trial_inf >= penalty_update_infeasibility_tol_);
    if (!increase) {
      info_alpha_primal_char='i';
    }
    if (increase) {
      Number max_step = Max(CGPenData().delta_cgpen()->x()->Amax(),
                            CGPenData().delta_cgpen()->s()->Amax());
      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
                     "Max norm of step = %8.2\n", max_step);
      increase = (max_step <= curr_eta_);
      if (!increase) {
        info_alpha_primal_char='d';
      }
    }
    // Lifeng: Should we use the new complementarity here?  If so, I
    // have to restructure BacktrackingLineSearch
    Number mu = IpData().curr_mu();
    if (increase) {
      Number min_compl = mu;
      Number max_compl = mu;
      if (IpNLP().x_L()->Dim()>0) {
        SmartPtr<const Vector> compl_x_L = IpCq().curr_compl_x_L();
        min_compl = Min(min_compl, compl_x_L->Min());
        max_compl = Max(max_compl, compl_x_L->Max());
      }
      if (IpNLP().x_U()->Dim()>0) {
        SmartPtr<const Vector> compl_x_U = IpCq().curr_compl_x_U();
        min_compl = Min(min_compl, compl_x_U->Min());
        max_compl = Max(max_compl, compl_x_U->Max());
      }
      if (IpNLP().d_L()->Dim()>0) {
        SmartPtr<const Vector> compl_s_L = IpCq().curr_compl_s_L();
        min_compl = Min(min_compl, compl_s_L->Min());
        max_compl = Max(max_compl, compl_s_L->Max());
      }
      if (IpNLP().d_U()->Dim()>0) {
        SmartPtr<const Vector> compl_s_U = IpCq().curr_compl_s_U();
        min_compl = Min(min_compl, compl_s_U->Min());
        max_compl = Max(max_compl, compl_s_U->Max());
      }
      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
                     "Minimal compl = %8.2\n", min_compl);
      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
                     "Maximal compl = %8.2\n", max_compl);
      increase = (min_compl >= mu*penalty_update_compl_tol_ &&
                  max_compl <= mu/penalty_update_compl_tol_);
      if (!increase) {
        info_alpha_primal_char='c';
      }
    }
    // Lifeng: Here I'm using the information from the current step
    // and the current infeasibility
    if (increase) {
      SmartPtr<Vector> vec = IpData().curr()->y_c()->MakeNewCopy();
      vec->AddTwoVectors(1., *CGPenData().delta_cgpen()->y_c(),
                         -1./CGPenCq().curr_cg_pert_fact(), *IpCq().curr_c(),
                         1.);
      Number omega_test = vec->Amax();
      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
                     "omega_test for c = %8.2\n", omega_test);
      increase = (omega_test < curr_eta_);
      if (increase) {
        SmartPtr<Vector> vec = IpData().curr()->y_d()->MakeNewCopy();
        vec->AddTwoVectors(1., *IpData().delta()->y_d(),
                           -1./CGPenCq().curr_cg_pert_fact(), *IpCq().curr_d_minus_s(),
                           1.);
        omega_test = vec->Amax();
        Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
                       "omega_test for d = %8.2\n", omega_test);
        increase = (omega_test < curr_eta_);
      }
      if (!increase) {
        info_alpha_primal_char='m';
      }
    }
    if (increase) {
      // Ok, now we should increase the penalty parameter
      counter_first_type_penalty_updates_++;
      // Update the eta tolerance
      curr_eta_ = Max(eta_min_, curr_eta_/2.);
      Jnlst().Printf(J_MOREDETAILED, J_LINE_SEARCH,
                     "Updating eta to = %8.2\n", curr_eta_);
      Number penalty = CGPenData().curr_kkt_penalty();
//.........这里部分代码省略.........

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

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

Number InexactLSAcceptor::ComputeAlphaForY(
   Number                    alpha_primal,
   Number                    alpha_dual,
   SmartPtr<IteratesVector>& delta
   )
{
   DBG_START_METH("InexactLSAcceptor::ComputeAlphaForY",
      dbg_verbosity);

   // Here, we choose as stepsize for y either alpha_primal, if the
   // conditions from the ineqaxt paper is satisfied for it, or we
   // compute the step size closest to alpha_primal but great than
   // it, that does give the same progress as the full step would
   // give.

   Number alpha_y = alpha_primal;

   SmartPtr<Vector> gx = IpCq().curr_grad_barrier_obj_x()->MakeNewCopy();
   gx->AddTwoVectors(1., *IpCq().curr_jac_cT_times_curr_y_c(), 1., *IpCq().curr_jac_dT_times_curr_y_d(), 1.);
   SmartPtr<Vector> Jxy = gx->MakeNew();
   IpCq().curr_jac_c()->TransMultVector(1., *delta->y_c(), 0., *Jxy);
   IpCq().curr_jac_d()->TransMultVector(1., *delta->y_d(), 1., *Jxy);

   SmartPtr<const Vector> curr_scaling_slacks = InexCq().curr_scaling_slacks();
   SmartPtr<Vector> gs = curr_scaling_slacks->MakeNew();
   gs->AddTwoVectors(1., *IpCq().curr_grad_barrier_obj_s(), -1., *IpData().curr()->y_d(), 0.);
   gs->ElementWiseMultiply(*curr_scaling_slacks);

   SmartPtr<Vector> Sdy = delta->y_d()->MakeNewCopy();
   Sdy->ElementWiseMultiply(*curr_scaling_slacks);

   // using the magic formula in my notebook
   Number a = pow(Jxy->Nrm2(), 2) + pow(Sdy->Nrm2(), 2);
   Number b = 2 * (gx->Dot(*Jxy) - gs->Dot(*Sdy));
   Number c = pow(gx->Nrm2(), 2) + pow(gs->Nrm2(), 2);

   // First we check if the primal step size is good enough:
   Number val_ap = alpha_primal * alpha_primal * a + alpha_primal * b + c;
   Number val_1 = a + b + c;

   if( val_ap <= val_1 )
   {
      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "  Step size for y: using alpha_primal\n.");
   }
   else
   {
      Number alpha_2 = -b / a - 1.;
      Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "  Step size for y candidate: %8.2e - ", alpha_2);
      if( alpha_2 > alpha_primal && alpha_2 < 1. )
      {
         alpha_y = alpha_2;
         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "using that one\n.");
      }
      else
      {
         alpha_y = 1.;
         Jnlst().Printf(J_DETAILED, J_LINE_SEARCH, "using 1 instead\n");
      }
   }

   return alpha_y;
}

  Number IpoptAlgorithm::correct_bound_multiplier(
    const Vector& trial_z,
    const Vector& trial_slack,
    const Vector& trial_compl,
    SmartPtr<const Vector>& new_trial_z)
  {
    DBG_START_METH("IpoptAlgorithm::CorrectBoundMultiplier",
                   dbg_verbosity);

    if (kappa_sigma_<1. || trial_z.Dim()==0) {
      new_trial_z = &trial_z;
      return 0.;
    }

    // We choose as barrier parameter to be used either the current
    // algorithmic barrier parameter (if we are not in the free mode),
    // or the average complementarity (at the trial point)
    Number mu;
    if (IpData().FreeMuMode()) {
      mu = IpCq().trial_avrg_compl();
      mu = Min(mu, 1e3);
    }
    else {
      mu = IpData().curr_mu();
    }
    DBG_PRINT((1,"mu = %8.2e\n", mu));
    DBG_PRINT_VECTOR(2, "trial_z", trial_z);

    // First check quickly if anything need to be corrected, using the
    // trial complementarity directly.  Here, Amax is the same as Max
    // (and we use Amax because that can be used later)
    if (trial_compl.Amax() <= kappa_sigma_*mu &&
        trial_compl.Min() >= 1./kappa_sigma_*mu) {
      new_trial_z = &trial_z;
      return 0.;
    }

    SmartPtr<Vector> one_over_s = trial_z.MakeNew();
    one_over_s->Copy(trial_slack);
    one_over_s->ElementWiseReciprocal();

    SmartPtr<Vector> step_z = trial_z.MakeNew();
    step_z->AddTwoVectors(kappa_sigma_*mu, *one_over_s, -1., trial_z, 0.);

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

    Number max_correction_up = Max(0., -step_z->Min());
    if (max_correction_up>0.) {
      SmartPtr<Vector> tmp = trial_z.MakeNew();
      tmp->Set(0.);
      step_z->ElementWiseMin(*tmp);
      tmp->AddTwoVectors(1., trial_z, 1., *step_z, 0.);
      new_trial_z = GetRawPtr(tmp);
    }
    else {
      new_trial_z = &trial_z;
    }

    step_z->AddTwoVectors(1./kappa_sigma_*mu, *one_over_s, -1., *new_trial_z, 0.);

    Number max_correction_low = Max(0., step_z->Max());
    if (max_correction_low>0.) {
      SmartPtr<Vector> tmp = trial_z.MakeNew();
      tmp->Set(0.);
      step_z->ElementWiseMax(*tmp);
      tmp->AddTwoVectors(1., *new_trial_z, 1., *step_z, 0.);
      new_trial_z = GetRawPtr(tmp);
    }

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

    return Max(max_correction_up, max_correction_low);
  }