NetThreads.h

#pragma once
 
#include <algorithm>
#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);
#ifdef _OPENMP
#   pragma omp parallel num_threads(chunks)
#endif
   {
#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
      for (int chunk = first; chunk < chunks; chunk += stride)
         updateNodes(ChunkBegin(numNodes, chunk, chunks), ChunkBegin(numNodes, chunk + 1, chunks));
#ifdef _OPENMP
#   pragma omp barrier
#endif
      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);
      }
   }
   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