NetThreads.h

#pragma once
 
#include <algorithm>
#include <exception>
#include <fhsim/simobject/ISimObjectCreator.h>
#include <stdexcept>
#include <string>
#include <vector>
 
#ifdef _OPENMP
#   include <omp.h>
#endif
 
namespace net_threads
{
inline int ReadNumThreads(ISimObjectCreator* creator)
{
   int numThreads = 1;
   creator->GetIntParam("NumThreads", &numThreads, numThreads);
   if (numThreads < 0) {
      const std::string message = "NumThreads = " + std::to_string(numThreads) + " is negative; expected 0 (OpenMP's default), 1 (serial) or more.";
      creator->ReportParameterError("NumThreads", message);
      throw std::invalid_argument(message);
   }
#ifdef _OPENMP
   if (numThreads == 0)
      numThreads = omp_get_max_threads();
#else
   numThreads = 1; // built without OpenMP
#endif
   return std::max(numThreads, 1);
}
 
inline int ChunkBegin(const int count, const int chunk, const int chunks)
{
   return count * chunk / chunks;
}
 
template <class UpdateNodes, class AddForces>
void SumForcesInChunks(const int chunks, const int numNodes, std::vector<double>& buffers, double* const nodeForce, UpdateNodes updateNodes, AddForces addForces)
{
   const size_t       bufferSize = 3 * static_cast<size_t>(numNodes);
   std::exception_ptr failure; // The first exception any thread met (MARE-0422).
#ifdef _OPENMP
#   pragma omp parallel num_threads(chunks)
#endif
   {
      // An exception that leaves an OpenMP region ends the process with std::terminate and no message. Each thread keeps
      // the first one and still reaches the barrier; it is rethrown on the calling thread after the region.
      const auto keepFailure = [&failure]() {
#ifdef _OPENMP
#   pragma omp critical(net_threads_failure)
#endif
         {
            if (!failure)
               failure = std::current_exception();
         }
      };
#ifdef _OPENMP
      const int first  = omp_get_thread_num();
      const int stride = omp_get_num_threads();
#else
      const int first  = 0;
      const int stride = 1;
#endif
      try {
         for (int chunk = first; chunk < chunks; chunk += stride)
            updateNodes(ChunkBegin(numNodes, chunk, chunks), ChunkBegin(numNodes, chunk + 1, chunks));
      } catch (...) {
         keepFailure();
      }
#ifdef _OPENMP
#   pragma omp barrier
#endif
      try {
         for (int chunk = first; chunk < chunks; chunk += stride) {
            double* const forces = buffers.data() + chunk * bufferSize;
            std::fill(forces, forces + bufferSize, 0.0);
            addForces(chunk, chunks, forces);
         }
      } catch (...) {
         keepFailure();
      }
   }
   if (failure)
      std::rethrow_exception(failure);
   for (int chunk = 0; chunk < chunks; chunk++) {
      const double* const forces = buffers.data() + chunk * bufferSize;
      for (size_t i = 0; i < bufferSize; i++)
         nodeForce[i] += forces[i];
   }
}
} // namespace net_threads