VectorBase ADFun<Base>::SparseHessian(
	const VectorBase& x, const VectorBase& w, const VectorSet& p
)
{	size_t i, j, k;

	size_t n = Domain();
	VectorBase hes(n * n);

	CPPAD_ASSERT_KNOWN(
		size_t(x.size()) == n,
		"SparseHessian: size of x not equal domain size for f."
	);

	typedef typename VectorSet::value_type Set_type;
	typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;

	// initialize the return value as zero
	Base zero(0);
	for(i = 0; i < n; i++)
		for(j = 0; j < n; j++)
			hes[i * n + j] = zero;

	// arguments to SparseHessianCompute
	Pattern_type          s;
	CppAD::vector<size_t> row;
	CppAD::vector<size_t> col; 
	sparse_hessian_work   work;
	bool transpose = false;
	sparsity_user2internal(s, p, n, n, transpose);
	k = 0;
	for(i = 0; i < n; i++)
	{	s.begin(i);
		j = s.next_element();
		while( j != s.end() )
		{	row.push_back(i);
			col.push_back(j);
			k++;
			j = s.next_element();
		}
	}
	size_t K = k;
	VectorBase H(K);

	// now we have folded this into the following case
	SparseHessianCompute(x, w, s, row, col, H, work);

	// now set the non-zero return values
	for(k = 0; k < K; k++)
		hes[ row[k] * n + col[k] ] = H[k];

	return hes;
}

 void print_int_vector(const VectorBase& vec, FILE* fp)
 {
   for(System::const_iterator it = vec.begin(); it != vec.end(); ++it) {
     fprintf(fp, "%d ", (int)*it);
   }
   fprintf(fp, "\n");
 }

VectorBase<scalar,index,SizeAtCompileTime> VectorBase<scalar,index,SizeAtCompileTime>::operator* (scalar s)
{
	VectorBase<scalar,index,Dynamic> result;
	result.setArray(this->size(),this->dataPtr(),this->interval());
	result.scale(s);
	return result;
}

  VectorBase subset(const VectorBase& x, size_t tapeid, int p=1){
    VectorBase y;
    y.resize(vecind(tapeid).size()*p);
    for(int i=0;i<y.size()/p;i++)
      for(int j=0;j<p;j++)
	{y(i*p+j)=x(vecind(tapeid)[i]*p+j);}
    return y;
  }

typename VectorBase<scalar,index,SizeAtCompileTime>::scalar VectorBase<scalar,index,SizeAtCompileTime>::dot(const VectorBase<scalar,index,S> &vec) const
{
	if(this->size()!=vec.size()) 
	{
		std::cerr<<"one size is "<<this->size()<<" and another size is "<<vec.size()<<std::endl;
		throw COException("VectorBase dot operation error: the length of two VectorBases do not equal to each other");
	}
	else{
		scalar sum = blas::copt_blas_dot(this->size(),this->dataPtr(),this->interval(),vec.dataPtr(),vec.interval());
		return sum;
	}
}

size_t ADFun<Base>::SparseJacobianReverse(
	const VectorBase&     x    ,
	const VectorSet&      p    ,
	const VectorSize&     row  ,
	const VectorSize&     col  ,
	VectorBase&           jac  ,
	sparse_jacobian_work& work )
{
	size_t m = Range();
	size_t n = Domain();
# ifndef NDEBUG
	size_t k, K = jac.size();
	CPPAD_ASSERT_KNOWN(
		size_t(x.size()) == n ,
		"SparseJacobianReverse: size of x not equal domain dimension for f."
	); 
	CPPAD_ASSERT_KNOWN(
		size_t(row.size()) == K && size_t(col.size()) == K ,
		"SparseJacobianReverse: either r or c does not have "
		"the same size as jac."
	); 
	CPPAD_ASSERT_KNOWN(
		work.color.size() == 0 || work.color.size() == m,
		"SparseJacobianReverse: invalid value in work."
	);
	for(k = 0; k < K; k++)
	{	CPPAD_ASSERT_KNOWN(
			row[k] < m,
			"SparseJacobianReverse: invalid value in r."
		);
		CPPAD_ASSERT_KNOWN(
			col[k] < n,
			"SparseJacobianReverse: invalid value in c."
		);
	}
	if( work.color.size() != 0 )
		for(size_t i = 0; i < m; i++) CPPAD_ASSERT_KNOWN(
			work.color[i] <= m,
			"SparseJacobianReverse: invalid value in work."
	);
# endif
 
	typedef typename VectorSet::value_type Set_type;
	typedef typename internal_sparsity<Set_type>::pattern_type Pattern_type;
	Pattern_type s;
	if( work.color.size() == 0 )
	{	bool transpose = false;
		sparsity_user2internal(s, p, m, n, transpose);
	}
	size_t n_sweep = SparseJacobianRev(x, s, row, col, jac, work);
	return n_sweep;
}

void VectorBase::Copy( const VectorBase& vector, SizeType elementSize )
{
  if( this != &vector )
  {
    // release old data
    Release();
    // reserve space based on source capacity
    const SizeType capacity = vector.Capacity();
    Reserve( capacity, elementSize );
    // copy over whole data
    const SizeType wholeAllocation = sizeof(SizeType) * 2 + capacity * elementSize;
    SizeType* srcData = reinterpret_cast< SizeType* >( vector.mData );
    SizeType* dstData = reinterpret_cast< SizeType* >( mData );
    memcpy( dstData - 2, srcData - 2, wholeAllocation );
  }
}

size_t ADFun<Base>::SparseHessianCompute(
	const VectorBase&           x           ,
	const VectorBase&           w           ,
	      VectorSet&            sparsity    ,
	const VectorSize&           row         ,
	const VectorSize&           col         ,
	      VectorBase&           hes         ,
	      sparse_hessian_work&  work        )
{
	using   CppAD::vectorBool;
	size_t i, k, ell;

	CppAD::vector<size_t>& color(work.color);
	CppAD::vector<size_t>& order(work.order);

	size_t n = Domain();

	// some values
	const Base zero(0);
	const Base one(1);

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

	CPPAD_ASSERT_UNKNOWN( size_t(x.size()) == n );
	CPPAD_ASSERT_UNKNOWN( color.size() == 0 || color.size() == n );

	// number of components of Hessian that are required
	size_t K = hes.size();
	CPPAD_ASSERT_UNKNOWN( row.size() == K );
	CPPAD_ASSERT_UNKNOWN( col.size() == K );

	// Point at which we are evaluating the Hessian
	Forward(0, x);

	// check for case where nothing (except Forward above) to do
	if( K == 0 )
		return 0;

	// Rows of the Hessian (i below) correspond to the forward mode index
	// and columns (j below) correspond to the reverse mode index.
	if( color.size() == 0 )
	{
		CPPAD_ASSERT_UNKNOWN( sparsity.n_set() ==  n );
		CPPAD_ASSERT_UNKNOWN( sparsity.end() ==  n );

		// execute coloring algorithm
		color.resize(n);
		color_general_cppad(sparsity, row, col, color);

		// put sorting indices in color order
		VectorSize key(K);
		order.resize(K);
		for(k = 0; k < K; k++)
			key[k] = color[ row[k] ];
		index_sort(key, order);

	}
	size_t n_color = 1;
	for(ell = 0; ell < n; ell++) if( color[ell] < n )
		n_color = std::max(n_color, color[ell] + 1);

	// direction vector for calls to forward (rows of the Hessian)
	VectorBase u(n);

	// location for return values from reverse (columns of the Hessian)
	VectorBase ddw(2 * n);

	// initialize the return value
	for(k = 0; k < K; k++)
		hes[k] = zero;

	// loop over colors
	k = 0;
	for(ell = 0; ell < n_color; ell++)
	{	CPPAD_ASSERT_UNKNOWN( color[ row[ order[k] ] ] == ell );

		// combine all rows with this color
		for(i = 0; i < n; i++)
		{	u[i] = zero;
			if( color[i] == ell )
				u[i] = one;
		}
		// call forward mode for all these rows at once
		Forward(1, u);

		// evaluate derivative of w^T * F'(x) * u
		ddw = Reverse(2, w);

		// set the corresponding components of the result
		while( k < K && color[ row[ order[k] ] ] == ell ) 
		{	hes[ order[k] ] = ddw[ col[ order[k] ] * 2 + 1 ];
			k++;
		}
	}
	return n_color;
}

VectorBase ADFun<Base>::Forward(
	size_t              q         , 
	const VectorBase&   xq        , 
	      std::ostream& s         )
{	// 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 Vector is Simple Vector class with Base type elements
	CheckSimpleVector<Base, VectorBase>();


	CPPAD_ASSERT_KNOWN(
		size_t(xq.size()) == n || size_t(xq.size()) == n*(q+1),
		"Forward(q, xq): xq.size() is not equal n or n*(q+1)"
	);

	// lowest order we are computing
	size_t p = q + 1 - size_t(xq.size()) / n;
	CPPAD_ASSERT_UNKNOWN( p == 0 || p == q );
	CPPAD_ASSERT_KNOWN(
		q <= num_order_taylor_ || p == 0,
		"Forward(q, xq): Number of Taylor coefficient orders stored in this"
		" ADFun\nis less than q and xq.size() != n*(q+1)."
	);  
	CPPAD_ASSERT_KNOWN(
		p <= 1 || num_direction_taylor_ == 1,
		"Forward(q, xq): computing order q >= 2"
		" and number of directions is not one."
		"\nMust use Forward(q, r, xq) for this case"
	);
	// does taylor_ need more orders or fewer directions
	if( (cap_order_taylor_ <= q) | (num_direction_taylor_ != 1) )
	{	if( p == 0 )
		{	// no need to copy old values during capacity_order
			num_order_taylor_ = 0;
		}
		else	num_order_taylor_ = q;
		size_t c = std::max(q + 1, cap_order_taylor_);
		size_t r = 1;
		capacity_order(c, r);
	}
	CPPAD_ASSERT_UNKNOWN( cap_order_taylor_ > q );
	CPPAD_ASSERT_UNKNOWN( num_direction_taylor_ == 1 );

	// short hand notation for order capacity
	size_t C = cap_order_taylor_;

	// set Taylor coefficients for independent variables
	for(j = 0; j < n; j++)
	{	CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_  );

		// ind_taddr_[j] is operator taddr for j-th independent variable
		CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == InvOp );

		if( p ==  q )
			taylor_[ C * ind_taddr_[j] + q] = xq[j];
		else
		{	for(k = 0; k <= q; k++)
				taylor_[ C * ind_taddr_[j] + k] = xq[ (q+1)*j + 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() );
	if( q == 0 )
	{	forward0sweep(s, true,
			n, num_var_tape_, &play_, C, 
			taylor_.data(), cskip_op_.data(), load_op_,
			compare_change_count_,
			compare_change_number_,
			compare_change_op_index_
		);
	}
	else
	{	forward1sweep(s, true, p, q, 
			n, num_var_tape_, &play_, C, 
			taylor_.data(), cskip_op_.data(), load_op_,
			compare_change_count_,
			compare_change_number_,
			compare_change_op_index_
		);
	}

	// return Taylor coefficients for dependent variables
	VectorBase yq;
	if( p == q )
	{	yq.resize(m);
		for(i = 0; i < m; i++)
		{	CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_  );
			yq[i] = taylor_[ C * dep_taddr_[i] + q];
		}
	}
	else
//.........这里部分代码省略.........

VectorBase ADFun<Base>::Reverse(size_t q, const VectorBase &w) 
{	// 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();

	pod_vector<Base> Partial;
	Partial.extend(num_var_tape_  * q);

	// update maximum memory requirement
	// memoryMax = std::max( memoryMax, 
	// 	Memory() + num_var_tape_  * q * sizeof(Base)
	// );

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

	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 for the corresponding ADFun."
	);
	CPPAD_ASSERT_KNOWN(
		q > 0,
		"The first argument to Reverse must be greater than zero."
	);  
	CPPAD_ASSERT_KNOWN(
		num_order_taylor_ >= q,
		"Less that 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
	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() );
	ReverseSweep(
		q - 1,
		n,
		num_var_tape_,
		&play_,
		cap_order_taylor_,
		taylor_.data(),
		q,
		Partial.data(),
		cskip_op_.data(),
		load_op_
	);

	// return the derivative values
	VectorBase 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] ) == 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];
//.........这里部分代码省略.........

 size_t num_rows(VectorBase<M,S> const& vec)
 {
   return vec.getSize();
 }

SEXP attribute_hidden do_subassign_dflt(SEXP call, SEXP op, SEXP argsarg,
					SEXP rho)
{
    GCStackRoot<PairList> args(SEXP_downcast<PairList*>(argsarg));

    SEXP ignored, x, y;
    PairList* subs;
    int nsubs = SubAssignArgs(args, &x, &subs, &y);
   
    /* If there are multiple references to an object we must */
    /* duplicate it so that only the local version is mutated. */
    /* This will duplicate more often than necessary, but saves */
    /* over always duplicating. */
    if (MAYBE_SHARED(CAR(args))) {
	x = SETCAR(args, shallow_duplicate(CAR(args)));
    }

    bool S4 = IS_S4_OBJECT(x);
    SEXPTYPE xorigtype = TYPEOF(x);
    if (xorigtype == LISTSXP || xorigtype == LANGSXP)
	x = PairToVectorList(x);

    /* bug PR#2590 coerce only if null */
    if (!x)
	x = coerceVector(x, TYPEOF(y));

    switch (TYPEOF(x)) {
    case LGLSXP:
    case INTSXP:
    case REALSXP:
    case CPLXSXP:
    case STRSXP:
    case EXPRSXP:
    case VECSXP:
    case RAWSXP:
	{
	    VectorBase* xv = static_cast<VectorBase*>(x);
	    if (xv->size() == 0 && Rf_length(y) == 0)
		return x;
	    size_t nsubs = listLength(subs);
	    switch (nsubs) {
	    case 0:
		x = VectorAssign(call, x, R_MissingArg, y);
		break;
	    case 1:
		x = VectorAssign(call, x, subs->car(), y);
		break;
	    default:
		x = ArrayAssign(call, x, subs, y);
		break;
	    }
	}
	break;
    default:
	error(R_MSG_ob_nonsub, TYPEOF(x));
	break;
    }

    if (xorigtype == LANGSXP) {
	if(Rf_length(x)) {
	    GCStackRoot<PairList> xlr(static_cast<PairList*>(VectorToPairList(x)));
	    GCStackRoot<Expression> xr(ConsCell::convert<Expression>(xlr));
	    x = xr;
	} else
	    error(_("result is zero-length and so cannot be a language object"));
    }

    /* Note the setting of NAMED(x) to zero here.  This means */
    /* that the following assignment will not duplicate the value. */
    /* This works because at this point, x is guaranteed to have */
    /* at most one symbol bound to it.  It does mean that there */
    /* will be multiple reference problems if "[<-" is used */
    /* in a naked fashion. */

    SET_NAMED(x, 0);
    if (S4)
	SET_S4_OBJECT(x);
    return x;
}

VectorBase ADFun<Base>::ForTwo(
	const VectorBase   &x, 
	const VectorSize_t &j,
	const VectorSize_t &k)
{	size_t i;
	size_t j1;
	size_t k1;
	size_t l;

	size_t n = Domain();
	size_t m = Range();
	size_t p = j.size();

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

	// check VectorSize_t is Simple Vector class with size_t elements
	CheckSimpleVector<size_t, VectorSize_t>();

	CPPAD_ASSERT_KNOWN(
		x.size() == n,
		"ForTwo: Length of x not equal domain dimension for f."
	); 
	CPPAD_ASSERT_KNOWN(
		j.size() == k.size(),
		"ForTwo: Lenght of the j and k vectors are not equal."
	);
	// point at which we are evaluating the second partials
	Forward(0, x);


	// dimension the return value
	VectorBase ddy(m * p);

	// allocate memory to hold all possible diagonal Taylor coefficients
	// (for large sparse cases, this is not efficient)
	VectorBase D(m * n);

	// boolean flag for which diagonal coefficients are computed
	CppAD::vector<bool> c(n);
	for(j1 = 0; j1 < n; j1++)
		c[j1] = false;

	// direction vector in argument space
	VectorBase dx(n);
	for(j1 = 0; j1 < n; j1++)
		dx[j1] = Base(0);

	// result vector in range space
	VectorBase dy(m);

	// compute the diagonal coefficients that are needed
	for(l = 0; l < p; l++)
	{	j1 = j[l];
		k1 = k[l];
		CPPAD_ASSERT_KNOWN(
		j1 < n,
		"ForTwo: an element of j not less than domain dimension for f."
		);
		CPPAD_ASSERT_KNOWN(
		k1 < n,
		"ForTwo: an element of k not less than domain dimension for f."
		);
		size_t count = 2;
		while(count)
		{	count--;
			if( ! c[j1] )
			{	// diagonal term in j1 direction
				c[j1]  = true;
				dx[j1] = Base(1);
				Forward(1, dx);

				dx[j1] = Base(0);
				dy     = Forward(2, dx);
				for(i = 0; i < m; i++)
					D[i * n + j1 ] = dy[i];
			} 
			j1 = k1;
		}
	}
	// compute all the requested cross partials
	for(l = 0; l < p; l++)
	{	j1 = j[l];
		k1 = k[l];
		if( j1 == k1 )
		{	for(i = 0; i < m; i++)
				ddy[i * p + l] = Base(2) * D[i * n + j1];
		}
		else
		{
			// cross term in j1 and k1 directions
			dx[j1] = Base(1);
			dx[k1] = Base(1);
			Forward(1, dx);

			dx[j1] = Base(0);
			dx[k1] = Base(0);
			dy = Forward(2, dx);

			// place result in return value
//.........这里部分代码省略.........

VectorBase ADFun<Base>::RevTwo(
    const VectorBase   &x,
    const VectorSize_t &i,
    const VectorSize_t &j)
{   size_t i1;
    size_t j1;
    size_t k;
    size_t l;

    size_t n = Domain();
    size_t m = Range();
    size_t p = i.size();

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

    // check VectorSize_t is Simple Vector class with size_t elements
    CheckSimpleVector<size_t, VectorSize_t>();

    CPPAD_ASSERT_KNOWN(
        x.size() == n,
        "RevTwo: Length of x not equal domain dimension for f."
    );
    CPPAD_ASSERT_KNOWN(
        i.size() == j.size(),
        "RevTwo: Lenght of the i and j vectors are not equal."
    );
    // point at which we are evaluating the second partials
    Forward(0, x);

    // dimension the return value
    VectorBase ddw(n * p);

    // direction vector in argument space
    VectorBase dx(n);
    for(j1 = 0; j1 < n; j1++)
        dx[j1] = Base(0);

    // direction vector in range space
    VectorBase w(m);
    for(i1 = 0; i1 < m; i1++)
        w[i1] = Base(0);

    // place to hold the results of a reverse calculation
    VectorBase r(n * 2);

    // check the indices in i and j
    for(l = 0; l < p; l++)
    {   i1 = i[l];
        j1 = j[l];
        CPPAD_ASSERT_KNOWN(
            i1 < m,
            "RevTwo: an eleemnt of i not less than range dimension for f."
        );
        CPPAD_ASSERT_KNOWN(
            j1 < n,
            "RevTwo: an element of j not less than domain dimension for f."
        );
    }

    // loop over all forward directions
    for(j1 = 0; j1 < n; j1++)
    {   // first order forward mode calculation done
        bool first_done = false;
        for(l = 0; l < p; l++) if( j[l] == j1 )
            {   if( ! first_done )
                {   first_done = true;

                    // first order forward mode in j1 direction
                    dx[j1] = Base(1);
                    Forward(1, dx);
                    dx[j1] = Base(0);
                }
                // execute a reverse in this component direction
                i1    = i[l];
                w[i1] = Base(1);
                r     = Reverse(2, w);
                w[i1] = Base(0);

                // place the reverse result in return value
                for(k = 0; k < n; k++)
                    ddw[k * p + l] = r[k * 2 + 1];
            }
    }
    return ddw;
}

VectorBase ADFun<Base>::Forward(
	size_t              q         ,
	const VectorBase&   xq        ,
	      std::ostream& s         )
{	// 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 Vector is Simple Vector class with Base type elements
	CheckSimpleVector<Base, VectorBase>();


	CPPAD_ASSERT_KNOWN(
		size_t(xq.size()) == n || size_t(xq.size()) == n*(q+1),
		"Forward(q, xq): xq.size() is not equal n or n*(q+1)"
	);

	// lowest order we are computing
	size_t p = q + 1 - size_t(xq.size()) / n;
	CPPAD_ASSERT_UNKNOWN( p == 0 || p == q );
	CPPAD_ASSERT_KNOWN(
		q <= num_order_taylor_ || p == 0,
		"Forward(q, xq): Number of Taylor coefficient orders stored in this"
		" ADFun\nis less than q and xq.size() != n*(q+1)."
	);
	CPPAD_ASSERT_KNOWN(
		p <= 1 || num_direction_taylor_ == 1,
		"Forward(q, xq): computing order q >= 2"
		" and number of directions is not one."
		"\nMust use Forward(q, r, xq) for this case"
	);
	// does taylor_ need more orders or fewer directions
	if( (cap_order_taylor_ <= q) | (num_direction_taylor_ != 1) )
	{	if( p == 0 )
		{	// no need to copy old values during capacity_order
			num_order_taylor_ = 0;
		}
		else	num_order_taylor_ = q;
		size_t c = std::max(q + 1, cap_order_taylor_);
		size_t r = 1;
		capacity_order(c, r);
	}
	CPPAD_ASSERT_UNKNOWN( cap_order_taylor_ > q );
	CPPAD_ASSERT_UNKNOWN( num_direction_taylor_ == 1 );

	// short hand notation for order capacity
	size_t C = cap_order_taylor_;

	// The optimizer may skip a step that does not affect dependent variables.
	// Initilaizing zero order coefficients avoids following valgrind warning:
	// "Conditional jump or move depends on uninitialised value(s)".
	for(j = 0; j < num_var_tape_; j++)
	{	for(k = p; k <= q; k++)
			taylor_[C * j + k] = CppAD::numeric_limits<Base>::quiet_NaN();
	}

	// set Taylor coefficients for independent variables
	for(j = 0; j < n; j++)
	{	CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < num_var_tape_  );

		// ind_taddr_[j] is operator taddr for j-th independent variable
		CPPAD_ASSERT_UNKNOWN( play_.GetOp( ind_taddr_[j] ) == local::InvOp );

		if( p == q )
			taylor_[ C * ind_taddr_[j] + q] = xq[j];
		else
		{	for(k = 0; k <= q; k++)
				taylor_[ C * ind_taddr_[j] + k] = xq[ (q+1)*j + 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() );
	if( q == 0 )
	{	local::forward0sweep(s, true,
			n, num_var_tape_, &play_, C,
			taylor_.data(), cskip_op_.data(), load_op_,
			compare_change_count_,
			compare_change_number_,
			compare_change_op_index_
		);
	}
	else
	{	local::forward1sweep(s, true, p, q,
			n, num_var_tape_, &play_, C,
			taylor_.data(), cskip_op_.data(), load_op_,
			compare_change_count_,
			compare_change_number_,
			compare_change_op_index_
		);
	}

	// return Taylor coefficients for dependent variables
	VectorBase yq;
//.........这里部分代码省略.........

bool VectorTester<Scalar>::check(
  const VectorBase<Scalar>       &v
  ,Teuchos::FancyOStream         *out_arg
  ) const
{

  using std::endl;
  using Teuchos::describe;
  using Teuchos::FancyOStream;
  using Teuchos::OSTab;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  //typedef typename ST::magnitudeType    ScalarMag;

  Teuchos::RCP<FancyOStream> out = Teuchos::rcp(out_arg,false);
  const Teuchos::EVerbosityLevel verbLevel = (dump_all()?Teuchos::VERB_EXTREME:Teuchos::VERB_MEDIUM);

  OSTab tab(out,1,"THYRA");

  bool result, success = true;

  if(out.get()) *out <<endl<< "*** Entering Thyra::VectorTester<"<<ST::name()<<">::check(v,...) ...\n";

  if(out.get()) *out <<endl<< "Testing a VectorBase object described as:\n" << describe(v,verbLevel);

  if(out.get()) *out <<endl<< "A) Creating temporary vector t1, t2, t3, and t4 from v.space() ...\n";
  Teuchos::RCP<const Thyra::VectorSpaceBase<Scalar> >
    vs = v.space();
  Teuchos::RCP<Thyra::VectorBase<Scalar> >
    t1 = createMember(vs), t2 = createMember(vs), t3 = createMember(vs), t4 = createMember(vs);

  if(out.get()) *out <<endl<< "B) Testing VectorBase::applyOp(...) by calling a few standard RTOp operations ... ";

  const Scalar
    one   = ST::one(),
    two   = Scalar(2)*one,
    three = Scalar(3)*one;

  {
    using Teuchos::inoutArg;

    TestResultsPrinter testResultsPrinter(out, show_all_tests());
    const RCP<FancyOStream> testOut = testResultsPrinter.getTestOStream();

    bool these_results = true;
    
    *testOut <<endl<< "assign(t1.ptr(),2.0) ...\n";
    Thyra::assign( t1.ptr(), two );
    if(dump_all()) *testOut <<endl<< "\nt1 =\n" << describe(*t1,verbLevel);
    
    result = Teuchos::testRelErr<Scalar>(
      "sum(t1)", sum(*t1), "2*vs->dim()", two*Scalar(vs->dim()),
      "error_tol()", error_tol(), "warning_tol()", warning_tol(),
      inoutArg(*testOut)
      );
    if(!result) these_results = false;
    
    *testOut <<endl<< "assign(t2.ptr(),3.0) ...\n";
    Thyra::assign( t2.ptr(), three );
    if(dump_all()) *testOut <<endl<< "t2 =\n" << *t1;
    
    result = Teuchos::testRelErr<Scalar>(
      "sum(t2)",sum(*t2),"3*vs->dim()",three*Scalar(vs->dim()),
      "error_tol()",error_tol(),"warning_tol()",warning_tol(),
      inoutArg(*testOut)
      );
    if(!result) these_results = false;
    
    result = Teuchos::testRelErr<Scalar>(
      "vs->scalarProd(*t1,*t2)",vs->scalarProd(*t1,*t2),"2*3*vs->dim()",two*three*Scalar(vs->dim()),
      "error_tol()",error_tol(),"warning_tol()",warning_tol(),
      inoutArg(*testOut)
      );
    if(!result) these_results = false;

    testResultsPrinter.printTestResults(these_results, inoutArg(success));

  }
    
  // ToDo: Test the rest of the specific VectorBase interface on v1

  if(out.get()) *out <<endl<< "C) Checking the MultiVectorBase interface of v ...\n";
  result = multiVectorTester_.check(v, out.ptr());
  if(!result) success = false;

  if(out.get()) *out <<endl<< "*** Leaving Thyra::VectorTester<"<<ST::name()<<">::check(v,...) ...\n";
  
  return success;

}

template <> void Swap<VectorBase>(VectorBase& first, VectorBase& second)
{
    first.Swap(second);
}

	VectorBase(const VectorBase<U, N>& vector)
	{
		for(std::size_t i=0; i!=N; ++i)
			_elem[i] = T(vector.At(i));
	}

 size_t size(VectorBase<M,S> const& vec)
 {
   return vec.getSize();
 }

VectorBase ADFun<Base>::Reverse(size_t p, const VectorBase &w) 
{	// 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();

	pod_vector<Base> Partial;
	Partial.extend(total_num_var_ * p);

	// update maximum memory requirement
	// memoryMax = std::max( memoryMax, 
	// 	Memory() + total_num_var_ * p * sizeof(Base)
	// );

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

	CPPAD_ASSERT_KNOWN(
		size_t(w.size()) == m || size_t(w.size()) == (m * p),
		"Argument w to Reverse does not have length equal to\n"
		"the dimension of the range for the corresponding ADFun."
	);
	CPPAD_ASSERT_KNOWN(
		p > 0,
		"The first argument to Reverse must be greater than zero."
	);  
	CPPAD_ASSERT_KNOWN(
		taylor_per_var_ >= p,
		"Less that p taylor_ coefficients are currently stored"
		" in this ADFun object."
	);  

	// initialize entire Partial matrix to zero
	for(i = 0; i < total_num_var_; i++)
		for(j = 0; j < p; j++)
			Partial[i * p + 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] < total_num_var_ );
		if( size_t(w.size()) == m )
			Partial[dep_taddr_[i] * p + p - 1] += w[i];
		else
		{	for(k = 0; k < p; k++)
				// ? should use += here, first make test to demonstrate bug
				Partial[ dep_taddr_[i] * p + k ] = w[i * p + k ];
		}
	}

	// evaluate the derivatives
	ReverseSweep(
		p - 1,
		n,
		total_num_var_,
		&play_,
		taylor_col_dim_,
		taylor_.data(),
		p,
		Partial.data(),
		cskip_op_
	);

	// return the derivative values
	VectorBase value(n * p);
	for(j = 0; j < n; j++)
	{	CPPAD_ASSERT_UNKNOWN( ind_taddr_[j] < total_num_var_ );

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

		// by the Reverse Identity Theorem 
		// partial of y^{(k)} w.r.t. u^{(0)} is equal to
		// partial of y^{(p-1)} w.r.t. u^{(p - 1 - k)}
		if( size_t(w.size()) == m )
		{	for(k = 0; k < p; k++)
				value[j * p + k ] = 
					Partial[ind_taddr_[j] * p + p - 1 - k];
		}
		else
		{	for(k = 0; k < p; k++)
				value[j * p + k ] =
					Partial[ind_taddr_[j] * p + k];
		}
	}
	CPPAD_ASSERT_KNOWN( ! ( hasnan(value) && check_for_nan_ ) ,
		"dw = f.Reverse(p, w): has a nan,\n"
		"but none of its Taylor coefficents are nan."
	);

	return value;
}