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

  SmartPtr<const Vector> AugRestoSystemSolver::Rhs_cR(const Vector& rhs_c,
      const SmartPtr<const Vector>& sigma_tilde_n_c_inv, const Vector& rhs_n_c,
      const SmartPtr<const Vector>& sigma_tilde_p_c_inv, const Vector& rhs_p_c)
  {
    DBG_START_METH("AugRestoSystemSolver::Rhs_cR",dbg_verbosity);
    SmartPtr<Vector> retVec;
    std::vector<const TaggedObject*> deps(5);
    std::vector<Number> scalar_deps;
    deps[0] = &rhs_c;
    deps[1] = GetRawPtr(sigma_tilde_n_c_inv);
    deps[2] = &rhs_n_c;
    deps[3] = GetRawPtr(sigma_tilde_p_c_inv);
    deps[4] = &rhs_p_c;
    if (!rhs_cR_cache_.GetCachedResult(retVec, deps, scalar_deps)) {
      DBG_PRINT((1,"Not found in cache\n"));
      retVec = rhs_c.MakeNew();
      retVec->Copy(rhs_c);

      SmartPtr<Vector> tmp = retVec->MakeNew();
      if (IsValid(sigma_tilde_n_c_inv)) {
        tmp->Copy(*sigma_tilde_n_c_inv);
        tmp->ElementWiseMultiply(rhs_n_c);
        retVec->Axpy(-1.0, *tmp);
      }

      if (IsValid(sigma_tilde_p_c_inv)) {
        tmp->Copy(*sigma_tilde_p_c_inv);
        tmp->ElementWiseMultiply(rhs_p_c);
        retVec->Axpy(1.0, *tmp);
      }
      rhs_cR_cache_.AddCachedResult(retVec, deps, scalar_deps);
    }
    return ConstPtr(retVec);
  }

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

  void PDFullSpaceSolver::ComputeResiduals(
    const SymMatrix& W,
    const Matrix& J_c,
    const Matrix& J_d,
    const Matrix& Px_L,
    const Matrix& Px_U,
    const Matrix& Pd_L,
    const Matrix& Pd_U,
    const Vector& z_L,
    const Vector& z_U,
    const Vector& v_L,
    const Vector& v_U,
    const Vector& slack_x_L,
    const Vector& slack_x_U,
    const Vector& slack_s_L,
    const Vector& slack_s_U,
    const Vector& sigma_x,
    const Vector& sigma_s,
    Number alpha,
    Number beta,
    const IteratesVector& rhs,
    const IteratesVector& res,
    IteratesVector& resid)
  {
    DBG_START_METH("PDFullSpaceSolver::ComputeResiduals", dbg_verbosity);

    DBG_PRINT_VECTOR(2, "res", res);
    IpData().TimingStats().ComputeResiduals().Start();

    // Get the current sizes of the perturbation factors
    Number delta_x;
    Number delta_s;
    Number delta_c;
    Number delta_d;
    perturbHandler_->CurrentPerturbation(delta_x, delta_s, delta_c, delta_d);

    SmartPtr<Vector> tmp;

    // x
    W.MultVector(1., *res.x(), 0., *resid.x_NonConst());
    J_c.TransMultVector(1., *res.y_c(), 1., *resid.x_NonConst());
    J_d.TransMultVector(1., *res.y_d(), 1., *resid.x_NonConst());
    Px_L.MultVector(-1., *res.z_L(), 1., *resid.x_NonConst());
    Px_U.MultVector(1., *res.z_U(), 1., *resid.x_NonConst());
    resid.x_NonConst()->AddTwoVectors(delta_x, *res.x(), -1., *rhs.x(), 1.);

    // s
    Pd_U.MultVector(1., *res.v_U(), 0., *resid.s_NonConst());
    Pd_L.MultVector(-1., *res.v_L(), 1., *resid.s_NonConst());
    resid.s_NonConst()->AddTwoVectors(-1., *res.y_d(), -1., *rhs.s(), 1.);
    if (delta_s!=0.) {
      resid.s_NonConst()->Axpy(delta_s, *res.s());
    }

    // c
    J_c.MultVector(1., *res.x(), 0., *resid.y_c_NonConst());
    resid.y_c_NonConst()->AddTwoVectors(-delta_c, *res.y_c(), -1., *rhs.y_c(), 1.);

    // d
    J_d.MultVector(1., *res.x(), 0., *resid.y_d_NonConst());
    resid.y_d_NonConst()->AddTwoVectors(-1., *res.s(), -1., *rhs.y_d(), 1.);
    if (delta_d!=0.) {
      resid.y_d_NonConst()->Axpy(-delta_d, *res.y_d());
    }

    // zL
    resid.z_L_NonConst()->Copy(*res.z_L());
    resid.z_L_NonConst()->ElementWiseMultiply(slack_x_L);
    tmp = z_L.MakeNew();
    Px_L.TransMultVector(1., *res.x(), 0., *tmp);
    tmp->ElementWiseMultiply(z_L);
    resid.z_L_NonConst()->AddTwoVectors(1., *tmp, -1., *rhs.z_L(), 1.);

    // zU
    resid.z_U_NonConst()->Copy(*res.z_U());
    resid.z_U_NonConst()->ElementWiseMultiply(slack_x_U);
    tmp = z_U.MakeNew();
    Px_U.TransMultVector(1., *res.x(), 0., *tmp);
    tmp->ElementWiseMultiply(z_U);
    resid.z_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.z_U(), 1.);

    // vL
    resid.v_L_NonConst()->Copy(*res.v_L());
    resid.v_L_NonConst()->ElementWiseMultiply(slack_s_L);
    tmp = v_L.MakeNew();
    Pd_L.TransMultVector(1., *res.s(), 0., *tmp);
    tmp->ElementWiseMultiply(v_L);
    resid.v_L_NonConst()->AddTwoVectors(1., *tmp, -1., *rhs.v_L(), 1.);

    // vU
    resid.v_U_NonConst()->Copy(*res.v_U());
    resid.v_U_NonConst()->ElementWiseMultiply(slack_s_U);
    tmp = v_U.MakeNew();
    Pd_U.TransMultVector(1., *res.s(), 0., *tmp);
    tmp->ElementWiseMultiply(v_U);
    resid.v_U_NonConst()->AddTwoVectors(-1., *tmp, -1., *rhs.v_U(), 1.);

    DBG_PRINT_VECTOR(2, "resid", resid);

    if (Jnlst().ProduceOutput(J_MOREVECTOR, J_LINEAR_ALGEBRA)) {
//.........这里部分代码省略.........

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