threads_pthread.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
 
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
 
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
 
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
#ifndef LIBMESH_THREADS_PTHREAD_H
#define LIBMESH_THREADS_PTHREAD_H
 
// Do not try to #include this header directly, it is designed to be
// #included directly by threads.h
#ifndef LIBMESH_SQUASH_HEADER_WARNING
# warning "This file is designed to be included through libmesh/threads.h"
#else
 
#ifdef LIBMESH_HAVE_PTHREAD
 
// C++ includes
#ifdef LIBMESH_HAVE_CXX11_THREAD
# include <thread>
#endif
 
#include <pthread.h>
#include <algorithm>
#include <vector>
#include <memory> // std::unique_ptr, std::make_unique
 
#ifdef __APPLE__
#  ifdef __MAC_10_12
#    include <os/lock.h>
#else
#    include <libkern/OSAtomic.h>
#  endif
#endif
 
// Thread-Local-Storage macros
#ifdef LIBMESH_HAVE_CXX11_THREAD
#  define LIBMESH_TLS_TYPE(type)  thread_local type
#  define LIBMESH_TLS_REF(value)  (value)
#else // Maybe support gcc __thread eventually?
#  define LIBMESH_TLS_TYPE(type)  type
#  define LIBMESH_TLS_REF(value)  (value)
#endif
 
namespace libMesh
{
 
namespace Threads
{
 
 
#ifdef LIBMESH_HAVE_CXX11_THREAD
typedef std::thread Thread;
 
#else
 
typedef NonConcurrentThread Thread;
 
#endif // LIBMESH_HAVE_CXX11_THREAD
 
 
#ifdef __APPLE__
#ifdef __MAC_10_12
class spin_mutex
{
public:
  spin_mutex() { ulock = OS_UNFAIR_LOCK_INIT; }
  ~spin_mutex() = default;
 
  void lock () { os_unfair_lock_lock(&ulock); }
  void unlock () { os_unfair_lock_unlock(&ulock); }
 
  class scoped_lock
  {
  public:
    scoped_lock () : smutex(nullptr) {}
    explicit scoped_lock ( spin_mutex & in_smutex ) : smutex(&in_smutex) { smutex->lock(); }
 
    ~scoped_lock () { release(); }
 
    void acquire ( spin_mutex & in_smutex ) { smutex = &in_smutex; smutex->lock(); }
    void release () { if (smutex) smutex->unlock(); smutex = nullptr; }
 
  private:
    spin_mutex * smutex;
  };
 
private:
  os_unfair_lock ulock;
};
#else
class spin_mutex
{
public:
  spin_mutex() : slock(0) {} // The convention is that the lock being zero is _unlocked_
  ~spin_mutex() = default;
 
  void lock () { OSSpinLockLock(&slock); }
  void unlock () { OSSpinLockUnlock(&slock); }
 
  class scoped_lock
  {
  public:
    scoped_lock () : smutex(nullptr) {}
    explicit scoped_lock ( spin_mutex & in_smutex ) : smutex(&in_smutex) { smutex->lock(); }
 
    ~scoped_lock () { release(); }
 
    void acquire ( spin_mutex & in_smutex ) { smutex = &in_smutex; smutex->lock(); }
    void release () { if (smutex) smutex->unlock(); smutex = nullptr; }
 
  private:
    spin_mutex * smutex;
  };
 
private:
  OSSpinLock slock;
};
#endif
#else
class spin_mutex
{
public:
  // Might want to use PTHREAD_MUTEX_ADAPTIVE_NP on Linux, but it's not available on OSX.
  spin_mutex() { pthread_spin_init(&slock, PTHREAD_PROCESS_PRIVATE); }
  ~spin_mutex() { pthread_spin_destroy(&slock); }
 
  void lock () { pthread_spin_lock(&slock); }
  void unlock () { pthread_spin_unlock(&slock); }
 
  class scoped_lock
  {
  public:
    scoped_lock () : smutex(nullptr) {}
    explicit scoped_lock ( spin_mutex & in_smutex ) : smutex(&in_smutex) { smutex->lock(); }
 
    ~scoped_lock () { release(); }
 
    void acquire ( spin_mutex & in_smutex ) { smutex = &in_smutex; smutex->lock(); }
    void release () { if (smutex) smutex->unlock(); smutex = nullptr; }
 
  private:
    spin_mutex * smutex;
  };
 
private:
  pthread_spinlock_t slock;
};
#endif // __APPLE__
 
 
 
class recursive_mutex
{
public:
  // Might want to use PTHREAD_MUTEX_ADAPTIVE_NP on Linux, but it's not available on OSX.
  recursive_mutex()
  {
    pthread_mutexattr_init(&attr);
    pthread_mutexattr_settype(&attr, PTHREAD_MUTEX_RECURSIVE);
    pthread_mutex_init(&mutex, &attr);
  }
  ~recursive_mutex() { pthread_mutex_destroy(&mutex); }
 
  void lock () { pthread_mutex_lock(&mutex); }
  void unlock () { pthread_mutex_unlock(&mutex); }
 
  class scoped_lock
  {
  public:
    scoped_lock () : rmutex(nullptr) {}
    explicit scoped_lock ( recursive_mutex & in_rmutex ) : rmutex(&in_rmutex) { rmutex->lock(); }
 
    ~scoped_lock () { release(); }
 
    void acquire ( recursive_mutex & in_rmutex ) { rmutex = &in_rmutex; rmutex->lock(); }
    void release () { if (rmutex) rmutex->unlock(); rmutex = nullptr; }
 
  private:
    recursive_mutex * rmutex;
  };
 
private:
  pthread_mutex_t mutex;
  pthread_mutexattr_t attr;
};
 
template <typename Range>
unsigned int num_pthreads(const Range & range,
                          unsigned int requested = libMesh::n_threads())
{
  std::size_t mn = std::min((std::size_t)requested, range.size()/range.grainsize());
  return mn > 0 ? cast_int<unsigned int>(mn) : 1;
}
 
template <typename Range, typename Body>
class RangeBody
{
public:
  Range * range;
  Body * body;
};
 
template <typename Range, typename Body>
void * run_body(void * args)
{
  RangeBody<Range, Body> * range_body = (RangeBody<Range, Body> *)args;
 
  Body & body = *range_body->body;
  Range & range = *range_body->range;
 
  body(range);
 
  return nullptr;
}
 
class task_scheduler_init
{
public:
  static const int automatic = -1;
  explicit task_scheduler_init (int = automatic) {}
  void initialize (int = automatic) {}
  void terminate () {}
};
 
//-------------------------------------------------------------------
class split {};
 
 
 
 
//-------------------------------------------------------------------
template <typename Range, typename Body>
inline
void parallel_for (const Range & range, const Body & body,
                   unsigned int n_threads = libMesh::n_threads())
{
  libmesh_error_msg_if(n_threads > libMesh::n_threads(),
                       "Requested n_threads (" << n_threads << ") exceeds the "
                       "global thread count (" << libMesh::n_threads() << ").");
  Threads::BoolAcquire set_in_threads(Threads::in_threads);
 
  unsigned int actual_threads = num_pthreads(range, n_threads);
 
  // If we're running in serial - just run!
  if (actual_threads == 1)
  {
    body(range);
    return;
  }
 
  Threads::RAIIAcquire<int> set_active_threads(Threads::active_threads, actual_threads);
 
  DisablePerfLogInScope disable_perf;
 
  std::vector<std::unique_ptr<Range>> ranges(actual_threads);
  std::vector<RangeBody<const Range, const Body>> range_bodies(actual_threads);
  std::vector<pthread_t> threads(actual_threads);
 
  // Create the ranges for each thread
  std::size_t range_size = range.size() / actual_threads;
 
  typename Range::const_iterator current_beginning = range.begin();
 
  for (unsigned int i=0; i<actual_threads; i++)
    {
      std::size_t this_range_size = range_size;
 
      if (i+1 == actual_threads)
        this_range_size += range.size() % actual_threads; // Give the last one the remaining work to do
 
      ranges[i] = std::make_unique<Range>(range, current_beginning, current_beginning + this_range_size);
 
      current_beginning = current_beginning + this_range_size;
    }
 
  // Create the RangeBody arguments
  for (unsigned int i=0; i<actual_threads; i++)
    {
      range_bodies[i].range = ranges[i].get();
      range_bodies[i].body = &body;
    }
 
  // Create the threads.  It may seem redundant to wrap a pragma in
  // #ifdefs... but GCC warns about an "unknown pragma" if it
  // encounters this line of code when -fopenmp is not passed to the
  // compiler.
#ifdef LIBMESH_HAVE_OPENMP
// Cap the OpenMP team to the actual number of threads this dispatch uses. Without this clause the
// team defaults to the process-wide libMesh::n_threads(), so omp_get_thread_num()
// could exceed a reduced per-application thread count and index past
// per-thread storage sized to that count.
#pragma omp parallel for schedule (static) num_threads(actual_threads)
#endif
  for (int i=0; i<static_cast<int>(actual_threads); i++)
    {
#if !LIBMESH_HAVE_OPENMP
      pthread_create(&threads[i], nullptr, &run_body<Range, Body>, (void *)&range_bodies[i]);
#else
      run_body<Range, Body>((void *)&range_bodies[i]);
#endif
    }
 
#if !LIBMESH_HAVE_OPENMP
  // Wait for them to finish
 
  // The use of 'int' instead of unsigned for the iteration variable
  // is deliberate here.  This is an OpenMP loop, and some older
  // compilers warn when you don't use int for the loop index.  The
  // reason has to do with signed vs. unsigned integer overflow
  // behavior and optimization.
  // http://blog.llvm.org/2011/05/what-every-c-programmer-should-know.html
  for (int i=0; i<static_cast<int>(actual_threads); i++)
    pthread_join(threads[i], nullptr);
#endif
}
 
template <typename Range, typename Body, typename Partitioner>
inline
void parallel_for (const Range & range, const Body & body, const Partitioner &,
                   unsigned int n_threads = libMesh::n_threads())
{
  parallel_for (range, body, n_threads);
}
 
template <typename Range, typename Body>
inline
void parallel_reduce (const Range & range, Body & body,
                      unsigned int n_threads = libMesh::n_threads())
{
  libmesh_error_msg_if(n_threads > libMesh::n_threads(),
                       "Requested n_threads (" << n_threads << ") exceeds the "
                       "global thread count (" << libMesh::n_threads() << ").");
  Threads::BoolAcquire set_in_threads(Threads::in_threads);
 
  unsigned int actual_threads = num_pthreads(range, n_threads);
 
  // If we're running in serial - just run!
  if (actual_threads == 1)
  {
    body(range);
    return;
  }
 
  Threads::RAIIAcquire<int> set_active_threads(Threads::active_threads, actual_threads);
 
  DisablePerfLogInScope disable_perf;
 
  std::vector<std::unique_ptr<Range>> ranges(actual_threads);
  std::vector<std::unique_ptr<Body>> managed_bodies(actual_threads); // bodies we are responsible for
  std::vector<Body *> bodies(actual_threads); // dumb pointers to managed_bodies
  std::vector<RangeBody<Range, Body>> range_bodies(actual_threads);
 
  // Create actual_threads-1 copies of "body". We manage the lifetime of
  // these copies with std::unique_ptrs.
  for (unsigned int i=1; i<actual_threads; i++)
    managed_bodies[i] = std::make_unique<Body>(body, Threads::split());
 
  // Set up the "bodies" vector. Use the passed in body for the first
  // one, point to managed_bodies entries for the others.
  bodies[0] = &body;
  for (unsigned int i=1; i<actual_threads; i++)
    bodies[i] = managed_bodies[i].get();
 
  // Create the ranges for each thread
  std::size_t range_size = range.size() / actual_threads;
 
  typename Range::const_iterator current_beginning = range.begin();
 
  for (unsigned int i=0; i<actual_threads; i++)
    {
      std::size_t this_range_size = range_size;
 
      if (i+1 == actual_threads)
        this_range_size += range.size() % actual_threads; // Give the last one the remaining work to do
 
      ranges[i] = std::make_unique<Range>(range, current_beginning, current_beginning + this_range_size);
 
      current_beginning = current_beginning + this_range_size;
    }
 
  // Create the RangeBody arguments
  for (unsigned int i=0; i<actual_threads; i++)
    {
      range_bodies[i].range = ranges[i].get();
      range_bodies[i].body = bodies[i];
    }
 
  // Create the threads
  std::vector<pthread_t> threads(actual_threads);
 
  // It may seem redundant to wrap a pragma in #ifdefs... but GCC
  // warns about an "unknown pragma" if it encounters this line of
  // code when -fopenmp is not passed to the compiler.
#ifdef LIBMESH_HAVE_OPENMP
// Cap the OpenMP team to the actual number of threads this dispatch uses. Without this clause the
// team defaults to the process-wide libMesh::n_threads(), so omp_get_thread_num()
// could exceed a reduced per-application thread count and index past
// per-thread storage sized to that count.
#pragma omp parallel for schedule (static) num_threads(actual_threads)
#endif
  // The use of 'int' instead of unsigned for the iteration variable
  // is deliberate here.  This is an OpenMP loop, and some older
  // compilers warn when you don't use int for the loop index.  The
  // reason has to do with signed vs. unsigned integer overflow
  // behavior and optimization.
  // http://blog.llvm.org/2011/05/what-every-c-programmer-should-know.html
  for (int i=0; i<static_cast<int>(actual_threads); i++)
    {
#if !LIBMESH_HAVE_OPENMP
      pthread_create(&threads[i], nullptr, &run_body<Range, Body>, (void *)&range_bodies[i]);
#else
      run_body<Range, Body>((void *)&range_bodies[i]);
#endif
    }
 
#if !LIBMESH_HAVE_OPENMP
  // Wait for them to finish
  for (unsigned int i=0; i<actual_threads; i++)
    pthread_join(threads[i], nullptr);
#endif
 
  // Join them all down to the original Body
  for (unsigned int i=actual_threads-1; i != 0; i--)
  {
    bodies[i-1]->join(*bodies[i]);
 
    // Free data ASAP lest we have O(actual_threads^2) at the end
    managed_bodies[i].reset();
  }
}
 
template <typename Range, typename Body, typename Partitioner>
inline
void parallel_reduce (const Range & range, Body & body, const Partitioner &,
                      unsigned int n_threads = libMesh::n_threads())
{
  parallel_reduce(range, body, n_threads);
}
 
 
template <typename T>
class atomic
{
public:
  atomic () : val(0) {}
  operator T () { return val; }
 
  T operator=( T value )
  {
    spin_mutex::scoped_lock lock(smutex);
    val = value;
    return val;
  }
 
  atomic<T> & operator=( const atomic<T> & value )
  {
    spin_mutex::scoped_lock lock(smutex);
    val = value;
    return *this;
  }
 
 
  T operator+=(T value)
  {
    spin_mutex::scoped_lock lock(smutex);
    val += value;
    return val;
  }
 
  T operator-=(T value)
  {
    spin_mutex::scoped_lock lock(smutex);
    val -= value;
    return val;
  }
 
  T operator++()
  {
    spin_mutex::scoped_lock lock(smutex);
    val++;
    return val;
  }
 
  T operator++(int)
  {
    spin_mutex::scoped_lock lock(smutex);
    val++;
    return val;
  }
 
  T operator--()
  {
    spin_mutex::scoped_lock lock(smutex);
    val--;
    return val;
  }
 
  T operator--(int)
  {
    spin_mutex::scoped_lock lock(smutex);
    val--;
    return val;
  }
 
private:
  T val;
  spin_mutex smutex;
};
 
} // namespace Threads
 
} // namespace libMesh
 
#endif // #ifdef LIBMESH_HAVE_PTHREAD
 
#endif // LIBMESH_SQUASH_HEADER_WARNING
 
#endif // LIBMESH_THREADS_PTHREAD_H