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

  CompoundVector::CompoundVector(const CompoundVectorSpace* owner_space, bool create_new)
      :
      Vector(owner_space),
      comps_(owner_space->NCompSpaces()),
      const_comps_(owner_space->NCompSpaces()),
      owner_space_(owner_space),
      vectors_valid_(false)
  {
    Index dim_check = 0;
    for (Index i=0; i<NComps(); i++) {
      SmartPtr<const VectorSpace> space = owner_space_->GetCompSpace(i);
      DBG_ASSERT(IsValid(space));
      dim_check += space->Dim();

      if (create_new) {
        comps_[i] = space->MakeNew();
      }
    }

    DBG_ASSERT(dim_check == Dim());

    if (create_new) {
      vectors_valid_ = VectorsValid();
    }
  }

  void TripletHelper::PutValuesInVector(Index dim, const Number* values, Vector& vector)
  {
    DBG_ASSERT(dim == vector.Dim());
    DenseVector* dv = dynamic_cast<DenseVector*>(&vector);
    if (dv) {
      Number* dv_vals = dv->Values();
      IpBlasDcopy(dim, values, 1, dv_vals, 1);
      return;
    }

    CompoundVector* cv = dynamic_cast<CompoundVector*>(&vector);
    if (cv) {
      Index ncomps = cv->NComps();
      Index total_dim = 0;
      for (Index i=0; i<ncomps; i++) {
        SmartPtr<Vector> comp = cv->GetCompNonConst(i);
        Index comp_dim = comp->Dim();
        PutValuesInVector(comp_dim, values, *comp);
        values += comp_dim;
        total_dim += comp_dim;
      }
      DBG_ASSERT(total_dim == dim);
      return;
    }

    THROW_EXCEPTION(UNKNOWN_VECTOR_TYPE,"Unknown vector type passed to TripletHelper::PutValuesInVector");
  }

  void TripletHelper::FillValuesFromVector(Index dim, const Vector& vector, Number* values)
  {
    DBG_ASSERT(dim == vector.Dim());
    const DenseVector* dv = dynamic_cast<const DenseVector*>(&vector);
    if (dv) {
      if (dv->IsHomogeneous()) {
        Number scalar = dv->Scalar();
        IpBlasDcopy(dim, &scalar, 0, values, 1);
      }
      else {
        const Number* dv_vals = dv->Values();
        IpBlasDcopy(dim, dv_vals, 1, values, 1);
      }
      return;
    }

    const CompoundVector* cv = dynamic_cast<const CompoundVector*>(&vector);
    if (cv) {
      Index ncomps = cv->NComps();
      Index total_dim = 0;
      for (Index i=0; i<ncomps; i++) {
        SmartPtr<const Vector> comp = cv->GetComp(i);
        Index comp_dim = comp->Dim();
        FillValuesFromVector(comp_dim, *comp, values);
        values += comp_dim;
        total_dim += comp_dim;
      }
      DBG_ASSERT(total_dim == dim);
      return;
    }

    THROW_EXCEPTION(UNKNOWN_VECTOR_TYPE,"Unknown vector type passed to TripletHelper::FillValues");
  }

 /** Constructor, given the number of row and columns blocks, as
  *  well as the totel number of rows and columns.
  */
 SymScaledMatrixSpace(const SmartPtr<const Vector>& row_col_scaling,
                      bool row_col_scaling_reciprocal,
                      const SmartPtr<const SymMatrixSpace>& unscaled_matrix_space)
     :
     SymMatrixSpace(unscaled_matrix_space->Dim()),
     unscaled_matrix_space_(unscaled_matrix_space)
 {
   scaling_ = row_col_scaling->MakeNewCopy();
   if (row_col_scaling_reciprocal) {
     scaling_->ElementWiseReciprocal();
   }
 }

  void EquilibrationScaling::DetermineScalingParametersImpl(
    const SmartPtr<const VectorSpace> x_space,
    const SmartPtr<const VectorSpace> c_space,
    const SmartPtr<const VectorSpace> d_space,
    const SmartPtr<const MatrixSpace> jac_c_space,
    const SmartPtr<const MatrixSpace> jac_d_space,
    const SmartPtr<const SymMatrixSpace> h_space,
    const Matrix& Px_L, const Vector& x_L,
    const Matrix& Px_U, const Vector& x_U,
    Number& df,
    SmartPtr<Vector>& dx,
    SmartPtr<Vector>& dc,
    SmartPtr<Vector>& dd)
  {
    DBG_ASSERT(IsValid(nlp_));

    SmartPtr<Vector> x0 = x_space->MakeNew();
    if (!nlp_->GetStartingPoint(GetRawPtr(x0), true,
                                NULL, false,
                                NULL, false,
                                NULL, false,
                                NULL, false)) {
      THROW_EXCEPTION(FAILED_INITIALIZATION,
                      "Error getting initial point from NLP in EquilibrationScaling.\n");
    }

    // We store the added absolute values of the Jacobian and
    // objective function gradient in an array of sufficient size

    SmartPtr<Matrix> jac_c = jac_c_space->MakeNew();
    SmartPtr<Matrix> jac_d = jac_d_space->MakeNew();
    SmartPtr<Vector> grad_f = x_space->MakeNew();
    const Index nnz_jac_c = TripletHelper::GetNumberEntries(*jac_c);
    const Index nnz_jac_d = TripletHelper::GetNumberEntries(*jac_d);
    const Index nc = jac_c_space->NRows();
    const Index nd = jac_d_space->NRows();
    const Index nx = x_space->Dim();
    Number* avrg_values = new Number[nnz_jac_c+nnz_jac_d+nx];
    Number* val_buffer = new Number[Max(nnz_jac_c,nnz_jac_d,nx)];

    SmartPtr<PointPerturber> perturber =
      new PointPerturber(*x0, point_perturbation_radius_,
                         Px_L, x_L, Px_U, x_U);

    const Index num_evals = 4;
    const Index max_num_eval_errors = 10;
    Index num_eval_errors = 0;
    for (Index ieval=0; ieval<num_evals; ieval++) {
      // Compute obj gradient and Jacobian at random perturbation point
      bool done = false;
      while (!done) {
        SmartPtr<Vector> xpert = perturber->MakeNewPerturbedPoint();
        done = (nlp_->Eval_grad_f(*xpert, *grad_f) &&
                nlp_->Eval_jac_c(*xpert, *jac_c) &&
                nlp_->Eval_jac_d(*xpert, *jac_d));
        if (!done) {
          Jnlst().Printf(J_WARNING, J_INITIALIZATION,
                         "Error evaluating first derivatives as at perturbed point for equilibration-based scaling.\n");
          num_eval_errors++;
        }
        if (num_eval_errors>max_num_eval_errors) {
          delete [] val_buffer;
          delete [] avrg_values;
          THROW_EXCEPTION(FAILED_INITIALIZATION,
                          "Too many evaluation failures during equilibiration-based scaling.");
        }
      }
      // Get the numbers out of the matrices and vectors, and add it
      // to avrg_values
      TripletHelper::FillValues(nnz_jac_c, *jac_c, val_buffer);
      if (ieval==0) {
        for (Index i=0; i<nnz_jac_c; i++) {
          avrg_values[i] = fabs(val_buffer[i]);
        }
      }
      else {
        for (Index i=0; i<nnz_jac_c; i++) {
          avrg_values[i] =+ fabs(val_buffer[i]);
        }
      }
      TripletHelper::FillValues(nnz_jac_d, *jac_d, val_buffer);
      if (ieval==0) {
        for (Index i=0; i<nnz_jac_d; i++) {
          avrg_values[nnz_jac_c+i] = fabs(val_buffer[i]);
        }
      }
      else {
        for (Index i=0; i<nnz_jac_d; i++) {
          avrg_values[nnz_jac_c+i] =+ fabs(val_buffer[i]);
        }
      }
      TripletHelper::FillValuesFromVector(nx, *grad_f, val_buffer);
      if (ieval==0) {
        for (Index i=0; i<nx; i++) {
          avrg_values[nnz_jac_c+nnz_jac_d+i] = fabs(val_buffer[i]);
        }
      }
      else {
        for (Index i=0; i<nx; i++) {
          avrg_values[nnz_jac_c+nnz_jac_d+i] =+ fabs(val_buffer[i]);
//.........这里部分代码省略.........

  void GradientScaling::DetermineScalingParametersImpl(
    const SmartPtr<const VectorSpace> x_space,
    const SmartPtr<const VectorSpace> p_space,
    const SmartPtr<const VectorSpace> c_space,
    const SmartPtr<const VectorSpace> d_space,
    const SmartPtr<const MatrixSpace> jac_c_space,
    const SmartPtr<const MatrixSpace> jac_d_space,
    const SmartPtr<const SymMatrixSpace> h_space,
    const Matrix& Px_L, const Vector& x_L,
    const Matrix& Px_U, const Vector& x_U,
    Number& df,
    SmartPtr<Vector>& dx,
    SmartPtr<Vector>& dc,
    SmartPtr<Vector>& dd)
  {
    DBG_ASSERT(IsValid(nlp_));

    SmartPtr<Vector> x = x_space->MakeNew();
    SmartPtr<Vector> p = p_space->MakeNew();
    if (!nlp_->GetStartingPoint(GetRawPtr(x), true,
				GetRawPtr(p), true,
                                NULL, false,
                                NULL, false,
                                NULL, false,
                                NULL, false)) {
      THROW_EXCEPTION(FAILED_INITIALIZATION,
                      "Error getting initial point from NLP in GradientScaling.\n");
    }

    //
    // Calculate grad_f scaling
    //
    SmartPtr<Vector> grad_f = x_space->MakeNew();
    if (nlp_->Eval_grad_f(*x, *p, *grad_f)) {
      double max_grad_f = grad_f->Amax();
      df = 1.;
      if (scaling_obj_target_gradient_ == 0.) {
        if (max_grad_f > scaling_max_gradient_) {
          df = scaling_max_gradient_ / max_grad_f;
        }
      }
      else {
        if (max_grad_f == 0.) {
          Jnlst().Printf(J_WARNING, J_INITIALIZATION,
                         "Gradient of objective function is zero at starting point.  Cannot determine scaling factor based on scaling_obj_target_gradient option.\n");
        }
        else {
          df = scaling_obj_target_gradient_ / max_grad_f;
        }
      }
      df = Max(df, scaling_min_value_);
      Jnlst().Printf(J_DETAILED, J_INITIALIZATION,
                     "Scaling parameter for objective function = %e\n", df);
    }
    else {
      Jnlst().Printf(J_WARNING, J_INITIALIZATION,
                     "Error evaluating objective gradient at user provided starting point.\n  No scaling factor for objective function computed!\n");
      df = 1.;
    }
    //
    // No x scaling
    //
    dx = NULL;

    dc = NULL;
    if (c_space->Dim()>0) {
      //
      // Calculate c scaling
      //
      SmartPtr<Matrix> jac_c = jac_c_space->MakeNew();
      if (nlp_->Eval_jac_c(*x, *p, *jac_c)) {
        dc = c_space->MakeNew();
        const double dbl_min = std::numeric_limits<double>::min();
        dc->Set(dbl_min);
        jac_c->ComputeRowAMax(*dc, false);
        Number arow_max = dc->Amax();
        if (scaling_constr_target_gradient_<=0.) {
          if (arow_max > scaling_max_gradient_) {
            dc->ElementWiseReciprocal();
            dc->Scal(scaling_max_gradient_);
            SmartPtr<Vector> dummy = dc->MakeNew();
            dummy->Set(1.);
            dc->ElementWiseMin(*dummy);
          }
          else {
            dc = NULL;
          }
        }
        else {
          dc->Set(scaling_constr_target_gradient_/arow_max);
        }
        if (IsValid(dc) && scaling_min_value_ > 0.) {
          SmartPtr<Vector> tmp = dc->MakeNew();
          tmp->Set(scaling_min_value_);
          dc->ElementWiseMax(*tmp);
        }
      }
      else {
        Jnlst().Printf(J_WARNING, J_INITIALIZATION,
                       "Error evaluating Jacobian of equality constraints at user provided starting point.\n  No scaling factors for equality constraints computed!\n");
//.........这里部分代码省略.........

 inline
 Index Vector::Dim() const
 {
   return owner_space_->Dim();
 }

  bool
  NLPBoundsRemover::GetSpaces(SmartPtr<const VectorSpace>& x_space,
                              SmartPtr<const VectorSpace>& c_space,
                              SmartPtr<const VectorSpace>& d_space,
                              SmartPtr<const VectorSpace>& x_l_space,
                              SmartPtr<const MatrixSpace>& px_l_space,
                              SmartPtr<const VectorSpace>& x_u_space,
                              SmartPtr<const MatrixSpace>& px_u_space,
                              SmartPtr<const VectorSpace>& d_l_space,
                              SmartPtr<const MatrixSpace>& pd_l_space,
                              SmartPtr<const VectorSpace>& d_u_space,
                              SmartPtr<const MatrixSpace>& pd_u_space,
                              SmartPtr<const MatrixSpace>& Jac_c_space,
                              SmartPtr<const MatrixSpace>& Jac_d_space,
                              SmartPtr<const SymMatrixSpace>& Hess_lagrangian_space,
							  bool bHostInit)
  {
    DBG_START_METH("NLPBoundsRemover::GetSpaces", dbg_verbosity);
    SmartPtr<const VectorSpace> d_space_orig;
    SmartPtr<const VectorSpace> x_l_space_orig;
    SmartPtr<const MatrixSpace> px_l_space_orig;
    SmartPtr<const VectorSpace> x_u_space_orig;
    SmartPtr<const MatrixSpace> px_u_space_orig;
    SmartPtr<const VectorSpace> d_l_space_orig;
    SmartPtr<const MatrixSpace> pd_l_space_orig;
    SmartPtr<const VectorSpace> d_u_space_orig;
    SmartPtr<const MatrixSpace> pd_u_space_orig;
    SmartPtr<const MatrixSpace> Jac_d_space_orig;

    bool retval = nlp_->GetSpaces(x_space, c_space, d_space_orig,
                                  x_l_space_orig, px_l_space_orig,
                                  x_u_space_orig, px_u_space_orig,
                                  d_l_space_orig, pd_l_space_orig,
                                  d_u_space_orig, pd_u_space_orig,
                                  Jac_c_space, Jac_d_space_orig,
                                  Hess_lagrangian_space);
    if (!retval) {
      return retval;
    }
    // Keep a copy of the expansion matrices for the x bounds
    Px_l_orig_ = px_l_space_orig->MakeNew();
    Px_u_orig_ = px_u_space_orig->MakeNew();

    // create the new d_space
    Index total_dim = d_space_orig->Dim() + x_l_space_orig->Dim() +
                      x_u_space_orig->Dim();
    SmartPtr<CompoundVectorSpace> d_space_new =
      new CompoundVectorSpace(3, total_dim);
    d_space_new->SetCompSpace(0, *d_space_orig);
    d_space_new->SetCompSpace(1, *x_l_space_orig);
    d_space_new->SetCompSpace(2, *x_u_space_orig);
    d_space = GetRawPtr(d_space_new);

    // create the new (emply) x_l and x_u spaces, and also the
    // corresponding projection matrix spaces
    x_l_space = new DenseVectorSpace(0);
    x_u_space = new DenseVectorSpace(0);
    px_l_space = new ZeroMatrixSpace(x_space->Dim(),0);
    px_u_space = new ZeroMatrixSpace(x_space->Dim(),0);

    // create the new d_l and d_u vector spaces
    total_dim = d_l_space_orig->Dim() + x_l_space_orig->Dim();
    SmartPtr<CompoundVectorSpace> d_l_space_new =
      new CompoundVectorSpace(2, total_dim);
    d_l_space_new->SetCompSpace(0, *d_l_space_orig);
    d_l_space_new->SetCompSpace(1, *x_l_space_orig);
    d_l_space = GetRawPtr(d_l_space_new);
    total_dim = d_u_space_orig->Dim() + x_u_space_orig->Dim();
    SmartPtr<CompoundVectorSpace> d_u_space_new =
      new CompoundVectorSpace(2, total_dim);
    d_u_space_new->SetCompSpace(0, *d_u_space_orig);
    d_u_space_new->SetCompSpace(1, *x_u_space_orig);
    d_u_space = GetRawPtr(d_u_space_new);

    // create the new d_l and d_u projection matrix spaces
    Index total_rows = d_space_orig->Dim() + x_l_space_orig->Dim() +
                       x_u_space_orig->Dim();
    Index total_cols = d_l_space_orig->Dim() + x_l_space_orig->Dim();
    SmartPtr<CompoundMatrixSpace> pd_l_space_new =
      new CompoundMatrixSpace(3, 2, total_rows, total_cols);
    pd_l_space_new->SetBlockRows(0, d_space_orig->Dim());
    pd_l_space_new->SetBlockRows(1, x_l_space_orig->Dim());
    pd_l_space_new->SetBlockRows(2, x_u_space_orig->Dim());
    pd_l_space_new->SetBlockCols(0, d_l_space_orig->Dim());
    pd_l_space_new->SetBlockCols(1, x_l_space_orig->Dim());
    pd_l_space_new->SetCompSpace(0, 0, *pd_l_space_orig, true);
    SmartPtr<const MatrixSpace> identity_space =
      new IdentityMatrixSpace(x_l_space_orig->Dim());
    pd_l_space_new->SetCompSpace(1, 1, *identity_space, true);
    pd_l_space = GetRawPtr(pd_l_space_new);

    total_cols = d_u_space_orig->Dim() + x_u_space_orig->Dim();
    SmartPtr<CompoundMatrixSpace> pd_u_space_new =
      new CompoundMatrixSpace(3, 2, total_rows, total_cols);
    pd_u_space_new->SetBlockRows(0, d_space_orig->Dim());
    pd_u_space_new->SetBlockRows(1, x_l_space_orig->Dim());
    pd_u_space_new->SetBlockRows(2, x_u_space_orig->Dim());
    pd_u_space_new->SetBlockCols(0, d_u_space_orig->Dim());
    pd_u_space_new->SetBlockCols(1, x_u_space_orig->Dim());
    pd_u_space_new->SetCompSpace(0, 0, *pd_u_space_orig, true);
//.........这里部分代码省略.........