C++ vector::push_back方法代码示例

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

// -----------------------------------------------------------------------
// get the result of the work
bool multi_newton_combine(CppAD::vector<double>& xout)
{	// number of threads in the calculation
	size_t num_threads  = std::max(num_threads_, size_t(1));

	// remove duplicates and points that are not solutions
	xout.resize(0);
	bool   ok = true;
	size_t thread_num;

	// initialize as more that sub_lenght_ / 2 from any possible solution
	double xlast = - sub_length_;
	for(thread_num = 0; thread_num < num_threads; thread_num++)
	{	vector<double>& x = work_all_[thread_num]->x;

		size_t i;
		for(i = 0; i < x.size(); i++)
		{	// check for case where this point is lower limit for this
			// thread and upper limit for previous thread
			if( fabs(x[i] - xlast) >= sub_length_ )
			{	xout.push_back( x[i] );
				xlast = x[i];
			}
			else
			{	double fcur, flast, df;
				fun_(x[i],   fcur, df);
				fun_(xlast, flast, df);
				if( fabs(fcur) < fabs(flast) )
				{	xout[ xout.size() - 1] = x[i];
					xlast                  = x[i];
				}
			}
		}
		ok &= work_all_[thread_num]->ok;
	}

	// go down so free memory for other threads before memory for master
	thread_num = num_threads;
	while(thread_num--)
	{
# if USE_THREAD_ALLOC_FOR_WORK_ALL
		// call the destructor for CppAD::vector destructor
		work_all_[thread_num]->x.~vector<double>();
		// delete the raw memory allocation
		void* v_ptr = static_cast<void*>( work_all_[thread_num] );
		thread_alloc::return_memory( v_ptr );
# else
		delete work_all_[thread_num];
# endif
		// Note that xout corresponds to memroy that is inuse by master
		// (so we can only chech have freed all their memory).
		if( thread_num > 0 )
		{	// check that there is no longer any memory inuse by this thread
			ok &= thread_alloc::inuse(thread_num) == 0;
			// return all memory being held for future use by this thread
			thread_alloc::free_available(thread_num);
		}
	}
	// now we are done with the work_all_ vector so free its memory
	// (becasue it is a static variable)
	work_all_.clear();

	return ok;
}

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

	size_t m = Range();
	size_t n = Domain();
	VectorBase jac(m * n);

	CPPAD_ASSERT_KNOWN(
		size_t(x.size()) == n,
		"SparseJacobian: size of x not equal domain size for f."
	);
	CheckSimpleVector<Base, VectorBase>();

	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 < m; i++)
		for(j = 0; j < n; j++)
			jac[i * n + j] = zero;

	sparse_jacobian_work work;
	CppAD::vector<size_t> row;
	CppAD::vector<size_t> col;
	if( n <= m )
	{
		// need an internal copy of sparsity pattern
		Pattern_type s_transpose;
		bool transpose = true;
		sparsity_user2internal(s_transpose, p, m, n, transpose);

		k = 0;
		for(j = 0; j < n; j++)
		{	s_transpose.begin(j);
			i = s_transpose.next_element();
			while( i != s_transpose.end() )
			{	row.push_back(i);
				col.push_back(j);
				k++;
				i = s_transpose.next_element();
			}
		} 
		size_t K = k;
		VectorBase J(K);
	
		// now we have folded this into the following case
		SparseJacobianFor(x, s_transpose, row, col, J, work);

		// now set the non-zero return values
		for(k = 0; k < K; k++)
			jac[ row[k] * n + col[k] ] = J[k];
	}
	else
	{
		// need an internal copy of sparsity pattern
		Pattern_type s;
		bool transpose = false;
		sparsity_user2internal(s, p, m, n, transpose);

		k = 0;
		for(i = 0; i < m; 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 J(K);

		// now we have folded this into the following case
		SparseJacobianRev(x, s, row, col, J, work);

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

	return jac;
}