parallel.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2009 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_parallel_h
#define dealii_parallel_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/multithread_info.h>
#include <deal.II/base/mutex.h>
#include <deal.II/base/std_cxx20/type_traits.h>
#include <deal.II/base/synchronous_iterator.h>
#include <deal.II/base/template_constraints.h>
 
#include <cstddef>
#include <functional>
#include <memory>
#include <tuple>
 
#ifdef DEAL_II_WITH_TASKFLOW
 
#  include <taskflow/algorithm/for_each.hpp>
#  include <taskflow/taskflow.hpp>
#endif
 
#ifdef DEAL_II_WITH_TBB
#  include <tbb/blocked_range.h>
#  include <tbb/parallel_for.h>
#  include <tbb/parallel_reduce.h>
#  include <tbb/partitioner.h>
#else
#  include <boost/range/iterator_range.hpp>
#endif
 
#ifdef DEAL_II_HAVE_CXX20
#  include <concepts>
#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
 
#ifdef DEAL_II_WITH_TASKFLOW
    template <typename Iterator, typename Functor>
    void
    taskflow_parallel_for(Iterator           x_begin,
                          Iterator           x_end,
                          const Functor     &functor,
                          const unsigned int grainsize)
    {
      tf::Executor &executor = MultithreadInfo::get_taskflow_executor();
      tf::Taskflow  taskflow;
 
      // TODO: We have several choices for the Partitioner and we should spend
      // some time benchmarking them:
      // 1. StaticPartitioner(grainsize): all work items have grainsize number
      // of items
      // 2. GuidedPartitioner(grainsize):
      // "The size of a partition is proportional to the number of unassigned
      // iterations divided by the number of workers, and the size will
      // gradually decrease to the given chunk size."
      // 3. GuidedPartitioner(0): The default.
      taskflow.for_each(
        x_begin,
        x_end,
        [&](const auto &item) { functor(item); },
        tf::StaticPartitioner(grainsize));
 
      // Schedule the work. If we are within a task, we are required to use
      // corun() without wait() to avoid a potential deadlock as described in
      // https://taskflow.github.io/taskflow/ExecuteTaskflow.html#ExecuteATaskflowFromAnInternalWorker:
      if (executor.this_worker_id() != -1)
        executor.corun(taskflow);
      else
        executor.run(taskflow).wait();
    }
#endif
 
#ifdef DEAL_II_WITH_TBB
    template <typename Iterator, typename Functor>
    void
    parallel_for(Iterator           x_begin,
                 Iterator           x_end,
                 const Functor     &functor,
                 const unsigned int grainsize)
    {
      tbb::parallel_for(tbb::blocked_range<Iterator>(x_begin, x_end, grainsize),
                        functor,
                        tbb::auto_partitioner());
    }
 
 

    template <typename Iterator, typename Functor>
    void
    parallel_for(Iterator                                          x_begin,
                 Iterator                                          x_end,
                 const Functor                                    &functor,
                 const unsigned int                                grainsize,
                 const std::shared_ptr<tbb::affinity_partitioner> &partitioner)
    {
      tbb::parallel_for(tbb::blocked_range<Iterator>(x_begin, x_end, grainsize),
                        functor,
                        *partitioner);
    }
 
#else

    template <typename Iterator, typename Functor>
    void
    parallel_for(Iterator       x_begin,
                 Iterator       x_end,
                 const Functor &functor,
                 const unsigned int)
    {
      functor(boost::iterator_range<Iterator>(x_begin, x_end));
    }
 
#endif
  } // namespace internal

  template <typename InputIterator, typename OutputIterator, typename Function>
  DEAL_II_CXX20_REQUIRES(
    (std::invocable<Function, decltype(*std::declval<InputIterator>())> &&
     std::assignable_from<
       decltype(*std::declval<OutputIterator>()),
       std::invoke_result_t<Function,
                            decltype(*std::declval<InputIterator>())>>))
  DEAL_II_DEPRECATED void transform(const InputIterator &begin_in,
                                    const InputIterator &end_in,
                                    OutputIterator       out,
                                    const Function      &function,
                                    const unsigned int   grainsize)
  {
#ifdef DEAL_II_WITH_TASKFLOW
    using Iterators     = std::tuple<InputIterator, OutputIterator>;
    using SyncIterators = SynchronousIterators<Iterators>;
    Iterators x_begin(begin_in, out);
    Iterators x_end(end_in, OutputIterator());
 
    internal::taskflow_parallel_for(
      SyncIterators(x_begin),
      SyncIterators(x_end),
      [function](const auto &it) {
        *std::get<1>(it) = function(*std::get<0>(it));
      },
      grainsize);
 
#elif defined(DEAL_II_WITH_TBB)
    using Iterators     = std::tuple<InputIterator, OutputIterator>;
    using SyncIterators = SynchronousIterators<Iterators>;
    Iterators x_begin(begin_in, out);
    Iterators x_end(end_in, OutputIterator());
    internal::parallel_for(
      SyncIterators(x_begin),
      SyncIterators(x_end),
      [function](const auto &range) {
        for (const auto &p : range)
          *std::get<1>(p) = function(*std::get<0>(p));
      },
      grainsize);
#else
    // make sure we don't get compiler
    // warnings about unused arguments
    (void)grainsize;
 
    for (InputIterator in = begin_in; in != end_in;)
      *out++ = function(*in++);
#endif
  }
 
 

  template <typename InputIterator1,
            typename InputIterator2,
            typename OutputIterator,
            typename Function>
  DEAL_II_CXX20_REQUIRES(
    (std::invocable<Function,
                    decltype(*std::declval<InputIterator1>()),
                    decltype(*std::declval<InputIterator2>())> &&
     std::assignable_from<
       decltype(*std::declval<OutputIterator>()),
       std::invoke_result_t<Function,
                            decltype(*std::declval<InputIterator1>()),
                            decltype(*std::declval<InputIterator2>())>>))
  DEAL_II_DEPRECATED void transform(const InputIterator1 &begin_in1,
                                    const InputIterator1 &end_in1,
                                    InputIterator2        in2,
                                    OutputIterator        out,
                                    const Function       &function,
                                    const unsigned int    grainsize)
  {
#ifdef DEAL_II_WITH_TASKFLOW
    using Iterators =
      std::tuple<InputIterator1, InputIterator2, OutputIterator>;
    using SyncIterators = SynchronousIterators<Iterators>;
    Iterators x_begin(begin_in1, in2, out);
    Iterators x_end(end_in1, InputIterator2(), OutputIterator());
 
    internal::taskflow_parallel_for(
      SyncIterators(x_begin),
      SyncIterators(x_end),
      [function](const auto &it) {
        *std::get<2>(it) = function(*std::get<0>(it), *std::get<1>(it));
      },
      grainsize);
 
#elif defined(DEAL_II_WITH_TBB)
    using Iterators =
      std::tuple<InputIterator1, InputIterator2, OutputIterator>;
    using SyncIterators = SynchronousIterators<Iterators>;
    Iterators x_begin(begin_in1, in2, out);
    Iterators x_end(end_in1, InputIterator2(), OutputIterator());
    internal::parallel_for(
      SyncIterators(x_begin),
      SyncIterators(x_end),
      [function](const auto &range) {
        for (const auto &p : range)
          *std::get<2>(p) = function(*std::get<0>(p), *std::get<1>(p));
      },
      grainsize);
 
#else
    // make sure we don't get compiler
    // warnings about unused arguments
    (void)grainsize;
 
    for (InputIterator1 in1 = begin_in1; in1 != end_in1;)
      *out++ = function(*in1++, *in2++);
#endif
  }
 
 

  template <typename InputIterator1,
            typename InputIterator2,
            typename InputIterator3,
            typename OutputIterator,
            typename Function>
  DEAL_II_CXX20_REQUIRES(
    (std::invocable<Function,
                    decltype(*std::declval<InputIterator1>()),
                    decltype(*std::declval<InputIterator2>()),
                    decltype(*std::declval<InputIterator3>())> &&
     std::assignable_from<
       decltype(*std::declval<OutputIterator>()),
       std::invoke_result_t<Function,
                            decltype(*std::declval<InputIterator1>()),
                            decltype(*std::declval<InputIterator2>()),
                            decltype(*std::declval<InputIterator3>())>>))
  DEAL_II_DEPRECATED void transform(const InputIterator1 &begin_in1,
                                    const InputIterator1 &end_in1,
                                    InputIterator2        in2,
                                    InputIterator3        in3,
                                    OutputIterator        out,
                                    const Function       &function,
                                    const unsigned int    grainsize)
  {
#ifdef DEAL_II_WITH_TASKFLOW
    using Iterators = std::
      tuple<InputIterator1, InputIterator2, InputIterator3, OutputIterator>;
    using SyncIterators = SynchronousIterators<Iterators>;
    Iterators x_begin(begin_in1, in2, in3, out);
    Iterators x_end(end_in1,
                    InputIterator2(),
                    InputIterator3(),
                    OutputIterator());
 
    internal::taskflow_parallel_for(
      SyncIterators(x_begin),
      SyncIterators(x_end),
      [function](const auto &it) {
        *std::get<3>(it) =
          function(*std::get<0>(it), *std::get<1>(it), *std::get<2>(it));
      },
      grainsize);
 
#elif defined(DEAL_II_WITH_TBB)
    using Iterators = std::
      tuple<InputIterator1, InputIterator2, InputIterator3, OutputIterator>;
    using SyncIterators = SynchronousIterators<Iterators>;
    Iterators x_begin(begin_in1, in2, in3, out);
    Iterators x_end(end_in1,
                    InputIterator2(),
                    InputIterator3(),
                    OutputIterator());
    internal::parallel_for(
      SyncIterators(x_begin),
      SyncIterators(x_end),
      [function](const auto &range) {
        for (const auto &p : range)
          *std::get<3>(p) =
            function(*std::get<0>(p), *std::get<1>(p), *std::get<2>(p));
      },
      grainsize);
#else
    // make sure we don't get compiler
    // warnings about unused arguments
    (void)grainsize;
 
    for (OutputIterator in1 = begin_in1; in1 != end_in1;)
      *out++ = function(*in1++, *in2++, *in3++);
#endif
  }
 
 
  namespace internal
  {
#ifdef DEAL_II_WITH_TBB
    template <typename Iterator, typename Function>
    DEAL_II_CXX20_REQUIRES((std::invocable<Function, Iterator, Iterator>))
    void apply_to_subranges(const tbb::blocked_range<Iterator> &range,
                            const Function                     &f)
    {
      f(range.begin(), range.end());
    }
#endif
  } // namespace internal
 

  template <typename Iterator, typename Function>
  DEAL_II_CXX20_REQUIRES((std::invocable<Function, Iterator, Iterator>))
  void apply_to_subranges(const Iterator                             &begin,
                          const std_cxx20::type_identity_t<Iterator> &end,
                          const Function                             &f,
                          const unsigned int                          grainsize)
  {
    Assert(grainsize > 0, ExcMessage("Grainsize must be a positive number."));
 
    // Count elements in the range
    const std::size_t n = [&]() -> auto {
      if constexpr (std::is_integral_v<Iterator>)
        return static_cast<std::size_t>(end - begin);
      else
        return static_cast<std::size_t>(std::distance(begin, end));
    }();
 
    const bool worth_doing_in_parallel =
      (MultithreadInfo::n_threads() > 1) && (n > grainsize);
 
#ifdef DEAL_II_WITH_TASKFLOW
    if (worth_doing_in_parallel)
      {
        tf::Executor &executor = MultithreadInfo::get_taskflow_executor();
        tf::Taskflow  taskflow;
 
        if constexpr (std::is_integral_v<Iterator>)
          {
            // Integral "iterator" types (e.g., int): use plain arithmetic
            using integral_type = Iterator;
 
            tf::IndexRange<integral_type> range(0, n, 1);
 
            taskflow.for_each_by_index(
              range,
              [&, begin](const tf::IndexRange<integral_type> &subrange) {
                integral_type subrange_begin =
                  begin + static_cast<integral_type>(subrange.begin());
 
                integral_type subrange_end =
                  subrange_begin + static_cast<integral_type>(subrange.size());
 
                f(subrange_begin, subrange_end);
              },
              tf::GuidedPartitioner(grainsize));
          }
        else
          {
            // Non-integral iterator types: advance from 'begin' using
            // std::distance/std::advance
            using diff_t =
              typename std::iterator_traits<Iterator>::difference_type;
 
            tf::IndexRange<std::size_t> range(0, n, 1);
 
            taskflow.for_each_by_index(
              range,
              [&, begin](const tf::IndexRange<std::size_t> &subrange) {
                Iterator subrange_begin = begin;
                std::advance(subrange_begin,
                             static_cast<diff_t>(subrange.begin()));
 
                Iterator subrange_end = subrange_begin;
                std::advance(subrange_end,
                             static_cast<diff_t>(subrange.size()));
 
                f(subrange_begin, subrange_end);
              },
              tf::GuidedPartitioner(grainsize));
          }
 
        // Schedule the work. If we are within a task, we are required to use
        // corun() without wait() to avoid a potential deadlock as described in
        // https://taskflow.github.io/taskflow/ExecuteTaskflow.html#ExecuteATaskflowFromAnInternalWorker:
        if (executor.this_worker_id() != -1)
          executor.corun(taskflow);
        else
          executor.run(taskflow).wait();
 
        return;
      }
#elif defined(DEAL_II_WITH_TBB)
    if (worth_doing_in_parallel)
      {
        internal::parallel_for(
          begin,
          end,
          [&f](const tbb::blocked_range<Iterator> &range) {
            internal::apply_to_subranges<Iterator, Function>(range, f);
          },
          grainsize);
 
        return;
      }
#else
    (void)worth_doing_in_parallel;
#endif
 
    // Without library support (taskflow or TBB) or if not worth running in
    // parallel, we just call the function directly with the whole range:
    f(begin, end);
  }
 
 

  struct ParallelForInteger
  {
    virtual ~ParallelForInteger() = default;

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

  template <typename ResultType, typename Iterator, typename Function>
  DEAL_II_CXX20_REQUIRES(
    (std::invocable<Function, Iterator, Iterator> &&
     std::convertible_to<std::invoke_result_t<Function, Iterator, Iterator>,
                         ResultType>))
  ResultType
    accumulate_from_subranges(const Function                             &f,
                              const Iterator                             &begin,
                              const std_cxx20::type_identity_t<Iterator> &end,
                              const unsigned int grainsize)
  {
#ifndef DEAL_II_WITH_TBB
    // make sure we don't get compiler
    // warnings about unused arguments
    (void)grainsize;
 
    return f(begin, end);
#else
    return tbb::parallel_reduce(
      tbb::blocked_range<Iterator>(begin, end, grainsize),
      ResultType(0),
      [f](const auto &range, const ResultType &starting_value) {
        ResultType value = starting_value;
        value += f(range.begin(), range.end());
        return value;
      },
      std::plus<ResultType>(),
      tbb::auto_partitioner());
#endif
  }
 
 
  // --------------------- for loop affinity partitioner -----------------------

  namespace internal
  {
    class TBBPartitioner
    {
    public:
      TBBPartitioner();
 
#ifdef DEAL_II_WITH_TBB
      ~TBBPartitioner();

      std::shared_ptr<tbb::affinity_partitioner>
      acquire_one_partitioner();

      void
      release_one_partitioner(
        const std::shared_ptr<tbb::affinity_partitioner> &p);
 
    private:
      std::shared_ptr<tbb::affinity_partitioner> my_partitioner;

      bool in_use;

      Threads::Mutex mutex;
#endif
    };
  } // namespace internal
} // namespace parallel
 
 
namespace internal
{
  namespace VectorImplementation
  {
    extern unsigned int minimum_parallel_grain_size;
  } // namespace VectorImplementation
 
 
  namespace SparseMatrixImplementation
  {
    extern unsigned int minimum_parallel_grain_size;
  } // namespace SparseMatrixImplementation
 
} // end of namespace internal
 
 
/* --------------------------- inline functions ------------------------- */
 
namespace parallel
{
  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_TBB
    // make sure we don't get compiler
    // warnings about unused arguments
    (void)minimum_parallel_grain_size;
 
    apply_to_subrange(begin, end);
#else
    internal::parallel_for(
      begin,
      end,
      [this](const tbb::blocked_range<std::size_t> &range) {
        apply_to_subrange(range.begin(), range.end());
      },
      minimum_parallel_grain_size);
#endif
  }
 
} // end of namespace parallel
 
DEAL_II_NAMESPACE_CLOSE
 
#endif