parallel.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 2016 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, 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.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii__parallel_h
#define dealii__parallel_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/synchronous_iterator.h>
#include <deal.II/base/thread_management.h>

#include <deal.II/base/std_cxx11/tuple.h>
#include <deal.II/base/std_cxx11/bind.h>
#include <deal.II/base/std_cxx11/function.h>
#include <deal.II/base/std_cxx11/shared_ptr.h>

#include <cstddef>

#ifdef DEAL_II_WITH_THREADS
#  include <tbb/parallel_for.h>
#  include <tbb/parallel_reduce.h>
#  include <tbb/partitioner.h>
#  include <tbb/blocked_range.h>
#endif


//TODO[WB]: allow calling functions to pass along a tbb::affinity_partitioner object to ensure that subsequent calls use the same cache lines

DEAL_II_NAMESPACE_OPEN

namespace parallel
{
  namespace internal
  {
    template <typename Number>
    struct EnableOpenMPSimdFor
    {
      static const bool value = true;
    };

#ifdef __INTEL_COMPILER
    // Disable long double SIMD instructions on ICC. This is to work around a bug
    // that generates wrong code at least up to intel 15 (see
    // tests/lac/vector-vector, tests/lac/intel-15-bug, and the discussion at
    // https://github.com/dealii/dealii/issues/598).
    template <>
    struct EnableOpenMPSimdFor<long double>
    {
      static const bool value = false;
    };
#endif



    template <typename F>
    struct Body
    {
      Body (const F &f)
        :
        f (f)
      {}

      template <typename Range>
      void
      operator () (const Range &range) const
      {
        for (typename Range::const_iterator p=range.begin();
             p != range.end(); ++p)
          apply (f, p.iterators);
      }

    private:
      const F f;

      template <typename I1, typename I2>
      static
      void
      apply (const F &f,
             const std_cxx11::tuple<I1,I2> &p)
      {
        *std_cxx11::get<1>(p) = f (*std_cxx11::get<0>(p));
      }

      template <typename I1, typename I2, typename I3>
      static
      void
      apply (const F &f,
             const std_cxx11::tuple<I1,I2,I3> &p)
      {
        *std_cxx11::get<2>(p) = f (*std_cxx11::get<0>(p),
                                   *std_cxx11::get<1>(p));
      }

      template <typename I1, typename I2,
                typename I3, typename I4>
      static
      void
      apply (const F &f,
             const std_cxx11::tuple<I1,I2,I3,I4> &p)
      {
        *std_cxx11::get<3>(p) = f (*std_cxx11::get<0>(p),
                                   *std_cxx11::get<1>(p),
                                   *std_cxx11::get<2>(p));
      }
    };


    template <typename F>
    Body<F> make_body(const F &f)
    {
      return Body<F>(f);
    }
  }

  template <typename InputIterator,
            typename OutputIterator,
            typename Predicate>
  void transform (const InputIterator &begin_in,
                  const InputIterator &end_in,
                  OutputIterator       out,
                  Predicate           &predicate,
                  const unsigned int   grainsize)
  {
#ifndef DEAL_II_WITH_THREADS
    // make sure we don't get compiler
    // warnings about unused arguments
    (void) grainsize;

    for (OutputIterator in = begin_in; in != end_in;)
      *out++ = predicate (*in++);
#else
    typedef std_cxx11::tuple<InputIterator,OutputIterator> Iterators;
    typedef SynchronousIterators<Iterators> SyncIterators;
    Iterators x_begin (begin_in, out);
    Iterators x_end (end_in, OutputIterator());
    tbb::parallel_for (tbb::blocked_range<SyncIterators>(x_begin,
                                                         x_end,
                                                         grainsize),
                       internal::make_body (predicate),
                       tbb::auto_partitioner());
#endif
  }



  template <typename InputIterator1,
            typename InputIterator2,
            typename OutputIterator,
            typename Predicate>
  void transform (const InputIterator1 &begin_in1,
                  const InputIterator1 &end_in1,
                  InputIterator2        in2,
                  OutputIterator        out,
                  Predicate            &predicate,
                  const unsigned int    grainsize)
  {
#ifndef DEAL_II_WITH_THREADS
    // make sure we don't get compiler
    // warnings about unused arguments
    (void) grainsize;

    for (OutputIterator in1 = begin_in1; in1 != end_in1;)
      *out++ = predicate (*in1++, *in2++);
#else
    typedef
    std_cxx11::tuple<InputIterator1,InputIterator2,OutputIterator>
    Iterators;
    typedef SynchronousIterators<Iterators> SyncIterators;
    Iterators x_begin (begin_in1, in2, out);
    Iterators x_end (end_in1, InputIterator2(), OutputIterator());
    tbb::parallel_for (tbb::blocked_range<SyncIterators>(x_begin,
                                                         x_end,
                                                         grainsize),
                       internal::make_body (predicate),
                       tbb::auto_partitioner());
#endif
  }



  template <typename InputIterator1,
            typename InputIterator2,
            typename InputIterator3,
            typename OutputIterator,
            typename Predicate>
  void transform (const InputIterator1 &begin_in1,
                  const InputIterator1 &end_in1,
                  InputIterator2        in2,
                  InputIterator3        in3,
                  OutputIterator        out,
                  Predicate            &predicate,
                  const unsigned int    grainsize)
  {
#ifndef DEAL_II_WITH_THREADS
    // make sure we don't get compiler
    // warnings about unused arguments
    (void) grainsize;

    for (OutputIterator in1 = begin_in1; in1 != end_in1;)
      *out++ = predicate (*in1++, *in2++, *in3++);
#else
    typedef
    std_cxx11::tuple<InputIterator1,InputIterator2,InputIterator3,OutputIterator>
    Iterators;
    typedef SynchronousIterators<Iterators> SyncIterators;
    Iterators x_begin (begin_in1, in2, in3, out);
    Iterators x_end (end_in1, InputIterator2(),
                     InputIterator3(), OutputIterator());
    tbb::parallel_for (tbb::blocked_range<SyncIterators>(x_begin,
                                                         x_end,
                                                         grainsize),
                       internal::make_body (predicate),
                       tbb::auto_partitioner());
#endif
  }


  namespace internal
  {
#ifdef DEAL_II_WITH_THREADS

    template <typename RangeType, typename Function>
    void apply_to_subranges (const tbb::blocked_range<RangeType> &range,
                             const Function  &f)
    {
      f (range.begin(), range.end());
    }
#endif
  }


  template <typename RangeType, typename Function>
  void apply_to_subranges (const RangeType                          &begin,
                           const typename identity<RangeType>::type &end,
                           const Function                           &f,
                           const unsigned int                        grainsize)
  {
#ifndef DEAL_II_WITH_THREADS
    // make sure we don't get compiler
    // warnings about unused arguments
    (void) grainsize;

#  ifndef DEAL_II_BIND_NO_CONST_OP_PARENTHESES
    f (begin, end);
#  else
    // work around a problem with MS VC++ where there is no const
    // operator() in 'Function' if 'Function' is the result of std::bind
    Function ff = f;
    ff (begin, end);
#  endif
#else
    tbb::parallel_for (tbb::blocked_range<RangeType>
                       (begin, end, grainsize),
                       std_cxx11::bind (&internal::apply_to_subranges<RangeType,Function>,
                                        std_cxx11::_1,
                                        std_cxx11::cref(f)),
                       tbb::auto_partitioner());
#endif
  }



  struct ParallelForInteger
  {
    virtual ~ParallelForInteger ();

    void apply_parallel (const std::size_t begin,
                         const std::size_t end,
                         const std::size_t minimum_parallel_grain_size) const;

    virtual void apply_to_subrange (const std::size_t,
                                    const std::size_t) const = 0;
  };



  namespace internal
  {
#ifdef DEAL_II_WITH_THREADS

    template <typename ResultType,
              typename Function>
    struct ReductionOnSubranges
    {
      ResultType result;

      template <typename Reductor>
      ReductionOnSubranges (const Function &f,
                            const Reductor &reductor,
                            const ResultType neutral_element = ResultType())
        :
        result (neutral_element),
        f (f),
        neutral_element (neutral_element),
        reductor (reductor)
      {}

      ReductionOnSubranges (const ReductionOnSubranges &r,
                            tbb::split)
        :
        result (r.neutral_element),
        f (r.f),
        neutral_element (r.neutral_element),
        reductor (r.reductor)
      {}

      void join (const ReductionOnSubranges &r)
      {
        result = reductor(result, r.result);
      }

      template <typename RangeType>
      void operator () (const tbb::blocked_range<RangeType> &range)
      {
        result = reductor(result,
                          f (range.begin(), range.end()));
      }

    private:
      const Function f;

      const ResultType neutral_element;

      const std_cxx11::function<ResultType (ResultType, ResultType)> reductor;
    };
#endif
  }


  template <typename ResultType, typename RangeType, typename Function>
  ResultType accumulate_from_subranges (const Function &f,
                                        const RangeType                          &begin,
                                        const typename identity<RangeType>::type &end,
                                        const unsigned int grainsize)
  {
#ifndef DEAL_II_WITH_THREADS
    // make sure we don't get compiler
    // warnings about unused arguments
    (void) grainsize;

#  ifndef DEAL_II_BIND_NO_CONST_OP_PARENTHESES
    return f (begin, end);
#  else
    // work around a problem with MS VC++ where there is no const
    // operator() in 'Function' if 'Function' is the result of std::bind
    Function ff = f;
    return ff (begin, end);
#  endif
#else
    internal::ReductionOnSubranges<ResultType,Function>
    reductor (f, std::plus<ResultType>(), 0);
    tbb::parallel_reduce (tbb::blocked_range<RangeType>(begin, end, grainsize),
                          reductor,
                          tbb::auto_partitioner());
    return reductor.result;
#endif
  }


// --------------------- for loop affinity partitioner -----------------------

  namespace internal
  {
    class TBBPartitioner
    {
    public:
      TBBPartitioner()
#ifdef DEAL_II_WITH_THREADS
        :
        my_partitioner(new tbb::affinity_partitioner()),
        in_use(false)
#endif
      {}

#ifdef DEAL_II_WITH_THREADS

      ~TBBPartitioner()
      {
        Assert(in_use == false,
               ExcInternalError("A vector partitioner goes out of scope, but "
                                "it appears to be still in use."));
      }

      std_cxx11::shared_ptr<tbb::affinity_partitioner>
      acquire_one_partitioner()
      {
        ::Threads::Mutex::ScopedLock lock(mutex);
        if (in_use)
          return std_cxx11::shared_ptr<tbb::affinity_partitioner>(new tbb::affinity_partitioner());

        in_use = true;
        return my_partitioner;
      }

      void release_one_partitioner(std_cxx11::shared_ptr<tbb::affinity_partitioner> &p)
      {
        if (p.get() == my_partitioner.get())
          {
            ::Threads::Mutex::ScopedLock lock(mutex);
            in_use = false;
          }
      }

    private:
      std_cxx11::shared_ptr<tbb::affinity_partitioner> my_partitioner;

      bool in_use;

      ::Threads::Mutex mutex;
#endif
    };
  }
}


namespace internal
{
  namespace Vector
  {
    extern unsigned int minimum_parallel_grain_size;
  }


  namespace SparseMatrix
  {
    extern unsigned int minimum_parallel_grain_size;
  }

} // end of namespace internal


/* --------------------------- inline functions ------------------------- */

namespace parallel
{

#ifdef DEAL_II_WITH_THREADS

  namespace internal
  {
    struct ParallelForWrapper
    {
      ParallelForWrapper (const parallel::ParallelForInteger &worker)
        :
        worker_ (worker)
      {}

      void operator() (const tbb::blocked_range<std::size_t> &range) const
      {
        worker_.apply_to_subrange (range.begin(), range.end());
      }

      const parallel::ParallelForInteger &worker_;
    };
  }

#endif


  inline
  ParallelForInteger::~ParallelForInteger ()
  {}


  inline
  void
  ParallelForInteger::apply_parallel (const std::size_t begin,
                                      const std::size_t end,
                                      const std::size_t minimum_parallel_grain_size) const
  {
#ifndef DEAL_II_WITH_THREADS
    // make sure we don't get compiler
    // warnings about unused arguments
    (void) minimum_parallel_grain_size;

    apply_to_subrange (begin, end);
#else
    internal::ParallelForWrapper worker(*this);
    tbb::parallel_for (tbb::blocked_range<std::size_t>
                       (begin, end, minimum_parallel_grain_size),
                       worker,
                       tbb::auto_partitioner());
#endif
  }

} // end of namespace parallel

DEAL_II_NAMESPACE_CLOSE

#endif