C++ BaseVector类代码示例

void Vector :: GetPart (int startpos, BaseVector & v2) const
{
  Vector & hv2 = v2.CastToVector();

  if (Length() >= startpos + v2.Length() - 1)
    {
      const double * p1 = &Get(startpos);
      double * p2 = &hv2.Elem(1);
      
      memcpy (p2, p1, hv2.Length() * sizeof(double));
    }
  else
    MyError ("Vector::GetPart: Vector to short");
}

bool HostMatrixCOO<ValueType>::Permute(const BaseVector<int> &permutation) {

  assert(&permutation != NULL);

  // symmetric permutation only
  assert( (permutation.get_size() == this->get_nrow()) &&
          (permutation.get_size() == this->get_ncol()) );

  const HostVector<int> *cast_perm = dynamic_cast<const HostVector<int>*> (&permutation) ;   
  assert(cast_perm != NULL);

  HostMatrixCOO<ValueType> src(this->local_backend_);       
  src.AllocateCOO(this->get_nnz(), this->get_nrow(), this->get_ncol());
  src.CopyFrom(*this);

  omp_set_num_threads(this->local_backend_.OpenMP_threads);  

#pragma omp parallel for      
  for (int i=0; i<this->get_nnz(); ++i) {
 
    this->mat_.row[i] = cast_perm->vec_[ src.mat_.row[i] ];
    this->mat_.col[i] = cast_perm->vec_[ src.mat_.col[i] ];

  }

  return true;

}

void HostVector<ValueType>::CopyFrom(const BaseVector<ValueType> &vec) {

  if (this != &vec)  {

    if (const HostVector<ValueType> *cast_vec = dynamic_cast<const HostVector<ValueType>*> (&vec)) {

      if (this->size_ == 0)
        this->Allocate(cast_vec->size_);

      assert(cast_vec->size_ == this->size_);

      _set_omp_backend_threads(this->local_backend_, this->size_);

#pragma omp parallel for
      for (int i=0; i<this->size_; ++i)
        this->vec_[i] = cast_vec->vec_[i];

    } else {

      // non-host type

    vec.CopyTo(this);

    }

  }

}

void Vector :: AddElementVector (const ARRAY<INDEX> & pnum,
				 const BaseVector & elvec)
{
  int i;
  const Vector & helvec = elvec.CastToVector();
  for (i = 1; i <= pnum.Size(); i++)
    Elem(pnum.Get(i)) += helvec.Get(i);
}

void Vector :: GetElementVector (const ARRAY<INDEX> & pnum,
				 BaseVector & elvec) const
{
  int i;
  Vector & helvec = elvec.CastToVector();
  for (i = 1; i <= pnum.Size(); i++)
    helvec.Elem(i) = Get(pnum.Get(i));
}

void Vector :: AddPart (int startpos, double val, const BaseVector & v2)
{
  const Vector & hv2 = v2.CastToVector();
  INDEX i;
  INDEX n = v2.Length();

  if (Length() >= startpos + n - 1)
    {
      double * p1 = &Elem(startpos);
      const double * p2 = &hv2.Get(1);

      for (i = 1; i <= n; i++)
	{
	  (*p1) += val * (*p2);
	  p1++;
	  p2++;
	}
    }
  else
    MyError ("Vector::AddPart: Vector to short");
}

void HostVector<ValueType>::CopyToHost(
                        BaseVector<ValueType> &otherVector) const
{

    DEBUGLOG(this, "HostVector::CopyToHost()", "Vec = " << &otherVector, 2);
    otherVector.Allocate(this->GetSize());
    const HostVector<ValueType> *cast_vec = 
        dynamic_cast<const HostVector<ValueType>*> (&otherVector);

    for (int i=0; i<this->mSize; i++)
    {
        cast_vec->mData[i] = mData[i];
    }
    DEBUGEND();

}

double Vector :: operator* (const BaseVector & v2) const
  {
  const Vector & hv2 = v2.CastToVector();

  double sum = 0;
  double * p1 = data;
  double * p2 = hv2.data;

  if (Length() == hv2.Length())
    {
    for (INDEX i = Length(); i > 0; i--)
      {
      sum += (*p1) * (*p2);
      p1++; p2++;
      }
    }
  else
    (*myerr) << "Scalarproduct illegal dimension" << endl;
  return sum;
  }

bool HostMatrixCOO<ValueType>::PermuteBackward(const BaseVector<int> &permutation) {

    assert(&permutation != NULL);

    // symmetric permutation only
    assert( (permutation.get_size() == this->nrow_) &&
            (permutation.get_size() == this->ncol_) );

    const HostVector<int> *cast_perm = dynamic_cast<const HostVector<int>*> (&permutation);
    assert(cast_perm != NULL);

    HostMatrixCOO<ValueType> src(this->local_backend_);
    src.AllocateCOO(this->nnz_, this->nrow_, this->ncol_);
    src.CopyFrom(*this);

    _set_omp_backend_threads(this->local_backend_, this->nnz_);

    // TODO
    // Is there a better way?
    int *pb = NULL;
    allocate_host(this->nrow_, &pb);

    #pragma omp parallel for
    for (int i=0; i<this->nrow_; ++i)
        pb [ cast_perm->vec_[i] ] = i;

    #pragma omp parallel for
    for (int i=0; i<this->nnz_; ++i) {

        this->mat_.row[i] = pb[src.mat_.row[i]];
        this->mat_.col[i] = pb[src.mat_.col[i]];

    }

    free_host(&pb);

    return true;

}

int opt_val_hes(
	const BaseVector&   x     , 
	const BaseVector&   y     , 
	Fun                 fun   , 
	BaseVector&         jac   ,
	BaseVector&         hes   )
{	// determine the base type
	typedef typename BaseVector::value_type Base;

	// check that BaseVector is a SimpleVector class with Base elements
	CheckSimpleVector<Base, BaseVector>();

	// determine the AD vector type
	typedef typename Fun::ad_vector ad_vector;

	// check that ad_vector is a SimpleVector class with AD<Base> elements
	CheckSimpleVector< AD<Base> , ad_vector >();

	// size of the x and y spaces
	size_t n = size_t(x.size());
	size_t m = size_t(y.size());

	// number of terms in the summation
	size_t ell = fun.ell();

	// check size of return values
	CPPAD_ASSERT_KNOWN(
		size_t(jac.size()) == n || jac.size() == 0,
		"opt_val_hes: size of the vector jac is not equal to n or zero"
	);
	CPPAD_ASSERT_KNOWN(
		size_t(hes.size()) == n * n || hes.size() == 0,
		"opt_val_hes: size of the vector hes is not equal to n * n or zero"
	);

	// some temporary indices
	size_t i, j, k;

	// AD version of S_k(x, y)
	ad_vector s_k(1);

	// ADFun version of S_k(x, y)
	ADFun<Base> S_k;

	// AD version of x 
	ad_vector a_x(n);

	// AD version of y
	ad_vector a_y(n);

	if( jac.size() > 0  )
	{	// this is the easy part, computing the V^{(1)} (x) which is equal 
		// to \partial_x F (x, y) (see Thoerem 2 of the reference).

		// copy x and y to AD version
		for(j = 0; j < n; j++)
			a_x[j] = x[j];
		for(j = 0; j < m; j++)
			a_y[j] = y[j];

		// initialize summation
		for(j = 0; j < n; j++)
			jac[j] = Base(0.);

		// add in \partial_x S_k (x, y)
		for(k = 0; k < ell; k++)
		{	// start recording
			Independent(a_x);
			// record
			s_k[0] = fun.s(k, a_x, a_y);
			// stop recording and store in S_k
			S_k.Dependent(a_x, s_k);
			// compute partial of S_k with respect to x
			BaseVector jac_k = S_k.Jacobian(x);
			// add \partial_x S_k (x, y) to jac
			for(j = 0; j < n; j++)
				jac[j] += jac_k[j];			 
		}
	}
	// check if we are done
	if( hes.size() == 0 )
		return 0;

	/*
	In this case, we need to compute the Hessian. Using Theorem 1 of the
	reference:
		Y^{(1)}(x) = - F_yy (x, y)^{-1} F_yx (x, y)
	Using Theorem 2 of the reference:
		V^{(2)}(x) = F_xx (x, y) + F_xy (x, y)  Y^{(1)}(x) 
	*/
	// Base and AD version of xy
	BaseVector xy(n + m);
	ad_vector a_xy(n + m);
	for(j = 0; j < n; j++)
		a_xy[j] = xy[j] = x[j]; 
	for(j = 0; j < m; j++)
		a_xy[n+j] = xy[n+j] = y[j]; 

	// Initialization summation for Hessian of F
	size_t nm_sq = (n + m) * (n + m);
//.........这里部分代码省略.........

BaseVector ADFun<Base,RecBase>::Reverse(size_t q, const BaseVector &w)
{   // used to identify the RecBase type in calls to sweeps
    RecBase not_used_rec_base;

    // constants
    const Base zero(0);

    // temporary indices
    size_t i, j, k;

    // number of independent variables
    size_t n = ind_taddr_.size();

    // number of dependent variables
    size_t m = dep_taddr_.size();

    // check BaseVector is Simple Vector class with Base type elements
    CheckSimpleVector<Base, BaseVector>();

    CPPAD_ASSERT_KNOWN(
        size_t(w.size()) == m || size_t(w.size()) == (m * q),
        "Argument w to Reverse does not have length equal to\n"
        "the dimension of the range or dimension of range times q."
    );
    CPPAD_ASSERT_KNOWN(
        q > 0,
        "The first argument to Reverse must be greater than zero."
    );
    CPPAD_ASSERT_KNOWN(
        num_order_taylor_ >= q,
        "Less than q Taylor coefficients are currently stored"
        " in this ADFun object."
    );
    // special case where multiple forward directions have been computed,
    // but we are only using the one direction zero order results
    if( (q == 1) & (num_direction_taylor_ > 1) )
    {   num_order_taylor_ = 1;        // number of orders to copy
        size_t c = cap_order_taylor_; // keep the same capacity setting
        size_t r = 1;                 // only keep one direction
        capacity_order(c, r);
    }
    CPPAD_ASSERT_KNOWN(
        num_direction_taylor_ == 1,
        "Reverse mode for Forward(q, r, xq) with more than one direction"
        "\n(r > 1) is not yet supported for q > 1."
    );

    // initialize entire Partial matrix to zero
    local::pod_vector_maybe<Base> Partial(num_var_tape_ * q);
    for(i = 0; i < num_var_tape_; i++)
        for(j = 0; j < q; j++)
            Partial[i * q + j] = zero;

    // set the dependent variable direction
    // (use += because two dependent variables can point to same location)
    for(i = 0; i < m; i++)
    {   CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_  );
        if( size_t(w.size()) == m )
            Partial[dep_taddr_[i] * q + q - 1] += w[i];
        else
        {   for(k = 0; k < q; k++)
                // ? should use += here, first make test to demonstrate bug
                Partial[ dep_taddr_[i] * q + k ] = w[i * q + k ];
        }
    }

    // evaluate the derivatives
    CPPAD_ASSERT_UNKNOWN( cskip_op_.size() == play_.num_op_rec() );
    CPPAD_ASSERT_UNKNOWN( load_op_.size()  == play_.num_load_op_rec() );
    local::play::const_sequential_iterator play_itr = play_.end();
    local::sweep::reverse(
        q - 1,
        n,
        num_var_tape_,
        &play_,
        cap_order_taylor_,
        taylor_.data(),
        q,
        Partial.data(),
        cskip_op_.data(),
        load_op_,
        play_itr,
        not_used_rec_base
    );

    // return the derivative values
    BaseVector value(n * q);
    for(j = 0; j < n; j++)
    {   CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_  );

        // independent variable taddr equals its operator taddr
        CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == local::InvOp );

        // by the Reverse Identity Theorem
        // partial of y^{(k)} w.r.t. u^{(0)} is equal to
        // partial of y^{(q-1)} w.r.t. u^{(q - 1 - k)}
        if( size_t(w.size()) == m )
        {   for(k = 0; k < q; k++)
                value[j * q + k ] =
                    Partial[ind_taddr_[j] * q + q - 1 - k];
//.........这里部分代码省略.........

void ADFun<Base,RecBase>::subgraph_jac_rev(
    const BoolVector&                   select_domain  ,
    const BoolVector&                   select_range   ,
    const BaseVector&                   x              ,
    sparse_rcv<SizeVector, BaseVector>& matrix_out     )
{   size_t m = Range();
    size_t n = Domain();
    //
    // point at which we are evaluating Jacobian
    Forward(0, x);
    //
    // nnz and row, column, and row_major vectors for subset
    local::pod_vector<size_t> row_out;
    local::pod_vector<size_t> col_out;
    local::pod_vector_maybe<Base>   val_out;
    //
    // initialize reverse mode computation on subgraphs
    subgraph_reverse(select_domain);
    //
    // memory used to hold subgraph_reverse results
    BaseVector dw;
    SizeVector col;
    //
    // loop through selected independent variables
    for(size_t i = 0; i < m; ++i) if( select_range[i] )
    {   // compute Jacobian and sparsity for this dependent variable
        size_t q   = 1;
        subgraph_reverse(q, i, col, dw);
        CPPAD_ASSERT_UNKNOWN( size_t( dw.size() ) == n );
        //
        // offset for this dependent variable
        size_t index = row_out.size();
        CPPAD_ASSERT_UNKNOWN( col_out.size() == index );
        CPPAD_ASSERT_UNKNOWN( val_out.size() == index );
        //
        // extend vectors to hold results for this dependent variable
        size_t col_size = size_t( col.size() );
        row_out.extend( col_size );
        col_out.extend( col_size );
        val_out.extend( col_size );
        //
        // store results for this dependent variable
        for(size_t c = 0; c < col_size; ++c)
        {   row_out[index + c] = i;
            col_out[index + c] = col[c];
            val_out[index + c] = dw[ col[c] ];
        }
    }
    //
    // create sparsity pattern corresponding to row_out, col_out
    size_t nr  = m;
    size_t nc  = n;
    size_t nnz = row_out.size();
    sparse_rc<SizeVector> pattern(nr, nc, nnz);
    for(size_t k = 0; k < nnz; ++k)
        pattern.set(k, row_out[k], col_out[k]);
    //
    // create sparse matrix
    sparse_rcv<SizeVector, BaseVector> matrix(pattern);
    for(size_t k = 0; k < nnz; ++k)
        matrix.set(k,  val_out[k]);
    //
    // return matrix
    matrix_out = matrix;
    //
    return;
}

void Vector :: Swap (BaseVector & v2)
  {
  Vector & hv2 = v2.CastToVector();
  swap (length, hv2.length);
  swap (data, hv2.data);
  }

TempVector operator/ (const BaseVector & v1, double scal) 
  {
  return (*v1.Copy()) /= scal;
  }