tpl_ca_parallel_engine.H

Go to the documentation of this file.

/*
                          Aleph_w
 
  Data structures & Algorithms
  version 2.0.0b
  https://github.com/lrleon/Aleph-w
 
  This file is part of Aleph-w library
 
  Copyright (c) 2002-2026 Leandro Rabindranath Leon
 
  Permission is hereby granted, free of charge, to any person obtaining a copy
  of this software and associated documentation files (the "Software"), to deal
  in the Software without restriction, including without limitation the rights
  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  copies of the Software, and to permit persons to whom the Software is
  furnished to do so, subject to the following conditions:
 
  The above copyright notice and this permission notice shall be included in all
  copies or substantial portions of the Software.
 
  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  SOFTWARE.
*/
 
#ifndef TPL_CA_PARALLEL_ENGINE_H
#define TPL_CA_PARALLEL_ENGINE_H
 
#include <array>
#include <concepts>
#include <cstddef>
#include <functional>
#include <future>
#include <span>
#include <thread>
#include <type_traits>
#include <utility>
#include <vector>
 
#include <thread_pool.H>
 
#include <ca-tiling.H>
#include <ca-traits.H>
#include <tpl_ca_bit_storage.H>
#include <tpl_ca_concepts.H>
#include <tpl_ca_neighborhood.H>
 
namespace Aleph {
namespace CA {
 
namespace ca_parallel_detail {
 
template <typename O>
struct is_tile : std::false_type
{
};
 
template <std::size_t W, std::size_t H>
struct is_tile<Tile<W, H>> : std::true_type
{
};
 
template <typename O>
inline constexpr bool is_tile_v = is_tile<O>::value;
 
template <typename Storage>
struct is_bit_cell_storage : std::false_type
{
};
 
template <std::size_t N>
struct is_bit_cell_storage<Bit_Cell_Storage<N>> : std::true_type
{
};
 
template <typename Lattice>
inline constexpr bool uses_bit_cell_storage_v = requires {
  typename Lattice::storage_type;
} and is_bit_cell_storage<typename Lattice::storage_type>::value;
 
}  // namespace ca_parallel_detail
 
struct Parallel_Engine_Config
{
  ThreadPool *pool = nullptr;
 
  std::size_t num_partitions = 0;
 
  std::size_t min_parallel_cells = 4096;
};
 
template <typename Lattice, typename Rule, typename Neighborhood, typename Order = RowMajor>
class Parallel_Synchronous_Engine
{
  static_assert(LatticeLike<Lattice>, "Parallel_Synchronous_Engine requires a LatticeLike lattice");
  static_assert(NeighborhoodLike<Neighborhood>,
                "Parallel_Synchronous_Engine requires a NeighborhoodLike neighborhood");
  static_assert(RuleLike<Rule, Lattice>, "Parallel_Synchronous_Engine requires a RuleLike rule");
  static_assert(Neighborhood::rank_v == Lattice::rank, "Neighborhood rank must match lattice rank");
 
public:
  using lattice_type = Lattice;
  using rule_type = Rule;
  using neighborhood_type = Neighborhood;
  using order_type = Order;
 
  using state_type = typename Lattice::state_type;
  using coord_type = typename Lattice::coord_type;
  using extents_type = typename Lattice::extents_type;
 
  static constexpr std::size_t rank = Lattice::rank;
  static constexpr std::size_t neighbour_count = Neighborhood::size_v;
 
  using hook_type = std::function<void(std::size_t, const Lattice &)>;
 
private:
  Lattice cur_buf_;
  Lattice nxt_buf_;
  Rule rule_;
  Neighborhood nh_;
  hook_type pre_hook_;
  hook_type post_hook_;
  std::size_t step_count_ = 0;
  Parallel_Engine_Config cfg_{};
 
  // -------------------------------------------------------------------
  // Per-cell update used by both the sequential fallback and the
  // worker tasks. Captured by value-of-reference; safe to invoke from
  // multiple workers as long as their `c` arguments are disjoint.
  // -------------------------------------------------------------------
  void compute_cell(const coord_type &c) const = delete;  // not used
 
  void update_cell(const coord_type &c, std::span<state_type> nbuf)
  {
    const std::size_t neigh_sz = std::min(nh_.size(), nbuf.size());
    gather_neighbors(nh_, cur_buf_, c, std::span<state_type>(nbuf.data(), neigh_sz));
    const Cell_Context<rank> ctx{step_count_, c};
    nxt_buf_.set(
      c, apply_rule(rule_, cur_buf_.at(c), Neighbor_View<state_type>(nbuf.data(), neigh_sz), ctx));
  }
 
  // -------------------------------------------------------------------
  // Iteration kernels for each (rank, order) combination. They take a
  // half-open row range `[r_begin, r_end)` so the parallel scheduler
  // can hand each worker a disjoint slab.
  // -------------------------------------------------------------------
 
  void update_slab_1d(const ca_size_t r_begin, const ca_size_t r_end)
  {
    std::array<state_type, neighbour_count> nbuf{};
    coord_type c{};
    for (ca_size_t i = r_begin; i < r_end; ++i)
      {
        c[0] = static_cast<ca_index_t>(i);
        update_cell(c, nbuf);
      }
  }
 
  void update_slab_2d_row_major(const ca_size_t r_begin, const ca_size_t r_end)
  {
    std::array<state_type, neighbour_count> nbuf{};
    const ca_size_t cols = cur_buf_.size(1);
    coord_type c{};
    for (ca_size_t i = r_begin; i < r_end; ++i)
      for (ca_size_t j = 0; j < cols; ++j)
        {
          c[0] = static_cast<ca_index_t>(i);
          c[1] = static_cast<ca_index_t>(j);
          update_cell(c, nbuf);
        }
  }
 
  template <std::size_t TileW, std::size_t TileH>
  void update_slab_2d_tile(const ca_size_t r_begin, const ca_size_t r_end)
  {
    std::array<state_type, neighbour_count> nbuf{};
    const ca_size_t cols = cur_buf_.size(1);
    coord_type c{};
    // Walk tiles strictly inside the assigned row slab; the schedule
    // therefore covers each cell exactly once across all workers.
    for (ca_size_t ti = r_begin; ti < r_end; ti += TileH)
      for (ca_size_t tj = 0; tj < cols; tj += TileW)
        {
          const ca_size_t i_end = ti + TileH < r_end ? ti + TileH : r_end;
          const ca_size_t j_end = tj + TileW < cols ? tj + TileW : cols;
          for (ca_size_t i = ti; i < i_end; ++i)
            for (ca_size_t j = tj; j < j_end; ++j)
              {
                c[0] = static_cast<ca_index_t>(i);
                c[1] = static_cast<ca_index_t>(j);
                update_cell(c, nbuf);
              }
        }
  }
 
  void update_slab_3d_row_major(const ca_size_t r_begin, const ca_size_t r_end)
  {
    std::array<state_type, neighbour_count> nbuf{};
    const ca_size_t n1 = cur_buf_.size(1);
    const ca_size_t n2 = cur_buf_.size(2);
    coord_type c{};
    for (ca_size_t i = r_begin; i < r_end; ++i)
      for (ca_size_t j = 0; j < n1; ++j)
        for (ca_size_t k = 0; k < n2; ++k)
          {
            c[0] = static_cast<ca_index_t>(i);
            c[1] = static_cast<ca_index_t>(j);
            c[2] = static_cast<ca_index_t>(k);
            update_cell(c, nbuf);
          }
  }
 
  void process_slab(const ca_size_t r_begin, const ca_size_t r_end)
  {
    if constexpr (rank == 1)
      {
        update_slab_1d(r_begin, r_end);
      }
    else if constexpr (rank == 2)
      {
        if constexpr (ca_parallel_detail::is_tile_v<Order>)
          update_slab_2d_tile<Order::width, Order::height>(r_begin, r_end);
        else
          update_slab_2d_row_major(r_begin, r_end);
      }
    else if constexpr (rank == 3)
      {
        update_slab_3d_row_major(r_begin, r_end);
      }
    else
      {
        static_assert(rank <= 3, "Parallel_Synchronous_Engine currently supports rank <= 3");
      }
  }
 
  [[nodiscard]] std::size_t effective_partitions(ThreadPool *pool) const noexcept
  {
    if (pool == nullptr)
      return 1;
    const std::size_t pool_workers = pool->num_threads();
    std::size_t parts = cfg_.num_partitions == 0 ? pool_workers : cfg_.num_partitions;
    if (parts == 0)
      parts = 1;
 
    // Cannot split into more pieces than there are rows along axis 0:
    // an empty range would still cost a task scheduling roundtrip.
    const ca_size_t outer = cur_buf_.size(0);
    if (parts > outer)
      parts = outer == 0 ? 1 : outer;
 
    if (should_run_sequential(cur_buf_.size(), parts, cfg_.min_parallel_cells))
      return 1;
    return parts;
  }
 
public:
  Parallel_Synchronous_Engine(Lattice initial, Rule r, Neighborhood n, Parallel_Engine_Config cfg = {})
    : cur_buf_(std::move(initial)), nxt_buf_(cur_buf_.extents()), rule_(std::move(r)),
      nh_(std::move(n)), cfg_(cfg)
  {}
 
  Parallel_Synchronous_Engine(Lattice initial, Rule r, Parallel_Engine_Config cfg = {})
    requires std::default_initializable<Neighborhood>
    : Parallel_Synchronous_Engine(std::move(initial), std::move(r), Neighborhood{}, cfg)
  {}
 
  [[nodiscard]] const Lattice &frame() const noexcept
  {
    return cur_buf_;
  }
 
  [[nodiscard]] std::size_t steps_run() const noexcept
  {
    return step_count_;
  }
 
  [[nodiscard]] const extents_type &extents() const noexcept
  {
    return cur_buf_.extents();
  }
 
  [[nodiscard]] const Parallel_Engine_Config &config() const noexcept
  {
    return cfg_;
  }
 
  void set_config(Parallel_Engine_Config cfg) noexcept
  {
    cfg_ = cfg;
  }
 
  template <typename F>
  void on_pre_step(F &&f)
  {
    pre_hook_ = std::forward<F>(f);
  }
 
  template <typename F>
  void on_post_step(F &&f)
  {
    post_hook_ = std::forward<F>(f);
  }
 
  void step()
  {
    if constexpr (HasRefreshHalo<Lattice>)
      cur_buf_.refresh_halo();
 
    if (pre_hook_)
      pre_hook_(step_count_, cur_buf_);
 
    ThreadPool *pool = cfg_.pool != nullptr ? cfg_.pool : &Aleph::default_pool();
    const std::size_t parts = effective_partitions(pool);
    const ca_size_t outer = cur_buf_.size(0);
 
    if (parts <= 1)
      {
        // Sequential fast path. Identical iteration to the synchronous
        // engine and to the parallel branch with parts == 1.
        process_slab(0, outer);
      }
    else
      {
        // Parallel branch: enqueue one task per partition and wait for
        // every future. We collect futures up-front so the destructor
        // semantics of `std::future` cannot block us in unexpected
        // orders if the rule throws.
        std::vector<std::future<void>> futures;
        futures.reserve(parts);
        const std::array<ca_size_t, rank> ext = cur_buf_.extents();
        for (std::size_t p = 0; p < parts; ++p)
          {
            const Range1D rng = Row_Partition<rank>::slab(ext, parts, p);
            if (rng.empty())
              continue;
            futures.emplace_back(pool->enqueue([this, b = rng.begin, e = rng.end]
            {
              this->process_slab(b, e);
            }));
          }
        // Phase 1: wait for every worker to finish without re-throwing,
        // so no worker is still writing to nxt_buf_ when we propagate.
        for (auto &f : futures)
          f.wait();
        // Phase 2: collect exceptions; re-throw the first one if any.
        for (auto &f : futures)
          f.get();
      }
 
    cur_buf_.swap(nxt_buf_);
    ++step_count_;
 
    if (post_hook_)
      post_hook_(step_count_, cur_buf_);
  }
 
  void run(const std::size_t steps)
  {
    for (std::size_t i = 0; i < steps; ++i)
      step();
  }
};
 
}  // namespace CA
}  // namespace Aleph
 
#endif  // TPL_CA_PARALLEL_ENGINE_H