bool link_sparse_hessian(
	size_t                           size     , 
	size_t                           repeat   , 
	CppAD::vector<double>&           x        ,
	const CppAD::vector<size_t>&     row      ,
	const CppAD::vector<size_t>&     col      ,
	CppAD::vector<double>&           hessian  )
{
	// -----------------------------------------------------
	// setup
	typedef vector<double>              DblVector;
	typedef vector< std::set<size_t> >  SetVector;
	typedef CppAD::AD<double>           ADScalar;
	typedef vector<ADScalar>            ADVector;

	size_t i, j, k;
	size_t order = 0;         // derivative order corresponding to function
	size_t m = 1;             // number of dependent variables
	size_t n = size;          // number of independent variables
	size_t K = row.size();    // number of non-zeros in lower triangle
	ADVector   a_x(n);        // AD domain space vector
	ADVector   a_y(m);        // AD range space vector
	DblVector  w(m);          // double range space vector
	DblVector hes(K);         // non-zeros in lower triangle
	CppAD::ADFun<double> f;   // AD function object

	// weights for hessian calculation (only one component of f)
	w[0] = 1.;

	// use the unspecified fact that size is non-decreasing between calls
	static size_t previous_size = 0;
	bool print    = (repeat > 1) & (previous_size != size);
	previous_size = size;

	// declare sparsity pattern
# if USE_SET_SPARSITY
	SetVector sparsity(n);
# else
	typedef vector<bool>                BoolVector;
	BoolVector sparsity(n * n);
# endif
	// initialize all entries as zero
	for(i = 0; i < n; i++)
	{	for(j = 0; j < n; j++)
			hessian[ i * n + j] = 0.;
	}
	// ------------------------------------------------------
	extern bool global_retape;
	if( global_retape) while(repeat--)
	{	// choose a value for x 
		CppAD::uniform_01(n, x);
		for(j = 0; j < n; j++)
			a_x[j] = x[j];

		// declare independent variables
		Independent(a_x);	

		// AD computation of f(x)
		CppAD::sparse_hes_fun<ADScalar>(n, a_x, row, col, order, a_y);

		// create function object f : X -> Y
		f.Dependent(a_x, a_y);

		extern bool global_optimize;
		if( global_optimize )
		{	print_optimize(f, print, "cppad_sparse_hessian_optimize", size);
			print = false;
		}

		// calculate the Hessian sparsity pattern for this function
		calc_sparsity(sparsity, f);

		// structure that holds some of work done by SparseHessian
		CppAD::sparse_hessian_work work;

		// calculate this Hessian at this x
		f.SparseHessian(x, w, sparsity, row, col, hes, work);
		for(k = 0; k < K; k++)
		{	hessian[ row[k] * n + col[k] ] = hes[k];
			hessian[ col[k] * n + row[k] ] = hes[k];
		}
	}
	else
	{	// choose a value for x 
		CppAD::uniform_01(n, x);
		for(j = 0; j < n; j++)
			a_x[j] = x[j];

		// declare independent variables
		Independent(a_x);	

		// AD computation of f(x)
		CppAD::sparse_hes_fun<ADScalar>(n, a_x, row, col, order, a_y);

		// create function object f : X -> Y
		f.Dependent(a_x, a_y);

		extern bool global_optimize;
		if( global_optimize )
		{	print_optimize(f, print, "cppad_sparse_hessian_optimize", size);
//.........这里部分代码省略.........

/* $$
$head Use Atomic Function$$
$codep */
bool get_started(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();
/* $$
$subhead Constructor$$
$codep */
	// Create the atomic get_started object
	atomic_get_started afun("atomic_get_started");
/* $$
$subhead Recording$$
$codep */
	// Create the function f(x)
	//
	// domain space vector
	size_t n  = 1;
	double  x0 = 0.5;
	vector< AD<double> > ax(n);
	ax[0]     = x0;

	// declare independent variables and start tape recording
	CppAD::Independent(ax);

	// range space vector 
	size_t m = 1;
	vector< AD<double> > ay(m);

	// call user function and store get_started(x) in au[0] 
	vector< AD<double> > au(m);
	afun(ax, au);        // u = 1 / x

	// now use AD division to invert to invert the operation
	ay[0] = 1.0 / au[0]; // y = 1 / u = x

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f;
	f.Dependent (ax, ay);  // f(x) = x
/* $$
$subhead forward$$
$codep */
	// check function value 
	double check = x0;
	ok &= NearEqual( Value(ay[0]) , check,  eps, eps);

	// check zero order forward mode
	size_t p;
	vector<double> x_p(n), y_p(m);
	p      = 0;
	x_p[0] = x0;
	y_p    = f.Forward(p, x_p);
	ok &= NearEqual(y_p[0] , check,  eps, eps);

	return ok;
}

bool LuRatio(void)
{	bool  ok = true;

	size_t  n = 2; // number rows in A 
	double  ratio;

	// values for independent and dependent variables
	CPPAD_TEST_VECTOR<double>  x(n*n), y(n*n+1);

	// pivot vectors
	CPPAD_TEST_VECTOR<size_t> ip(n), jp(n);

	// set x equal to the identity matrix
	x[0] = 1.; x[1] = 0;
	x[2] = 0.; x[3] = 1.;

	// create a fnction object corresponding to this value of x
	CppAD::ADFun<double> *FunPtr = NewFactor(n, x, ok, ip, jp);

	// use function object to factor matrix
	y     = FunPtr->Forward(0, x);
	ratio = y[n*n];
	ok   &= (ratio == 1.);
	ok   &= CheckLuFactor(n, x, y, ip, jp);

	// set x so that the pivot ratio will be infinite
	x[0] = 0. ; x[1] = 1.;
	x[2] = 1. ; x[3] = 0.;

	// try to use old function pointer to factor matrix
	y     = FunPtr->Forward(0, x);
	ratio = y[n*n];

	// check to see if we need to refactor matrix
	ok &= (ratio > 10.);
	if( ratio > 10. )
	{	delete FunPtr; // to avoid a memory leak	
		FunPtr = NewFactor(n, x, ok, ip, jp);
	}

	//  now we can use the function object to factor matrix
	y     = FunPtr->Forward(0, x);
	ratio = y[n*n];
	ok    &= (ratio == 1.);
	ok    &= CheckLuFactor(n, x, y, ip, jp);

	delete FunPtr;  // avoid memory leak
	return ok;
}

/* %$$
$head Use Atomic Function$$
$srccode%cpp% */
bool norm_sq(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();
/* %$$
$subhead Constructor$$
$srccode%cpp% */
	// --------------------------------------------------------------------
	// Create the atomic reciprocal object
	atomic_norm_sq afun("atomic_norm_sq");
/* %$$
$subhead Recording$$
$srccode%cpp% */
	// Create the function f(x)
	//
	// domain space vector
	size_t  n  = 2;
	double  x0 = 0.25;
	double  x1 = 0.75;
	vector< AD<double> > ax(n);
	ax[0]      = x0;
	ax[1]      = x1;

	// declare independent variables and start tape recording
	CppAD::Independent(ax);

	// range space vector
	size_t m = 1;
	vector< AD<double> > ay(m);

	// call user function and store norm_sq(x) in au[0]
	afun(ax, ay);        // y_0 = x_0 * x_0 + x_1 * x_1

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f;
	f.Dependent (ax, ay);
/* %$$
$subhead forward$$
$srccode%cpp% */
	// check function value
	double check = x0 * x0 + x1 * x1;
	ok &= NearEqual( Value(ay[0]) , check,  eps, eps);

	// check zero order forward mode
	size_t q;
	vector<double> x_q(n), y_q(m);
	q      = 0;
	x_q[0] = x0;
	x_q[1] = x1;
	y_q    = f.Forward(q, x_q);
	ok &= NearEqual(y_q[0] , check,  eps, eps);

	// check first order forward mode
	q      = 1;
	x_q[0] = 0.3;
	x_q[1] = 0.7;
	y_q    = f.Forward(q, x_q);
	check  = 2.0 * x0 * x_q[0] + 2.0 * x1 * x_q[1];
	ok &= NearEqual(y_q[0] , check,  eps, eps);

/* %$$
$subhead reverse$$
$srccode%cpp% */
	// first order reverse mode
	q     = 1;
	vector<double> w(m), dw(n * q);
	w[0]  = 1.;
	dw    = f.Reverse(q, w);
	check = 2.0 * x0;
	ok &= NearEqual(dw[0] , check,  eps, eps);
	check = 2.0 * x1;
	ok &= NearEqual(dw[1] , check,  eps, eps);
/* %$$
$subhead for_sparse_jac$$
$srccode%cpp% */
	// forward mode sparstiy pattern
	size_t p = n;
	CppAD::vectorBool r1(n * p), s1(m * p);
	r1[0] = true;  r1[1] = false; // sparsity pattern identity
	r1[2] = false; r1[3] = true;
	//
	s1    = f.ForSparseJac(p, r1);
	ok  &= s1[0] == true;  // f[0] depends on x[0]
	ok  &= s1[1] == true;  // f[0] depends on x[1]
/* %$$
$subhead rev_sparse_jac$$
$srccode%cpp% */
	// reverse mode sparstiy pattern
	q = m;
	CppAD::vectorBool s2(q * m), r2(q * n);
	s2[0] = true;          // compute sparsity pattern for f[0]
	//
	r2    = f.RevSparseJac(q, s2);
	ok  &= r2[0] == true;  // f[0] depends on x[0]
	ok  &= r2[1] == true;  // f[0] depends on x[1]
/* %$$
//.........这里部分代码省略.........

bool sparse_hessian(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	size_t i, j, k, ell;
	typedef CPPAD_TESTVECTOR(AD<double>)               a_vector;
	typedef CPPAD_TESTVECTOR(double)                     d_vector;
	typedef CPPAD_TESTVECTOR(size_t)                     i_vector;
	typedef CPPAD_TESTVECTOR(bool)                       b_vector;
	typedef CPPAD_TESTVECTOR(std::set<size_t>)         s_vector;
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();

	// domain space vector
	size_t n = 12;  // must be greater than or equal 3; see n_sweep below
	a_vector a_x(n);
	for(j = 0; j < n; j++)
		a_x[j] = AD<double> (0);

	// declare independent variables and starting recording
	CppAD::Independent(a_x);

	// range space vector
	size_t m = 1;
	a_vector a_y(m);
	a_y[0] = a_x[0]*a_x[1];
	for(j = 0; j < n; j++)
		a_y[0] += a_x[j] * a_x[j] * a_x[j];

	// create f: x -> y and stop tape recording
	// (without executing zero order forward calculation)
	CppAD::ADFun<double> f;
	f.Dependent(a_x, a_y);

	// new value for the independent variable vector, and weighting vector
	d_vector w(m), x(n);
	for(j = 0; j < n; j++)
		x[j] = double(j);
	w[0] = 1.0;

	// vector used to check the value of the hessian
	d_vector check(n * n);
	for(ell = 0; ell < n * n; ell++)
		check[ell] = 0.0;
	ell        = 0 * n + 1;
	check[ell] = 1.0;
	ell        = 1 * n + 0;
	check[ell] = 1.0 ;
	for(j = 0; j < n; j++)
	{	ell = j * n + j;
		check[ell] = 6.0 * x[j];
	}

	// -------------------------------------------------------------------
	// second derivative of y[0] w.r.t x
	d_vector hes(n * n);
	hes = f.SparseHessian(x, w);
	for(ell = 0; ell < n * n; ell++)
		ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );

	// --------------------------------------------------------------------
	// example using vectors of bools to compute sparsity pattern for Hessian
	b_vector r_bool(n * n);
	for(i = 0; i < n; i++)
	{	for(j = 0; j < n; j++)
			r_bool[i * n + j] = false;
		r_bool[i * n + i] = true;
	}
	f.ForSparseJac(n, r_bool);
	//
	b_vector s_bool(m);
	for(i = 0; i < m; i++)
		s_bool[i] = w[i] != 0;
	b_vector p_bool = f.RevSparseHes(n, s_bool);

	hes = f.SparseHessian(x, w, p_bool);
	for(ell = 0; ell < n * n; ell++)
		ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );

	// --------------------------------------------------------------------
	// example using vectors of sets to compute sparsity pattern for Hessian
	s_vector r_set(n);
	for(i = 0; i < n; i++)
		r_set[i].insert(i);
	f.ForSparseJac(n, r_set);
	//
	s_vector s_set(m);
	for(i = 0; i < m; i++)
		if( w[i] != 0. )
			s_set[0].insert(i);
	s_vector p_set = f.RevSparseHes(n, s_set);

	// example passing sparsity pattern to SparseHessian
	hes = f.SparseHessian(x, w, p_set);
	for(ell = 0; ell < n * n; ell++)
		ok &=  NearEqual(w[0] * check[ell], hes[ell], eps, eps );

	// --------------------------------------------------------------------
	// use row and column indices to specify upper triangle of
	// non-zero elements of Hessian
	size_t K = n + 1;
//.........这里部分代码省略.........

bool sub_sparse_hes(void)
{	bool ok = true;
	using CppAD::AD;
	typedef AD<double>   adouble;
	typedef AD<adouble> a2double;
	typedef vector< std::set<size_t> > pattern;
	double eps = 10. * std::numeric_limits<double>::epsilon();
	size_t i, j;

	// start recording with x = (u , v)
	size_t nu = 10;
	size_t nv = 5;
	size_t n  = nu + nv;
	vector<adouble> ax(n);
	for(j = 0; j < n; j++)
		ax[j] = adouble(j + 2);
	CppAD::Independent(ax);

	// extract u as independent variables
	vector<a2double> a2u(nu);
	for(j = 0; j < nu; j++)
		a2u[j] = a2double(j + 2);
	CppAD::Independent(a2u);

	// extract v as parameters
	vector<a2double> a2v(nv);
	for(j = 0; j < nv; j++)
		a2v[j] = ax[nu+j];

	// record g(u)
	vector<a2double> a2y(1);
	a2y[0] = f(a2u, a2v);
	CppAD::ADFun<adouble> g;
	g.Dependent(a2u, a2y);

	// compue sparsity pattern for Hessian of g(u)
	pattern r(nu), s(1);
	for(j = 0; j < nu; j++)
		r[j].insert(j);
	g.ForSparseJac(nu, r);
	s[0].insert(0);
	pattern p = g.RevSparseHes(nu, s);

	// Row and column indices for non-zeros in lower triangle of Hessian
	vector<size_t> row, col;
	for(i = 0; i < nu; i++)
	{	std::set<size_t>::const_iterator itr;
		for(itr = p[i].begin(); itr != p[i].end(); itr++)
		{	j = *itr;
			if( j <= i )
			{	row.push_back(i);
				col.push_back(j);
			}
		}
	}
	size_t K = row.size();
	CppAD::sparse_hessian_work work;
	vector<adouble> au(nu), ahes(K), aw(1);
	aw[0] = 1.0;
	for(j = 0; j < nu; j++)
		au[j] = ax[j];
	size_t n_sweep = g.SparseHessian(au, aw, p, row, col, ahes, work);

	// The Hessian w.r.t u is diagonal
	ok &= n_sweep == 1;

	// record H(u, v) = Hessian of f w.r.t u
	CppAD::ADFun<double> H(ax, ahes);

	// remove unecessary operations
	H.optimize();

	// Now evaluate the Hessian at a particular value for u, v
	vector<double> u(nu), v(nv), x(n);
	for(j = 0; j < n; j++)
		x[j] = double(j + 2);
	vector<double> hes = H.Forward(0, x);

	// Now check the Hessian
	double sum_v = 0.0;
	for(j = 0; j < nv; j++)
		sum_v += x[nu + j];
	for(size_t k = 0; k < K; k++)
	{	i     = row[k];
		j     = col[k];
		ok   &= i == j;
		double check = sum_v * x[i];
		ok &= CppAD::NearEqual(hes[k], check, eps, eps);
	}
	return ok;
}

/* %$$
$head Test Atomic Function$$
$srccode%cpp% */
bool set_sparsity(void)
{   bool ok = true;
    using CppAD::AD;
    using CppAD::NearEqual;
    double eps = 10. * std::numeric_limits<double>::epsilon();
/* %$$
$subhead Constructor$$
$srccode%cpp% */
    // Create the atomic get_started object
    atomic_set_sparsity afun("atomic_set_sparsity");
/* %$$
$subhead Recording$$
$srccode%cpp% */
    size_t n = 3;
    size_t m = 2;
    vector< AD<double> > ax(n), ay(m);
    for(size_t j = 0; j < n; j++)
        ax[j] = double(j + 1);

    // declare independent variables and start tape recording
    CppAD::Independent(ax);

    // call atomic function
    afun(ax, ay);

    // create f: x -> y and stop tape recording
    CppAD::ADFun<double> f;
    f.Dependent (ax, ay);  // f(x) = x

    // check function value
    ok &= NearEqual(ay[0] , ax[2],  eps, eps);
    ok &= NearEqual(ay[1] , ax[0] * ax[1],  eps, eps);

/* %$$
$subhead for_sparse_jac$$
$srccode%cpp% */
    // correct Jacobian result
    set_vector check_s(m);
    check_s[0].insert(2);
    check_s[1].insert(0);
    check_s[1].insert(1);
    // compute and test forward mode
    {   set_vector r(n), s(m);
        for(size_t i = 0; i < n; i++)
            r[i].insert(i);
        s = f.ForSparseJac(n, r);
        for(size_t i = 0; i < m; i++)
            ok &= s[i] == check_s[i];
    }
/* %$$
$subhead rev_sparse_jac$$
$srccode%cpp% */
    // compute and test reverse mode
    {   set_vector r(m), s(m);
        for(size_t i = 0; i < m; i++)
            r[i].insert(i);
        s = f.RevSparseJac(m, r);
        for(size_t i = 0; i < m; i++)
            ok &= s[i] == check_s[i];
    }
/* %$$
$subhead for_sparse_hes$$
$srccode%cpp% */
    // correct Hessian result
    set_vector check_h(n);
    check_h[0].insert(1);
    check_h[1].insert(0);
    // compute and test forward mode
    {   set_vector r(1), s(1), h(n);
        for(size_t i = 0; i < m; i++)
            s[0].insert(i);
        for(size_t j = 0; j < n; j++)
            r[0].insert(j);
        h = f.ForSparseHes(r, s);
        for(size_t i = 0; i < n; i++)
            ok &= h[i] == check_h[i];
    }
/* %$$
$subhead rev_sparse_hes$$
Note the previous call to $code ForSparseJac$$ above
stored its results in $icode f$$.
$srccode%cpp% */
    // compute and test reverse mode
    {   set_vector s(1), h(n);
        for(size_t i = 0; i < m; i++)
            s[0].insert(i);
        h = f.RevSparseHes(n, s);
        for(size_t i = 0; i < n; i++)
            ok &= h[i] == check_h[i];
    }
/* %$$
$subhead Test Result$$
$srccode%cpp% */
    return ok;
}

/* %$$
$head Use Atomic Function$$
$srccode%cpp% */
bool forward(void)
{   bool ok = true;
    using CppAD::AD;
    using CppAD::NearEqual;
    double eps = 10. * CppAD::numeric_limits<double>::epsilon();
    //
    // Create the atomic_forward object corresponding to g(x)
    atomic_forward afun("atomic_forward");
    //
    // Create the function f(u) = g(u) for this example.
    //
    // domain space vector
    size_t n  = 3;
    double u_0 = 1.00;
    double u_1 = 2.00;
    double u_2 = 3.00;
    vector< AD<double> > au(n);
    au[0] = u_0;
    au[1] = u_1;
    au[2] = u_2;

    // declare independent variables and start tape recording
    CppAD::Independent(au);

    // range space vector
    size_t m = 2;
    vector< AD<double> > ay(m);

    // call atomic function
    vector< AD<double> > ax = au;
    afun(ax, ay);

    // create f: u -> y and stop tape recording
    CppAD::ADFun<double> f;
    f.Dependent (au, ay);  // y = f(u)
    //
    // check function value
    double check = u_2 * u_2;
    ok &= NearEqual( Value(ay[0]) , check,  eps, eps);
    check = u_0 * u_1;
    ok &= NearEqual( Value(ay[1]) , check,  eps, eps);

    // --------------------------------------------------------------------
    // zero order forward
    //
    vector<double> u0(n), y0(m);
    u0[0] = u_0;
    u0[1] = u_1;
    u0[2] = u_2;
    y0   = f.Forward(0, u0);
    check = u_2 * u_2;
    ok &= NearEqual(y0[0] , check,  eps, eps);
    check = u_0 * u_1;
    ok &= NearEqual(y0[1] , check,  eps, eps);
    // --------------------------------------------------------------------
    // first order forward
    //
    // value of Jacobian of f
    double check_jac[] = {
        0.0, 0.0, 2.0 * u_2,
        u_1, u_0,       0.0
    };
    vector<double> u1(n), y1(m);
    // check first order forward mode
    for(size_t j = 0; j < n; j++)
        u1[j] = 0.0;
    for(size_t j = 0; j < n; j++)
    {   // compute partial in j-th component direction
        u1[j] = 1.0;
        y1    = f.Forward(1, u1);
        u1[j] = 0.0;
        // check this direction
        for(size_t i = 0; i < m; i++)
            ok &= NearEqual(y1[i], check_jac[i * n + j], eps, eps);
    }
    // --------------------------------------------------------------------
    // second order forward
    //
    // value of Hessian of g_0
    double check_hes_0[] = {
        0.0, 0.0, 0.0,
        0.0, 0.0, 0.0,
        0.0, 0.0, 2.0
    };
    //
    // value of Hessian of g_1
    double check_hes_1[] = {
        0.0, 1.0, 0.0,
        1.0, 0.0, 0.0,
        0.0, 0.0, 0.0
    };
    vector<double> u2(n), y2(m);
    for(size_t j = 0; j < n; j++)
        u2[j] = 0.0;
    // compute diagonal elements of the Hessian
    for(size_t j = 0; j < n; j++)
    {   // first order forward in j-th direction
//.........这里部分代码省略.........

bool change_param(void)
{   bool ok = true;                     // initialize test result

    typedef CppAD::AD<double> a1type;   // for first level of taping
    typedef CppAD::AD<a1type>  a2type;  // for second level of taping

    size_t nu = 3;       // number components in u
    size_t nx = 2;       // number components in x
    size_t ny = 2;       // num components in f(x)
    size_t nJ = ny * nx; // number components in Jacobian of f(x)

    // temporary indices
    size_t j;

    // declare first level of independent variables
    // (Start taping now so can record dependency of a1f on a1p.)
    CPPAD_TESTVECTOR(a1type) a1u(nu);
    for(j = 0; j < nu; j++)
        a1u[j] = 0.;
    CppAD::Independent(a1u);

    // parameter in computation of Jacobian
    a1type a1p = a1u[2];

    // declare second level of independent variables
    CPPAD_TESTVECTOR(a2type) a2x(nx);
    for(j = 0; j < nx; j++)
        a2x[j] = 0.;
    CppAD::Independent(a2x);

    // compute dependent variables at second level
    CPPAD_TESTVECTOR(a2type) a2y(ny);
    a2y[0] = sin( a2x[0] ) * a1p;
    a2y[1] = sin( a2x[1] ) * a1p;

    // declare function object that computes values at the first level
    // (make sure we do not run zero order forward during constructor)
    CppAD::ADFun<a1type> a1f;
    a1f.Dependent(a2x, a2y);

    // compute the Jacobian of a1f at a1u[0], a1u[1]
    CPPAD_TESTVECTOR(a1type) a1x(nx);
    a1x[0] = a1u[0];
    a1x[1] = a1u[1];
    CPPAD_TESTVECTOR(a1type) a1J(nJ);
    a1J = a1f.Jacobian( a1x );

    // declare function object that maps u = (x, p) to Jacobian of f
    // (make sure we do not run zero order forward during constructor)
    CppAD::ADFun<double> g;
    g.Dependent(a1u, a1J);

    // remove extra variables used during the reconding of a1f,
    // but not needed any more.
    g.optimize();

    // compute the Jacobian of f using zero order forward
    // sweep with double values
    CPPAD_TESTVECTOR(double) J(nJ), u(nu);
    for(j = 0; j < nu; j++)
        u[j] = double(j+1);
    J = g.Forward(0, u);

    // accuracy for tests
    double eps = 100. * CppAD::numeric_limits<double>::epsilon();

    // y[0] = sin( x[0] ) * p
    // y[1] = sin( x[1] ) * p
    CPPAD_TESTVECTOR(double) x(nx);
    x[0]      = u[0];
    x[1]      = u[1];
    double p  = u[2];

    // J[0] = partial y[0] w.r.t x[0] = cos( x[0] ) * p
    double check = cos( x[0] ) * p;
    ok   &= fabs( check - J[0] ) <= eps;

    // J[1] = partial y[0] w.r.t x[1] = 0.;
    check = 0.;
    ok   &= fabs( check - J[1] ) <= eps;

    // J[2] = partial y[1] w.r.t. x[0] = 0.
    check = 0.;
    ok   &= fabs( check - J[2] ) <= eps;

    // J[3] = partial y[1] w.r.t x[1] = cos( x[1] ) * p
    check = cos( x[1] ) * p;
    ok   &= fabs( check - J[3] ) <= eps;

    return ok;
}

bool link_poly(
	size_t                     size     ,
	size_t                     repeat   ,
	CppAD::vector<double>     &a        ,  // coefficients of polynomial
	CppAD::vector<double>     &z        ,  // polynomial argument value
	CppAD::vector<double>     &ddp      )  // second derivative w.r.t z
{
	// speed test global option values
	if( global_atomic )
		return false;

	// -----------------------------------------------------
	// setup
	typedef CppAD::AD<double>     ADScalar;
	typedef CppAD::vector<ADScalar> ADVector;

	size_t i;      // temporary index
	size_t m = 1;  // number of dependent variables
	size_t n = 1;  // number of independent variables
	ADVector Z(n); // AD domain space vector
	ADVector P(m); // AD range space vector

	// choose the polynomial coefficients
	CppAD::uniform_01(size, a);

	// AD copy of the polynomial coefficients
	ADVector A(size);
	for(i = 0; i < size; i++)
		A[i] = a[i];

	// forward mode first and second differentials
	CppAD::vector<double> p(1), dp(1), dz(1), ddz(1);
	dz[0]  = 1.;
	ddz[0] = 0.;

	// AD function object
	CppAD::ADFun<double> f;

	// --------------------------------------------------------------------
	if( ! global_onetape ) while(repeat--)
	{
		// choose an argument value
		CppAD::uniform_01(1, z);
		Z[0] = z[0];

		// declare independent variables
		Independent(Z);

		// AD computation of the function value
		P[0] = CppAD::Poly(0, A, Z[0]);

		// create function object f : A -> detA
		f.Dependent(Z, P);

		if( global_optimize )
			f.optimize();

		// skip comparison operators
		f.compare_change_count(0);

		// pre-allocate memory for three forward mode calculations
		f.capacity_order(3);

		// evaluate the polynomial
		p = f.Forward(0, z);

		// evaluate first order Taylor coefficient
		dp = f.Forward(1, dz);

		// second derivative is twice second order Taylor coef
		ddp     = f.Forward(2, ddz);
		ddp[0] *= 2.;
	}
	else
	{
		// choose an argument value
		CppAD::uniform_01(1, z);
		Z[0] = z[0];

		// declare independent variables
		Independent(Z);

		// AD computation of the function value
		P[0] = CppAD::Poly(0, A, Z[0]);

		// create function object f : A -> detA
		f.Dependent(Z, P);

		if( global_optimize )
			f.optimize();

		// skip comparison operators
		f.compare_change_count(0);

		while(repeat--)
		{	// sufficient memory is allocated by second repetition

			// get the next argument value
			CppAD::uniform_01(1, z);

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

/* $$
$head Use Atomic Function$$
$codep */
bool reciprocal(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();
/* $$
$subhead Constructor$$
$codep */
	// --------------------------------------------------------------------
	// Create the atomic reciprocal object
	atomic_reciprocal afun("atomic_reciprocal");
/* $$
$subhead Recording$$
$codep */
	// Create the function f(x)
	//
	// domain space vector
	size_t n  = 1;
	double  x0 = 0.5;
	vector< AD<double> > ax(n);
	ax[0]     = x0;

	// declare independent variables and start tape recording
	CppAD::Independent(ax);

	// range space vector 
	size_t m = 1;
	vector< AD<double> > ay(m);

	// call user function and store reciprocal(x) in au[0] 
	vector< AD<double> > au(m);
	afun(ax, au);        // u = 1 / x

	// now use AD division to invert to invert the operation
	ay[0] = 1.0 / au[0]; // y = 1 / u = x

	// create f: x -> y and stop tape recording
	CppAD::ADFun<double> f;
	f.Dependent (ax, ay);  // f(x) = x
/* $$
$subhead forward$$
$codep */
	// check function value 
	double check = x0;
	ok &= NearEqual( Value(ay[0]) , check,  eps, eps);

	// check zero order forward mode
	size_t p;
	vector<double> x_p(n), y_p(m);
	p      = 0;
	x_p[0] = x0;
	y_p    = f.Forward(p, x_p);
	ok &= NearEqual(y_p[0] , check,  eps, eps);

	// check first order forward mode
	p      = 1;
	x_p[0] = 1;
	y_p    = f.Forward(p, x_p);
	check  = 1.;
	ok &= NearEqual(y_p[0] , check,  eps, eps);

	// check second order forward mode
	p      = 2;
	x_p[0] = 0;
	y_p    = f.Forward(p, x_p);
	check  = 0.;
	ok &= NearEqual(y_p[0] , check,  eps, eps);
/* $$
$subhead reverse$$
$codep */
	// third order reverse mode 
	p     = 3;
	vector<double> w(m), dw(n * p);
	w[0]  = 1.;
	dw    = f.Reverse(p, w);
	check = 1.;
	ok &= NearEqual(dw[0] , check,  eps, eps);
	check = 0.;
	ok &= NearEqual(dw[1] , check,  eps, eps);
	ok &= NearEqual(dw[2] , check,  eps, eps);
/* $$
$subhead for_sparse_jac$$
$codep */
	// forward mode sparstiy pattern
	size_t q = n;
	CppAD::vectorBool r1(n * q), s1(m * q);
	r1[0] = true;          // compute sparsity pattern for x[0]
	//
	afun.option( CppAD::atomic_base<double>::bool_sparsity_enum );
	s1    = f.ForSparseJac(q, r1);
	ok  &= s1[0] == true;  // f[0] depends on x[0]  
	//
	afun.option( CppAD::atomic_base<double>::set_sparsity_enum );
	s1    = f.ForSparseJac(q, r1);
	ok  &= s1[0] == true;  // f[0] depends on x[0]  
/* $$
$subhead rev_sparse_jac$$
//.........这里部分代码省略.........

bool mul_cond_rev(void)
{
	bool ok = true;
	using CppAD::vector;
	using CppAD::NearEqual;
	double eps = 10. * std::numeric_limits<double>::epsilon();
	//
	typedef CppAD::AD<double>   a1double;
	typedef CppAD::AD<a1double> a2double;
	//
	a1double a1zero = 0.0;
	a2double a2zero = a1zero;
	a1double a1one  = 1.0;
	a2double a2one  = a1one;
	//
	// --------------------------------------------------------------------
	// create a1f = f(x)
	size_t n = 1;
	size_t m = 25;
	//
	vector<a2double> a2x(n), a2y(m);
	a2x[0] = a2double( 5.0 );
	Independent(a2x);
	//
	size_t i = 0;
	// variable that is greater than one when x[0] is zero
	// and less than one when x[0] is 1.0 or greater
	a2double a2switch  = a2one / (a2x[0] + a2double(0.5));
	// variable that is infinity when x[0] is zero
	// and a normal number when x[0] is 1.0 or greater
	a2double a2inf_var = a2one / a2x[0];
	// variable that is nan when x[0] is zero
	// and a normal number when x[0] is 1.0 or greater
	a2double a2nan_var = ( a2one / a2inf_var ) / a2x[0];
	// variable that is one when x[0] is zero
	// and less then one when x[0] is 1.0 or greater
	a2double a2one_var = a2one / ( a2one + a2x[0] );
	// div
	a2y[i++]  = CondExpGt(a2x[0], a2zero, a2nan_var, a2zero);
	// abs
	a2y[i++]  = CondExpGt(a2x[0], a2zero, abs( a2y[0] ), a2zero);
	// add
	a2y[i++]  = CondExpGt(a2x[0], a2zero, a2nan_var + a2nan_var, a2zero);
	// acos
	a2y[i++]  = CondExpGt(a2x[0], a2zero, acos(a2switch), a2zero);
	// asin
	a2y[i++]  = CondExpGt(a2x[0], a2zero, asin(a2switch), a2zero);
	// atan
	a2y[i++]  = CondExpGt(a2x[0], a2zero, atan(a2nan_var), a2zero);
	// cos
	a2y[i++]  = CondExpGt(a2x[0], a2zero, cos(a2nan_var), a2zero);
	// cosh
	a2y[i++]  = CondExpGt(a2x[0], a2zero, cosh(a2nan_var), a2zero);
	// exp
	a2y[i++]  = CondExpGt(a2x[0], a2zero, exp(a2nan_var), a2zero);
	// log
	a2y[i++]  = CondExpGt(a2x[0], a2zero, log(a2x[0]), a2zero);
	// mul
	a2y[i++]  = CondExpGt(a2x[0], a2zero, a2x[0] * a2inf_var, a2zero);
	// pow
	a2y[i++]  = CondExpGt(a2x[0], a2zero, pow(a2inf_var, a2x[0]), a2zero);
	// sin
	a2y[i++]  = CondExpGt(a2x[0], a2zero, sin(a2nan_var), a2zero);
	// sinh
	a2y[i++]  = CondExpGt(a2x[0], a2zero, sinh(a2nan_var), a2zero);
	// sqrt
	a2y[i++]  = CondExpGt(a2x[0], a2zero, sqrt(a2x[0]), a2zero);
	// sub
	a2y[i++]  = CondExpGt(a2x[0], a2zero, a2inf_var - a2nan_var, a2zero);
	// tan
	a2y[i++]  = CondExpGt(a2x[0], a2zero, tan(a2nan_var), a2zero);
	// tanh
	a2y[i++]  = CondExpGt(a2x[0], a2zero, tanh(a2nan_var), a2zero);
	// azmul
	a2y[i++]  = CondExpGt(a2x[0], a2zero, azmul(a2x[0], a2inf_var), a2zero);
	//
	// Operations that are C+11 atomic
	//
	// acosh
	a2y[i++]  = CondExpGt(a2x[0], a2zero, acosh( a2x[0] ), a2zero);
	// asinh
	a2y[i++]  = CondExpGt(a2x[0], a2zero, asinh( a2nan_var ), a2zero);
	// atanh
	a2y[i++]  = CondExpGt(a2x[0], a2zero, atanh( a2one_var ), a2zero);
	// erf
	a2y[i++]  = CondExpGt(a2x[0], a2zero, erf( a2nan_var ), a2zero);
	// expm1
	a2y[i++]  = CondExpGt(a2x[0], a2zero, expm1(a2nan_var), a2zero);
	// log1p
	a2y[i++]  = CondExpGt(a2x[0], a2zero, log1p(- a2one_var ), a2zero);
	//
	ok &= i == m;
	CppAD::ADFun<a1double> a1f;
	a1f.Dependent(a2x, a2y);
	// --------------------------------------------------------------------
	// create h = f(x)
	vector<a1double> a1x(n), a1y(m);
	a1x[0] = 5.0;
	//
	Independent(a1x);
//.........这里部分代码省略.........

bool link_ode(
	size_t                     size       ,
	size_t                     repeat     ,
	CppAD::vector<double>      &x         ,
	CppAD::vector<double>      &jacobian
)
{
	// speed test global option values
	if( global_option["atomic"] )
		return false;

	// optimization options: no conditional skips or compare operators
	std::string options="no_compare_op";
	// --------------------------------------------------------------------
	// setup
	assert( x.size() == size );
	assert( jacobian.size() == size * size );

	typedef CppAD::AD<double>       ADScalar;
	typedef CppAD::vector<ADScalar> ADVector;

	size_t j;
	size_t p = 0;              // use ode to calculate function values
	size_t n = size;           // number of independent variables
	size_t m = n;              // number of dependent variables
	ADVector  X(n), Y(m);      // independent and dependent variables
	CppAD::ADFun<double>  f;   // AD function

	// -------------------------------------------------------------
	if( ! global_option["onetape"] ) while(repeat--)
	{	// choose next x value
		uniform_01(n, x);
		for(j = 0; j < n; j++)
			X[j] = x[j];

		// declare the independent variable vector
		Independent(X);

		// evaluate function
		CppAD::ode_evaluate(X, p, Y);

		// create function object f : X -> Y
		f.Dependent(X, Y);

		if( global_option["optimize"] )
			f.optimize(options);

		// skip comparison operators
		f.compare_change_count(0);

		jacobian = f.Jacobian(x);
	}
	else
	{	// an x value
		uniform_01(n, x);
		for(j = 0; j < n; j++)
			X[j] = x[j];

		// declare the independent variable vector
		Independent(X);

		// evaluate function
		CppAD::ode_evaluate(X, p, Y);

		// create function object f : X -> Y
		f.Dependent(X, Y);

		if( global_option["optimize"] )
			f.optimize(options);

		// skip comparison operators
		f.compare_change_count(0);

		while(repeat--)
		{	// get next argument value
			uniform_01(n, x);

			// evaluate jacobian
			jacobian = f.Jacobian(x);
		}
	}
	return true;
}

bool old_mat_mul(void)
{	bool ok = true;
	using CppAD::AD;

	// matrix sizes for this test
	size_t nr_result = 2;
	size_t n_middle  = 2;
	size_t nc_result = 2;
	
	// declare the AD<double> vectors ax and ay and X 
	size_t n = nr_result * n_middle + n_middle * nc_result;
	size_t m = nr_result * nc_result;
	CppAD::vector< AD<double> > X(4), ax(n), ay(m);
	size_t i, j;
	for(j = 0; j < X.size(); j++)
		X[j] = (j + 1);

	// X is the vector of independent variables
	CppAD::Independent(X);
	// left matrix
	ax[0]  = X[0];  // left[0,0]   = x[0] = 1
	ax[1]  = X[1];  // left[0,1]   = x[1] = 2
	ax[2]  = 5.;    // left[1,0]   = 5
	ax[3]  = 6.;    // left[1,1]   = 6
	// right matrix
	ax[4]  = X[2];  // right[0,0]  = x[2] = 3
	ax[5]  = 7.;    // right[0,1]  = 7
	ax[6]  = X[3];  // right[1,0]  = x[3] = 4 
	ax[7]  = 8.;    // right[1,1]  = 8
	/*
	[ x0 , x1 ] * [ x2 , 7 ] = [ x0*x2 + x1*x3 , x0*7 + x1*8 ]
	[ 5  , 6 ]    [ x3 , 8 ]   [ 5*x2  + 6*x3  , 5*7 + 6*8 ]
	*/

	// The call back routines need to know the dimensions of the matrices.
	// Store information about the matrix multiply for this call to mat_mul.
	call_info info;
	info.nr_result = nr_result;
	info.n_middle  = n_middle;
	info.nc_result = nc_result;
	// info.vx gets set by forward during call to mat_mul below
	assert( info.vx.size() == 0 ); 
	size_t id      = info_.size();
	info_.push_back(info);

	// user defined AD<double> version of matrix multiply
	mat_mul(id, ax, ay);
	//----------------------------------------------------------------------
	// check AD<double>  results
	ok &= ay[0] == (1*3 + 2*4); ok &= Variable( ay[0] );
	ok &= ay[1] == (1*7 + 2*8); ok &= Variable( ay[1] );
	ok &= ay[2] == (5*3 + 6*4); ok &= Variable( ay[2] );
	ok &= ay[3] == (5*7 + 6*8); ok &= Parameter( ay[3] );
	//----------------------------------------------------------------------
	// use mat_mul to define a function g : X -> ay
	CppAD::ADFun<double> G;
	G.Dependent(X, ay);
	// g(x) = [ x0*x2 + x1*x3 , x0*7 + x1*8 , 5*x2  + 6*x3  , 5*7 + 6*8 ]^T
	//----------------------------------------------------------------------
	// Test zero order forward mode evaluation of g(x)
	CppAD::vector<double> x( X.size() ), y(m);
	for(j = 0; j <  X.size() ; j++)
		x[j] = j + 2;
	y = G.Forward(0, x);
	ok &= y[0] == x[0] * x[2] + x[1] * x[3];
	ok &= y[1] == x[0] * 7.   + x[1] * 8.;
	ok &= y[2] == 5. * x[2]   + 6. * x[3];
	ok &= y[3] == 5. * 7.     + 6. * 8.;

	//----------------------------------------------------------------------
	// Test first order forward mode evaluation of g'(x) * [1, 2, 3, 4]^T 
	// g'(x) = [ x2, x3, x0, x1 ]
	//         [ 7 ,  8,  0, 0  ]
	//         [ 0 ,  0,  5, 6  ]
	//         [ 0 ,  0,  0, 0  ] 
	CppAD::vector<double> dx( X.size() ), dy(m);
	for(j = 0; j <  X.size() ; j++)
		dx[j] = j + 1;
	dy = G.Forward(1, dx);
	ok &= dy[0] == 1. * x[2] + 2. * x[3] + 3. * x[0] + 4. * x[1];
	ok &= dy[1] == 1. * 7.   + 2. * 8.   + 3. * 0.   + 4. * 0.;
	ok &= dy[2] == 1. * 0.   + 2. * 0.   + 3. * 5.   + 4. * 6.;
	ok &= dy[3] == 1. * 0.   + 2. * 0.   + 3. * 0.   + 4. * 0.;

	//----------------------------------------------------------------------
	// Test second order forward mode 
	// g_0^2 (x) = [ 0, 0, 1, 0 ], g_0^2 (x) * [1] = [3]
	//             [ 0, 0, 0, 1 ]              [2]   [4]
	//             [ 1, 0, 0, 0 ]              [3]   [1]
	//             [ 0, 1, 0, 0 ]              [4]   [2]
	CppAD::vector<double> ddx( X.size() ), ddy(m);
	for(j = 0; j <  X.size() ; j++)
		ddx[j] = 0.;
	ddy = G.Forward(2, ddx);
	// [1, 2, 3, 4] * g_0^2 (x) * [1, 2, 3, 4]^T = 1*3 + 2*4 + 3*1 + 4*2
	ok &= 2. * ddy[0] == 1. * 3. + 2. * 4. + 3. * 1. + 4. * 2.; 
	// for i > 0, [1, 2, 3, 4] * g_i^2 (x) * [1, 2, 3, 4]^T = 0
	ok &= ddy[1] == 0.;
	ok &= ddy[2] == 0.;
	ok &= ddy[3] == 0.;
//.........这里部分代码省略.........

bool link_det_minor(
	size_t                     size     , 
	size_t                     repeat   , 
	CppAD::vector<double>     &matrix   ,
	CppAD::vector<double>     &gradient )
{
	// -----------------------------------------------------
	// setup

	// object for computing determinant
	typedef CppAD::AD<double>       ADScalar; 
	typedef CppAD::vector<ADScalar> ADVector; 
	CppAD::det_by_minor<ADScalar>   Det(size);

	size_t i;               // temporary index
	size_t m = 1;           // number of dependent variables
	size_t n = size * size; // number of independent variables
	ADVector   A(n);        // AD domain space vector
	ADVector   detA(m);     // AD range space vector
	
	// vectors of reverse mode weights 
	CppAD::vector<double> w(1);
	w[0] = 1.;

	// the AD function object
	CppAD::ADFun<double> f;

	static bool printed = false;
	bool print_this_time = (! printed) & (repeat > 1) & (size >= 3);

	extern bool global_retape;
	if( global_retape ) while(repeat--)
	{
		// choose a matrix
		CppAD::uniform_01(n, matrix);
		for( i = 0; i < size * size; i++)
			A[i] = matrix[i];
	
		// declare independent variables
		Independent(A);
	
		// AD computation of the determinant
		detA[0] = Det(A);
	
		// create function object f : A -> detA
		f.Dependent(A, detA);

		extern bool global_optimize;
		if( global_optimize )
		{	size_t before, after;
			before = f.size_var();
			f.optimize();
			if( print_this_time ) 
			{	after = f.size_var();
				std::cout << "cppad_det_minor_optimize_size_" 
				          << int(size) << " = [ " << int(before) 
				          << ", " << int(after) << "]" << std::endl;
				printed         = true;
				print_this_time = false;
			}
		}
	
		// get the next matrix
		CppAD::uniform_01(n, matrix);
	
		// evaluate the determinant at the new matrix value
		f.Forward(0, matrix);
	
		// evaluate and return gradient using reverse mode
		gradient = f.Reverse(1, w);
	}
	else
	{
		// choose a matrix
		CppAD::uniform_01(n, matrix);
		for( i = 0; i < size * size; i++)
			A[i] = matrix[i];
	
		// declare independent variables
		Independent(A);
	
		// AD computation of the determinant
		detA[0] = Det(A);
	
		// create function object f : A -> detA
		CppAD::ADFun<double> f;
		f.Dependent(A, detA);

		extern bool global_optimize;
		if( global_optimize )
		{	size_t before, after;
			before = f.size_var();
			f.optimize();
			if( print_this_time ) 
			{	after = f.size_var();
				std::cout << "optimize: size = " << size
				          << ": size_var() = "
				          << before << "(before) " 
				          << after << "(after) " 
				          << std::endl;
//.........这里部分代码省略.........

bool fun_assign(void)
{	bool ok = true;
	using CppAD::AD;
	using CppAD::NearEqual;
	size_t i, j;

	// ten times machine percision
	double eps = 10. * CppAD::numeric_limits<double>::epsilon();

	// two ADFun<double> objects
	CppAD::ADFun<double> g;

	// domain space vector
	size_t n  = 3;
	CPPAD_TESTVECTOR(AD<double>) x(n);
	for(j = 0; j < n; j++)
		x[j] = AD<double>(j + 2);

	// declare independent variables and start tape recording
	CppAD::Independent(x);

	// range space vector
	size_t m = 2;
	CPPAD_TESTVECTOR(AD<double>) y(m);
	y[0] = x[0] + x[0] * x[1];
	y[1] = x[1] * x[2] + x[2];

	// Store operation sequence, and order zero forward results, in f.
	CppAD::ADFun<double> f(x, y);

	// sparsity pattern for the identity matrix
	CPPAD_TESTVECTOR(std::set<size_t>) r(n);
	for(j = 0; j < n; j++)
		r[j].insert(j);

	// Store forward mode sparsity pattern in f
	f.ForSparseJac(n, r);

	// make a copy in g
	g = f;

	// check values that should be equal
	ok &= ( g.size_order()       == f.size_order() );
	ok &= ( g.size_forward_bool() == f.size_forward_bool() );
	ok &= ( g.size_forward_set()  == f.size_forward_set() );

	// Use zero order Taylor coefficient from f for first order
	// calculation using g.
	CPPAD_TESTVECTOR(double) dx(n), dy(m);
	for(i = 0; i < n; i++)
		dx[i] = 0.;
	dx[1] = 1;
	dy    = g.Forward(1, dx);
	ok &= NearEqual(dy[0], x[0], eps, eps); // partial y[0] w.r.t x[1]
	ok &= NearEqual(dy[1], x[2], eps, eps); // partial y[1] w.r.t x[1]

	// Use forward Jacobian sparsity pattern from f to calculate
	// Hessian sparsity pattern using g.
	CPPAD_TESTVECTOR(std::set<size_t>) s(1), h(n);
	s[0].insert(0); // Compute sparsity pattern for Hessian of y[0]
	h =  f.RevSparseHes(n, s);

	// check sparsity pattern for Hessian of y[0] = x[0] + x[0] * x[1]
	ok  &= ( h[0].find(0) == h[0].end() ); // zero     w.r.t x[0], x[0]
	ok  &= ( h[0].find(1) != h[0].end() ); // non-zero w.r.t x[0], x[1]
	ok  &= ( h[0].find(2) == h[0].end() ); // zero     w.r.t x[0], x[2]

	ok  &= ( h[1].find(0) != h[1].end() ); // non-zero w.r.t x[1], x[0]
	ok  &= ( h[1].find(1) == h[1].end() ); //