C++ VectorBase::size方法代码示例

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 VectorBase<scalar,index,SizeAtCompileTime>::stCombine(const VectorBase& v1,const VectorBase& v2,VectorBase& v)
{
	v.resize(v1.size()+v2.size());
	index n = v1.size();
	for ( index i = 0 ; i < n ; ++ i )
		v[i] = v1[i];
	for ( index i = 0 ; i < v2.size() ; ++ i )
		v[i+n] = v2[i];
}

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

size_t ADFun<Base>::SparseHessian(
	const VectorBase&     x    ,
	const VectorBase&     w    ,
	const VectorSet&      p    ,
	const VectorSize&     row  ,
	const VectorSize&     col  ,
	VectorBase&           hes  ,
	sparse_hessian_work&  work )
{
	size_t n    = Domain();
# ifndef NDEBUG
	size_t k, K = hes.size();
	CPPAD_ASSERT_KNOWN(
		size_t(x.size()) == n ,
		"SparseHessian: size of x not equal domain dimension for f."
	); 
	CPPAD_ASSERT_KNOWN(
		size_t(row.size()) == K && size_t(col.size()) == K ,
		"SparseHessian: either r or c does not have the same size as ehs."
	); 
	CPPAD_ASSERT_KNOWN(
		work.color.size() == 0 || work.color.size() == n,
		"SparseHessian: invalid value in work."
	);
	for(k = 0; k < K; k++)
	{	CPPAD_ASSERT_KNOWN(
			row[k] < n,
			"SparseHessian: invalid value in r."
		);
		CPPAD_ASSERT_KNOWN(
			col[k] < n,
			"SparseHessian: invalid value in c."
		);
	}
	if( work.color.size() != 0 )
		for(size_t j = 0; j < n; j++) CPPAD_ASSERT_KNOWN(
			work.color[j] <= n,
			"SparseHessian: 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, n, n, transpose);
	}
	size_t n_sweep = SparseHessianCompute(x, w, s, row, col, hes, work);
	return n_sweep;
}

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

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

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

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         , 
	size_t              r         , 
	const VectorBase&   xq        )
{	// temporary indices
	size_t i, j, ell;

	// 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( q > 0, "Forward(q, r, xq): q == 0" );
	CPPAD_ASSERT_KNOWN(
		size_t(xq.size()) == r * n,
		"Forward(q, r, xq): xq.size() is not equal r * n"
	);
	CPPAD_ASSERT_KNOWN(
		q <= num_order_taylor_ ,
		"Forward(q, r, xq): Number of Taylor coefficient orders stored in"
		" this ADFun is less than q"
	);  
	CPPAD_ASSERT_KNOWN(
		q == 1 || num_direction_taylor_ == r ,
		"Forward(q, r, xq): q > 1 and number of Taylor directions r"
		" is not same as previous Forward(1, r, xq)"
	);
		
	// does taylor_ need more orders or new number of directions
	if( cap_order_taylor_ <= q || num_direction_taylor_ != r )
	{	if( num_direction_taylor_ != r )
			num_order_taylor_ = 1;

		size_t c = std::max(q + 1, cap_order_taylor_);
		capacity_order(c, r);
	}
	CPPAD_ASSERT_UNKNOWN( cap_order_taylor_ > q );
	CPPAD_ASSERT_UNKNOWN( num_direction_taylor_ == r )

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

		for(ell = 0; ell < r; ell++)
		{	size_t index = ((c-1)*r + 1)*ind_taddr_[j] + (q-1)*r + ell + 1;
			taylor_[ index ] = xq[ r * j + ell ];
		}
	}

	// evaluate the derivatives
	CPPAD_ASSERT_UNKNOWN( cskip_op_.size() == play_.num_op_rec() );
	CPPAD_ASSERT_UNKNOWN( load_op_.size()  == play_.num_load_op_rec() );
	forward2sweep(
		q, 
		r, 
		n, 
		num_var_tape_, 
		&play_, 
		c, 
		taylor_.data(), 
		cskip_op_.data(), 
		load_op_
	);

	// return Taylor coefficients for dependent variables
	VectorBase yq;
	yq.resize(r * m);
	for(i = 0; i < m; i++)
	{	CPPAD_ASSERT_UNKNOWN( dep_taddr_[i] < num_var_tape_  );
		for(ell = 0; ell < r; ell++)
		{	size_t index = ((c-1)*r + 1)*dep_taddr_[i] + (q-1)*r + ell + 1;
			yq[ r * i + ell ] = taylor_[ index ];
		}
	}
# ifndef NDEBUG
	if( check_for_nan_ )
	{	bool ok = true;
		for(i = 0; i < m; i++)
		{	for(ell = 0; ell < r; ell++)
			{	// Studio 2012, CppAD required in front of isnan ?
				ok &= ! CppAD::isnan( yq[ r * i + ell ] );
			}
		}
		CPPAD_ASSERT_KNOWN(ok,
		"yq = f.Forward(q, r, xq): has a non-zero order Taylor coefficient\n"
		"with the value nan (but zero order coefficients are not nan)."
		);
	}
# endif

//.........这里部分代码省略.........

size_t ADFun<Base>::SparseJacobianFor(
	const VectorBase&            x           ,
	      VectorSet&             p_transpose ,
	const VectorSize&            row         ,
	const VectorSize&            col         ,
	      VectorBase&            jac         ,
	       sparse_jacobian_work& work        )
{
	size_t j, k, ell;

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

	size_t m = Range();
	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 Jacobian that are required
	size_t K = size_t(jac.size());
	CPPAD_ASSERT_UNKNOWN( row.size() == K );
	CPPAD_ASSERT_UNKNOWN( col.size() == K );

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

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

	if( color.size() == 0 )
	{
		CPPAD_ASSERT_UNKNOWN( p_transpose.n_set() ==  n );
		CPPAD_ASSERT_UNKNOWN( p_transpose.end() ==  m );

		// execute coloring algorithm
		color.resize(n);
		if(	work.color_method == "cppad" )
			color_general_cppad(p_transpose, col, row, color);
		else if( work.color_method == "colpack" )
		{
# if CPPAD_HAS_COLPACK
			color_general_colpack(p_transpose, col, row, color);
# else
			CPPAD_ASSERT_KNOWN(
				false,
				"SparseJacobianForward: work.color_method = colpack "
				"and colpack_prefix missing from cmake command line."
			);
# endif
		}
		else CPPAD_ASSERT_KNOWN(
			false,
			"SparseJacobianForward: work.color_method is not valid."
		);

		// put sorting indices in color order
		VectorSize key(K);
		order.resize(K);
		for(k = 0; k < K; k++)
			key[k] = color[ col[k] ];
		index_sort(key, order);
	}
	size_t n_color = 1;
	for(j = 0; j < n; j++) if( color[j] < n )
		n_color = std::max(n_color, color[j] + 1);

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

# if CPPAD_SPARSE_JACOBIAN_MAX_MULTIPLE_DIRECTION == 1
	// direction vector and return values for calls to forward
	VectorBase dx(n), dy(m);

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

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

		// set the corresponding components of the result
		while( k < K && color[ col[order[k]] ] == ell ) 
		{	jac[ order[k] ] = dy[row[order[k]]];
//.........这里部分代码省略.........

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();
	size_t k, K = jac.size();
	if( work.user_row.size() == 0 )
	{	// create version of (row, col, k) sorted by row row value
		work.user_col.resize(K+1);
		work.user_row.resize(K+1);
		work.sort_row.resize(K+1);

		// put sorted indices in user_row and user_col
		for(k = 0; k < K; k++)
		{	work.user_row[k] = row[k];
			work.user_col[k] = col[k];
		}
		work.user_row[K] = m;
		work.user_col[K] = n;

		// put sorting indices in sort_row
		index_sort(work.user_row, work.sort_row);
	}
# ifndef NDEBUG
	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.user_row.size() == K+1 &&
		work.user_col.size() == K+1 &&
		work.sort_row.size() == K+1 ,
		"SparseJacobianReverse: invalid value in work."
	);
	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."
		);
		CPPAD_ASSERT_KNOWN(
			work.sort_row[k] < K,
			"SparseJacobianReverse: invalid value in work."
		);
		CPPAD_ASSERT_KNOWN(
			work.user_row[k] == row[k],
			"SparseJacobianReverse: invalid value in work."
		);
		CPPAD_ASSERT_KNOWN(
			work.user_col[k] == col[k],
			"SparseJacobianReverse: invalid value in work."
		);
	}
	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, jac, work);
	return n_sweep;
}

size_t ADFun<Base>::SparseJacobianRev(
	const VectorBase&     x           ,
	VectorSet&            p           ,
	VectorBase&           jac         ,
	sparse_jacobian_work& work        )
{
	using   CppAD::vectorBool;
	size_t i, j, k, ell;

	CppAD::vector<size_t>& row(work.user_row);
	CppAD::vector<size_t>& col(work.user_col);
	CppAD::vector<size_t>& sort_row(work.sort_row);
	CppAD::vector<size_t>& color(work.color);

	size_t m = Range();
	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() == m || color.size() == 0 );

	// number of components of Jacobian that are required
	size_t K = size_t(jac.size());
	CPPAD_ASSERT_UNKNOWN( row.size() == K+1 );
	CPPAD_ASSERT_UNKNOWN( col.size() == K+1 );
	CPPAD_ASSERT_UNKNOWN( row[K] == m );
	CPPAD_ASSERT_UNKNOWN( col[K] == n );

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

	if( color.size() == 0 )
	{	CPPAD_ASSERT_UNKNOWN( p.n_set() ==  m );
		CPPAD_ASSERT_UNKNOWN( p.end() ==  n );

		// rows and columns that are in the returned jacobian
		VectorSet r_used, c_used;
		r_used.resize(n, m);
		c_used.resize(m, n);
		k = 0;
		while( k < K )
		{	CPPAD_ASSERT_UNKNOWN( 
				row[sort_row[k]] < m && col[sort_row[k]] < n 
			);
			CPPAD_ASSERT_UNKNOWN( 
				k == 0 || row[sort_row[k-1]] <= row[sort_row[k]] 
			);
			CPPAD_ASSERT_KNOWN(
				p.is_element(row[sort_row[k]], col[sort_row[k]]) ,
				"SparseJacobianReverse: "
				"an (row, col) pair is not in sparsity pattern."
			);
			r_used.add_element(col[sort_row[k]], row[sort_row[k]]);
			c_used.add_element(row[sort_row[k]], col[sort_row[k]]);
			k++;
		}
	
		// given a column index, which rows are non-zero and not used
		VectorSet not_used;
		not_used.resize(n, m);
		for(i = 0; i < m; i++)
		{	p.begin(i);
			j = p.next_element();
			while( j != p.end() )
			{	if( ! r_used.is_element(j , i) )
					not_used.add_element(j, i);
				j = p.next_element();
			}
		}
	
		// initial coloring
		color.resize(m);
		for(i = 0; i < m; i++)
			color[i] = i;
	
		// See GreedyPartialD2Coloring Algorithm Section 3.6.2 of
		// Graph Coloring in Optimization Revisited by
		// Assefaw Gebremedhin, Fredrik Maane, Alex Pothen
		vectorBool forbidden(m);
		for(i = 1; i < m; i++)
		{
			// initial all colors as ok for this row
			// (value of forbidden for ell > i does not matter)
			for(ell = 0; ell <= i; ell++)
				forbidden[ell] = false;
	
			// -----------------------------------------------------
			// Forbid colors for which this row would destroy results 
			// for each column that is non-zero for this row
			p.begin(i);
			j = p.next_element();
			while( j != p.end() )
			{	// for each row that this column uses
				r_used.begin(j);
//.........这里部分代码省略.........

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

	Base *Partial = CPPAD_NULL;
	Partial       = CPPAD_TRACK_NEW_VEC(total_num_var_ * p, Partial);

	// 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(
		w.size() == m,
		"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] = Base(0);

	// 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_ );
		Partial[dep_taddr_[i] * p + p - 1] += w[i];
	}

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

	// 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)}
		for(k = 0; k < p; k++)
			value[j * p + k ] = 
				Partial[ind_taddr_[j] * p + p - 1 - k];
	}

	// done with the Partial array
	CPPAD_TRACK_DEL_VEC(Partial);

	return value;
}