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

   1 // -*- C++ -*-
   2
   3 // Copyright (C) 2007 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 2, 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 // You should have received a copy of the GNU General Public License
  17 // along with this library; see the file COPYING.  If not, write to
  18 // the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
  19 // MA 02111-1307, USA.
  20
  21 // As a special exception, you may use this file as part of a free
  22 // software library without restriction.  Specifically, if other files
  23 // instantiate templates or use macros or inline functions from this
  24 // file, or you compile this file and link it with other files to
  25 // produce an executable, this file does not by itself cause the
  26 // resulting executable to be covered by the GNU General Public
  27 // License.  This exception does not however invalidate any other
  28 // reasons why the executable file might be covered by the GNU General
  29 // Public License.
  30
  31 /** @file parallel/partial_sum.h
  32  *  @brief Parallel implementation of std::partial_sum(), i. e. prefix
  33  *  sums.
  34  *  This file is a GNU parallel extension to the Standard C++ Library.
  35  */
  36
  37 // Written by Johannes Singler.
  38
  39 #ifndef _GLIBCXX_PARALLEL_PARTIAL_SUM_H
  40 #define _GLIBCXX_PARALLEL_PARTIAL_SUM_H 1
  41
  42
  43 #include <omp.h>
  44 #include <bits/stl_algobase.h>
  45 #include <parallel/parallel.h>
  46 #include <parallel/numericfwd.h>
  47
  48 namespace __gnu_parallel
  49 {
  50   // Problem: there is no 0-element given.
  51
  52   /** @brief Base case prefix sum routine.
  53    *  @param begin Begin iterator of input sequence.
  54    *  @param end End iterator of input sequence.
  55    *  @param result Begin iterator of output sequence.
  56    *  @param bin_op Associative binary function.
  57    *  @param value Start value. Must be passed since the neutral
  58    *  element is unknown in general.
  59    *  @return End iterator of output sequence. */
  60   template<typename InputIterator, typename OutputIterator, typename BinaryOperation>
  61   inline OutputIterator
  62   parallel_partial_sum_basecase(InputIterator begin, InputIterator end,
  63                                 OutputIterator result, BinaryOperation bin_op,
  64                                 typename std::iterator_traits<InputIterator>::value_type value)
  65   {
  66     if (begin == end)
  67       return result;
  68
  69     while (begin != end)
  70       {
  71         value = bin_op(value, *begin);
  72         *result = value;
  73         result++;
  74         begin++;
  75       }
  76     return result;
  77   }
  78
  79   /** @brief Parallel partial sum implmenetation, two-phase approach,
  80       no recursion.
  81       *  @param begin Begin iterator of input sequence.
  82       *  @param end End iterator of input sequence.
  83       *  @param result Begin iterator of output sequence.
  84       *  @param bin_op Associative binary function.
  85       *  @param n Length of sequence.
  86       *  @param num_threads Number of threads to use.
  87       *  @return End iterator of output sequence.
  88       */
  89   template<typename InputIterator, typename OutputIterator, typename BinaryOperation>
  90   OutputIterator
  91   parallel_partial_sum_linear(InputIterator begin, InputIterator end,
  92                               OutputIterator result, BinaryOperation bin_op,
  93                               typename std::iterator_traits<InputIterator>::difference_type n, int num_threads)
  94   {
  95     typedef std::iterator_traits<InputIterator> traits_type;
  96     typedef typename traits_type::value_type value_type;
  97     typedef typename traits_type::difference_type difference_type;
  98
  99     if (num_threads > (n - 1))
 100       num_threads = static_cast<thread_index_t>(n - 1);
 101     if (num_threads < 2)
 102       {
 103         *result = *begin;
 104         return parallel_partial_sum_basecase(begin + 1, end, result + 1, bin_op, *begin);
 105       }
 106
 107     difference_type borders[num_threads + 2];
 108
 109     if (Settings::partial_sum_dilatation == 1.0f)
 110       equally_split(n, num_threads + 1, borders);
 111     else
 112       {
 113         difference_type chunk_length = (int)((double)n / ((double)num_threads + Settings::partial_sum_dilatation)), borderstart = n - num_threads * chunk_length;
 114         borders[0] = 0;
 115         for (int i = 1; i < (num_threads + 1); i++)
 116           {
 117             borders[i] = borderstart;
 118             borderstart += chunk_length;
 119           }
 120         borders[num_threads + 1] = n;
 121       }
 122
 123     value_type* sums = new value_type[num_threads];
 124     OutputIterator target_end;
 125
 126 #pragma omp parallel num_threads(num_threads)
 127     {
 128       int id = omp_get_thread_num();
 129       if (id == 0)
 130         {
 131           *result = *begin;
 132           parallel_partial_sum_basecase(begin + 1, begin + borders[1], result + 1, bin_op, *begin);
 133           sums[0] = *(result + borders[1] - 1);
 134         }
 135       else
 136         {
 137           sums[id] = std::accumulate(begin + borders[id] + 1, begin + borders[id + 1], *(begin + borders[id]), bin_op, __gnu_parallel::sequential_tag());
 138         }
 139
 140 #pragma omp barrier
 141
 142 #pragma omp single
 143       parallel_partial_sum_basecase(sums + 1, sums + num_threads, sums + 1, bin_op, sums[0]);
 144
 145 #pragma omp barrier
 146
 147       // Still same team.
 148       parallel_partial_sum_basecase(begin + borders[id + 1], begin + borders[id + 2], result + borders[id + 1], bin_op, sums[id]);
 149     }
 150
 151     delete[] sums;
 152
 153     return result + n;
 154   }
 155
 156   /** @brief Parallel partial sum front-end.
 157    *  @param begin Begin iterator of input sequence.
 158    *  @param end End iterator of input sequence.
 159    *  @param result Begin iterator of output sequence.
 160    *  @param bin_op Associative binary function.
 161    *  @return End iterator of output sequence. */
 162   template<typename InputIterator, typename OutputIterator, typename BinaryOperation>
 163   OutputIterator
 164   parallel_partial_sum(InputIterator begin, InputIterator end,
 165                        OutputIterator result, BinaryOperation bin_op)
 166   {
 167     _GLIBCXX_CALL(begin - end);
 168
 169     typedef std::iterator_traits<InputIterator> traits_type;
 170     typedef typename traits_type::value_type value_type;
 171     typedef typename traits_type::difference_type difference_type;
 172
 173     difference_type n = end - begin;
 174
 175     int num_threads = get_max_threads();
 176
 177     switch (Settings::partial_sum_algorithm)
 178       {
 179       case Settings::LINEAR:
 180         // Need an initial offset.
 181         return parallel_partial_sum_linear(begin, end, result, bin_op,
 182                                            n, num_threads);
 183       default:
 184         // Partial_sum algorithm not implemented.
 185         _GLIBCXX_PARALLEL_ASSERT(0);
 186         return end;
 187       }
 188   }
 189 }
 190
 191 #endif