omptl_numeric.h

omptl使用openmp多核库优化过的STL库

// Copyright (C) 2006 Fokko Beekhof
// Email contact: [e-mail]

// The OMPTL library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA

#ifndef OMPTL_NUMERIC_H
#define OMPTL_NUMERIC_H

#include <utility>
#include <numeric>

#include <omptl/omptl_algorithm>
#include <omptl/omptl_tools.h>

namespace omptl
{

// accumulate()

template <class InputIterator, class T, class BinaryFunction>
T _par_accumulate(InputIterator first, InputIterator last, T init, T par_init,
		BinaryFunction binary_op, const unsigned P)
{
	assert(P > 0u);
	if (_linear_serial_is_faster(first, last, P))
		return ::std::accumulate(first, last, init);
	assert(2*(int)P <= std::distance(first, last));

	::std::pair<InputIterator, InputIterator> partitions[P];
	::omptl::_partition_range(first, last, partitions, P);

	T results[P];
	int t;
	#pragma omp parallel for // default(shared), private(t)
	for (t = 0; t < int(P); ++t)
		results[t] = ::std::accumulate( partitions[t].first,
						partitions[t].second, par_init);

	return ::std::accumulate(results, results + P, init, binary_op);
}

template <class InputIterator, class T>
T _accumulate(InputIterator first, InputIterator last, T init,
	::std::plus<typename ::std::iterator_traits<InputIterator>::value_type>
		binary_op, const unsigned P)
{
	assert(P > 0u);
	return _par_accumulate(first, last, init, T(0), binary_op, P);
}

template <class InputIterator, class T>
T _accumulate(InputIterator first, InputIterator last, T init,
   ::std::multiplies<typename ::std::iterator_traits<InputIterator>::value_type>
	binary_op, const unsigned P)
{
	assert(P > 0u);
	return _par_accumulate(first, last, init, T(1), binary_op, P);
}

template <class InputIterator, class T, class BinaryFunction>
T _accumulate(InputIterator first, InputIterator last, T init,
             BinaryFunction binary_op, const unsigned P)
{
	assert(P > 0u);
	return ::std::accumulate(first, last, init, binary_op);
}

template <class InputIterator, class T, class BinaryFunction>
T accumulate(InputIterator first, InputIterator last, T init,
             BinaryFunction binary_op, const unsigned P)
{
	assert(P > 0u);
	return _accumulate(first, last, init, binary_op, P);
}

template <class InputIterator, class T>
T accumulate(InputIterator first, InputIterator last, T init, const unsigned P)
{
	assert(P > 0u);
	typedef typename std::iterator_traits<InputIterator>::value_type VT;

	return _accumulate(first, last, init, std::plus<VT>(), P);
}

// transform_accumulate
template <class InputIterator,class T, class UnaryFunction,class BinaryFunction>
T _ser_transform_accumulate(InputIterator first, InputIterator last, T init,
		UnaryFunction unary_op, BinaryFunction binary_op)
{
	// serial version
	while (first != last)
	{
		init = binary_op(unary_op(*first), init);
		++first;
	}

	return init;
}

template <class InputIterator, class UnaryFunction,class BinaryFunction>
typename BinaryFunction::result_type
 _par_accumulate(InputIterator first, InputIterator last,
		typename BinaryFunction::result_type init,
		typename BinaryFunction::result_type par_init,
		UnaryFunction unary_op, BinaryFunction binary_op, const unsigned P)
{
	assert(P > 0u);
	if (_linear_serial_is_faster(first, last, P))
		return _ser_transform_accumulate(first, last, init,
						 unary_op, binary_op);
	assert(2*(int)P <= std::distance(first, last));

	::std::pair<InputIterator, InputIterator> partitions[P];
	::omptl::_partition_range(first, last, partitions, P);

	typename BinaryFunction::result_type results[P];
	int t;
	#pragma omp parallel for // default(shared), private(t)
	for (t = 0; t < int(P); ++t)
		results[t] = _ser_transform_accumulate(partitions[t].first,
					partitions[t].second, par_init,
					unary_op, binary_op);

	return ::std::accumulate(results, results + P, init, binary_op);
}

template <class InputIterator, class UnaryFunction>
typename UnaryFunction::result_type
 _transform_accumulate(InputIterator first, InputIterator last,
		typename UnaryFunction::result_type init,
		UnaryFunction unary_op,
		::std::plus<typename UnaryFunction::result_type> binary_op,
		unsigned P)
{
	return _par_transform_accumulate(first, last, init,
				typename UnaryFunction::result_type(0),
					 unary_op, binary_op, P);
}

template <class InputIterator, class UnaryFunction>
typename UnaryFunction::result_type
 _transform_accumulate(InputIterator first, InputIterator last,
		typename UnaryFunction::result_type init,
		UnaryFunction unary_op,
		::std::multiplies<typename UnaryFunction::result_type>binary_op,
		unsigned P)
{
	return _par_transform_accumulate(first, last, init,
				typename UnaryFunction::result_type(1),
					 unary_op, binary_op, P);
}

template <class InputIterator, class UnaryFunction, class BinaryFunction>
typename BinaryFunction::result_type
 _transform_accumulate(InputIterator first, InputIterator last,
		typename BinaryFunction::result_type init,
		UnaryFunction unary_op, BinaryFunction binary_op, const unsigned P)
{
	return _ser_transform_accumulate(first, last, init, unary_op,binary_op);
}

template <class InputIterator, class UnaryFunction, class BinaryFunction>
typename BinaryFunction::result_type
transform_accumulate(InputIterator first, InputIterator last,
		typename BinaryFunction::result_type init,
		UnaryFunction unary_op, BinaryFunction binary_op,
		const unsigned P)
{
	return _transform_accumulate(first, last, init, unary_op, binary_op, P);
}

template <class InputIterator, class UnaryFunction>
typename UnaryFunction::result_type
transform_accumulate(InputIterator first, InputIterator last,
		typename UnaryFunction::result_type init,
		UnaryFunction unary_op,	unsigned P)
{
	typedef typename UnaryFunction::result_type RT;
	return ::omptl::transform_accumulate(first, last, init, unary_op,
					     ::std::plus<RT>(), P);
}

template <class InputIterator, class OutputIterator, class BinaryFunction>
OutputIterator adjacent_difference(InputIterator first, InputIterator last,
                                   OutputIterator result,
				   BinaryFunction binary_op, const unsigned P)
{
	return ::std::adjacent_difference(first, last, result, binary_op);
}

template <class InputIterator, class OutputIterator>
OutputIterator adjacent_difference(InputIterator first, InputIterator last,
                                   OutputIterator result, const unsigned P)
{
// 	::std::minus<typename ::std::iterator_traits<InputIterator>::value_type>()
	return ::std::adjacent_difference(first, last, result);
}

// TODO
template <class InputIterator1, class InputIterator2, class T,
          class BinaryFunction1, class BinaryFunction2>
T inner_product(InputIterator1 first1, InputIterator1 last1,
                InputIterator2 first2, T init,
		BinaryFunction1 binary_op1, BinaryFunction2 binary_op2,
		const unsigned P)
{
	return ::std::inner_product(first1, last1, first2, init,
				    binary_op1, binary_op2);
}

template <class InputIterator1, class InputIterator2, class T>
T inner_product(InputIterator1 first1, InputIterator1 last1,
                InputIterator2 first2, T init, const unsigned P)
{
// ::std::plus<typename ::std::iterator_traits<InputIterator1>::value_type>(),
// ::std::multiplies<typename::std::iterator_traits<InputIterator2>::value_type>(),
	return ::std::inner_product(first1, last1, first2, init);
}

template <class InputIterator, class OutputIterator, class BinaryOperation>
OutputIterator partial_sum(InputIterator first, InputIterator last,
                           OutputIterator result, BinaryOperation binary_op,
			   const unsigned P)
{
	return ::std::partial_sum(first, last, result, binary_op);
}

template <class InputIterator, class OutputIterator>
OutputIterator partial_sum(InputIterator first, InputIterator last,
                           OutputIterator result, const unsigned P)
{
// ::std::plus<typename ::std::iterator_traits<InputIterator>::value_type>(),
	return ::std::partial_sum(first, last, result);
}

template <class InputIterator1, class InputIterator2, class T>
typename ::std::iterator_traits<InputIterator1>::value_type
L1(InputIterator1 first1, InputIterator1 last1,
   InputIterator2 first2, const unsigned P)
{
	typedef typename ::std::iterator_traits<InputIterator1>::value_type VT;
	return ::omptl::inner_product(	first1, last1, first2,
					std::minus<VT>(), std::plus<VT>());
}

template <typename T>
struct _MinusSq
{
	T operator()(const T &lhs, const T &rhs) const
	{
		const T d = lhs - rhs;
		return d*d;
	}
};

template <class InputIterator1, class InputIterator2, class T>
typename ::std::iterator_traits<InputIterator1>::value_type
L2(InputIterator1 first1, InputIterator1 last1,
   InputIterator2 first2, const unsigned P)
{
	typedef typename ::std::iterator_traits<InputIterator1>::value_type VT;
	return ::omptl::inner_product(	first1, last1, first2,
					_MinusSq<VT>(), std::plus<VT>());
}

} // namespace omptl

#endif /* OMPTL_NUMERIC_H */