void ode_evaluate(
	CppAD::vector<Float> &x  , 
	size_t m                 , 
	CppAD::vector<Float> &fm )
{
	typedef CppAD::vector<Float> Vector;

	size_t n = x.size();
	size_t ell;
	CPPAD_ASSERT_KNOWN( m == 0 || m == 1,
		"ode_evaluate: m is not zero or one"
	);
	CPPAD_ASSERT_KNOWN( 
		((m==0) & (fm.size()==n)) || ((m==1) & (fm.size()==n*n)),
		"ode_evaluate: the size of fm is not correct"
	);
	if( m == 0 )
		ell = n;
	else	ell = n + n * n;

	// set up the case we are integrating
	Float  ti   = 0.;
	Float  tf   = 1.;
	Float  smin = 1e-5;
	Float smax  = 1.;
	Float scur  = 1.;
	Float erel  = 0.;
	vector<Float> yi(ell), eabs(ell);
	size_t i, j;
	for(i = 0; i < ell; i++)
	{	eabs[i] = 1e-10;
		if( i < n )
			yi[i] = 1.;
		else	yi[i]  = 0.;
	}

	// return values
	Vector yf(ell), ef(ell), maxabs(ell);
	size_t nstep;

	// construct ode method for taking one step
	ode_evaluate_method<Float> method(m, x);

	// solve differential equation
	yf = OdeErrControl(method, 
		ti, tf, yi, smin, smax, scur, eabs, erel, ef, maxabs, nstep);

	if( m == 0 )
	{	for(i = 0; i < n; i++)
			fm[i] = yf[i];
	}
	else
	{	for(i = 0; i < n; i++)
			for(j = 0; j < n; j++)
				fm[i * n + j] = yf[n + i * n + j];
	}
	return;
}

void put_check_for_nan(const CppAD::vector<Base>& vec, std::string& file_name)
{
	size_t char_size       = sizeof(Base) * vec.size();
	const char* char_ptr   = reinterpret_cast<const char*>( vec.data() );
# if CPPAD_HAS_MKSTEMP
	char pattern[] = "/tmp/fileXXXXXX";
	int fd = mkstemp(pattern);
	file_name = pattern;
	write(fd, char_ptr, char_size);
	close(fd);
# else
# if CPPAD_HAS_TMPNAM_S
		std::vector<char> name(L_tmpnam_s);
		if( tmpnam_s( name.data(), L_tmpnam_s ) != 0 )
		{	CPPAD_ASSERT_KNOWN(
				false,
				"Cannot create a temporary file name"
			);
		}
		file_name = name.data();
# else
		file_name = tmpnam( CPPAD_NULL );
# endif
	std::fstream file_out(file_name.c_str(), std::ios::out|std::ios::binary );
	file_out.write(char_ptr, char_size);
	file_out.close();
# endif
	return;
}

	void sparse_jac_fun(
		size_t                       m    ,
		size_t                       n    ,
		const FloatVector&           x    ,
		const CppAD::vector<size_t>& row  , 
		const CppAD::vector<size_t>& col  , 
		size_t                       p    ,
		FloatVector&                 fp   )
	{
		// check numeric type specifications
		CheckNumericType<Float>();
		// check value of p
		CPPAD_ASSERT_KNOWN(
			p == 0 || p == 1,
			"sparse_jac_fun: p != 0 and p != 1"
		);
		size_t K = row.size();
		CPPAD_ASSERT_KNOWN(
			K >= m,
			"sparse_jac_fun: row.size() < m"
		);
		size_t i, j, k;

		if( p == 0 )
			for(i = 0; i < m; i++)
				fp[i] = Float(0);

		Float t;
		for(k = 0; k < K; k++)
		{	i    = row[k];
			j    = col[k];
			t    = exp( x[j] * x[j] / 2.0 );	
			switch(p)
			{
				case 0:
				fp[i] += t;
				break;

				case 1:
				fp[k] = t * x[j];
				break;
			}
		}
	}

	void sparse_hes_fun(
		size_t                       n    ,
		const FloatVector&           x    ,
		const CppAD::vector<size_t>& row  , 
		const CppAD::vector<size_t>& col  , 
		size_t                       p    ,
		FloatVector&                fp    )
	{
		// check numeric type specifications
		CheckNumericType<Float>();

		// check value of p
		CPPAD_ASSERT_KNOWN(
			p < 3,
			"sparse_hes_fun: p > 2"
		);

		size_t i, j, k;
		size_t size = 1;
		for(k = 0; k < p; k++)
			size *= n;
		for(k = 0; k < size; k++)
			fp[k] = Float(0);

		size_t K = row.size();
		Float t;
		Float dt_i;
		Float dt_j;
		for(k = 0; k < K; k++)
		{	i    = row[k];
			j    = col[k];
			t    = exp( x[i] * x[j] );	
			dt_i = t * x[j];
			dt_j = t * x[i];
			switch(p)
			{
				case 0:
				fp[0] += t;
				break;

				case 1:
				fp[i] += dt_i;
				fp[j] += dt_j;
				break;

				case 2:
				fp[i * n + i] += dt_i * x[j];
				fp[i * n + j] += t + dt_j * x[j];
				//
				fp[j * n + i] += t + dt_i * x[i];
				fp[j * n + j] += dt_j * x[i];
				break;
			}
		}
			
	}

Vector ADFun<Base>::ForOne(const Vector &x, size_t j)
{	size_t j1;

	size_t n = Domain();
	size_t m = Range();

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

	CPPAD_ASSERT_KNOWN(
		x.size() == n,
		"ForOne: Length of x not equal domain dimension for f"
	);
	CPPAD_ASSERT_KNOWN(
		j < n,
		"ForOne: the index j is not less than domain dimension for f"
	);

	// point at which we are evaluating the second partials
	Forward(0, x);

	// direction in which are are taking the derivative
	Vector dx(n);
	for(j1 = 0; j1 < n; j1++)
		dx[j1] = Base(0);
	dx[j] = Base(1);

	// dimension the return value
	Vector dy(m);

	// compute the return value
	dy = Forward(1, dx);

	return dy;
}

void parallel_ad(void)
{   CPPAD_ASSERT_KNOWN(
        ! thread_alloc::in_parallel() ,
        "parallel_ad must be called before entering parallel execution mode."
    );
    CPPAD_ASSERT_KNOWN(
        AD<Base>::tape_ptr() == CPPAD_NULL ,
        "parallel_ad cannot be called while a tape recording is in progress"
    );

    // ensure statics in following functions are initialized
    elapsed_seconds();
    ErrorHandler::Current();
    local::NumArg(local::BeginOp);
    local::NumRes(local::BeginOp);
    local::one_element_std_set<size_t>();
    local::two_element_std_set<size_t>();

    // the sparse_pack class has member functions with static data
    local::sparse_pack sp;
    sp.resize(1, 1);       // so can call add_element
    sp.add_element(0, 0);  // has static data
    sp.clear(0);           // has static data
    sp.is_element(0, 0);   // has static data
    local::sparse_pack::const_iterator itr(sp, 0); // has static data
    ++itr;                                  // has static data

    // statics that depend on the value of Base
    AD<Base>::tape_id_ptr(0);
    AD<Base>::tape_handle(0);
    discrete<Base>::List();
    CheckSimpleVector< Base, CppAD::vector<Base> >();
    CheckSimpleVector< AD<Base>, CppAD::vector< AD<Base> > >();

}

void  AD<Base>::tape_delete(size_t id_old)
{
	size_t thread = id_old % CPPAD_MAX_NUM_THREADS;
	CPPAD_ASSERT_KNOWN(
		thread == thread_alloc::thread_num(),
		"AD tape recording must stop in same thread as it started in."
	);
	size_t        *id   = id_handle(thread);
	ADTape<Base> **tape = tape_handle(thread);

	CPPAD_ASSERT_UNKNOWN( *id   == id_old     );
	CPPAD_ASSERT_UNKNOWN( *tape != CPPAD_NULL );

	// increase the id for this thread in a way such that 
	// thread = id % CPPAD_MAX_NUM_THREADS
	CPPAD_ASSERT_KNOWN(
	size_t(*id) + CPPAD_MAX_NUM_THREADS <= 
		std::numeric_limits<CPPAD_TAPE_ID_TYPE>::max(),
	"To many different tapes given the type used for CPPAD_TAPE_ID_TYPE"
	);
	*id  += CPPAD_MAX_NUM_THREADS;

	// delete the old tape for this thread
	delete ( *tape );
	*tape = CPPAD_NULL;

	return;
}

Vector ADFun<Base>::RevOne(const Vector  &x, size_t i)
{	size_t i1;

	size_t n = Domain();
	size_t m = Range();

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

	CPPAD_ASSERT_KNOWN(
		x.size() == n,
		"RevOne: Length of x not equal domain dimension for f"
	);
	CPPAD_ASSERT_KNOWN(
		i < m,
		"RevOne: the index i is not less than range dimension for f"
	);

	// point at which we are evaluating the derivative
	Forward(0, x);

	// component which are are taking the derivative of
	Vector w(m);
	for(i1 = 0; i1 < m; i1++)
		w[i1] = 0.;
	w[i] = Base(1);

	// dimension the return value
	Vector dw(n);

	// compute the return value
	dw = Reverse(1, w);

	return dw;
}

void color_general_colpack(
	const VectorSet&        pattern ,
	const VectorSize&       row     ,
	const VectorSize&       col     ,
	CppAD::vector<size_t>&  color   )
{	size_t i, j, k;
	size_t m = pattern.n_set();
	size_t n = pattern.end();

	// Determine number of non-zero entries in each row
	CppAD::vector<size_t> n_nonzero(m);
	size_t n_nonzero_total = 0;
	for(i = 0; i < m; i++)
	{	n_nonzero[i] = 0;
		typename VectorSet::const_iterator pattern_itr(pattern, i);
		j = *pattern_itr;
		while( j != pattern.end() )
		{	n_nonzero[i]++;
			j = *(++pattern_itr);
		}
		n_nonzero_total += n_nonzero[i];
	}

	// Allocate memory and fill in Adolc sparsity pattern
	CppAD::vector<unsigned int*> adolc_pattern(m);
	CppAD::vector<unsigned int>  adolc_memory(m + n_nonzero_total);
	size_t i_memory = 0;
	for(i = 0; i < m; i++)
	{	adolc_pattern[i]    = adolc_memory.data() + i_memory;
		CPPAD_ASSERT_KNOWN(
			std::numeric_limits<unsigned int>::max() >= n_nonzero[i],
			"Matrix is too large for colpack"
		);
		adolc_pattern[i][0] = static_cast<unsigned int>( n_nonzero[i] );
		typename VectorSet::const_iterator pattern_itr(pattern, i);
		j = *pattern_itr;
		k = 1;
		while(j != pattern.end() )
		{
			CPPAD_ASSERT_KNOWN(
				std::numeric_limits<unsigned int>::max() >= j,
				"Matrix is too large for colpack"
			);
			adolc_pattern[i][k++] = static_cast<unsigned int>( j );
			j = *(++pattern_itr);
		}
		CPPAD_ASSERT_UNKNOWN( k == 1 + n_nonzero[i] );
		i_memory += k;
	}
	CPPAD_ASSERT_UNKNOWN( i_memory == m + n_nonzero_total );

	// Must use an external routine for this part of the calculation because
	// ColPack/ColPackHeaders.h has as 'using namespace std' at global level.
	cppad_colpack_general(color, m, n, adolc_pattern);

	return;
}

Vector ADFun<Base>::Hessian(const Vector &x, const Vector &w)
{	size_t j;
	size_t k;

	size_t n = Domain();

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

	CPPAD_ASSERT_KNOWN(
		size_t(x.size()) == n,
		"Hessian: length of x not equal domain dimension for f"
	);
	CPPAD_ASSERT_KNOWN(
		size_t(w.size()) == Range(),
		"Hessian: length of w not equal range dimension for f"
	);

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

	// define the return value
	Vector hes(n * n);

	// direction vector for calls to forward
	Vector u(n);
	for(j = 0; j < n; j++)
		u[j] = Base(0);


	// location for return values from Reverse
	Vector ddw(n * 2);

	// loop over forward directions
	for(j = 0; j < n; j++)
	{	// evaluate partials of entire function w.r.t. j-th coordinate
		u[j] = Base(1);
		Forward(1, u);
		u[j] = Base(0);

		// evaluate derivative of partial corresponding to F_i
		ddw = Reverse(2, w);

		// return desired components
		for(k = 0; k < n; k++)
			hes[k * n + j] = ddw[k * 2 + 1];
	}

	return hes;
}

void ADFun<Base>::RevSparseHesCase(
	bool              set_type         ,
	bool              transpose        ,  
	size_t            q                ,  
	const VectorSet&  s                ,
	VectorSet&        h                )
{	size_t n = Domain(); 	
	h.resize(q * n );

	CPPAD_ASSERT_KNOWN( 
		for_jac_sparse_pack_.n_set() > 0,
		"RevSparseHes: previous stored call to ForSparseJac did not "
		"use bool for the elements of r."
	);
	CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == 0 );
	CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == total_num_var_ );
	
	// use sparse_pack for the calculation
	CppAD::RevSparseHesBool( 
		transpose                ,
		q                        ,
		s                        ,
		h                        ,
		total_num_var_           ,
		dep_taddr_               ,
		ind_taddr_               ,
		play_                    ,
		for_jac_sparse_pack_ 
	);
}

void fun_record(
	cppad_ipopt_fg_info*                              fg_info ,
	size_t                                            k       ,
	const SizeVector&                                 p       ,
	const SizeVector&                                 q       ,
	size_t                                            n       ,
	const NumVector&                                  x       ,
	const SizeVector&                                 J       ,
	CppAD::vector< CppAD::ADFun<Ipopt::Number> >&     r_fun   )
{	size_t j;

	// extract u from x
	ADVector u(q[k]);
	for(j = 0; j < q[k]; j++)
	{	// when NDEBUG is not defined, this error should be caught
		// during the cppad_ipopt_nlp constructor.
		CPPAD_ASSERT_UNKNOWN( J[j] < n );
		u[j] = x[ J[j] ];
	}

	// start the recording
	CppAD::Independent(u);

	// record the evaulation of r_k (u)
	ADVector r_k = fg_info->eval_r(k, u);
	CPPAD_ASSERT_KNOWN( r_k.size() == p[k] ,
	"cppad_ipopt_nlp: eval_r return value size not equal to p[k]."
	);

	// stop the recording and store operation sequence in
	r_fun[k].Dependent(u, r_k);
}

 /// not a number
 static Float quiet_NaN(void)
 {   CPPAD_ASSERT_KNOWN(
         false,
         "numeric_limits<Float>::quiet_NaN() is not specialized for this Float"
     );
     return Float(0);
 }

	/*!
 	Link from forward mode sweep to users routine.

	\param index
	index for this function in the list of all user_atomic objects

	\param id
	extra information vector that is just passed through by CppAD,
	and possibly used by user's routines.

	\param k
	order for this forward mode calculation.

	\param n
	domain space size for this calcualtion.

	\param m
	range space size for this calculation.

	\param tx
	Taylor coefficients corresponding to \c x for this calculation.

	\param ty
	Taylor coefficient corresponding to \c y for this calculation

	See the forward mode in user's documentation for user_atomic 
 	*/
	static void forward(
		size_t                index , 
		size_t                   id ,
		size_t                    k ,
		size_t                    n , 
		size_t                    m , 
		const vector<Base>&      tx ,
		vector<Base>&            ty )
	{	
# ifdef _OPENMP
		vector<bool> empty(0);
# else
		static vector<bool> empty(0);
# endif
		
		CPPAD_ASSERT_UNKNOWN( tx.size() >= n * k );
		CPPAD_ASSERT_UNKNOWN( ty.size() >= m * k );

		CPPAD_ASSERT_UNKNOWN(index < List().size() );
		user_atomic* op = List()[index];

		bool ok = op->f_(id, k, n, m, empty, empty, tx, ty);
		if( ! ok )
		{	std::stringstream ss;
			ss << k;
			std::string msg = "user_atomic: ";
			msg = msg + op->name_ + ": ok returned false from " + ss.str()
			    + " order forward mode calculation";
			CPPAD_ASSERT_KNOWN(false, msg.c_str());
		}
	}

inline void forward_sqrt_op_dir(
	size_t q           ,
	size_t r           ,
	size_t i_z         ,
	size_t i_x         ,
	size_t cap_order   ,
	Base*  taylor      )
{
	// check assumptions
	CPPAD_ASSERT_UNKNOWN( NumArg(SqrtOp) == 1 );
	CPPAD_ASSERT_UNKNOWN( NumRes(SqrtOp) == 1 );
	CPPAD_ASSERT_UNKNOWN( 0 < q );
	CPPAD_ASSERT_UNKNOWN( q < cap_order );

	// Taylor coefficients corresponding to argument and result
	size_t num_taylor_per_var = (cap_order-1) * r + 1;
	Base* z = taylor + i_z * num_taylor_per_var;
	Base* x = taylor + i_x * num_taylor_per_var;
	CPPAD_ASSERT_KNOWN(
		x[0] != Base(0),
		"Forward: attempt to take derivatve of square root of zero"
	)

	size_t m = (q-1) * r + 1;
	for(size_t ell = 0; ell < r; ell++)
	{	z[m+ell] = Base(0);
		for(size_t k = 1; k < q; k++)
			z[m+ell] -= Base(k) * z[(k-1)*r+1+ell] * z[(q-k-1)*r+1+ell];
		z[m+ell] /= Base(q);
		z[m+ell] += x[m+ell] / Base(2);
		z[m+ell] /= z[0];
	}
}

	void operator()(const ADVector& ax, ADVector& ay, size_t id = 0)
	{	CPPAD_ASSERT_KNOWN(
			id == 0,
			"checkpoint: id is non-zero in afun(ax, ay, id)"
		);
		this->atomic_base<Base>::operator()(ax, ay, id);
	}

	/*!
 	Link from reverse mode sweep to users routine.

	\param index
	index in the list of all user_atomic objects
	corresponding to this function.


	\param id
	extra information vector that is just passed through by CppAD,
	and possibly used by user's routines.

	\param k
	order for this forward mode calculation.

	\param n
	domain space size for this calcualtion.

	\param m
	range space size for this calculation.

	\param tx
	Taylor coefficients corresponding to \c x for this calculation.

	\param ty
	Taylor coefficient corresponding to \c y for this calculation

	\param px
	Partials w.r.t. the \c x Taylor coefficients.

	\param py
	Partials w.r.t. the \c y Taylor coefficients.

	See reverse mode documentation for user_atomic 
 	*/
	static void reverse(
		size_t               index , 
		size_t                  id ,
		size_t                   k ,
		size_t                   n , 
		size_t                   m , 
		const vector<Base>&     tx ,
		const vector<Base>&     ty ,
		vector<Base>&           px ,
		const vector<Base>&     py )
	{
		CPPAD_ASSERT_UNKNOWN(index < List().size() );
		CPPAD_ASSERT_UNKNOWN( tx.size() >= n * k );
		CPPAD_ASSERT_UNKNOWN( px.size() >= n * k );
		CPPAD_ASSERT_UNKNOWN( ty.size() >= m * k );
		CPPAD_ASSERT_UNKNOWN( py.size() >= m * k );
		user_atomic* op = List()[index];

		bool ok = op->r_(id, k, n, m, tx, ty, px, py);
		if( ! ok )
		{	std::stringstream ss;
			ss << k;
			std::string msg = "user_atomic: ";
			msg = op->name_ + ": ok returned false from " + ss.str() 
			    + " order reverse mode calculation";
			CPPAD_ASSERT_KNOWN(false, msg.c_str());
		}
	}

 /// machine epsilon
 static Float epsilon(void)
 {   CPPAD_ASSERT_KNOWN(
         false,
         "numeric_limits<Float>::epsilon() is not specialized for this Float"
     );
     return Float(0);
 }

	/*!
 	Link from reverse Hessian sparsity sweep to users routine.

	\param index
	index in the list of all user_atomic objects
	corresponding to this function.

	\param id
	extra information vector that is just passed through by CppAD,
	and possibly used by user's routines.

	\param n
	domain space size for this calcualtion.

	\param m
	range space size for this calculation.

	\param q
	is the column dimension for the sparsity partterns.

	\param r
	is the forward Jacobian sparsity pattern w.r.t the argument vector x

	\param s
	is the reverse Jacobian sparsity pattern w.r.t the result vector y.

	\param t
	is the reverse Jacobian sparsity pattern w.r.t the argument vector x.

	\param u
	is the Hessian sparsity pattern w.r.t the result vector y.

	\param v
	is the Hessian sparsity pattern w.r.t the argument vector x.
	*/
	static void rev_hes_sparse(
		size_t                            index ,
		size_t                               id ,
		size_t                                n , 
		size_t                                m , 
		size_t                                q ,
		vector< std::set<size_t> >&           r ,
		const vector<bool>&                   s ,
		vector<bool>&                         t ,
		const vector< std::set<size_t> >&     u ,
		vector< std::set<size_t> >&           v )
	{
		CPPAD_ASSERT_UNKNOWN(index < List().size() );
		CPPAD_ASSERT_UNKNOWN( r.size() >= n );
		CPPAD_ASSERT_UNKNOWN( s.size() >= m );
		CPPAD_ASSERT_UNKNOWN( t.size() >= n );
		CPPAD_ASSERT_UNKNOWN( u.size() >= m );
		CPPAD_ASSERT_UNKNOWN( v.size() >= n );
		user_atomic* op = List()[index];

		bool ok = op->rhs_(id, n, m, q, r, s, t, u, v);
		if( ! ok )
		{	std::string msg = "user_atomic: ";
			msg = msg + op->name_ 
			    + ": ok returned false from rev_jac_sparse calculation";
			CPPAD_ASSERT_KNOWN(false, msg.c_str());
		}
	}

void ADFun<Base>::RevSparseHesCase(
	const std::set<size_t>&   set_type         ,
	bool                      transpose        ,  
	size_t                    q                ,  
	const VectorSet&          s                ,
	VectorSet&                h                )
{	size_t n = Domain();
	if( transpose )
		h.resize(n);
	else	h.resize(q);

	CPPAD_ASSERT_KNOWN( 
		for_jac_sparse_set_.n_set() > 0,
		"RevSparseHes: previous stored call to ForSparseJac did not "
		"use std::set<size_t> for the elements of r."
	);
	CPPAD_ASSERT_UNKNOWN( for_jac_sparse_pack_.n_set() == 0 );
	CPPAD_ASSERT_UNKNOWN( for_jac_sparse_set_.n_set() == num_var_tape_  );
	
	// use sparse_pack for the calculation
	CppAD::RevSparseHesSet( 
		transpose                ,
		q                        ,
		s                        ,
		h                        ,
		num_var_tape_            ,
		dep_taddr_               ,
		ind_taddr_               ,
		play_                    ,
		for_jac_sparse_set_ 
	);
}