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

  void IpoptAlgorithm::calc_number_of_bounds(
    const Vector& x,
    const Vector& x_L,
    const Vector& x_U,
    const Matrix& Px_L,
    const Matrix& Px_U,
    Index& n_tot,
    Index& n_only_lower,
    Index& n_both,
    Index& n_only_upper)
  {
    DBG_START_METH("IpoptAlgorithm::calc_number_of_bounds",
                   dbg_verbosity);

    n_tot = x.Dim();

    SmartPtr<Vector> tmpx = x.MakeNew();
    SmartPtr<Vector> tmpxL = x_L.MakeNew();
    SmartPtr<Vector> tmpxU = x_U.MakeNew();

    tmpxL->Set(-1.);
    tmpxU->Set(2.);
    Px_L.MultVector(1.0, *tmpxL, 0.0, *tmpx);
    Px_U.MultVector(1.0, *tmpxU, 1.0, *tmpx);
    // Now, x has elements
    //  -1 : if component has only lower bound
    //   0 : if component has no bound
    //   1 : if component has both lower and upper bound
    //   2 : if component has only upper bound
    DBG_PRINT_VECTOR(2, "x-indicator", *tmpx);

    SmartPtr<Vector> tmpx0 = x.MakeNew();
    tmpx0->Set(0.);

    SmartPtr<Vector> tmpx2 = x.MakeNew();
    tmpx2->Set(-1.0);
    tmpx2->Axpy(1.0, *tmpx);
    tmpx2->ElementWiseMax(*tmpx0); // tmpx2 is now 1 in those
    // components with only upper bounds
    n_only_upper = (Index)tmpx2->Asum();

    tmpx->Axpy(-2., *tmpx2);       // now make all those entries for
    // only upper bounds zero in tmpx

    tmpx2->Copy(*tmpx);
    tmpx2->ElementWiseMax(*tmpx0); // tmpx2 is now 1 in those
    // components with both bounds
    n_both = (Index)tmpx2->Asum();

    tmpx->Axpy(-1., *tmpx2);
    tmpx->ElementWiseMin(*tmpx);   // tmpx is now -1 in those with only
    // lower bounds
    n_only_lower = (Index)tmpx->Asum();

  }

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

  bool
  RestoRestorationPhase::PerformRestoration()
  {
    DBG_START_METH("RestoRestorationPhase::PerformRestoration",
                   dbg_verbosity);
    Jnlst().Printf(J_DETAILED, J_MAIN,
                   "Performing second level restoration phase for current constriant violation %8.2e\n", IpCq().curr_constraint_violation());

    DBG_ASSERT(IpCq().curr_constraint_violation()>0.);

    // Get a grip on the restoration phase NLP and obtain the pointers
    // to the original NLP data
    SmartPtr<RestoIpoptNLP> resto_ip_nlp =
      static_cast<RestoIpoptNLP*> (&IpNLP());
    DBG_ASSERT(dynamic_cast<RestoIpoptNLP*> (&IpNLP()));
    SmartPtr<IpoptNLP> orig_ip_nlp =
      static_cast<IpoptNLP*> (&resto_ip_nlp->OrigIpNLP());
    DBG_ASSERT(dynamic_cast<IpoptNLP*> (&resto_ip_nlp->OrigIpNLP()));

    // Get the current point and create a new vector for the result
    SmartPtr<const CompoundVector> Ccurr_x =
      static_cast<const CompoundVector*> (GetRawPtr(IpData().curr()->x()));
    SmartPtr<Vector> new_x = IpData().curr()->x()->MakeNew();
    SmartPtr<CompoundVector> Cnew_x =
      static_cast<CompoundVector*> (GetRawPtr(new_x));

    // The x values remain unchanged
    SmartPtr<Vector> x = Cnew_x->GetCompNonConst(0);
    x->Copy(*Ccurr_x->GetComp(0));

    // ToDo in free mu mode - what to do here?
    Number mu = IpData().curr_mu();

    // Compute the initial values for the n and p variables for the
    // equality constraints
    Number rho = resto_ip_nlp->Rho();
    SmartPtr<Vector> nc = Cnew_x->GetCompNonConst(1);
    SmartPtr<Vector> pc = Cnew_x->GetCompNonConst(2);
    SmartPtr<const Vector> cvec = orig_ip_nlp->c(*Ccurr_x->GetComp(0));
    SmartPtr<Vector> a = nc->MakeNew();
    SmartPtr<Vector> b = nc->MakeNew();
    a->Set(mu/(2.*rho));
    a->Axpy(-0.5, *cvec);
    b->Copy(*cvec);
    b->Scal(mu/(2.*rho));
    solve_quadratic(*a, *b, *nc);
    pc->Copy(*cvec);
    pc->Axpy(1., *nc);
    DBG_PRINT_VECTOR(2, "nc", *nc);
    DBG_PRINT_VECTOR(2, "pc", *pc);

    // initial values for the n and p variables for the inequality
    // constraints
    SmartPtr<Vector> nd = Cnew_x->GetCompNonConst(3);
    SmartPtr<Vector> pd = Cnew_x->GetCompNonConst(4);
    SmartPtr<Vector> dvec = pd->MakeNew();
    dvec->Copy(*orig_ip_nlp->d(*Ccurr_x->GetComp(0)));
    dvec->Axpy(-1., *IpData().curr()->s());
    a = nd->MakeNew();
    b = nd->MakeNew();
    a->Set(mu/(2.*rho));
    a->Axpy(-0.5, *dvec);
    b->Copy(*dvec);
    b->Scal(mu/(2.*rho));
    solve_quadratic(*a, *b, *nd);
    pd->Copy(*dvec);
    pd->Axpy(1., *nd);
    DBG_PRINT_VECTOR(2, "nd", *nd);
    DBG_PRINT_VECTOR(2, "pd", *pd);

    // Now set the trial point to the solution of the restoration phase
    // s and all multipliers remain unchanged
    SmartPtr<IteratesVector> new_trial = IpData().curr()->MakeNewContainer();
    new_trial->Set_x(*new_x);
    IpData().set_trial(new_trial);

    IpData().Append_info_string("R");

    return true;
  }

  bool RestoIterateInitializer::SetInitialIterates()
  {
    DBG_START_METH("RestoIterateInitializer::SetInitialIterates",
                   dbg_verbosity);

    // Get a grip on the restoration phase NLP and obtain the pointers
    // to the original NLP data
    SmartPtr<RestoIpoptNLP> resto_ip_nlp =
      static_cast<RestoIpoptNLP*> (&IpNLP());
    SmartPtr<IpoptNLP> orig_ip_nlp =
      static_cast<IpoptNLP*> (&resto_ip_nlp->OrigIpNLP());
    SmartPtr<IpoptData> orig_ip_data =
      static_cast<IpoptData*> (&resto_ip_nlp->OrigIpData());
    SmartPtr<IpoptCalculatedQuantities> orig_ip_cq =
      static_cast<IpoptCalculatedQuantities*> (&resto_ip_nlp->OrigIpCq());

    // Set the value of the barrier parameter
    Number resto_mu;
    resto_mu = Max(orig_ip_data->curr_mu(),
                   orig_ip_cq->curr_c()->Amax(),
                   orig_ip_cq->curr_d_minus_s()->Amax());
    IpData().Set_mu(resto_mu);
    Jnlst().Printf(J_DETAILED, J_INITIALIZATION,
                   "Initial barrier parameter resto_mu = %e\n", resto_mu);

    /////////////////////////////////////////////////////////////////////
    //                   Initialize primal varialbes                   //
    /////////////////////////////////////////////////////////////////////

    // initialize the data structures in the restoration phase NLP
    IpData().InitializeDataStructures(IpNLP(), false, false, false,
                                      false, false);

    SmartPtr<Vector> new_x = IpData().curr()->x()->MakeNew();
    SmartPtr<CompoundVector> Cnew_x =
      static_cast<CompoundVector*> (GetRawPtr(new_x));

    // Set the trial x variables from the original NLP
    Cnew_x->GetCompNonConst(0)->Copy(*orig_ip_data->curr()->x());

    // Compute the initial values for the n and p variables for the
    // equality constraints
    Number rho = resto_ip_nlp->Rho();
    DBG_PRINT((1,"rho = %e\n", rho));
    SmartPtr<Vector> nc = Cnew_x->GetCompNonConst(1);
    SmartPtr<Vector> pc = Cnew_x->GetCompNonConst(2);
    SmartPtr<const Vector> cvec = orig_ip_cq->curr_c();
    DBG_PRINT_VECTOR(2, "cvec", *cvec);
    SmartPtr<Vector> a = nc->MakeNew();
    SmartPtr<Vector> b = nc->MakeNew();
    a->Set(resto_mu/(2.*rho));
    a->Axpy(-0.5, *cvec);
    b->Copy(*cvec);
    b->Scal(resto_mu/(2.*rho));
    DBG_PRINT_VECTOR(2, "a", *a);
    DBG_PRINT_VECTOR(2, "b", *b);
    solve_quadratic(*a, *b, *nc);
    pc->Copy(*cvec);
    pc->Axpy(1., *nc);
    DBG_PRINT_VECTOR(2, "nc", *nc);
    DBG_PRINT_VECTOR(2, "pc", *pc);

    // initial values for the n and p variables for the inequality
    // constraints
    SmartPtr<Vector> nd = Cnew_x->GetCompNonConst(3);
    SmartPtr<Vector> pd = Cnew_x->GetCompNonConst(4);
    cvec = orig_ip_cq->curr_d_minus_s();
    a = nd->MakeNew();
    b = nd->MakeNew();
    a->Set(resto_mu/(2.*rho));
    a->Axpy(-0.5, *cvec);
    b->Copy(*cvec);
    b->Scal(resto_mu/(2.*rho));
    solve_quadratic(*a, *b, *nd);
    pd->Copy(*cvec);
    pd->Axpy(1., *nd);
    DBG_PRINT_VECTOR(2, "nd", *nd);
    DBG_PRINT_VECTOR(2, "pd", *pd);

    // Leave the slacks unchanged
    SmartPtr<const Vector> new_s = orig_ip_data->curr()->s();

    // Now set the primal trial variables
    DBG_PRINT_VECTOR(2,"new_s",*new_s);
    DBG_PRINT_VECTOR(2,"new_x",*new_x);
    SmartPtr<IteratesVector> trial = IpData().curr()->MakeNewContainer();
    trial->Set_primal(*new_x, *new_s);
    IpData().set_trial(trial);

    DBG_PRINT_VECTOR(2, "resto_c", *IpCq().trial_c());
    DBG_PRINT_VECTOR(2, "resto_d_minus_s", *IpCq().trial_d_minus_s());

    /////////////////////////////////////////////////////////////////////
    //                   Initialize bound multipliers                  //
    /////////////////////////////////////////////////////////////////////

    SmartPtr<Vector> new_z_L = IpData().curr()->z_L()->MakeNew();
    SmartPtr<CompoundVector> Cnew_z_L =
      static_cast<CompoundVector*> (GetRawPtr(new_z_L));
    DBG_ASSERT(IsValid(Cnew_z_L));
//.........这里部分代码省略.........

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

    bool improve_solution = false;
    if (IpData().HaveDeltas()) {
      improve_solution = true;
    }

    bool retval;
    if (improve_solution && fast_step_computation_) {
      retval = true;
    }
    else {
      SmartPtr<IteratesVector> rhs = IpData().curr()->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(*IpCq().curr_c());
      rhs->Set_y_d(*IpCq().curr_d_minus_s());
      Index nbounds = IpNLP().x_L()->Dim()+ IpNLP().x_U()->Dim() +
                      IpNLP().d_L()->Dim()+ IpNLP().d_U()->Dim();
      if (nbounds>0 && mehrotra_algorithm_) {
        // set up the right hand side a la Mehrotra
        DBG_ASSERT(IpData().HaveAffineDeltas());
        DBG_ASSERT(!IpData().HaveDeltas());
        const SmartPtr<const IteratesVector> delta_aff = IpData().delta_aff();

        SmartPtr<Vector> tmpvec = delta_aff->z_L()->MakeNew();
        IpNLP().Px_L()->TransMultVector(1., *delta_aff->x(), 0., *tmpvec);
        tmpvec->ElementWiseMultiply(*delta_aff->z_L());
        tmpvec->Axpy(1., *IpCq().curr_relaxed_compl_x_L());
        rhs->Set_z_L(*tmpvec);

        tmpvec = delta_aff->z_U()->MakeNew();
        IpNLP().Px_U()->TransMultVector(-1., *delta_aff->x(), 0., *tmpvec);
        tmpvec->ElementWiseMultiply(*delta_aff->z_U());
        tmpvec->Axpy(1., *IpCq().curr_relaxed_compl_x_U());
        rhs->Set_z_U(*tmpvec);

        tmpvec = delta_aff->v_L()->MakeNew();
        IpNLP().Pd_L()->TransMultVector(1., *delta_aff->s(), 0., *tmpvec);
        tmpvec->ElementWiseMultiply(*delta_aff->v_L());
        tmpvec->Axpy(1., *IpCq().curr_relaxed_compl_s_L());
        rhs->Set_v_L(*tmpvec);

        tmpvec = delta_aff->v_U()->MakeNew();
        IpNLP().Pd_U()->TransMultVector(-1., *delta_aff->s(), 0., *tmpvec);
        tmpvec->ElementWiseMultiply(*delta_aff->v_U());
        tmpvec->Axpy(1., *IpCq().curr_relaxed_compl_s_U());
        rhs->Set_v_U(*tmpvec);
      }
      else {
        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());
      }

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

      // Get space for the search direction
      SmartPtr<IteratesVector> delta =
        IpData().curr()->MakeNewIteratesVector(true);

      if (improve_solution) {
        // We can probably avoid copying and scaling...
        delta->AddOneVector(-1., *IpData().delta(), 0.);
      }

      bool& allow_inexact = fast_step_computation_;
      retval = pd_solver_->Solve(-1.0, 0.0, *rhs, *delta, allow_inexact,
                                 improve_solution);
      if (retval) {
        // Store the search directions in the IpData object
        IpData().set_delta(delta);
      }
    }
    return retval;
  }