parallel.hpp

Go to the documentation of this file.

 
#ifndef NUMERICAL_PARALLEL_HPP
#define NUMERICAL_PARALLEL_HPP
 
#include <algorithm>
#include <atomic>
#include <chrono>
#include <condition_variable>
#include <cstddef>
#include <exception>
#include <limits>
#include <memory>
#include <mutex>
#include <span>
#include <thread>
#include <vector>
 
#if defined(_WIN32)
#  include <process.h>
#else
#  include <unistd.h>
#endif
#if defined(__linux__)
#  include <fstream>
#  include <map>
#  include <string>
#  include <utility>
#  include <sched.h>
#elif defined(__APPLE__)
#  include <sys/sysctl.h>
#endif
 
namespace numerical::detail {
 
  [[nodiscard]] inline long current_process() noexcept
  {
#if defined(_WIN32)
    return static_cast<long>(_getpid());
#else
    return static_cast<long>(getpid());
#endif
  }
 
  inline constexpr std::size_t no_core {std::numeric_limits<std::size_t>::max()};
 
  [[nodiscard]] inline std::vector<std::size_t> worker_cores(std::size_t cores,
                                                             std::size_t threads)
  {
    std::vector<std::size_t> out(cores != 0 ? cores : threads - 1, no_core);
    for (std::size_t c {0}; c < cores; ++c) {
      out[c] = c;
    }
    return out;
  }
 
  [[nodiscard]] constexpr bool sits_out(std::size_t worker,
                                        std::size_t caller) noexcept
  {
    return worker != no_core && worker == caller;
  }
 
#if defined(__linux__)
  [[nodiscard]] inline std::string first_line(const std::string& path)
  {
    std::ifstream file {path};
    std::string   line;
    std::getline(file, line);
    return line;
  }
 
  [[nodiscard]] inline std::vector<cpu_set_t> core_sets_from(const cpu_set_t&   allowed,
                                                             const std::string& cpu_root) noexcept
  {
    try {
      std::map<std::pair<long, long>, std::size_t> index; // (package, core) → its set
      std::vector<cpu_set_t>                       sets;
      for (int cpu {0}; cpu < CPU_SETSIZE; ++cpu) {
        if (CPU_ISSET(cpu, &allowed) == 0) { // NOLINT(hicpp-signed-bitwise): the system's macro
          continue;
        }
        const std::string topology {cpu_root + "cpu" + std::to_string(cpu) + "/topology/"};
        const long        core     {std::stol(first_line(topology + "core_id"))};
        const long        package  {std::stol(first_line(topology + "physical_package_id"))};
        const auto [where, added] {index.try_emplace({package, core}, sets.size())};
        if (added) {
          sets.emplace_back();
          CPU_ZERO(&sets.back()); // NOLINT(hicpp-signed-bitwise): the system's macro
        }
        CPU_SET(cpu, &sets[where->second]); // NOLINT(hicpp-signed-bitwise): the system's macro
      }
      return sets;
    }
    catch (...) {
      return {};
    }
  }
 
  [[nodiscard]] inline std::vector<cpu_set_t> core_sets() noexcept
  {
    cpu_set_t allowed;
    CPU_ZERO(&allowed);                                                 // NOLINT(hicpp-signed-bitwise): the system's macro
    static_cast<void>(sched_getaffinity(0, sizeof(allowed), &allowed)); // on failure the set stays empty
    return core_sets_from(allowed, "/sys/devices/system/cpu/");
  }
 
  [[nodiscard]] inline std::vector<std::size_t> processor_cores(const std::vector<cpu_set_t>& cores)
  {
    std::vector<std::size_t> table(CPU_SETSIZE, no_core);
    for (std::size_t c {0}; c < cores.size(); ++c) {
      for (std::size_t cpu {0}; cpu < table.size(); ++cpu) {
        if (CPU_ISSET(cpu, &cores[c]) != 0) { // NOLINT(hicpp-signed-bitwise): the system's macro
          table[cpu] = c;
        }
      }
    }
    return table;
  }
 
  [[nodiscard]] inline std::size_t core_of(const std::vector<std::size_t>& table,
                                           int                             cpu) noexcept
  {
    // Through the subspan from `cpu`: GCC 15 at -O3 reads `table[cpu]` behind an index check (a conditional or an
    // early return) before the check, and refuses the subscript past the end (-Werror=array-bounds).
    const std::span<const std::size_t> rest {std::span {table}.subspan(cpu < 0 ? table.size() : std::min(table.size(), static_cast<std::size_t>(cpu)))};
    return rest.empty() ? no_core : rest.front();
  }
 
  inline void pin_to_core(const std::vector<cpu_set_t>& cores,
                          std::size_t                   index) noexcept
  {
    if (index < cores.size()) {
      (void)sched_setaffinity(0, sizeof(cpu_set_t), &cores[index]);
    }
  }
 
#endif
 
  [[nodiscard]] inline std::size_t physical_cores() noexcept
  {
#if defined(__linux__)
    return core_sets().size();
#elif defined(__APPLE__)
    int         cores {0};
    std::size_t size  {sizeof(cores)};
    (void)sysctlbyname("hw.physicalcpu", &cores, &size, nullptr, 0); // on failure cores stays 0
    return static_cast<std::size_t>(std::max(cores, 0));
#else
    return 0;
#endif
  }
 
  [[nodiscard]] constexpr std::size_t pool_threads(std::size_t logical,
                                                   std::size_t physical) noexcept
  {
    const std::size_t hardware {std::max<std::size_t>(1, logical)};
    return physical == 0 ? hardware : std::min(physical, hardware);
  }
 
  [[nodiscard]] inline std::size_t hardware_threads() noexcept
  {
    static const std::size_t threads {pool_threads(std::thread::hardware_concurrency(), physical_cores())};
    return threads;
  }
 
  inline constexpr std::chrono::microseconds spin_time {500};
 
  template <typename Ready>
  void spin_until(const Ready& ready)
  {
    const std::chrono::steady_clock::time_point until {std::chrono::steady_clock::now() + spin_time};
    while (!ready() && std::chrono::steady_clock::now() < until) {
      std::this_thread::yield();
    }
  }
 
  using chunk_runner = void (*)(const void* body, std::size_t begin, std::size_t end);
 
  using side_runner = void (*)(const void* side);
 
  struct pool_job
  {
    chunk_runner                    run      {};        
    const void*                     body     {};        
    std::size_t                     count    {};        
    std::size_t                     chunk    {};        
    std::size_t                     chunks   {};        
    std::size_t                     helpers  {};        
    std::size_t                     core     {no_core}; 
    std::atomic<std::size_t>        joined   {0};       
    std::atomic<std::size_t>        next     {0};       
    std::atomic<std::size_t>        finished {0};       
    std::mutex                      lock;               
    std::condition_variable         done;               
    std::vector<std::exception_ptr> failed;             
 
    void work() noexcept
    {
      for (std::size_t c {next.fetch_add(1)}; c < chunks; c = next.fetch_add(1)) {
        try {
          run(body, c * chunk, std::min(count, (c + 1) * chunk));
        }
        catch (...) {
          failed[c] = std::current_exception();
        }
        if (finished.fetch_add(1) + 1 == chunks) {
          const std::lock_guard<std::mutex> held {lock};
          done.notify_all();
        }
      }
    }
  };
 
  struct worker_bell
  {
    explicit worker_bell(std::size_t on) : core {on}
    {}
 
    const std::size_t       core;             
    std::atomic<bool>       asleep {false};   
    std::mutex              lock;             
    std::condition_variable rung;             
  };
 
  [[nodiscard]] constexpr bool wants(std::size_t core,
                                     std::size_t caller,
                                     std::size_t helpers) noexcept
  {
    if (core == no_core) {
      return true;
    }
    if (core == caller) {
      return false;
    }
    const std::size_t rank {caller == no_core || core < caller ? core : core - 1};
    return rank < helpers;
  }
 
  class thread_pool
  {
  public:
 
    thread_pool(const thread_pool&)            = delete;
    thread_pool& operator=(const thread_pool&) = delete;
    thread_pool(thread_pool&&)                 = delete;
    thread_pool& operator=(thread_pool&&)      = delete;
    ~thread_pool()                             = delete;
 
    static thread_pool& shared()
    {
      static thread_pool* const pool {[] {
                                        auto* const started {new thread_pool {}};
                                        started->start();
                                        return started;
                                      }()};
      return *pool;
    }
 
    bool try_run(chunk_runner run,
                 const void*  body,
                 std::size_t  count,
                 std::size_t  chunk,
                 std::size_t  chunks,
                 std::size_t  helpers,
                 side_runner  side      = nullptr,
                 const void*  side_body = nullptr)
    {
      bool idle {false};
      // A forked child has none of the parent's workers, and their locks may be held for good: it runs alone.
      if (owner_ != current_process() || !busy_.compare_exchange_strong(idle, true)) {
        return false;
      }
      struct release
      {
        std::atomic<bool>& busy;
 
        explicit release(std::atomic<bool>& flag) : busy {flag}
        {}
 
        release(const release&)            = delete;
        release& operator=(const release&) = delete;
        release(release&&)                 = delete;
        release& operator=(release&&)      = delete;
        ~release()
        {
          busy.store(false);
        }
      } const freed  {busy_};
      const auto job {std::make_shared<pool_job>()};
      job->run     = run;
      job->body    = body;
      job->count   = count;
      job->chunk   = chunk;
      job->chunks  = chunks;
      job->helpers = helpers;
#if defined(__linux__)
      job->core = core_of(processor_cores_, sched_getcpu());
#endif
      job->failed.resize(chunks);
      {
        const std::lock_guard<std::mutex> held {lock_};
        current_ = job;
        generation_.fetch_add(1);
      }
      for (std::size_t w {0}; w < bells_.size(); ++w) { // the workers the run wants, asleep: not the one on the caller's core
        worker_bell& bell {*bells_[w]};
#if defined(__linux__)                                  // elsewhere no worker is pinned, and every one may join
        if (!wants(bell.core, job->core, helpers)) {
          continue;
        }
#endif
        if (bell.asleep.load()) {
          { const std::lock_guard<std::mutex> held {bell.lock}; }
          bell.rung.notify_one();
        }
      }
      std::exception_ptr side_failed;
      if (side != nullptr) {
        try {
          side(side_body);
        }
        catch (...) {
          side_failed = std::current_exception();
        }
      }
      job->work();
      const auto all_done {[&job] { return job->finished.load() == job->chunks; }};
      spin_until(all_done);
      {
        std::unique_lock<std::mutex> held {job->lock};
        job->done.wait(held, all_done);
      }
      if (side_failed) {
        std::rethrow_exception(side_failed);
      }
      for (const std::exception_ptr& error : job->failed) {
        if (error) {
          std::rethrow_exception(error);
        }
      }
      return true;
    }
 
  private:
 
    thread_pool() = default;
 
    void start()
    {
#if defined(__linux__)
      const std::vector<cpu_set_t> cores {core_sets()};
      processor_cores_ = processor_cores(cores);
      for (const std::size_t core : worker_cores(cores.size(), hardware_threads())) {
        bells_.push_back(std::make_unique<worker_bell>(core));
        worker_bell* const bell {bells_.back().get()};
        std::thread {[this, cores, core, bell] {
                       pin_to_core(cores, core);
                       serve(*bell);
                     }}.detach();
        ++workers_;
      }
#else
      for (const std::size_t core : worker_cores(0, hardware_threads())) {
        bells_.push_back(std::make_unique<worker_bell>(core));
        worker_bell* const bell {bells_.back().get()};
        std::thread {[this, bell] { serve(*bell); }}.detach();
        ++workers_;
      }
#endif
    }
 
    [[noreturn]] void serve(worker_bell& bell)
    {
      [[maybe_unused]] const std::size_t core {bell.core};
      std::size_t                        seen {0};
#if defined(__linux__)
      bool spin                               {true}; // the worker that sat a run out sleeps till the next: polling, it slowed its core's caller
#endif
      for (;;) {
        const auto started {[this, &seen] { return generation_.load() != seen; }};
#if defined(__linux__)
        if (spin) {
          spin_until(started);
        }
#else
        spin_until(started);
#endif
        if (!started()) { // asleep till a run that wants it rings: one that leaves it out never wakes it
          bell.asleep.store(true);
          std::unique_lock<std::mutex> held {bell.lock};
          bell.rung.wait(held, started);
          bell.asleep.store(false);
        }
        std::shared_ptr<pool_job> job;
        {
          const std::lock_guard<std::mutex> held {lock_};
          seen = generation_.load();
          job  = current_;
        }
#if defined(__linux__) // elsewhere no worker is pinned, so none sits out
        spin = !sits_out(core, job->core);
        if (!spin) {
          continue;
        }
#endif
#if defined(__linux__)
        // By rank among the cores, the caller's aside, not by arrival: a run of few threads then takes the same cores
        // from one run to the next, their caches still warm (ten AVX-512 cores, an LU of order 1024 on two threads:
        // 21.5 → 14.6 ms, on four 18 → 11.4; six AVX2 cores unchanged).
        const bool        wanted {wants(core, job->core, job->helpers)};
#else
        const bool wanted        {true};
#endif
        if (wanted && job->joined.fetch_add(1) < job->helpers) { // as many workers as the run's threads, the caller aside
          job->work();
        }
      }
    }
 
    std::mutex                                lock_;               
    std::vector<std::unique_ptr<worker_bell>> bells_;              
    std::shared_ptr<pool_job>                 current_;            
    std::atomic<std::size_t>                  generation_ {0};     
    std::atomic<bool>                         busy_       {false}; 
    std::size_t                               workers_    {0};     
#if defined(__linux__)
    std::vector<std::size_t> processor_cores_;                     
#endif
    long                      owner_      {current_process()};     
  };
 
  [[nodiscard]] inline std::atomic<std::size_t>& thread_cap() noexcept
  {
    static std::atomic<std::size_t> cap {0};
    return cap;
  }
 
  inline constexpr std::size_t work_per_thread {std::size_t {1} << 16};
 
  inline constexpr std::size_t elementwise_chunk {work_per_thread / 4};
 
  inline constexpr std::size_t vector_speedup {4};
 
  [[nodiscard]] inline std::size_t threads_for(std::size_t count,
                                               std::size_t work_per_index,
                                               std::size_t hardware)
  {
    const std::size_t total {work_per_index != 0 && count > std::numeric_limits<std::size_t>::max() / work_per_index
                             ? std::numeric_limits<std::size_t>::max()
                             : count * work_per_index};
    return std::min({total / work_per_thread, count, std::max<std::size_t>(1, hardware)});
  }
 
  template <typename Body>
  void parallel_for(std::size_t count,
                    std::size_t work_per_index,
                    const Body& body)
  {
    const std::size_t cap     {thread_cap().load()};
    const std::size_t threads {threads_for(count, work_per_index, cap == 0 ? hardware_threads() : std::min(cap, hardware_threads()))};
    const auto        runner  {[](const void* erased, std::size_t begin, std::size_t end) {
                                 const Body& work {*static_cast<const Body*>(erased)};
                                 for (std::size_t i {begin}; i < end; ++i) {
                                   work(i);
                                 }
                               }};
    if (threads < 2) {
      runner(&body, 0, count);
      return;
    }
    // Up to four chunks a thread, taken as threads free up, so a slower core holds back no one.
    const std::size_t chunks {std::max<std::size_t>(1, std::min(count, threads * 4))};
    const std::size_t chunk  {((count + chunks) - 1) / chunks}; // a last chunk past the count runs nothing
    if (!thread_pool::shared().try_run(runner, &body, count, chunk, chunks, threads - 1)) {
      runner(&body, 0, count);
    }
  }
 
  template <typename Body,
            typename Side>
  void parallel_for_beside(std::size_t count,
                           std::size_t work_per_index,
                           const Body& body,
                           const Side& side)
  {
    const std::size_t cap     {thread_cap().load()};
    const std::size_t threads {threads_for(count, work_per_index, cap == 0 ? hardware_threads() : std::min(cap, hardware_threads()))};
    const auto        runner  {[](const void* erased, std::size_t begin, std::size_t end) {
                                 const Body& work {*static_cast<const Body*>(erased)};
                                 for (std::size_t i {begin}; i < end; ++i) {
                                   work(i);
                                 }
                               }};
    const auto beside {[](const void* erased) { (*static_cast<const Side*>(erased))(); }};
    if (threads >= 2) {
      const std::size_t chunks {std::max<std::size_t>(1, std::min(count, threads * 4))};
      const std::size_t chunk  {((count + chunks) - 1) / chunks};
      if (thread_pool::shared().try_run(runner, &body, count, chunk, chunks, threads - 1, beside, &side)) {
        return;
      }
    }
    side();
    runner(&body, 0, count);
  }
 
  inline constexpr std::size_t zero_chunk {std::size_t {1} << 12};
 
  template <typename T>
  [[nodiscard]] constexpr std::size_t fault_chunk(std::size_t count) noexcept
  {
#if defined(__linux__)
    constexpr std::size_t huge_page  {std::size_t {1} << 21U};
    constexpr std::size_t huge_least {std::size_t {1} << 24U};
    if (count >= huge_least / sizeof(T)) {
      return huge_page / sizeof(T);
    }
#else
    (void)count;
#endif
    return zero_chunk;
  }
 
  inline constexpr std::size_t zero_work {1};
 
  template <typename T>
  void parallel_zero(T*          out,
                     std::size_t count)
  {
    if (count < 2 * (work_per_thread / zero_work)) { // too little to pay for a second thread
      std::fill(out, out + count, T {});
      return;
    }
    const std::size_t piece {fault_chunk<T>(count)};
    parallel_for((count + piece - 1) / piece, piece * zero_work, [out, count, piece](std::size_t chunk) {
                   const std::size_t begin {chunk * piece};
                   std::fill(out + begin, out + std::min(count, begin + piece), T {});
                 });
  }
 
  template <typename T>
  void parallel_copy(const T*    in,
                     std::size_t count,
                     T*          out)
  {
    if (count < 2 * (work_per_thread / zero_work)) { // too little to pay for a second thread
      std::copy(in, in + count, out);
      return;
    }
    const std::size_t piece {fault_chunk<T>(count)};
    parallel_for((count + piece - 1) / piece, piece * zero_work, [in, out, count, piece](std::size_t chunk) {
                   const std::size_t begin {chunk * piece};
                   std::copy(in + begin, in + std::min(count, begin + piece), out + begin);
                 });
  }
} // namespace numerical::detail
 
namespace numerical {
 
  inline void set_max_threads(std::size_t threads) noexcept
  {
    detail::thread_cap().store(threads);
  }
 
  [[nodiscard]] inline std::size_t max_threads()
  {
    const std::size_t cap      {detail::thread_cap().load()};
    const std::size_t hardware {detail::hardware_threads()};
    return cap == 0 ? hardware : std::min(cap, hardware);
  }
} // namespace numerical
 
#endif // NUMERICAL_PARALLEL_HPP