threads.h

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2014 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_H
#define LIBMESH_THREADS_H

// Local includes
#include "libmesh/libmesh_config.h"
#include "libmesh/libmesh_common.h"  // for libmesh_assert

// Threading building blocks includes
#ifdef LIBMESH_HAVE_TBB_API
#  include "libmesh/libmesh_logging.h" // only mess with the perflog if we are really multithreaded
#  include "tbb/tbb_stddef.h"
#  include "tbb/blocked_range.h"
#  include "tbb/parallel_for.h"
#  include "tbb/parallel_reduce.h"
#  include "tbb/task_scheduler_init.h"
#  include "tbb/partitioner.h"
#  include "tbb/spin_mutex.h"
#  include "tbb/recursive_mutex.h"
//  include "boost/thread/recursive_mutex.hpp"
#  include "tbb/atomic.h"
#endif

// C++ includes
#ifdef LIBMESH_HAVE_STD_THREAD
#  include <thread>
#elif LIBMESH_HAVE_TBB_CXX_THREAD
#  include "tbb/tbb_thread.h"
#endif

#ifdef LIBMESH_HAVE_PTHREAD
#  include "libmesh/libmesh_logging.h" // only mess with the perflog if we are really multithreaded
#  include <pthread.h>
#  include <algorithm>
#  include <vector>

#ifdef __APPLE__
#include <libkern/OSAtomic.h>
#endif

#endif

// Thread-Local-Storage macros

#ifdef LIBMESH_HAVE_STD_THREAD
#  define LIBMESH_TLS_TYPE(type)  thread_local type
#  define LIBMESH_TLS_REF(value)  (value)
#else
#  ifdef LIBMESH_HAVE_TBB_API
#    include "tbb/enumerable_thread_specific.h"
#    define LIBMESH_TLS_TYPE(type)  tbb::enumerable_thread_specific<type>
#    define LIBMESH_TLS_REF(value)  (value).local()
#  else // Maybe support gcc __thread eventually?
#    define LIBMESH_TLS_TYPE(type)  type
#    define LIBMESH_TLS_REF(value)  (value)
#  endif
#endif



namespace libMesh
{


namespace Threads
{
extern bool in_threads;

class BoolAcquire {
public:
  explicit
  BoolAcquire(bool& b) : _b(b) { libmesh_assert(!_b); _b = true; }

  ~BoolAcquire() { libmesh_exceptionless_assert(_b); _b = false; }
private:
  bool& _b;
};


#ifdef LIBMESH_HAVE_STD_THREAD
//--------------------------------------------------------------------
typedef std::thread Thread;

#elif LIBMESH_HAVE_TBB_CXX_THREAD
//--------------------------------------------------------------------
typedef tbb::tbb_thread Thread;

#else
//--------------------------------------------------------------------
class Thread
{
public:
  template <typename Callable>
  Thread (Callable f) { f(); }

  void join() {}

  bool joinable() const { return true; }
};

#endif



#ifdef LIBMESH_HAVE_TBB_API
//-------------------------------------------------------------------
typedef tbb::task_scheduler_init task_scheduler_init;



//-------------------------------------------------------------------
typedef tbb::split split;



//-------------------------------------------------------------------
template <typename Range, typename Body>
inline
void parallel_for (const Range &range, const Body &body)
{
  BoolAcquire b(in_threads);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  const bool logging_was_enabled = libMesh::perflog.logging_enabled();

  if (libMesh::n_threads() > 1)
    libMesh::perflog.disable_logging();
#endif

  if (libMesh::n_threads() > 1)
    tbb::parallel_for (range, body, tbb::auto_partitioner());

  else
    body(range);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  if (libMesh::n_threads() > 1 && logging_was_enabled)
    libMesh::perflog.enable_logging();
#endif
}



//-------------------------------------------------------------------
template <typename Range, typename Body, typename Partitioner>
inline
void parallel_for (const Range &range, const Body &body, const Partitioner &partitioner)
{
  BoolAcquire b(in_threads);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  const bool logging_was_enabled = libMesh::perflog.logging_enabled();

  if (libMesh::n_threads() > 1)
    libMesh::perflog.disable_logging();
#endif

  if (libMesh::n_threads() > 1)
    tbb::parallel_for (range, body, partitioner);

  else
    body(range);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  if (libMesh::n_threads() > 1 && logging_was_enabled)
    libMesh::perflog.enable_logging();
#endif
}



//-------------------------------------------------------------------
template <typename Range, typename Body>
inline
void parallel_reduce (const Range &range, Body &body)
{
  BoolAcquire b(in_threads);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  const bool logging_was_enabled = libMesh::perflog.logging_enabled();

  if (libMesh::n_threads() > 1)
    libMesh::perflog.disable_logging();
#endif

  if (libMesh::n_threads() > 1)
    tbb::parallel_reduce (range, body, tbb::auto_partitioner());

  else
    body(range);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  if (libMesh::n_threads() > 1 && logging_was_enabled)
    libMesh::perflog.enable_logging();
#endif
}



//-------------------------------------------------------------------
template <typename Range, typename Body, typename Partitioner>
inline
void parallel_reduce (const Range &range, Body &body, const Partitioner &partitioner)
{
  BoolAcquire b(in_threads);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  const bool logging_was_enabled = libMesh::perflog.logging_enabled();

  if (libMesh::n_threads() > 1)
    libMesh::perflog.disable_logging();
#endif

  if (libMesh::n_threads() > 1)
    tbb::parallel_reduce (range, body, partitioner);

  else
    body(range);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  if (libMesh::n_threads() > 1 && logging_was_enabled)
    libMesh::perflog.enable_logging();
#endif
}



//-------------------------------------------------------------------
typedef tbb::spin_mutex spin_mutex;

//-------------------------------------------------------------------
typedef tbb::recursive_mutex recursive_mutex;

//-------------------------------------------------------------------
template <typename T>
class atomic : public tbb::atomic<T> {};

#else //LIBMESH_HAVE_TBB_API
#ifdef LIBMESH_HAVE_PTHREAD

//-------------------------------------------------------------------
#ifdef __APPLE__
class spin_mutex
{
public:
  spin_mutex() : slock(0) {} // The convention is that the lock being zero is _unlocked_
  ~spin_mutex() {}

  void lock () { OSSpinLockLock(&slock); }
  void unlock () { OSSpinLockUnlock(&slock); }

  class scoped_lock
  {
  public:
    scoped_lock () : smutex(NULL) {}
    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 = NULL; }

  private:
    spin_mutex * smutex;
  };

private:
  OSSpinLock slock;
};
#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(NULL) {}
    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 = NULL; }

  private:
    spin_mutex * smutex;
  };

private:
  pthread_spinlock_t slock;
};
#endif

//-------------------------------------------------------------------
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(NULL) {}
    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 = NULL; }

  private:
    recursive_mutex * rmutex;
  };

private:
  pthread_mutex_t mutex;
  pthread_mutexattr_t attr;
};

extern std::map<pthread_t, unsigned int> _pthread_unique_ids;
extern spin_mutex _pthread_unique_id_mutex;

unsigned int pthread_unique_id();

template <typename Range>
unsigned int num_pthreads(Range & range)
{
  unsigned int min = std::min((std::size_t)libMesh::n_threads(), range.size());
  return min > 0 ? min : 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 NULL;
}

//-------------------------------------------------------------------
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)
{
  Threads::BoolAcquire b(Threads::in_threads);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  const bool logging_was_enabled = libMesh::perflog.logging_enabled();

  if (libMesh::n_threads() > 1)
    libMesh::perflog.disable_logging();
#endif

  unsigned int n_threads = num_pthreads(range);

  std::vector<Range *> ranges(n_threads);
  std::vector<RangeBody<const Range, const Body> > range_bodies(n_threads);
  std::vector<pthread_t> threads(n_threads);

  // Create the ranges for each thread
  unsigned int range_size = range.size() / n_threads;

  typename Range::const_iterator current_beginning = range.begin();

  for(unsigned int i=0; i<n_threads; i++)
    {
      unsigned int this_range_size = range_size;

      if(i+1 == n_threads)
        this_range_size += range.size() % n_threads; // Give the last one the remaining work to do

      ranges[i] = new 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<n_threads; i++)
    {
      range_bodies[i].range = ranges[i];
      range_bodies[i].body = &body;
    }

  // Create the threads
#pragma omp parallel for schedule (static)
  for(unsigned int i=0; i<n_threads; i++)
    {
#if LIBMESH_HAVE_OPENMP
      run_body<Range, Body>((void*)&range_bodies[i]);
#else // Just use Pthreads
      spin_mutex::scoped_lock lock(_pthread_unique_id_mutex);
      pthread_create(&threads[i], NULL, &run_body<Range, Body>, (void*)&range_bodies[i]);
      _pthread_unique_ids[threads[i]] = 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>(n_threads); i++)
    {
      pthread_join(threads[i], NULL);
      spin_mutex::scoped_lock lock(_pthread_unique_id_mutex);
      _pthread_unique_ids.erase(threads[i]);
    }
#endif

  // Clean up
  for(unsigned int i=0; i<n_threads; i++)
    delete ranges[i];

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  if (libMesh::n_threads() > 1 && logging_was_enabled)
    libMesh::perflog.enable_logging();
#endif
}

//-------------------------------------------------------------------
template <typename Range, typename Body, typename Partitioner>
inline
void parallel_for (const Range &range, const Body &body, const Partitioner &)
{
  parallel_for(range, body);
}

//-------------------------------------------------------------------
template <typename Range, typename Body>
inline
void parallel_reduce (const Range &range, Body &body)
{
  Threads::BoolAcquire b(Threads::in_threads);

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  const bool logging_was_enabled = libMesh::perflog.logging_enabled();

  if (libMesh::n_threads() > 1)
    libMesh::perflog.disable_logging();
#endif

  unsigned int n_threads = num_pthreads(range);

  std::vector<Range *> ranges(n_threads);
  std::vector<Body *> bodies(n_threads);
  std::vector<RangeBody<Range, Body> > range_bodies(n_threads);

  // Create copies of the body for each thread
  bodies[0] = &body; // Use the original body for the first one
  for(unsigned int i=1; i<n_threads; i++)
    bodies[i] = new Body(body, Threads::split());

  // Create the ranges for each thread
  unsigned int range_size = range.size() / n_threads;

  typename Range::const_iterator current_beginning = range.begin();

  for(unsigned int i=0; i<n_threads; i++)
    {
      unsigned int this_range_size = range_size;

      if(i+1 == n_threads)
        this_range_size += range.size() % n_threads; // Give the last one the remaining work to do

      ranges[i] = new 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<n_threads; i++)
    {
      range_bodies[i].range = ranges[i];
      range_bodies[i].body = bodies[i];
    }

  // Create the threads
  std::vector<pthread_t> threads(n_threads);

#pragma omp parallel for schedule (static)
  // 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>(n_threads); i++)
    {
#if LIBMESH_HAVE_OPENMP
      run_body<Range, Body>((void*)&range_bodies[i]);
#else // Just use Pthreads
      spin_mutex::scoped_lock lock(_pthread_unique_id_mutex);
      pthread_create(&threads[i], NULL, &run_body<Range, Body>, (void*)&range_bodies[i]);
      _pthread_unique_ids[threads[i]] = i;
#endif
    }

#if !LIBMESH_HAVE_OPENMP
  // Wait for them to finish
  for(unsigned int i=0; i<n_threads; i++)
    {
      pthread_join(threads[i], NULL);
      spin_mutex::scoped_lock lock(_pthread_unique_id_mutex);
      _pthread_unique_ids.erase(threads[i]);
    }
#endif

  // Join them all down to the original Body
  for(unsigned int i=n_threads-1; i != 0; i--)
    bodies[i-1]->join(*bodies[i]);

  // Clean up
  for(unsigned int i=1; i<n_threads; i++)
    delete bodies[i];
  for(unsigned int i=0; i<n_threads; i++)
    delete ranges[i];

#ifdef LIBMESH_ENABLE_PERFORMANCE_LOGGING
  if (libMesh::n_threads() > 1 && logging_was_enabled)
    libMesh::perflog.enable_logging();
#endif
}

//-------------------------------------------------------------------
template <typename Range, typename Body, typename Partitioner>
inline
void parallel_reduce (const Range &range, Body &body, const Partitioner &)
{
  parallel_reduce(range, body);
}


//-------------------------------------------------------------------
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;
};

#else //LIBMESH_HAVE_PTHREAD

//-------------------------------------------------------------------
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)
{
  BoolAcquire b(in_threads);
  body(range);
}

//-------------------------------------------------------------------
template <typename Range, typename Body, typename Partitioner>
inline
void parallel_for (const Range &range, const Body &body, const Partitioner &)
{
  BoolAcquire b(in_threads);
  body(range);
}

//-------------------------------------------------------------------
template <typename Range, typename Body>
inline
void parallel_reduce (const Range &range, Body &body)
{
  BoolAcquire b(in_threads);
  body(range);
}

//-------------------------------------------------------------------
template <typename Range, typename Body, typename Partitioner>
inline
void parallel_reduce (const Range &range, Body &body, const Partitioner &)
{
  BoolAcquire b(in_threads);
  body(range);
}

//-------------------------------------------------------------------
class spin_mutex
{
public:
  spin_mutex() {}
  void lock () {}
  void unlock () {}

  class scoped_lock
  {
  public:
    scoped_lock () {}
    explicit scoped_lock ( spin_mutex&  ) {}
    void acquire ( spin_mutex& ) {}
    void release () {}
  };
};

//-------------------------------------------------------------------
class recursive_mutex
{
public:
  recursive_mutex() {}

  class scoped_lock
  {
  public:
    scoped_lock () {}
    explicit scoped_lock ( recursive_mutex&  ) {}
    void acquire ( recursive_mutex& ) {}
    void release () {}
  };
};

//-------------------------------------------------------------------
template <typename T>
class atomic
{
public:
  atomic () : _val(0) {}
  operator T& () { return _val; }
private:
  T _val;
};

#endif // LIBMESH_HAVE_PTHREAD
#endif // #ifdef LIBMESH_HAVE_TBB_API



template <typename T>
class BlockedRange
{
public:
  typedef T const_iterator;

  explicit BlockedRange (const unsigned int new_grainsize = 1000) :
    _grainsize(new_grainsize)
  {}

  BlockedRange (const const_iterator first,
                const const_iterator last,
                const unsigned int new_grainsize = 1000) :
    _grainsize(new_grainsize)
  {
    this->reset(first, last);
  }

  BlockedRange (const BlockedRange<T> &r):
    _end(r._end),
    _begin(r._begin),
    _grainsize(r._grainsize)
  {}

  BlockedRange (BlockedRange<T> &r, Threads::split ) :
    _end(r._end),
    _begin(r._begin),
    _grainsize(r._grainsize)
  {
    const_iterator
      beginning = r._begin,
      ending    = r._end,
      middle    = beginning + (ending - beginning)/2u;

    r._end = _begin = middle;
  }

  void reset (const const_iterator first,
              const const_iterator last)
  {
    _begin = first;
    _end   = last;
  }

  const_iterator begin () const { return _begin; }

  const_iterator end () const { return _end; }

  unsigned int grainsize () const {return _grainsize;}

  void grainsize (const unsigned int &gs) {_grainsize = gs;}

  int size () const { return (_end -_begin); }

  //------------------------------------------------------------------------
  // Methods that implement Range concept
  //------------------------------------------------------------------------

  bool empty() const { return (_begin == _end); }

  bool is_divisible() const { return ((_begin + this->grainsize()) < _end); }

private:

  const_iterator _end;
  const_iterator _begin;
  unsigned int _grainsize;
};



extern spin_mutex spin_mtx;

extern recursive_mutex recursive_mtx;

} // namespace Threads

} // namespace libMesh

#endif // LIBMESH_THREADS_H