work_stream.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_work_stream_h
#  define dealii_work_stream_h
 
 
#  include <deal.II/base/config.h>
 
#  include <deal.II/base/array_view.h>
#  include <deal.II/base/graph_coloring.h>
#  include <deal.II/base/iterator_range.h>
#  include <deal.II/base/multithread_info.h>
#  include <deal.II/base/parallel.h>
#  include <deal.II/base/std_cxx20/type_traits.h>
#  include <deal.II/base/template_constraints.h>
#  include <deal.II/base/thread_local_storage.h>
#  include <deal.II/base/thread_management.h>
 
#  ifdef DEAL_II_WITH_TBB
#    ifdef DEAL_II_TBB_WITH_ONEAPI
#      include <tbb/parallel_pipeline.h>
#    else
#      include <tbb/pipeline.h>
#    endif
#    include <tbb/blocked_range.h>
#  endif
 
#  ifdef DEAL_II_WITH_TASKFLOW
#    include <taskflow/taskflow.hpp>
#  endif
 
#  include <functional>
#  include <iterator>
#  include <list>
#  include <memory>
#  include <utility>
#  include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
 
 
namespace WorkStream
{
  namespace internal
  {
    template <typename Iterator, typename ScratchData, typename CopyData>
    struct ScratchAndCopyDataObjects
    {
      std::unique_ptr<ScratchData> scratch_data;
      std::unique_ptr<CopyData>    copy_data;
      bool                         currently_in_use;
 
      ScratchAndCopyDataObjects()
        : currently_in_use(false)
      {}
 
      ScratchAndCopyDataObjects(std::unique_ptr<ScratchData> &&p,
                                std::unique_ptr<CopyData>    &&q,
                                const bool                     in_use)
        : scratch_data(std::move(p))
        , copy_data(std::move(q))
        , currently_in_use(in_use)
      {}
 
      // Provide a copy constructor that actually doesn't copy the
      // internal state. This makes handling ScratchAndCopyDataObjects
      // easier to handle with STL containers.
      ScratchAndCopyDataObjects(const ScratchAndCopyDataObjects &)
        : currently_in_use(false)
      {}
    };
 
    template <typename ScratchData>
    struct ScratchDataObject
    {
      std::unique_ptr<ScratchData> scratch_data;
      bool                         currently_in_use;
 
      ScratchDataObject()
        : currently_in_use(false)
      {}
 
      ScratchDataObject(std::unique_ptr<ScratchData> &&p, const bool in_use)
        : scratch_data(std::move(p))
        , currently_in_use(in_use)
      {}
 
      ScratchDataObject(ScratchData *p, const bool in_use)
        : scratch_data(p)
        , currently_in_use(in_use)
      {}
 
      // Provide a copy constructor that actually doesn't copy the
      // internal state. This makes handling ScratchAndCopyDataObjects
      // easier to handle with STL containers.
      ScratchDataObject(const ScratchDataObject &)
        : currently_in_use(false)
      {}
 
      ScratchDataObject(ScratchDataObject &&o) noexcept = default;
    };
 
#  ifdef DEAL_II_WITH_TBB
    namespace tbb_no_coloring
    {
      template <typename Iterator, typename ScratchData, typename CopyData>
      class IteratorRangeToItemStream
      {
      public:
        struct ItemType
        {
          using ScratchDataList = std::list<ScratchDataObject<ScratchData>>;
 
          std::vector<Iterator> iterators;
 
          std::vector<CopyData> copy_datas;
 
          unsigned int n_iterators;
 
          Threads::ThreadLocalStorage<ScratchDataList> *scratch_data;
 
          const ScratchData *sample_scratch_data;
 
          bool currently_in_use;
 
 
          ItemType()
            : n_iterators(0)
            , scratch_data(nullptr)
            , sample_scratch_data(nullptr)
            , currently_in_use(false)
          {}
        };
 
 
        IteratorRangeToItemStream(const Iterator    &begin,
                                  const Iterator    &end,
                                  const unsigned int buffer_size,
                                  const unsigned int chunk_size,
                                  const ScratchData &sample_scratch_data,
                                  const CopyData    &sample_copy_data)
          : remaining_iterator_range(begin, end)
          , item_buffer(buffer_size)
          , sample_scratch_data(sample_scratch_data)
          , chunk_size(chunk_size)
        {
          // initialize the elements of the ring buffer
          for (auto &item : item_buffer)
            {
              Assert(item.n_iterators == 0, ExcInternalError());
 
              item.iterators.resize(chunk_size,
                                    remaining_iterator_range.second);
              item.scratch_data        = &thread_local_scratch;
              item.sample_scratch_data = &sample_scratch_data;
              item.copy_datas.resize(chunk_size, sample_copy_data);
              item.currently_in_use = false;
            }
        }
 
 
        ItemType *
        get_item()
        {
          // find first unused item. we know that there must be one
          // because we have set the maximal number of tokens in flight
          // and have set the ring buffer to have exactly this size. so
          // if this function is called, we know that less than the
          // maximal number of items in currently in flight
          //
          // note that we need not lock access to this array since
          // the current stage is run sequentially and we can therefore
          // enter the following block only once at any given time.
          // thus, there can be no race condition between checking that
          // a flag is false and setting it to true. (there may be
          // another thread where we release items and set 'false'
          // flags to 'true', but that too does not produce any
          // problems)
          ItemType *current_item = nullptr;
          for (unsigned int i = 0; i < item_buffer.size(); ++i)
            if (item_buffer[i].currently_in_use == false)
              {
                item_buffer[i].currently_in_use = true;
                current_item                    = &item_buffer[i];
                break;
              }
          Assert(current_item != nullptr,
                 ExcMessage("This can't be. There must be a free item!"));
 
          // initialize the next item. it may
          // consist of at most chunk_size
          // elements
          current_item->n_iterators = 0;
          while ((remaining_iterator_range.first !=
                  remaining_iterator_range.second) &&
                 (current_item->n_iterators < chunk_size))
            {
              current_item->iterators[current_item->n_iterators] =
                remaining_iterator_range.first;
 
              ++remaining_iterator_range.first;
              ++current_item->n_iterators;
            }
 
          if (current_item->n_iterators == 0)
            // there were no items
            // left. terminate the pipeline
            return nullptr;
          else
            return current_item;
        }
 
      private:
        std::pair<Iterator, Iterator> remaining_iterator_range;
 
        std::vector<ItemType> item_buffer;
 
        Threads::ThreadLocalStorage<typename ItemType::ScratchDataList>
          thread_local_scratch;
 
        const ScratchData &sample_scratch_data;
 
        const unsigned int chunk_size;
      };
 
 
 
      template <typename Worker,
                typename Copier,
                typename Iterator,
                typename ScratchData,
                typename CopyData>
      void
      run(const Iterator                             &begin,
          const std_cxx20::type_identity_t<Iterator> &end,
          Worker                                      worker,
          Copier                                      copier,
          const ScratchData                          &sample_scratch_data,
          const CopyData                             &sample_copy_data,
          const unsigned int                          queue_length,
          const unsigned int                          chunk_size)
      {
        using ItemType = typename IteratorRangeToItemStream<Iterator,
                                                            ScratchData,
                                                            CopyData>::ItemType;
 
        // Define the three stages of the pipeline:
 
        //
        // ----- Stage 1 -----
        //
        // The first stage is the one that provides us with chunks of data
        // to work on (the stream of "items"). This stage will run sequentially.
        IteratorRangeToItemStream<Iterator, ScratchData, CopyData>
             iterator_range_to_item_stream(begin,
                                        end,
                                        queue_length,
                                        chunk_size,
                                        sample_scratch_data,
                                        sample_copy_data);
        auto item_generator = [&](tbb::flow_control &fc) -> ItemType * {
          if (const auto item = iterator_range_to_item_stream.get_item())
            return item;
          else
            {
              fc.stop();
              return nullptr;
            }
        };
 
        //
        // ----- Stage 2 -----
        //
        // The second stage is the one that does the actual work. This is the
        // stage that runs in parallel
        auto item_worker =
          [worker =
             std::function<void(const Iterator &, ScratchData &, CopyData &)>(
               worker),
           copier_exists =
             static_cast<bool>(std::function<void(const CopyData &)>(copier))](
            ItemType *current_item) {
            // we need to find an unused scratch data object in the list that
            // corresponds to the current thread and then mark it as used. if
            // we can't find one, create one
            //
            // as discussed in the discussion of the documentation of the
            // IteratorRangeToItemStream::scratch_data variable, there is no
            // need to synchronize access to this variable using a mutex
            // as long as we have no yield-point in between. this means that
            // we can't take an iterator into the list now and expect it to
            // still be valid after calling the worker, but we at least do
            // not have to lock the following section
            ScratchData *scratch_data = nullptr;
            {
              // see if there is an unused object. if so, grab it and mark
              // it as used
              for (auto &p : current_item->scratch_data->get())
                if (p.currently_in_use == false)
                  {
                    scratch_data       = p.scratch_data.get();
                    p.currently_in_use = true;
 
                    break;
                  }
 
              // if no object was found, create one and mark it as used
              if (scratch_data == nullptr)
                {
                  scratch_data =
                    new ScratchData(*current_item->sample_scratch_data);
                  current_item->scratch_data->get().emplace_back(scratch_data,
                                                                 true);
                }
            };
 
            // then call the worker function on each element of the chunk we
            // were given. since these worker functions are called on separate
            // threads, nothing good can happen if they throw an exception and
            // we are best off catching it and showing an error message
            for (unsigned int i = 0; i < current_item->n_iterators; ++i)
              {
                try
                  {
                    if (worker)
                      worker(current_item->iterators[i],
                             *scratch_data,
                             current_item->copy_datas[i]);
                  }
                catch (const std::exception &exc)
                  {
                    Threads::internal::handle_std_exception(exc);
                  }
                catch (...)
                  {
                    Threads::internal::handle_unknown_exception();
                  }
              }
 
            // finally mark the scratch object as unused again. as above, there
            // is no need to lock anything here since the object we work on
            // is thread-local
            for (auto &p : current_item->scratch_data->get())
              if (p.scratch_data.get() == scratch_data)
                {
                  Assert(p.currently_in_use == true, ExcInternalError());
                  p.currently_in_use = false;
 
                  break;
                }
 
            // if there is no copier, mark current item as usable again
            if (copier_exists == false)
              current_item->currently_in_use = false;
 
 
            // Then return the original pointer
            // to the now modified object. The copier will work on it next.
            return current_item;
          };
 
        //
        // ----- Stage 3 -----
        //
        // The last stage is the one that copies data from the CopyData objects
        // to the final destination. This stage runs sequentially again.
        auto item_copier = [copier = std::function<void(const CopyData &)>(
                              copier)](ItemType *current_item) {
          if (copier)
            {
              // Initiate copying data. For the same reasons as in the worker
              // class above, catch exceptions rather than letting them
              // propagate into unknown territories:
              for (unsigned int i = 0; i < current_item->n_iterators; ++i)
                {
                  try
                    {
                      copier(current_item->copy_datas[i]);
                    }
                  catch (const std::exception &exc)
                    {
                      Threads::internal::handle_std_exception(exc);
                    }
                  catch (...)
                    {
                      Threads::internal::handle_unknown_exception();
                    }
                }
            }
          // mark current item as usable again
          current_item->currently_in_use = false;
        };
 
 
        // Now we just have to set up the pipeline and run it:
        auto tbb_item_stream_filter = tbb::make_filter<void, ItemType *>(
#    ifdef DEAL_II_TBB_WITH_ONEAPI
          tbb::filter_mode::serial_in_order,
#    else
          tbb::filter::serial,
#    endif
          item_generator);
 
        auto tbb_worker_filter = tbb::make_filter<ItemType *, ItemType *>(
#    ifdef DEAL_II_TBB_WITH_ONEAPI
          tbb::filter_mode::parallel,
#    else
          tbb::filter::parallel,
#    endif
          item_worker);
 
        auto tbb_copier_filter = tbb::make_filter<ItemType *, void>(
#    ifdef DEAL_II_TBB_WITH_ONEAPI
          tbb::filter_mode::serial_in_order,
#    else
          tbb::filter::serial,
#    endif
          item_copier);
 
        tbb::parallel_pipeline(queue_length,
                               tbb_item_stream_filter & tbb_worker_filter &
                                 tbb_copier_filter);
      }
 
    }    // namespace tbb_no_coloring
#  endif // DEAL_II_WITH_TBB
 
 
 
#  ifdef DEAL_II_WITH_TASKFLOW
    namespace taskflow_no_coloring
    {
      template <typename Worker,
                typename Copier,
                typename Iterator,
                typename ScratchData,
                typename CopyData>
      void
      run(const Iterator                             &begin,
          const std_cxx20::type_identity_t<Iterator> &end,
          Worker                                      worker,
          Copier                                      copier,
          const ScratchData                          &sample_scratch_data,
          const CopyData                             &sample_copy_data,
          const unsigned int /*queue_length*/ = 2 *
                                                MultithreadInfo::n_threads(),
          const unsigned int chunk_size = 8)
 
      {
        tf::Executor &executor = MultithreadInfo::get_taskflow_executor();
        tf::Taskflow  taskflow;
 
        Threads::ThreadLocalStorage<std::deque<ScratchDataObject<ScratchData>>>
          thread_safe_scratch_datas;
 
        tf::Task last_copier;
 
 
        // A collection of chunk_size copy objects which each represent the
        // contribution of a single worker. Chunk_size workers are thus chunked
        // together by having their outputs stored in this object in a single
        // task.
        struct CopyChunk
        {
          std::vector<CopyData> copy_datas;
 
          CopyChunk(const unsigned int chunk_size,
                    const CopyData    &sample_copy_data)
            : copy_datas(chunk_size, sample_copy_data)
          {}
        };
 
        // CopyData objects are processed in the order they are created: record
        // that with a unique chunk_no for each one
        std::size_t chunk_no = 0;
 
        // A collection of handles to a CopyChunk object for each chunk. The
        // actual data will be allocated only when a worker arrives at this
        // chunk and will be freed as soon as the result has been copied.
        std::vector<std::unique_ptr<CopyChunk>> copy_chunks;
 
        // Generate a static task graph. Here we generate a task for each cell
        // that will be worked on. The tasks are not executed until all of them
        // are created, this code runs sequentially.
        for (Iterator it = begin; it != end;)
          {
            // Since Iterator can be an iterator, or an int, etc. (anything
            // incrementable) we can't rely on anything but operator++ being
            // defined (and it being copyable). We can also assume that it isn't
            // an input iterator since we can run workers in an arbitrary order.
            // Hence, we can't use std::distance(), so instead just calculate
            // each chunk size directly.
            std::size_t current_chunk_size  = 0;
            Iterator    current_chunk_start = it;
            for (unsigned int i = 0; (i < chunk_size) && (it != end); ++it, ++i)
              ++current_chunk_size;
            // Create a worker task.
            auto worker_task =
              taskflow
                .emplace([current_chunk_start,
                          current_chunk_size,
                          chunk_no,
                          &thread_safe_scratch_datas,
                          &sample_scratch_data,
                          &sample_copy_data,
                          &copy_chunks,
                          &worker]() {
                  ScratchData *scratch_data = nullptr;
 
                  auto &scratch_datas = thread_safe_scratch_datas.get();
                  // We need to find an unused scratch data object in the list
                  // that corresponds to the current thread and then mark it as
                  // used. if we can't find one, create one. There is no need to
                  // synchronize access to this variable using a mutex since
                  // each object is local to its own thread.
                  for (auto &p : scratch_datas)
                    {
                      if (p.currently_in_use == false)
                        {
                          scratch_data       = p.scratch_data.get();
                          p.currently_in_use = true;
                          break;
                        }
                    }
                  // If no element in the list was found, create
                  // one and mark it as used.
                  if (scratch_data == nullptr)
                    {
                      scratch_datas.emplace_back(std::make_unique<ScratchData>(
                                                   sample_scratch_data),
                                                 true);
                      scratch_data = scratch_datas.back().scratch_data.get();
                    }
 
                  // Create a unique copy chunk object where this
                  // worker's work will be stored.
                  copy_chunks[chunk_no] =
                    std::make_unique<CopyChunk>(current_chunk_size,
                                                sample_copy_data);
                  auto &copy_chunk = copy_chunks[chunk_no];
                  auto  it         = current_chunk_start;
                  for (std::size_t i = 0; i < current_chunk_size; ++i)
                    {
                      worker(it, *scratch_data, (copy_chunk->copy_datas)[i]);
                      ++it;
                    }
 
                  // Find our currently used scratch data and
                  // mark it as unused.
                  for (auto &p : scratch_datas)
                    {
                      if (p.scratch_data.get() == scratch_data)
                        {
                          Assert(p.currently_in_use == true,
                                 ExcInternalError());
                          p.currently_in_use = false;
                        }
                    }
                })
                .name("worker");
 
            // Create a copier task. This task is a separate object from the
            // worker task.
            tf::Task copier_task =
              taskflow
                .emplace([chunk_no, &copy_chunks, &copier]() {
                  auto copy_chunk = copy_chunks[chunk_no].get();
                  for (auto &copy_data : copy_chunk->copy_datas)
                    {
                      copier(copy_data);
                    }
                  // Finally free the memory.
                  copy_chunks[chunk_no].reset();
                })
                .name("copy");
 
            // Ensure the copy task runs after the worker task.
            worker_task.precede(copier_task);
 
            // Ensure that only one copy task can run at a time. The code
            // below makes each copy task wait until the previous one has
            // finished before it can start
            if (!last_copier.empty())
              last_copier.precede(copier_task);
 
            // Keep a handle to the last copier. Tasks in taskflow are
            // basically handles to internally stored data, so this does not
            // perform a copy:
            last_copier = copier_task;
            ++chunk_no;
          }
        copy_chunks.resize(chunk_no);
        // Now we run all the tasks in the task graph. They will be run in
        // parallel and are eligible to run when their dependencies established
        // above are met.
 
        // 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();
      }
    } // namespace taskflow_no_coloring
#  endif
 
    namespace sequential
    {
      template <typename Worker,
                typename Copier,
                typename Iterator,
                typename ScratchData,
                typename CopyData>
      void
      run(const Iterator                             &begin,
          const std_cxx20::type_identity_t<Iterator> &end,
          Worker                                      worker,
          Copier                                      copier,
          const ScratchData                          &sample_scratch_data,
          const CopyData                             &sample_copy_data)
      {
        // need to copy the sample since it is marked const
        ScratchData scratch_data = sample_scratch_data;
        CopyData    copy_data    = sample_copy_data; // NOLINT
 
        // Optimization: Check if the functions are not the zero function. To
        // check zero-ness, create a C++ function out of it:
        const bool have_worker =
          (static_cast<const std::function<
             void(const Iterator &, ScratchData &, CopyData &)> &>(worker)) !=
          nullptr;
        const bool have_copier =
          (static_cast<const std::function<void(const CopyData &)> &>(
            copier)) != nullptr;
 
        // Finally loop over all items and perform the necessary work:
        for (Iterator i = begin; i != end; ++i)
          {
            if (have_worker)
              worker(i, scratch_data, copy_data);
            if (have_copier)
              copier(copy_data);
          }
      }
 
 
 
      template <typename Worker,
                typename Copier,
                typename Iterator,
                typename ScratchData,
                typename CopyData>
      void
      run(const std::vector<std::vector<Iterator>> &colored_iterators,
          Worker                                    worker,
          Copier                                    copier,
          const ScratchData                        &sample_scratch_data,
          const CopyData                           &sample_copy_data)
      {
        // need to copy the sample since it is marked const
        ScratchData scratch_data = sample_scratch_data;
        CopyData    copy_data    = sample_copy_data; // NOLINT
 
        // Optimization: Check if the functions are not the zero function. To
        // check zero-ness, create a C++ function out of it:
        const bool have_worker =
          (static_cast<const std::function<
             void(const Iterator &, ScratchData &, CopyData &)> &>(worker)) !=
          nullptr;
        const bool have_copier =
          (static_cast<const std::function<void(const CopyData &)> &>(
            copier)) != nullptr;
 
        // Finally loop over all items and perform the necessary work:
        for (unsigned int color = 0; color < colored_iterators.size(); ++color)
          if (colored_iterators[color].size() > 0)
            for (auto &it : colored_iterators[color])
              {
                if (have_worker)
                  worker(it, scratch_data, copy_data);
                if (have_copier)
                  copier(copy_data);
              }
      }
 
    } // namespace sequential
 
 
 
#  ifdef DEAL_II_WITH_TBB
    namespace tbb_colored
    {
      template <typename Iterator, typename ScratchData, typename CopyData>
      class WorkerAndCopier
      {
      public:
        WorkerAndCopier(
          const std::function<void(const Iterator &, ScratchData &, CopyData &)>
                                                      &worker,
          const std::function<void(const CopyData &)> &copier,
          const ScratchData                           &sample_scratch_data,
          const CopyData                              &sample_copy_data)
          : worker(worker)
          , copier(copier)
          , sample_scratch_data(sample_scratch_data)
          , sample_copy_data(sample_copy_data)
        {}
 
 
        void
        operator()(const tbb::blocked_range<
                   typename std::vector<Iterator>::const_iterator> &range)
        {
          // we need to find an unused scratch and corresponding copy
          // data object in the list that corresponds to the current
          // thread and then mark it as used. If we can't find one,
          // create one as discussed in the discussion of the documentation
          // of the IteratorRangeToItemStream::scratch_data variable,
          // there is no need to synchronize access to this variable
          // using a mutex as long as we have no yield-point in between.
          // This means that we can't take an iterator into the list
          // now and expect it to still be valid after calling the worker,
          // but we at least do not have to lock the following section.
          ScratchData *scratch_data = nullptr;
          CopyData    *copy_data    = nullptr;
          {
            ScratchAndCopyDataList &scratch_and_copy_data_list = data.get();
 
            // see if there is an unused object. if so, grab it and mark
            // it as used
            for (typename ScratchAndCopyDataList::iterator p =
                   scratch_and_copy_data_list.begin();
                 p != scratch_and_copy_data_list.end();
                 ++p)
              if (p->currently_in_use == false)
                {
                  scratch_data        = p->scratch_data.get();
                  copy_data           = p->copy_data.get();
                  p->currently_in_use = true;
                  break;
                }
 
            // if no element in the list was found, create one and mark it as
            // used
            if (scratch_data == nullptr)
              {
                Assert(copy_data == nullptr, ExcInternalError());
 
                scratch_and_copy_data_list.emplace_back(
                  std::make_unique<ScratchData>(sample_scratch_data),
                  std::make_unique<CopyData>(sample_copy_data),
                  true);
                scratch_data =
                  scratch_and_copy_data_list.back().scratch_data.get();
                copy_data = scratch_and_copy_data_list.back().copy_data.get();
              }
          }
 
          // then call the worker and copier functions on each
          // element of the chunk we were given.
          for (typename std::vector<Iterator>::const_iterator p = range.begin();
               p != range.end();
               ++p)
            {
              try
                {
                  if (worker)
                    worker(*p, *scratch_data, *copy_data);
                  if (copier)
                    copier(*copy_data);
                }
              catch (const std::exception &exc)
                {
                  Threads::internal::handle_std_exception(exc);
                }
              catch (...)
                {
                  Threads::internal::handle_unknown_exception();
                }
            }
 
          // finally mark the scratch object as unused again. as above, there
          // is no need to lock anything here since the object we work on
          // is thread-local
          {
            ScratchAndCopyDataList &scratch_and_copy_data_list = data.get();
 
            for (typename ScratchAndCopyDataList::iterator p =
                   scratch_and_copy_data_list.begin();
                 p != scratch_and_copy_data_list.end();
                 ++p)
              if (p->scratch_data.get() == scratch_data)
                {
                  Assert(p->currently_in_use == true, ExcInternalError());
                  p->currently_in_use = false;
                }
          }
        }
 
      private:
        using ScratchAndCopyDataObjects = typename internal::
          ScratchAndCopyDataObjects<Iterator, ScratchData, CopyData>;
 
        using ScratchAndCopyDataList = std::list<ScratchAndCopyDataObjects>;
 
        Threads::ThreadLocalStorage<ScratchAndCopyDataList> data;
 
        const std::function<void(const Iterator &, ScratchData &, CopyData &)>
          worker;
 
        const std::function<void(const CopyData &)> copier;
 
        const ScratchData &sample_scratch_data;
        const CopyData    &sample_copy_data;
      };
 
      template <typename Worker,
                typename Copier,
                typename Iterator,
                typename ScratchData,
                typename CopyData>
      void
      run(const std::vector<std::vector<Iterator>> &colored_iterators,
          Worker                                    worker,
          Copier                                    copier,
          const ScratchData                        &sample_scratch_data,
          const CopyData                           &sample_copy_data,
          const unsigned int                        chunk_size)
      {
        // loop over the various colors of what we're given
        for (unsigned int color = 0; color < colored_iterators.size(); ++color)
          if (colored_iterators[color].size() > 0)
            {
              using WorkerAndCopier = internal::tbb_colored::
                WorkerAndCopier<Iterator, ScratchData, CopyData>;
 
              WorkerAndCopier worker_and_copier(worker,
                                                copier,
                                                sample_scratch_data,
                                                sample_copy_data);
 
              parallel::internal::parallel_for(
                colored_iterators[color].begin(),
                colored_iterators[color].end(),
                [&worker_and_copier](
                  const tbb::blocked_range<
                    typename std::vector<Iterator>::const_iterator> &range) {
                  worker_and_copier(range);
                },
                chunk_size);
            }
      }
 
    }    // namespace tbb_colored
#  endif // DEAL_II_WITH_TBB
 
 
 
#  ifdef DEAL_II_WITH_TASKFLOW
    namespace taskflow_colored
    {
      template <typename Iterator, typename ScratchData, typename CopyData>
      class WorkerAndCopier
      {
      public:
        WorkerAndCopier(
          const std::function<void(const Iterator &, ScratchData &, CopyData &)>
                                                      &worker,
          const std::function<void(const CopyData &)> &copier,
          const ScratchData                           &sample_scratch_data,
          const CopyData                              &sample_copy_data)
          : worker(worker)
          , copier(copier)
          , sample_scratch_data(sample_scratch_data)
          , sample_copy_data(sample_copy_data)
        {}
 
 
        void
        operator()(const ArrayView<const Iterator> &range)
        {
          // we need to find an unused scratch and corresponding copy
          // data object in the list that corresponds to the current
          // thread and then mark it as used. If we can't find one,
          // create one as discussed in the discussion of the documentation
          // of the IteratorRangeToItemStream::scratch_data variable,
          // there is no need to synchronize access to this variable
          // using a mutex as long as we have no yield-point in between.
          // This means that we can't take an iterator into the list
          // now and expect it to still be valid after calling the worker,
          // but we at least do not have to lock the following section.
          ScratchData *scratch_data = nullptr;
          CopyData    *copy_data    = nullptr;
          {
            ScratchAndCopyDataList &scratch_and_copy_data_list = data.get();
 
            // see if there is an unused object. if so, grab it and mark
            // it as used
            for (typename ScratchAndCopyDataList::iterator p =
                   scratch_and_copy_data_list.begin();
                 p != scratch_and_copy_data_list.end();
                 ++p)
              if (p->currently_in_use == false)
                {
                  scratch_data        = p->scratch_data.get();
                  copy_data           = p->copy_data.get();
                  p->currently_in_use = true;
                  break;
                }
 
            // if no element in the list was found, create one and mark it as
            // used
            if (scratch_data == nullptr)
              {
                Assert(copy_data == nullptr, ExcInternalError());
 
                scratch_and_copy_data_list.emplace_back(
                  std::make_unique<ScratchData>(sample_scratch_data),
                  std::make_unique<CopyData>(sample_copy_data),
                  true);
                scratch_data =
                  scratch_and_copy_data_list.back().scratch_data.get();
                copy_data = scratch_and_copy_data_list.back().copy_data.get();
              }
          }
 
          // then call the worker and copier functions on each
          // element of the chunk we were given.
 
          try
            {
              for (const auto &it : range)
                {
                  if (worker)
                    worker(it, *scratch_data, *copy_data);
                  if (copier)
                    copier(*copy_data);
                }
            }
          catch (const std::exception &exc)
            {
              Threads::internal::handle_std_exception(exc);
            }
          catch (...)
            {
              Threads::internal::handle_unknown_exception();
            }
 
 
 
          // finally mark the scratch object as unused again. as above, there
          // is no need to lock anything here since the object we work on
          // is thread-local
          {
            ScratchAndCopyDataList &scratch_and_copy_data_list = data.get();
 
            for (typename ScratchAndCopyDataList::iterator p =
                   scratch_and_copy_data_list.begin();
                 p != scratch_and_copy_data_list.end();
                 ++p)
              if (p->scratch_data.get() == scratch_data)
                {
                  Assert(p->currently_in_use == true, ExcInternalError());
                  p->currently_in_use = false;
                }
          }
        }
 
      private:
        using ScratchAndCopyDataObjects = typename internal::
          ScratchAndCopyDataObjects<Iterator, ScratchData, CopyData>;
 
        using ScratchAndCopyDataList = std::list<ScratchAndCopyDataObjects>;
 
        Threads::ThreadLocalStorage<ScratchAndCopyDataList> data;
 
        const std::function<void(const Iterator &, ScratchData &, CopyData &)>
          worker;
 
        const std::function<void(const CopyData &)> copier;
 
        const ScratchData &sample_scratch_data;
        const CopyData    &sample_copy_data;
      };
 
      template <typename Worker,
                typename Copier,
                typename Iterator,
                typename ScratchData,
                typename CopyData>
      void
      run(const std::vector<std::vector<Iterator>> &colored_iterators,
          const Worker                             &worker,
          const Copier                             &copier,
          const ScratchData                        &sample_scratch_data,
          const CopyData                           &sample_copy_data,
          const unsigned int                        chunk_size)
      {
        using WorkerAndCopier = internal::taskflow_colored::
          WorkerAndCopier<Iterator, ScratchData, CopyData>;
 
        WorkerAndCopier worker_and_copier(worker,
                                          copier,
                                          sample_scratch_data,
                                          sample_copy_data);
 
        tf::Executor &executor = MultithreadInfo::get_taskflow_executor();
        tf::Taskflow  taskflow;
 
        tf::Task last_task;
        for (auto &color : colored_iterators)
          if (color.size() > 0)
            {
              const tf::IndexRange<int> range(0, color.size(), 1);
              auto                      task = taskflow.for_each_by_index(
                range,
                [&color,
                 &worker_and_copier](const tf::IndexRange<int> &subrange) {
                  const ArrayView<const Iterator> chunk(
                    &color[subrange.begin()], subrange.size());
                  worker_and_copier(chunk);
                },
                tf::GuidedPartitioner(chunk_size));
 
              if (!last_task.empty())
                last_task.precede(task);
 
              last_task = task;
            }
 
        // 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();
      }
 
    }    // namespace taskflow_colored
#  endif // DEAL_II_WITH_TASKFLOW
 
 
  } // namespace internal
 
 
 
  template <typename Worker,
            typename Copier,
            typename Iterator,
            typename ScratchData,
            typename CopyData>
  void
  run(const std::vector<std::vector<Iterator>> &colored_iterators,
      Worker                                    worker,
      Copier                                    copier,
      const ScratchData                        &sample_scratch_data,
      const CopyData                           &sample_copy_data,
      const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
      const unsigned int chunk_size   = 8);
 
 
  template <typename Worker,
            typename Copier,
            typename Iterator,
            typename ScratchData,
            typename CopyData>
  void
  run(const Iterator                             &begin,
      const std_cxx20::type_identity_t<Iterator> &end,
      Worker                                      worker,
      Copier                                      copier,
      const ScratchData                          &sample_scratch_data,
      const CopyData                             &sample_copy_data,
      const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
      const unsigned int chunk_size   = 8)
  {
    Assert(queue_length > 0,
           ExcMessage("The queue length must be at least one, and preferably "
                      "larger than the number of processors on this system."));
    (void)queue_length; // removes -Wunused-parameter warning in optimized mode
    Assert(chunk_size > 0, ExcMessage("The chunk_size must be at least one."));
    (void)chunk_size; // removes -Wunused-parameter warning in optimized mode
 
    // If no work then skip. (only use operator!= for iterators since we may
    // not have an equality comparison operator)
    if (!(begin != end))
      return;
 
    if (MultithreadInfo::n_threads() > 1)
      {
#  if defined(DEAL_II_WITH_TBB) || defined(DEAL_II_WITH_TASKFLOW)
        if (static_cast<const std::function<void(const CopyData &)> &>(copier))
          {
            // If we have a copier, run the algorithm:
#    if defined(DEAL_II_WITH_TASKFLOW)
            internal::taskflow_no_coloring::run(begin,
                                                end,
                                                worker,
                                                copier,
                                                sample_scratch_data,
                                                sample_copy_data,
                                                queue_length,
                                                chunk_size);
#    elif defined(DEAL_II_WITH_TBB)
            internal::tbb_no_coloring::run(begin,
                                           end,
                                           worker,
                                           copier,
                                           sample_scratch_data,
                                           sample_copy_data,
                                           queue_length,
                                           chunk_size);
#    endif
          }
        else
          {
            // There is no copier function. in this case, we have an
            // embarrassingly parallel problem where we can
            // essentially apply parallel_for. because parallel_for
            // requires subdividing the range for which operator- is
            // necessary between iterators, it is often inefficient to
            // apply it directly to cell ranges and similar iterator
            // types for which operator- is expensive or, in fact,
            // nonexistent. rather, in that case, we simply copy the
            // iterators into a large array and use operator- on
            // iterators to this array of iterators.
            //
            // instead of duplicating code, this is essentially the
            // same situation we have in the colored implementation below, so we
            // just defer to that place
            std::vector<std::vector<Iterator>> all_iterators(1);
            for (Iterator p = begin; p != end; ++p)
              all_iterators[0].push_back(p);
 
            run(all_iterators,
                worker,
                copier,
                sample_scratch_data,
                sample_copy_data,
                queue_length,
                chunk_size);
          }
 
        // exit this function to not run the sequential version below:
        return;
#  endif
      }
 
    // no TBB or Taskflow installed or we are requested to run sequentially:
    internal::sequential::run(
      begin, end, worker, copier, sample_scratch_data, sample_copy_data);
  }
 
 
 
  template <
    typename Worker,
    typename Copier,
    typename IteratorRangeType,
    typename ScratchData,
    typename CopyData,
    typename = std::enable_if_t<
      has_begin_and_end<IteratorRangeType> &&
      !std::is_same_v<IteratorRangeType,
                      IteratorRange<typename IteratorRangeType::iterator>>>>
  void
  run(IteratorRangeType  iterator_range,
      Worker             worker,
      Copier             copier,
      const ScratchData &sample_scratch_data,
      const CopyData    &sample_copy_data,
      const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
      const unsigned int chunk_size   = 8)
  {
    // Call the function above
    run(iterator_range.begin(),
        iterator_range.end(),
        worker,
        copier,
        sample_scratch_data,
        sample_copy_data,
        queue_length,
        chunk_size);
  }
 
 
 
  template <typename Worker,
            typename Copier,
            typename Iterator,
            typename ScratchData,
            typename CopyData>
  void
  run(const IteratorRange<Iterator> &iterator_range,
      Worker                         worker,
      Copier                         copier,
      const ScratchData             &sample_scratch_data,
      const CopyData                &sample_copy_data,
      const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
      const unsigned int chunk_size   = 8)
  {
    // Call the function above
    run(iterator_range.begin(),
        iterator_range.end(),
        worker,
        copier,
        sample_scratch_data,
        sample_copy_data,
        queue_length,
        chunk_size);
  }
 
 
 
  template <typename Worker,
            typename Copier,
            typename Iterator,
            typename ScratchData,
            typename CopyData>
  void
  run(const std::vector<std::vector<Iterator>> &colored_iterators,
      Worker                                    worker,
      Copier                                    copier,
      const ScratchData                        &sample_scratch_data,
      const CopyData                           &sample_copy_data,
      const unsigned int                        queue_length,
      const unsigned int                        chunk_size)
  {
    Assert(queue_length > 0,
           ExcMessage("The queue length must be at least one, and preferably "
                      "larger than the number of processors on this system."));
    (void)queue_length; // removes -Wunused-parameter warning in optimized mode
    Assert(chunk_size > 0, ExcMessage("The chunk_size must be at least one."));
    (void)chunk_size; // removes -Wunused-parameter warning in optimized mode
 
 
    if (MultithreadInfo::n_threads() > 1)
      {
#  ifdef DEAL_II_WITH_TASKFLOW
        internal::taskflow_colored::run(colored_iterators,
                                        worker,
                                        copier,
                                        sample_scratch_data,
                                        sample_copy_data,
                                        chunk_size);
 
        // exit this function to not run the sequential version below:
        return;
#  elif defined(DEAL_II_WITH_TBB)
        internal::tbb_colored::run(colored_iterators,
                                   worker,
                                   copier,
                                   sample_scratch_data,
                                   sample_copy_data,
                                   chunk_size);
 
        // exit this function to not run the sequential version below:
        return;
#  endif
      }
 
    // run all colors sequentially:
    {
      internal::sequential::run(colored_iterators,
                                worker,
                                copier,
                                sample_scratch_data,
                                sample_copy_data);
    }
  }
 
 
 
  template <typename MainClass,
            typename Iterator,
            typename ScratchData,
            typename CopyData>
  void
  run(const Iterator                             &begin,
      const std_cxx20::type_identity_t<Iterator> &end,
      MainClass                                  &main_object,
      void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &),
      void (MainClass::*copier)(const CopyData &),
      const ScratchData &sample_scratch_data,
      const CopyData    &sample_copy_data,
      const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
      const unsigned int chunk_size   = 8)
  {
    // forward to the other function
    run(
      begin,
      end,
      [&main_object, worker](const Iterator &iterator,
                             ScratchData    &scratch_data,
                             CopyData       &copy_data) {
        (main_object.*worker)(iterator, scratch_data, copy_data);
      },
      [&main_object, copier](const CopyData &copy_data) {
        (main_object.*copier)(copy_data);
      },
      sample_scratch_data,
      sample_copy_data,
      queue_length,
      chunk_size);
  }
 
 
  template <typename MainClass,
            typename Iterator,
            typename ScratchData,
            typename CopyData>
  void
  run(const IteratorOverIterators<Iterator>                             &begin,
      const IteratorOverIterators<std_cxx20::type_identity_t<Iterator>> &end,
      MainClass &main_object,
      void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &),
      void (MainClass::*copier)(const CopyData &),
      const ScratchData &sample_scratch_data,
      const CopyData    &sample_copy_data,
      const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
      const unsigned int chunk_size   = 8)
  {
    // forward to the other function
    run(
      begin,
      end,
      [&main_object, worker](const Iterator &iterator,
                             ScratchData    &scratch_data,
                             CopyData       &copy_data) {
        (main_object.*worker)(iterator, scratch_data, copy_data);
      },
      [&main_object, copier](const CopyData &copy_data) {
        (main_object.*copier)(copy_data);
      },
      sample_scratch_data,
      sample_copy_data,
      queue_length,
      chunk_size);
  }
 
 
 
  template <
    typename MainClass,
    typename IteratorRangeType,
    typename ScratchData,
    typename CopyData,
    typename = std::enable_if_t<
      has_begin_and_end<IteratorRangeType> &&
      !std::is_same_v<IteratorRangeType,
                      IteratorRange<typename IteratorRangeType::iterator>>>>
  void
  run(
    IteratorRangeType iterator_range,
    MainClass        &main_object,
    void (MainClass::*worker)(
      const typename std_cxx20::type_identity_t<IteratorRangeType>::iterator &,
      ScratchData &,
      CopyData &),
    void (MainClass::*copier)(const CopyData &),
    const ScratchData &sample_scratch_data,
    const CopyData    &sample_copy_data,
    const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
    const unsigned int chunk_size   = 8)
  {
    // Call the function above
    run(std::begin(iterator_range),
        std::end(iterator_range),
        main_object,
        worker,
        copier,
        sample_scratch_data,
        sample_copy_data,
        queue_length,
        chunk_size);
  }
 
 
 
  template <typename MainClass,
            typename Iterator,
            typename ScratchData,
            typename CopyData>
  void
  run(IteratorRange<Iterator> iterator_range,
      MainClass              &main_object,
      void (MainClass::*worker)(const Iterator &, ScratchData &, CopyData &),
      void (MainClass::*copier)(const CopyData &),
      const ScratchData &sample_scratch_data,
      const CopyData    &sample_copy_data,
      const unsigned int queue_length = 2 * MultithreadInfo::n_threads(),
      const unsigned int chunk_size   = 8)
  {
    // Call the function above
    run(std::begin(iterator_range),
        std::end(iterator_range),
        main_object,
        worker,
        copier,
        sample_scratch_data,
        sample_copy_data,
        queue_length,
        chunk_size);
  }
 
} // namespace WorkStream
 
 
 
DEAL_II_NAMESPACE_CLOSE
 
 
 
//----------------------------   work_stream.h     ---------------------------
// end of #ifndef dealii_work_stream_h
#endif
//----------------------------   work_stream.h     ---------------------------