<commit_before>// $Id$

// The libMesh Finite Element Library.
// Copyright (C) 2002-2008 <NAME>, <NAME>, <NAME>
  
// 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 __parallel_h__
#define __parallel_h__

// System includes
#include <string>
#include <vector>

// Local includes
#include "libmesh_common.h" // for Real
#include "libmesh_logging.h"


// Macro to identify and debug functions which should only be called in
// parallel on every processor at once

#undef parallel_only
#ifndef NDEBUG
  #define parallel_only() do { libmesh_assert(Parallel::verify(std::string(__FILE__))); libmesh_assert(Parallel::verify(__LINE__)); } while (0)
#else
  #define parallel_only()
#endif

/**
 * The Parallel namespace is for wrapper functions
 * for common general parallel synchronization tasks.
 *
 * For MPI 1.1 compatibility, temporary buffers are used
 * instead of MPI 2's MPI_IN_PLACE
 */
namespace Parallel
{
  //-------------------------------------------------------------------
  /**
   * Forward declarations of classes we will define later.
   */
  class DataType;
  class Request;
  class Status;

#ifdef LIBMESH_HAVE_MPI
  //-------------------------------------------------------------------
  /**
   * Data types for communication
   */
  typedef MPI_Datatype data_type;

  /**
   * Request object for non-blocking I/O
   */
  typedef MPI_Request request;

  /**
   * Status object for querying messages
   */
  typedef MPI_Status status;
 
  /**
   * Templated function to return the appropriate MPI datatype
   * for use with built-in C types
   */
  template <typename T>
  inline data_type datatype();

  /**
   * Default message tag id
   */
  const int any_tag=MPI_ANY_TAG;

  /**
   * Accept from any source
   */
  const int any_source=MPI_ANY_SOURCE;
  
#else

  // These shouldn't actually be needed, but must be 
  // unique types for function overloading to work 
  // properly.
  struct data_type { /* unsigned int t; */ };
  struct request   { /* unsigned int r; */ };
  struct status    { /* unsigned int s; */ };

  template <typename T>
  inline data_type datatype() { return data_type(); }

  const int any_tag=-1;
  const int any_source=0;
#endif // LIBMESH_HAVE_MPI



  //-------------------------------------------------------------------
  /**
   * Encapsulates the MPI_Datatype.
   */
  class DataType
  {
  public:
    DataType () {}

    DataType (const DataType &other) :
      _datatype(other._datatype)
    {}

    DataType (const data_type &type) :
      _datatype(type)
    {}

    DataType & operator = (const DataType &other) 
    { _datatype = other._datatype; return *this; }

    DataType & operator = (const data_type &type)
    { _datatype = type; return *this; }

    operator const data_type & () const
    { return _datatype; }

    operator data_type & ()
    { return _datatype; }
    
//     operator data_type const * () const
//     { return &_datatype; }

//     operator data_type * ()
//     { return &_datatype; }

    void commit ()
    {
#ifdef LIBMESH_HAVE_MPI
      MPI_Type_commit (&_datatype);
#endif
    }

    void free ()
    {
#ifdef LIBMESH_HAVE_MPI
      MPI_Type_free (&_datatype);
#endif
    }
    
  private:

    data_type _datatype;    
  };



  //-------------------------------------------------------------------
  /**
   * Encapsulates the MPI_Status struct.  Allows the source and size 
   * of the message to be determined.
   */
  class Status
  {
  public:
    Status () {}
    
    Status (const data_type &type) :
    _datatype(type)           
    {}

    Status (const status &status) :
      _status(status)
    {}

    Status (const status    &status,
	    const data_type &type) :
      _status(status),
      _datatype(type)           
    {}

    Status (const Status &status) :
      _status(status._status),
      _datatype(status._datatype)
    {}

    Status (const Status    &status,
	    const data_type &type) :
      _status(status._status),
      _datatype(type)
    {}
    
    operator status * () 
    { return &_status; }

    operator status const * () const
    { return &_status; }
    
//     operator status & ()
//     { return _status; }

//     operator const status & () const
//     { return _status; }

    int source () const
    {
#ifdef LIBMESH_HAVE_MPI 
      return _status.MPI_SOURCE; 
#else
      return 0;
#endif
    }    

    int tag () const
    {
#ifdef LIBMESH_HAVE_MPI 
      return _status.MPI_TAG; 
#else
      libmesh_error();
      return 0;
#endif
    }

    data_type& datatype () 
    { return _datatype; }

    const data_type& datatype () const
    { return _datatype; }
  
#ifdef LIBMESH_HAVE_MPI
    unsigned int size (const data_type &type) const
    {
      int msg_size;
      MPI_Get_count (const_cast<MPI_Status*>(&_status), type, &msg_size);
      libmesh_assert (msg_size >= 0);
      return msg_size;
    }
#else
    unsigned int size (const data_type &) const
    {
      libmesh_error();
      return 0;
    }
#endif

    unsigned int size () const
    { return this->size (this->datatype()); }

  private:

    status    _status;
    data_type _datatype;
  };



  //-------------------------------------------------------------------
  /**
   * Encapsulates the MPI_Request
   */
  class Request
  {
  public:
    Request () 
    {
#ifdef LIBMESH_HAVE_MPI
      _request = MPI_REQUEST_NULL;
#endif
    }

    Request (const request &r) :
      _request(r)
    {}

    Request & operator = (const Request &other) 
    { _request = other._request; return *this; }

    Request & operator = (const request &r)
    { _request = r; return *this; }

    ~Request ()
    {
#ifdef LIBMESH_HAVE_MPI
      // explicitly free this request if not 
      // done so already, otherwise this would
      // be a memory leak!
      if (_request != MPI_REQUEST_NULL)
	MPI_Request_free (&_request);
#endif
    }

    operator const request & () const
    { return _request; }

    operator request & () 
    { return _request; }

    status wait ()
    {
      status status;
#ifdef LIBMESH_HAVE_MPI
      MPI_Wait (&_request, &status);
#endif
      return status;
    }

    bool test ()
    {
#ifdef LIBMESH_HAVE_MPI
      int val=0;

      MPI_Test (&_request,
		&val,
		MPI_STATUS_IGNORE);
      if (val)
	{
	  libmesh_assert (_request == MPI_REQUEST_NULL);
	  libmesh_assert (val == 1);
	}

      return val;
#else
      return true;
#endif
    }

#ifdef LIBMESH_HAVE_MPI
    bool test (status &status)
    {
      int val=0;

      MPI_Test (&_request,
		&val,
		&status);

      return val;
#else
    bool test (status &)
    {
      return true;
#endif
    }


  private:

    request _request;    
  };


  
  //-------------------------------------------------------------------
  /**
   * Pause execution until all processors reach a certain point.
   */
  inline void barrier ()
  {
#ifdef LIBMESH_HAVE_MPI
    MPI_Barrier (libMesh::COMM_WORLD);
#endif
  }    
  
  //-------------------------------------------------------------------
  /**
   * Verify that a local variable has the same value on all processors
   */
  template <typename T>
  inline bool verify(const T &r);

  //-------------------------------------------------------------------
  /**
   * Take a local variable and replace it with the minimum of it's values
   * on all processors
   */
  template <typename T>
  inline void min(T &r);

  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and replace each entry with the minimum
   * of it's values on all processors
   */
  template <typename T>
  inline void min(std::vector<T> &r);

  //-------------------------------------------------------------------
  /**
   * Take a local variable and replace it with the maximum of it's values
   * on all processors
   */
  template <typename T>
  inline void max(T &r);

  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and replace each entry with the maximum
   * of it's values on all processors
   */
  template <typename T>
  inline void max(std::vector<T> &r);

  //-------------------------------------------------------------------
  /**
   * Take a local variable and replace it with the sum of it's values
   * on all processors
   */
  template <typename T>
  inline void sum(T &r);

  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and replace each entry with the sum of
   * it's values on all processors
   */
  template <typename T>
  inline void sum(std::vector<T> &r);

  //-------------------------------------------------------------------
  /**
   * Blocking message probe.  Allows information about a message to be 
   * examined before the message is actually received.
   */
  inline status probe (const int src_processor_id,
		       const int tag=any_tag);

  //-------------------------------------------------------------------
  /**
   * Blocking-send vector to one processor with user-defined type.
   */
  template <typename T>
  inline void send (const unsigned int dest_processor_id,
		    std::vector<T> &buf,
		    const DataType &type,
		    const int tag=0);

  //-------------------------------------------------------------------
  /**
   * Nonblocking-send vector to one processor with user-defined type.
   */
  template <typename T>
  inline void send (const unsigned int dest_processor_id,
		    std::vector<T> &buf,
		    const DataType &type,
		    request &req,
		    const int tag=0);

  //-------------------------------------------------------------------
  /**
   * Blocking-send vector to one processor where the communication type 
   * is inferred from the template argument.
   */
  template <typename T>
  inline void send (const unsigned int dest_processor_id,
		    std::vector<T> &buf,
		    const int tag=0)
  {
    send (dest_processor_id,
	  buf,
	  datatype<T>(),
	  tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-send vector to one processor where the communication type 
   * is inferred from the template argument.
   */
  template <typename T>
  inline void send (const unsigned int dest_processor_id,
		    std::vector<T> &buf,
		    request &req,
		    const int tag=0)
  {
    send (dest_processor_id,
	  buf,
	  datatype<T>(),
	  req,
	  tag);
  }


  //-------------------------------------------------------------------
  /**
   * Nonblocking-send vector to one processor with user-defined type.
   */
  template <typename T>
  inline void nonblocking_send (const unsigned int dest_processor_id,
		                std::vector<T> &buf,
		                const DataType &type,
		                request &r,
		                const int tag=0)
  {
    send (dest_processor_id,
	  buf,
	  type,
	  r,
	  tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-send vector to one processor.
   */
  template <typename T>
  inline void nonblocking_send (const unsigned int dest_processor_id,
		                std::vector<T> &buf,
		                request &r,
		                const int tag=0)
  {
    send (dest_processor_id,
	  buf,
	  datatype<T>(),
	  r,
	  tag);
  }

  //-------------------------------------------------------------------
  /**
   * Blocking-receive vector from one processor with user-defined type.
   */
  template <typename T>
  inline Status receive (const int src_processor_id,
		         std::vector<T> &buf,
		         const DataType &type,
		         const int tag=any_tag);

  //-------------------------------------------------------------------
  /**
   * Nonblocking-receive vector from one processor with user-defined type.
   */
  template <typename T>
  inline void receive (const int src_processor_id,
		       std::vector<T> &buf,
		       const DataType &type,
		       request &req,
		       const int tag=any_tag);

  //-------------------------------------------------------------------
  /**
   * Blocking-receive vector from one processor where the communication type 
   * is inferred from the template argument.
   */
  template <typename T>
  inline Status receive (const int src_processor_id,
		         std::vector<T> &buf,
		         const int tag=any_tag)
  {
    return receive (src_processor_id,
		    buf,
		    datatype<T>(),
		    tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-receive vector from one processor where the communication type 
   * is inferred from the template argument.
   */
  template <typename T>
  inline void receive (const int src_processor_id,
		       std::vector<T> &buf,
		       request &req,
		       const int tag=any_tag)
  {
    receive (src_processor_id,
	     buf,
	     datatype<T>(),
	     req,
	     tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-receive vector from one processor with user-defined type
   */
  template <typename T>
  inline void nonblocking_receive (const int src_processor_id,
		                   std::vector<T> &buf,
				   const DataType &type,
		                   request &r,
		                   const int tag=any_tag)
  {
    receive (src_processor_id,
	     buf,
	     type,
	     r,
	     tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-receive vector from one processor.
   */
  template <typename T>
  inline void nonblocking_receive (const int src_processor_id,
		                   std::vector<T> &buf,
		                   request &r,
		                   const int tag=any_tag)
  {
    receive (src_processor_id,
	     buf,
	     datatype<T>(),
	     r,
	     tag);
  }

  
  //-------------------------------------------------------------------
  /**
   * Wait for a non-blocking send or receive to finish
   */
  inline status wait (request &r);
  
  //-------------------------------------------------------------------
  /**
   * Wait for a non-blocking send or receive to finish
   */
  inline void wait (std::vector<request> &r);
  
  //-------------------------------------------------------------------
  /**
   * Wait for a non-blocking send or receive to finish
   */
  inline void wait (std::vector<Request> &r)
  { for (unsigned int i=0; i<r.size(); i++) r[i].wait(); }
  
  //-------------------------------------------------------------------
  /**
   * Send vector \p send to one processor while simultaneously receiving
   * another vector \p recv from a (potentially different) processor.
   */
  template <typename T>
  inline void send_receive(const unsigned int dest_processor_id,
                           T &send,
			   const unsigned int source_processor_id,
                           T &recv);

  //-------------------------------------------------------------------
  /**
   * Send vector \p send to one processor while simultaneously receiving
   * another vector \p recv from a (potentially different) processor using
   * a user-specified MPI Dataype.
   */
  template <typename T>
  inline void send_receive(const unsigned int dest_processor_id,
                           T &send,
			   const unsigned int source_processor_id,
                           T &recv,
			   const DataType &type);

  //-------------------------------------------------------------------
  /**
   * Take a vector of length n_processors, and on processor root_id fill in
   * recv[processor_id] = the value of send on processor processor_id
   */
  template <typename T>
  inline void gather(const unsigned int root_id,
		     T send,
		     std::vector<T> &recv);

  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and expand it on processor root_id 
   * to include values from all processors
   */
  template <typename T>
  inline void gather(const unsigned int root_id,
		     std::vector<T> &r);

  //-------------------------------------------------------------------
  /**
   * Take a vector of length n_processors, and fill in 
   * \p recv[processor_id] = the value of \p send on that processor
   */
  template <typename T>
  inline void allgather(T send,
			std::vector<T> &recv);


  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and expand it to include 
   * values from all processors. By default, each processor is 
   * allowed to have its own unique input buffer length. If 
   * it is known that all processors have the same input sizes
   * additional communication can be avoided.
   */
  template <typename T>
  inline void allgather(std::vector<T> &r,
			const bool identical_buffer_sizes = false);



  //-------------------------------------------------------------------
  /**
   * Effectively transposes the input vector across all processors.  
   * The jth entry on processor i is replaced with the ith entry 
   * from processor j.
   */
  template <typename T>
  inline void alltoall(std::vector<T> &r);



  //-------------------------------------------------------------------
  /**
   * Take a local value and broadcast it to all processors.
   * Optionally takes the \p root_id processor, which specifies
   * the processor intiating the broadcast.
   */
  template <typename T>
  inline void broadcast(T &data, const unsigned int root_id=0);



  //-------------------------------------------------------------------
  /**
   * Take a local vector and broadcast it to all processors.
   * Optionally takes the \p root_id processor, which specifies
   * the processor intiating the broadcast.  The user is responsible
   * for appropriately sizing the input buffer on all processors.
   */
  template <typename T>
  inline void broadcast(std::vector<T> &data, const unsigned int root_id=0);

  // gcc appears to need an additional declaration to make sure it
  // uses the right definition below
  template <typename T>
  inline void broadcast(std::vector<std::complex<T> > &data, const unsigned int root_id=0);


  //-----------------------------------------------------------------------
  // Parallel members

  // Internal helper function to create vector<something_useable> from
  // vector<bool> for compatibility with MPI bitwise operations
  template <typename T>
  inline void pack_vector_bool(const std::vector<bool> &in,
			       std::vector<T> &out)
  {
    unsigned int data_bits = 8*sizeof(T);
    unsigned int in_size = in.size();
    unsigned int out_size = in_size/data_bits + (in_size%data_bits?1:0);
    out.clear();
    out.resize(out_size);
    for (unsigned int i=0; i != in_size; ++i)
      {
        unsigned int index = i/data_bits;
        unsigned int offset = i%data_bits;
        out[index] += (in[i]?1:0) << offset;
      }
  }

  // Internal helper function to create vector<something_useable> from
  // vector<bool> for compatibility with MPI byte operations
  template <typename T>
  inline void unpack_vector_bool(const std::vector<T> &in,
				 std::vector<bool> &out)
  {
    unsigned int data_bits = 8*sizeof(T);
    // We need the output vector to already be properly sized
    unsigned int out_size = out.size();
    libmesh_assert(out_size/data_bits + (out_size%data_bits?1:0) == in.size());

    for (unsigned int i=0; i != out_size; ++i)
      {
        unsigned int index = i/data_bits;
        unsigned int offset = i%data_bits;
        out[i] = in[index] << (data_bits-1-offset) >> (data_bits-1);
      }
  }

#ifdef LIBMESH_HAVE_MPI
 template<>
 inline data_type datatype<char>() { return MPI_CHAR; }

 template<>
 inline data_type datatype<unsigned char>() { return MPI_UNSIGNED_CHAR; }

  template<>
  inline data_type datatype<short int>() { return MPI_SHORT; }

  template<>
  inline data_type datatype<unsigned short int>() { return MPI_UNSIGNED_SHORT; }

  template<>
  inline data_type datatype<int>() { return MPI_INT; }

  template<>
  inline data_type datatype<unsigned int>() { return MPI_UNSIGNED; }

  template<>
  inline data_type datatype<long>() { return MPI_LONG; }

  template<>
  inline data_type datatype<unsigned long>() { return MPI_UNSIGNED_LONG; }

  template<>
  inline data_type datatype<float>() { return MPI_FLOAT; }

  template<>
  inline data_type datatype<double>() { return MPI_DOUBLE; }

  template<>
  inline data_type datatype<long double>() { return MPI_LONG_DOUBLE; }

  template <typename T>
  inline bool verify(const T &r)
  {
    if (libMesh::n_processors() > 1)
      {
	T tempmin = r, tempmax = r;
	Parallel::min(tempmin);
	Parallel::max(tempmax);
	bool verified = (r == tempmin) &&
	                (r == tempmax);
	Parallel::min(verified);
	return verified;
      }
    return true;
  }


  template <>
  inline bool verify(const std::string & r)
  {
    if (libMesh::n_processors() > 1)
      {
	// Cannot use <char> since MPI_MIN is not
	// strictly defined for chars!
	std::vector<short int> temp; temp.reserve(r.size());
	for (unsigned int i=0; i != r.size(); ++i)
	  temp.push_back(r[i]);
	return Parallel::verify(temp);
      }
    return true;
  }


  template <typename T>
  inline void min(T &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("min()", "Parallel");
    
	T temp = r;
	MPI_Allreduce (&temp,
		       &r,
		       1,
		       datatype<T>(),
		       MPI_MIN,
		       libMesh::COMM_WORLD);

	STOP_LOG("min()", "Parallel");
      }
  }


  template <>
  inline void min(bool &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("min()", "Parallel");
    
	unsigned int tempsend = r;
	unsigned int temp;
	MPI_Allreduce (&tempsend,
		       &temp,
		       1,
		       datatype<unsigned int>(),
		       MPI_MIN,
		       libMesh::COMM_WORLD);
	r = temp;

	STOP_LOG("min()", "Parallel");
      }
  }


  template <typename T>
  inline void min(std::vector<T> &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("min()", "Parallel");
    
	std::vector<T> temp(r);
	MPI_Allreduce (&temp[0],
		       &r[0],
		       r.size(),
		       datatype<T>(),
		       MPI_MIN,
		       libMesh::COMM_WORLD);

	STOP_LOG("min()", "Parallel");
      }
  }


  template <>
  inline void min(std::vector<bool> &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("min()", "Parallel");
    
        std::vector<unsigned int> ruint;
        pack_vector_bool(r, ruint);
	std::vector<unsigned int> temp(ruint.size());
	MPI_Allreduce (&ruint[0],
		       &temp[0],
		       ruint.size(),
		       datatype<unsigned int>(),
		       MPI_BAND,
		       libMesh::COMM_WORLD);
        unpack_vector_bool(temp, r);

	STOP_LOG("min()", "Parallel");
      }
  }


  template <typename T>
  inline void max(T &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("max()", "Parallel");
    
	T temp;
	MPI_Allreduce (&r,
		       &temp,
		       1,
		       datatype<T>(),
		       MPI_MAX,
		       libMesh::COMM_WORLD);
	r = temp;

	STOP_LOG("max()", "Parallel");
      }
  }


  template <>
  inline void max(bool &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("max()", "Parallel");
    
	unsigned int tempsend = r;
	unsigned int temp;
	MPI_Allreduce (&tempsend,
		       &temp,
		       1,
		       datatype<unsigned int>(),
		       MPI_MAX,
		       libMesh::COMM_WORLD);
	r = temp;

	STOP_LOG("max()", "Parallel");
      }
  }


  template <typename T>
  inline void max(std::vector<T> &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("max()", "Parallel");
    
	std::vector<T> temp(r);
	MPI_Allreduce (&temp[0],
		       &r[0],
		       r.size(),
		       datatype<T>(),
		       MPI_MAX,
		       libMesh::COMM_WORLD);

	STOP_LOG("max()", "Parallel");
      }
  }


  template <>
  inline void max(std::vector<bool> &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("max()", "Parallel");
    
        std::vector<unsigned int> ruint;
        pack_vector_bool(r, ruint);
	std::vector<unsigned int> temp(ruint.size());
	MPI_Allreduce (&ruint[0],
		       &temp[0],
		       ruint.size(),
		       datatype<unsigned int>(),
		       MPI_BOR,
		       libMesh::COMM_WORLD);
        unpack_vector_bool(temp, r);

	STOP_LOG("max()", "Parallel");
      }
  }


  template <typename T>
  inline void sum(T &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("sum()", "Parallel");
    
	T temp = r;
	MPI_Allreduce (&temp,<commit_msg>Fixes for Parallel:: operations on empty vectors<commit_after>// $Id$

// The libMesh Finite Element Library.
// Copyright (C) 2002-2008 <NAME>, <NAME>, <NAME>
  
// 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 __parallel_h__
#define __parallel_h__

// System includes
#include <string>
#include <vector>

// Local includes
#include "libmesh_common.h" // for Real
#include "libmesh_logging.h"


// Macro to identify and debug functions which should only be called in
// parallel on every processor at once

#undef parallel_only
#ifndef NDEBUG
  #define parallel_only() do { libmesh_assert(Parallel::verify(std::string(__FILE__))); libmesh_assert(Parallel::verify(__LINE__)); } while (0)
#else
  #define parallel_only()
#endif

/**
 * The Parallel namespace is for wrapper functions
 * for common general parallel synchronization tasks.
 *
 * For MPI 1.1 compatibility, temporary buffers are used
 * instead of MPI 2's MPI_IN_PLACE
 */
namespace Parallel
{
  //-------------------------------------------------------------------
  /**
   * Forward declarations of classes we will define later.
   */
  class DataType;
  class Request;
  class Status;

#ifdef LIBMESH_HAVE_MPI
  //-------------------------------------------------------------------
  /**
   * Data types for communication
   */
  typedef MPI_Datatype data_type;

  /**
   * Request object for non-blocking I/O
   */
  typedef MPI_Request request;

  /**
   * Status object for querying messages
   */
  typedef MPI_Status status;
 
  /**
   * Templated function to return the appropriate MPI datatype
   * for use with built-in C types
   */
  template <typename T>
  inline data_type datatype();

  /**
   * Default message tag id
   */
  const int any_tag=MPI_ANY_TAG;

  /**
   * Accept from any source
   */
  const int any_source=MPI_ANY_SOURCE;
  
#else

  // These shouldn't actually be needed, but must be 
  // unique types for function overloading to work 
  // properly.
  struct data_type { /* unsigned int t; */ };
  struct request   { /* unsigned int r; */ };
  struct status    { /* unsigned int s; */ };

  template <typename T>
  inline data_type datatype() { return data_type(); }

  const int any_tag=-1;
  const int any_source=0;
#endif // LIBMESH_HAVE_MPI



  //-------------------------------------------------------------------
  /**
   * Encapsulates the MPI_Datatype.
   */
  class DataType
  {
  public:
    DataType () {}

    DataType (const DataType &other) :
      _datatype(other._datatype)
    {}

    DataType (const data_type &type) :
      _datatype(type)
    {}

    DataType & operator = (const DataType &other) 
    { _datatype = other._datatype; return *this; }

    DataType & operator = (const data_type &type)
    { _datatype = type; return *this; }

    operator const data_type & () const
    { return _datatype; }

    operator data_type & ()
    { return _datatype; }
    
//     operator data_type const * () const
//     { return &_datatype; }

//     operator data_type * ()
//     { return &_datatype; }

    void commit ()
    {
#ifdef LIBMESH_HAVE_MPI
      MPI_Type_commit (&_datatype);
#endif
    }

    void free ()
    {
#ifdef LIBMESH_HAVE_MPI
      MPI_Type_free (&_datatype);
#endif
    }
    
  private:

    data_type _datatype;    
  };



  //-------------------------------------------------------------------
  /**
   * Encapsulates the MPI_Status struct.  Allows the source and size 
   * of the message to be determined.
   */
  class Status
  {
  public:
    Status () {}
    
    Status (const data_type &type) :
    _datatype(type)           
    {}

    Status (const status &status) :
      _status(status)
    {}

    Status (const status    &status,
	    const data_type &type) :
      _status(status),
      _datatype(type)           
    {}

    Status (const Status &status) :
      _status(status._status),
      _datatype(status._datatype)
    {}

    Status (const Status    &status,
	    const data_type &type) :
      _status(status._status),
      _datatype(type)
    {}
    
    operator status * () 
    { return &_status; }

    operator status const * () const
    { return &_status; }
    
//     operator status & ()
//     { return _status; }

//     operator const status & () const
//     { return _status; }

    int source () const
    {
#ifdef LIBMESH_HAVE_MPI 
      return _status.MPI_SOURCE; 
#else
      return 0;
#endif
    }    

    int tag () const
    {
#ifdef LIBMESH_HAVE_MPI 
      return _status.MPI_TAG; 
#else
      libmesh_error();
      return 0;
#endif
    }

    data_type& datatype () 
    { return _datatype; }

    const data_type& datatype () const
    { return _datatype; }
  
#ifdef LIBMESH_HAVE_MPI
    unsigned int size (const data_type &type) const
    {
      int msg_size;
      MPI_Get_count (const_cast<MPI_Status*>(&_status), type, &msg_size);
      libmesh_assert (msg_size >= 0);
      return msg_size;
    }
#else
    unsigned int size (const data_type &) const
    {
      libmesh_error();
      return 0;
    }
#endif

    unsigned int size () const
    { return this->size (this->datatype()); }

  private:

    status    _status;
    data_type _datatype;
  };



  //-------------------------------------------------------------------
  /**
   * Encapsulates the MPI_Request
   */
  class Request
  {
  public:
    Request () 
    {
#ifdef LIBMESH_HAVE_MPI
      _request = MPI_REQUEST_NULL;
#endif
    }

    Request (const request &r) :
      _request(r)
    {}

    Request & operator = (const Request &other) 
    { _request = other._request; return *this; }

    Request & operator = (const request &r)
    { _request = r; return *this; }

    ~Request ()
    {
#ifdef LIBMESH_HAVE_MPI
      // explicitly free this request if not 
      // done so already, otherwise this would
      // be a memory leak!
      if (_request != MPI_REQUEST_NULL)
	MPI_Request_free (&_request);
#endif
    }

    operator const request & () const
    { return _request; }

    operator request & () 
    { return _request; }

    status wait ()
    {
      status status;
#ifdef LIBMESH_HAVE_MPI
      MPI_Wait (&_request, &status);
#endif
      return status;
    }

    bool test ()
    {
#ifdef LIBMESH_HAVE_MPI
      int val=0;

      MPI_Test (&_request,
		&val,
		MPI_STATUS_IGNORE);
      if (val)
	{
	  libmesh_assert (_request == MPI_REQUEST_NULL);
	  libmesh_assert (val == 1);
	}

      return val;
#else
      return true;
#endif
    }

#ifdef LIBMESH_HAVE_MPI
    bool test (status &status)
    {
      int val=0;

      MPI_Test (&_request,
		&val,
		&status);

      return val;
#else
    bool test (status &)
    {
      return true;
#endif
    }


  private:

    request _request;    
  };


  
  //-------------------------------------------------------------------
  /**
   * Pause execution until all processors reach a certain point.
   */
  inline void barrier ()
  {
#ifdef LIBMESH_HAVE_MPI
    MPI_Barrier (libMesh::COMM_WORLD);
#endif
  }    
  
  //-------------------------------------------------------------------
  /**
   * Verify that a local variable has the same value on all processors
   */
  template <typename T>
  inline bool verify(const T &r);

  //-------------------------------------------------------------------
  /**
   * Take a local variable and replace it with the minimum of it's values
   * on all processors
   */
  template <typename T>
  inline void min(T &r);

  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and replace each entry with the minimum
   * of it's values on all processors
   */
  template <typename T>
  inline void min(std::vector<T> &r);

  //-------------------------------------------------------------------
  /**
   * Take a local variable and replace it with the maximum of it's values
   * on all processors
   */
  template <typename T>
  inline void max(T &r);

  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and replace each entry with the maximum
   * of it's values on all processors
   */
  template <typename T>
  inline void max(std::vector<T> &r);

  //-------------------------------------------------------------------
  /**
   * Take a local variable and replace it with the sum of it's values
   * on all processors
   */
  template <typename T>
  inline void sum(T &r);

  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and replace each entry with the sum of
   * it's values on all processors
   */
  template <typename T>
  inline void sum(std::vector<T> &r);

  //-------------------------------------------------------------------
  /**
   * Blocking message probe.  Allows information about a message to be 
   * examined before the message is actually received.
   */
  inline status probe (const int src_processor_id,
		       const int tag=any_tag);

  //-------------------------------------------------------------------
  /**
   * Blocking-send vector to one processor with user-defined type.
   */
  template <typename T>
  inline void send (const unsigned int dest_processor_id,
		    std::vector<T> &buf,
		    const DataType &type,
		    const int tag=0);

  //-------------------------------------------------------------------
  /**
   * Nonblocking-send vector to one processor with user-defined type.
   */
  template <typename T>
  inline void send (const unsigned int dest_processor_id,
		    std::vector<T> &buf,
		    const DataType &type,
		    request &req,
		    const int tag=0);

  //-------------------------------------------------------------------
  /**
   * Blocking-send vector to one processor where the communication type 
   * is inferred from the template argument.
   */
  template <typename T>
  inline void send (const unsigned int dest_processor_id,
		    std::vector<T> &buf,
		    const int tag=0)
  {
    send (dest_processor_id,
	  buf,
	  datatype<T>(),
	  tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-send vector to one processor where the communication type 
   * is inferred from the template argument.
   */
  template <typename T>
  inline void send (const unsigned int dest_processor_id,
		    std::vector<T> &buf,
		    request &req,
		    const int tag=0)
  {
    send (dest_processor_id,
	  buf,
	  datatype<T>(),
	  req,
	  tag);
  }


  //-------------------------------------------------------------------
  /**
   * Nonblocking-send vector to one processor with user-defined type.
   */
  template <typename T>
  inline void nonblocking_send (const unsigned int dest_processor_id,
		                std::vector<T> &buf,
		                const DataType &type,
		                request &r,
		                const int tag=0)
  {
    send (dest_processor_id,
	  buf,
	  type,
	  r,
	  tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-send vector to one processor.
   */
  template <typename T>
  inline void nonblocking_send (const unsigned int dest_processor_id,
		                std::vector<T> &buf,
		                request &r,
		                const int tag=0)
  {
    send (dest_processor_id,
	  buf,
	  datatype<T>(),
	  r,
	  tag);
  }

  //-------------------------------------------------------------------
  /**
   * Blocking-receive vector from one processor with user-defined type.
   */
  template <typename T>
  inline Status receive (const int src_processor_id,
		         std::vector<T> &buf,
		         const DataType &type,
		         const int tag=any_tag);

  //-------------------------------------------------------------------
  /**
   * Nonblocking-receive vector from one processor with user-defined type.
   */
  template <typename T>
  inline void receive (const int src_processor_id,
		       std::vector<T> &buf,
		       const DataType &type,
		       request &req,
		       const int tag=any_tag);

  //-------------------------------------------------------------------
  /**
   * Blocking-receive vector from one processor where the communication type 
   * is inferred from the template argument.
   */
  template <typename T>
  inline Status receive (const int src_processor_id,
		         std::vector<T> &buf,
		         const int tag=any_tag)
  {
    return receive (src_processor_id,
		    buf,
		    datatype<T>(),
		    tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-receive vector from one processor where the communication type 
   * is inferred from the template argument.
   */
  template <typename T>
  inline void receive (const int src_processor_id,
		       std::vector<T> &buf,
		       request &req,
		       const int tag=any_tag)
  {
    receive (src_processor_id,
	     buf,
	     datatype<T>(),
	     req,
	     tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-receive vector from one processor with user-defined type
   */
  template <typename T>
  inline void nonblocking_receive (const int src_processor_id,
		                   std::vector<T> &buf,
				   const DataType &type,
		                   request &r,
		                   const int tag=any_tag)
  {
    receive (src_processor_id,
	     buf,
	     type,
	     r,
	     tag);
  }

  //-------------------------------------------------------------------
  /**
   * Nonblocking-receive vector from one processor.
   */
  template <typename T>
  inline void nonblocking_receive (const int src_processor_id,
		                   std::vector<T> &buf,
		                   request &r,
		                   const int tag=any_tag)
  {
    receive (src_processor_id,
	     buf,
	     datatype<T>(),
	     r,
	     tag);
  }

  
  //-------------------------------------------------------------------
  /**
   * Wait for a non-blocking send or receive to finish
   */
  inline status wait (request &r);
  
  //-------------------------------------------------------------------
  /**
   * Wait for a non-blocking send or receive to finish
   */
  inline void wait (std::vector<request> &r);
  
  //-------------------------------------------------------------------
  /**
   * Wait for a non-blocking send or receive to finish
   */
  inline void wait (std::vector<Request> &r)
  { for (unsigned int i=0; i<r.size(); i++) r[i].wait(); }
  
  //-------------------------------------------------------------------
  /**
   * Send vector \p send to one processor while simultaneously receiving
   * another vector \p recv from a (potentially different) processor.
   */
  template <typename T>
  inline void send_receive(const unsigned int dest_processor_id,
                           T &send,
			   const unsigned int source_processor_id,
                           T &recv);

  //-------------------------------------------------------------------
  /**
   * Send vector \p send to one processor while simultaneously receiving
   * another vector \p recv from a (potentially different) processor using
   * a user-specified MPI Dataype.
   */
  template <typename T>
  inline void send_receive(const unsigned int dest_processor_id,
                           T &send,
			   const unsigned int source_processor_id,
                           T &recv,
			   const DataType &type);

  //-------------------------------------------------------------------
  /**
   * Take a vector of length n_processors, and on processor root_id fill in
   * recv[processor_id] = the value of send on processor processor_id
   */
  template <typename T>
  inline void gather(const unsigned int root_id,
		     T send,
		     std::vector<T> &recv);

  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and expand it on processor root_id 
   * to include values from all processors
   */
  template <typename T>
  inline void gather(const unsigned int root_id,
		     std::vector<T> &r);

  //-------------------------------------------------------------------
  /**
   * Take a vector of length n_processors, and fill in 
   * \p recv[processor_id] = the value of \p send on that processor
   */
  template <typename T>
  inline void allgather(T send,
			std::vector<T> &recv);


  //-------------------------------------------------------------------
  /**
   * Take a vector of local variables and expand it to include 
   * values from all processors. By default, each processor is 
   * allowed to have its own unique input buffer length. If 
   * it is known that all processors have the same input sizes
   * additional communication can be avoided.
   */
  template <typename T>
  inline void allgather(std::vector<T> &r,
			const bool identical_buffer_sizes = false);



  //-------------------------------------------------------------------
  /**
   * Effectively transposes the input vector across all processors.  
   * The jth entry on processor i is replaced with the ith entry 
   * from processor j.
   */
  template <typename T>
  inline void alltoall(std::vector<T> &r);



  //-------------------------------------------------------------------
  /**
   * Take a local value and broadcast it to all processors.
   * Optionally takes the \p root_id processor, which specifies
   * the processor intiating the broadcast.
   */
  template <typename T>
  inline void broadcast(T &data, const unsigned int root_id=0);



  //-------------------------------------------------------------------
  /**
   * Take a local vector and broadcast it to all processors.
   * Optionally takes the \p root_id processor, which specifies
   * the processor intiating the broadcast.  The user is responsible
   * for appropriately sizing the input buffer on all processors.
   */
  template <typename T>
  inline void broadcast(std::vector<T> &data, const unsigned int root_id=0);

  // gcc appears to need an additional declaration to make sure it
  // uses the right definition below
  template <typename T>
  inline void broadcast(std::vector<std::complex<T> > &data, const unsigned int root_id=0);


  //-----------------------------------------------------------------------
  // Parallel members

  // Internal helper function to create vector<something_useable> from
  // vector<bool> for compatibility with MPI bitwise operations
  template <typename T>
  inline void pack_vector_bool(const std::vector<bool> &in,
			       std::vector<T> &out)
  {
    unsigned int data_bits = 8*sizeof(T);
    unsigned int in_size = in.size();
    unsigned int out_size = in_size/data_bits + (in_size%data_bits?1:0);
    out.clear();
    out.resize(out_size);
    for (unsigned int i=0; i != in_size; ++i)
      {
        unsigned int index = i/data_bits;
        unsigned int offset = i%data_bits;
        out[index] += (in[i]?1:0) << offset;
      }
  }

  // Internal helper function to create vector<something_useable> from
  // vector<bool> for compatibility with MPI byte operations
  template <typename T>
  inline void unpack_vector_bool(const std::vector<T> &in,
				 std::vector<bool> &out)
  {
    unsigned int data_bits = 8*sizeof(T);
    // We need the output vector to already be properly sized
    unsigned int out_size = out.size();
    libmesh_assert(out_size/data_bits + (out_size%data_bits?1:0) == in.size());

    for (unsigned int i=0; i != out_size; ++i)
      {
        unsigned int index = i/data_bits;
        unsigned int offset = i%data_bits;
        out[i] = in[index] << (data_bits-1-offset) >> (data_bits-1);
      }
  }

#ifdef LIBMESH_HAVE_MPI
 template<>
 inline data_type datatype<char>() { return MPI_CHAR; }

 template<>
 inline data_type datatype<unsigned char>() { return MPI_UNSIGNED_CHAR; }

  template<>
  inline data_type datatype<short int>() { return MPI_SHORT; }

  template<>
  inline data_type datatype<unsigned short int>() { return MPI_UNSIGNED_SHORT; }

  template<>
  inline data_type datatype<int>() { return MPI_INT; }

  template<>
  inline data_type datatype<unsigned int>() { return MPI_UNSIGNED; }

  template<>
  inline data_type datatype<long>() { return MPI_LONG; }

  template<>
  inline data_type datatype<unsigned long>() { return MPI_UNSIGNED_LONG; }

  template<>
  inline data_type datatype<float>() { return MPI_FLOAT; }

  template<>
  inline data_type datatype<double>() { return MPI_DOUBLE; }

  template<>
  inline data_type datatype<long double>() { return MPI_LONG_DOUBLE; }

  template <typename T>
  inline bool verify(const T &r)
  {
    if (libMesh::n_processors() > 1)
      {
	T tempmin = r, tempmax = r;
	Parallel::min(tempmin);
	Parallel::max(tempmax);
	bool verified = (r == tempmin) &&
	                (r == tempmax);
	Parallel::min(verified);
	return verified;
      }
    return true;
  }


  template <>
  inline bool verify(const std::string & r)
  {
    if (libMesh::n_processors() > 1)
      {
	// Cannot use <char> since MPI_MIN is not
	// strictly defined for chars!
	std::vector<short int> temp; temp.reserve(r.size());
	for (unsigned int i=0; i != r.size(); ++i)
	  temp.push_back(r[i]);
	return Parallel::verify(temp);
      }
    return true;
  }


  template <typename T>
  inline void min(T &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("min()", "Parallel");
    
	T temp = r;
	MPI_Allreduce (&temp,
		       &r,
		       1,
		       datatype<T>(),
		       MPI_MIN,
		       libMesh::COMM_WORLD);

	STOP_LOG("min()", "Parallel");
      }
  }


  template <>
  inline void min(bool &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("min()", "Parallel");
    
	unsigned int tempsend = r;
	unsigned int temp;
	MPI_Allreduce (&tempsend,
		       &temp,
		       1,
		       datatype<unsigned int>(),
		       MPI_MIN,
		       libMesh::COMM_WORLD);
	r = temp;

	STOP_LOG("min()", "Parallel");
      }
  }


  template <typename T>
  inline void min(std::vector<T> &r)
  {
    if (libMesh::n_processors() > 1 && !r.empty())
      {
	START_LOG("min()", "Parallel");
    
	std::vector<T> temp(r);
	MPI_Allreduce (&temp[0],
		       &r[0],
		       r.size(),
		       datatype<T>(),
		       MPI_MIN,
		       libMesh::COMM_WORLD);

	STOP_LOG("min()", "Parallel");
      }
  }


  template <>
  inline void min(std::vector<bool> &r)
  {
    if (libMesh::n_processors() > 1 && !r.empty())
      {
	START_LOG("min()", "Parallel");
    
        std::vector<unsigned int> ruint;
        pack_vector_bool(r, ruint);
	std::vector<unsigned int> temp(ruint.size());
	MPI_Allreduce (&ruint[0],
		       &temp[0],
		       ruint.size(),
		       datatype<unsigned int>(),
		       MPI_BAND,
		       libMesh::COMM_WORLD);
        unpack_vector_bool(temp, r);

	STOP_LOG("min()", "Parallel");
      }
  }


  template <typename T>
  inline void max(T &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("max()", "Parallel");
    
	T temp;
	MPI_Allreduce (&r,
		       &temp,
		       1,
		       datatype<T>(),
		       MPI_MAX,
		       libMesh::COMM_WORLD);
	r = temp;

	STOP_LOG("max()", "Parallel");
      }
  }


  template <>
  inline void max(bool &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("max()", "Parallel");
    
	unsigned int tempsend = r;
	unsigned int temp;
	MPI_Allreduce (&tempsend,
		       &temp,
		       1,
		       datatype<unsigned int>(),
		       MPI_MAX,
		       libMesh::COMM_WORLD);
	r = temp;

	STOP_LOG("max()", "Parallel");
      }
  }


  template <typename T>
  inline void max(std::vector<T> &r)
  {
    if (libMesh::n_processors() > 1 && !r.empty())
      {
	START_LOG("max()", "Parallel");
    
	std::vector<T> temp(r);
	MPI_Allreduce (&temp[0],
		       &r[0],
		       r.size(),
		       datatype<T>(),
		       MPI_MAX,
		       libMesh::COMM_WORLD);

	STOP_LOG("max()", "Parallel");
      }
  }


  template <>
  inline void max(std::vector<bool> &r)
  {
    if (libMesh::n_processors() > 1 && !r.empty())
      {
	START_LOG("max()", "Parallel");
    
        std::vector<unsigned int> ruint;
        pack_vector_bool(r, ruint);
	std::vector<unsigned int> temp(ruint.size());
	MPI_Allreduce (&ruint[0],
		       &temp[0],
		       ruint.size(),
		       datatype<unsigned int>(),
		       MPI_BOR,
		       libMesh::COMM_WORLD);
        unpack_vector_bool(temp, r);

	STOP_LOG("max()", "Parallel");
      }
  }


  template <typename T>
  inline void sum(T &r)
  {
    if (libMesh::n_processors() > 1)
      {
	START_LOG("sum()", "Parallel");
    
	T temp = r;
	MPI_Allreduce (&temp,
{'fraction_non_alphanumeric': 0.2054208841251218, 'fraction_numerical': 0.0033010519617396148, 'mean_word_length': 1.74183523253912, 'pattern_counts': {'":': 0, '<': 275, '<?xml version=': 0, '>': 297, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1178435', 'n_tokens_mistral': 15406, 'n_tokens_neox': 14293, 'n_words': 4351}