libstdc++-v3/include/parallel/partial_sum.h

   1 // -*- C++ -*-
   2
   3 // Copyright (C) 2007, 2008, 2009 Free Software Foundation, Inc.
   4 //
   5 // This file is part of the GNU ISO C++ Library.  This library is free
   6 // software; you can redistribute it and/or modify it under the terms
   7 // of the GNU General Public License as published by the Free Software
   8 // Foundation; either version 3, or (at your option) any later
   9 // version.
  10
  11 // This library is distributed in the hope that it will be useful, but
  12 // WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // General Public License for more details.
  15
  16 // Under Section 7 of GPL version 3, you are granted additional
  17 // permissions described in the GCC Runtime Library Exception, version
  18 // 3.1, as published by the Free Software Foundation.
  19
  20 // You should have received a copy of the GNU General Public License and
  21 // a copy of the GCC Runtime Library Exception along with this program;
  22 // see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  23 // <http://www.gnu.org/licenses/>.
  24
  25 /** @file parallel/partial_sum.h
  26  *  @brief Parallel implementation of std::partial_sum(), i.e. prefix
  27 *  sums.
  28  *  This file is a GNU parallel extension to the Standard C++ Library.
  29  */
  30
  31 // Written by Johannes Singler.
  32
  33 #ifndef _GLIBCXX_PARALLEL_PARTIAL_SUM_H
  34 #define _GLIBCXX_PARALLEL_PARTIAL_SUM_H 1
  35
  36 #include <omp.h>
  37 #include <new>
  38 #include <bits/stl_algobase.h>
  39 #include <parallel/parallel.h>
  40 #include <parallel/numericfwd.h>
  41
  42 namespace __gnu_parallel
  43 {
  44   // Problem: there is no 0-element given.
  45
  46   /** @brief Base case prefix sum routine.
  47    *  @param __begin Begin iterator of input sequence.
  48    *  @param __end End iterator of input sequence.
  49    *  @param __result Begin iterator of output sequence.
  50    *  @param __bin_op Associative binary function.
  51    *  @param __value Start value. Must be passed since the neutral
  52    *  element is unknown in general.
  53    *  @return End iterator of output sequence. */
  54   template<typename _IIter,
  55            typename _OutputIterator,
  56            typename _BinaryOperation>
  57     _OutputIterator
  58     __parallel_partial_sum_basecase(_IIter __begin, _IIter __end,
  59                                     _OutputIterator __result,
  60                                     _BinaryOperation __bin_op,
  61       typename std::iterator_traits <_IIter>::value_type __value)
  62     {
  63       if (__begin == __end)
  64         return __result;
  65
  66       while (__begin != __end)
  67         {
  68           __value = __bin_op(__value, *__begin);
  69           *__result = __value;
  70           ++__result;
  71           ++__begin;
  72         }
  73       return __result;
  74     }
  75
  76   /** @brief Parallel partial sum implementation, two-phase approach,
  77       no recursion.
  78       *  @param __begin Begin iterator of input sequence.
  79       *  @param __end End iterator of input sequence.
  80       *  @param __result Begin iterator of output sequence.
  81       *  @param __bin_op Associative binary function.
  82       *  @param __n Length of sequence.
  83       *  @param __num_threads Number of threads to use.
  84       *  @return End iterator of output sequence.
  85       */
  86   template<typename _IIter,
  87            typename _OutputIterator,
  88            typename _BinaryOperation>
  89     _OutputIterator
  90     __parallel_partial_sum_linear(_IIter __begin, _IIter __end,
  91                                   _OutputIterator __result,
  92                                   _BinaryOperation __bin_op,
  93       typename std::iterator_traits<_IIter>::difference_type __n)
  94     {
  95       typedef std::iterator_traits<_IIter> _TraitsType;
  96       typedef typename _TraitsType::value_type _ValueType;
  97       typedef typename _TraitsType::difference_type _DifferenceType;
  98
  99       if (__begin == __end)
 100         return __result;
 101
 102       _ThreadIndex __num_threads =
 103         std::min<_DifferenceType>(__get_max_threads(), __n - 1);
 104
 105       if (__num_threads < 2)
 106         {
 107           *__result = *__begin;
 108           return __parallel_partial_sum_basecase(__begin + 1, __end,
 109                                                  __result + 1, __bin_op,
 110                                                  *__begin);
 111         }
 112
 113       _DifferenceType* __borders;
 114       _ValueType* __sums;
 115
 116       const _Settings& __s = _Settings::get();
 117
 118 #     pragma omp parallel num_threads(__num_threads)
 119       {
 120 #       pragma omp single
 121         {
 122           __num_threads = omp_get_num_threads();
 123
 124           __borders = new _DifferenceType[__num_threads + 2];
 125
 126           if (__s.partial_sum_dilation == 1.0f)
 127             equally_split(__n, __num_threads + 1, __borders);
 128           else
 129             {
 130               _DifferenceType __first_part_length =
 131                   std::max<_DifferenceType>(1,
 132                     __n / (1.0f + __s.partial_sum_dilation * __num_threads));
 133               _DifferenceType __chunk_length =
 134                   (__n - __first_part_length) / __num_threads;
 135               _DifferenceType __borderstart =
 136                   __n - __num_threads * __chunk_length;
 137               __borders[0] = 0;
 138               for (_ThreadIndex __i = 1; __i < (__num_threads + 1); ++__i)
 139                 {
 140                   __borders[__i] = __borderstart;
 141                   __borderstart += __chunk_length;
 142                 }
 143               __borders[__num_threads + 1] = __n;
 144             }
 145
 146           __sums = static_cast<_ValueType*>(::operator new(sizeof(_ValueType)
 147                                                            * __num_threads));
 148           _OutputIterator __target_end;
 149         } //single
 150
 151         _ThreadIndex __iam = omp_get_thread_num();
 152         if (__iam == 0)
 153           {
 154             *__result = *__begin;
 155             __parallel_partial_sum_basecase(__begin + 1,
 156                                             __begin + __borders[1],
 157                                             __result + 1,
 158                                             __bin_op, *__begin);
 159             ::new(&(__sums[__iam])) _ValueType(*(__result + __borders[1] - 1));
 160           }
 161         else
 162           {
 163             ::new(&(__sums[__iam]))
 164               _ValueType(__gnu_parallel::accumulate(
 165                                          __begin + __borders[__iam] + 1,
 166                                          __begin + __borders[__iam + 1],
 167                                          *(__begin + __borders[__iam]),
 168                                          __bin_op,
 169                                          __gnu_parallel::sequential_tag()));
 170           }
 171
 172 #       pragma omp barrier
 173
 174 #       pragma omp single
 175         __parallel_partial_sum_basecase(__sums + 1, __sums + __num_threads,
 176                                         __sums + 1, __bin_op, __sums[0]);
 177
 178 #       pragma omp barrier
 179
 180         // Still same team.
 181         __parallel_partial_sum_basecase(__begin + __borders[__iam + 1],
 182                                         __begin + __borders[__iam + 2],
 183                                         __result + __borders[__iam + 1],
 184                                         __bin_op, __sums[__iam]);
 185       } //parallel
 186
 187       ::operator delete(__sums);
 188       delete[] __borders;
 189
 190       return __result + __n;
 191     }
 192
 193   /** @brief Parallel partial sum front-__end.
 194    *  @param __begin Begin iterator of input sequence.
 195    *  @param __end End iterator of input sequence.
 196    *  @param __result Begin iterator of output sequence.
 197    *  @param __bin_op Associative binary function.
 198    *  @return End iterator of output sequence. */
 199   template<typename _IIter,
 200            typename _OutputIterator,
 201            typename _BinaryOperation>
 202     _OutputIterator
 203     __parallel_partial_sum(_IIter __begin, _IIter __end,
 204                            _OutputIterator __result, _BinaryOperation __bin_op)
 205     {
 206       _GLIBCXX_CALL(__begin - __end)
 207
 208       typedef std::iterator_traits<_IIter> _TraitsType;
 209       typedef typename _TraitsType::value_type _ValueType;
 210       typedef typename _TraitsType::difference_type _DifferenceType;
 211
 212       _DifferenceType __n = __end - __begin;
 213
 214       switch (_Settings::get().partial_sum_algorithm)
 215         {
 216         case LINEAR:
 217           // Need an initial offset.
 218           return __parallel_partial_sum_linear(__begin, __end, __result,
 219                                                __bin_op, __n);
 220         default:
 221           // Partial_sum algorithm not implemented.
 222           _GLIBCXX_PARALLEL_ASSERT(0);
 223           return __result + __n;
 224         }
 225     }
 226 }
 227
 228 #endif /* _GLIBCXX_PARALLEL_PARTIAL_SUM_H */