WriteHDF5Parallel.cpp

Go to the documentation of this file.

#undef DEBUG
#undef TIME_DEBUG
 
#include <cstdarg>
#include <ctime>
#include <cstdlib>
 
#include <cstring>
#include <cassert>
 
#include <vector>
#include <set>
#include <map>
#include <utility>
#include <iostream>
#include <sstream>
#include <string>
 
#include "moab/Interface.hpp"
#include "Internals.hpp"
#include "MBTagConventions.hpp"
#include "MBParallelConventions.h"
#include "moab/ParallelComm.hpp"
#include "moab/CN.hpp"
#include "moab/Range.hpp"
#include "moab/CpuTimer.hpp"
 
#include "WriteHDF5Parallel.hpp"
 
#ifndef MOAB_HAVE_HDF5
#error Attempt to compile WriteHDF5Parallel with HDF5 support disabled
#endif
 
#include <H5Tpublic.h>
#include <H5Ppublic.h>
#include <H5FDmpi.h>
#include <H5FDmpio.h>
 
#include "mhdf.h"
 
#include "IODebugTrack.hpp"
#include "moab/FileOptions.hpp"
 
namespace
{
template < bool Condition >
struct STATIC_ASSERTION;
template <>
struct STATIC_ASSERTION< true >
{
};
} // namespace
 
#define PP_CAT_( a, b ) a##b
#define PP_CAT( a, b ) PP_CAT_( a, b )
#define STATIC_ASSERT( Condition ) \
 enum \
 { \
 PP_CAT( dummy, __LINE__ ) = sizeof( ::STATIC_ASSERTION< (bool)( Condition ) > ) \
 }
 
namespace moab
{
 
#ifndef _WIN32 // problematic for windows
// Need an MPI type that we can put handles in
STATIC_ASSERT( sizeof( unsigned long ) >= sizeof( EntityHandle ) );
 
// Need an MPI type that we can put file IDs in
STATIC_ASSERT( sizeof( unsigned long ) >= sizeof( WriteHDF5::wid_t ) );
#endif
 
// This function doesn't do anything useful. It's just a nice
// place to set a break point to determine why the reader fails.
static inline ErrorCode error( ErrorCode rval )
{
 return rval;
}
 
const char* mpi_err_str( int errorcode )
{
 static char buffer[2048];
 int len = sizeof( buffer );
 MPI_Error_string( errorcode, buffer, &len );
 buffer[std::min( (size_t)len, sizeof( buffer ) - 1 )] = '\0';
 return buffer;
}
 
#define MPI_FAILURE_MSG( A ) \
 "MPI Failure at " __FILE__ ":%d : (Code %d) %s\n", __LINE__, (int)( A ), mpi_err_str( ( A ) )
 
#define CHECK_MPI( A ) \
 do \
 { \
 if( MPI_SUCCESS != ( A ) ) \
 { \
 MB_SET_ERR_CONT( "MPI Failure : (Code " << (int)( A ) << ") " << mpi_err_str( ( A ) ) ); \
 dbgOut.printf( 1, MPI_FAILURE_MSG( ( A ) ) ); \
 return error( MB_FAILURE ); \
 } \
 } while( false )
 
#define MB_FAILURE_MSG( A ) "MOAB_Failure at " __FILE__ ":%d : %s (%d)\n", __LINE__, ErrorCodeStr[( A )], (int)( A )
 
#define CHECK_MB( A ) \
 do \
 { \
 if( MB_SUCCESS != ( A ) ) \
 { \
 MB_SET_ERR_CONT( "MOAB Failure : " << ErrorCodeStr[( A )] ); \
 dbgOut.printf( 1, MB_FAILURE_MSG( ( A ) ) ); \
 return error( A ); \
 } \
 } while( false )
 
#define HDF_FAILURE_MSG( A ) "MHDF Failure at " __FILE__ ":%d : %s\n", __LINE__, mhdf_message( &( A ) )
 
#define CHECK_HDF( A ) \
 do \
 { \
 if( mhdf_isError( &( A ) ) ) \
 { \
 MB_SET_ERR_CONT( "MHDF Failure : " << mhdf_message( &( A ) ) ); \
 dbgOut.printf( 1, HDF_FAILURE_MSG( ( A ) ) ); \
 return error( MB_FAILURE ); \
 } \
 } while( false )
 
#define CHECK_HDFN( A ) \
 do \
 { \
 if( mhdf_isError( &( A ) ) ) \
 { \
 MB_SET_ERR_CONT( "MHDF Failure : " << mhdf_message( &( A ) ) ); \
 return error( MB_FAILURE ); \
 } \
 } while( false )
 
#ifdef VALGRIND
#include <valgrind/memcheck.h>
 
template < typename T >
inline void VALGRIND_MAKE_VEC_UNDEFINED( std::vector< T >& v )
{
 if( v.size() )
 {
 }
 (void)VALGRIND_MAKE_MEM_UNDEFINED( &v[0], v.size() * sizeof( T ) );
}
 
#else
#ifndef VALGRIND_CHECK_MEM_IS_DEFINED
#define VALGRIND_CHECK_MEM_IS_DEFINED( a, b ) ( (void)0 )
#endif
#ifndef VALGRIND_CHECK_MEM_IS_ADDRESSABLE
#define VALGRIND_CHECK_MEM_IS_ADDRESSABLE( a, b ) ( (void)0 )
#endif
#ifndef VALGRIND_MAKE_MEM_UNDEFINED
#define VALGRIND_MAKE_MEM_UNDEFINED( a, b ) ( (void)0 )
#endif
 
template < typename T >
inline void VALGRIND_MAKE_VEC_UNDEFINED( std::vector< T >& )
{
 /* Nothing to do */
}
 
#endif
 
#ifndef NDEBUG
#define START_SERIAL \
 for( unsigned _x = 0; _x < myPcomm->proc_config().proc_size(); ++_x ) \
 { \
 MPI_Barrier( myPcomm->proc_config().proc_comm() ); \
 if( _x != myPcomm->proc_config().proc_rank() ) continue
#define END_SERIAL \
 } \
 MPI_Barrier( myPcomm->proc_config().proc_comm() )
#else
#define START_SERIAL
#define END_SERIAL
#endif
 
static int my_Gatherv( void* sendbuf,
 int sendcount,
 MPI_Datatype sendtype,
 std::vector< unsigned char >& recvbuf,
 std::vector< int >& recvcounts,
 int root,
 MPI_Comm comm )
{
 int nproc, rank, bytes, err;
 MPI_Comm_size( comm, &nproc );
 MPI_Comm_rank( comm, &rank );
 MPI_Type_size( sendtype, &bytes );
 
 recvcounts.resize( rank == root ? nproc : 0 );
 err = MPI_Gather( &sendcount, 1, MPI_INT, &recvcounts[0], 1, MPI_INT, root, comm );
 if( MPI_SUCCESS != err ) return err;
 
 std::vector< int > disp( recvcounts.size() );
 if( root == rank )
 {
 disp[0] = 0;
 for( int i = 1; i < nproc; ++i )
 disp[i] = disp[i - 1] + recvcounts[i - 1];
 recvbuf.resize( bytes * ( disp.back() + recvcounts.back() ) );
 }
 
 return MPI_Gatherv( sendbuf, sendcount, sendtype, &recvbuf[0], &recvcounts[0], &disp[0], sendtype, root, comm );
}
 
static void print_type_sets( Interface* iFace, DebugOutput* str, Range& sets )
{
 const unsigned VB = 2;
 if( str->get_verbosity() < VB ) return;
 
 Tag gid, did, bid, sid, nid;
 gid = iFace->globalId_tag();
 iFace->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, did );
 iFace->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, bid );
 iFace->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, nid );
 iFace->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, sid );
 Range typesets[10];
 const char* typenames[] = { "Block ", "Sideset ", "NodeSet", "Vertex", "Curve",
 "Surface", "Volume", "Body", "Other" };
 for( Range::iterator riter = sets.begin(); riter != sets.end(); ++riter )
 {
 unsigned dim, id; //, oldsize;
 if( MB_SUCCESS == iFace->tag_get_data( bid, &*riter, 1, &id ) )
 dim = 0;
 else if( MB_SUCCESS == iFace->tag_get_data( sid, &*riter, 1, &id ) )
 dim = 1;
 else if( MB_SUCCESS == iFace->tag_get_data( nid, &*riter, 1, &id ) )
 dim = 2;
 else if( MB_SUCCESS == iFace->tag_get_data( did, &*riter, 1, &dim ) )
 {
 id = 0;
 iFace->tag_get_data( gid, &*riter, 1, &id );
 dim += 3;
 }
 else
 {
 id = *riter;
 dim = 9;
 }
 
 // oldsize = typesets[dim].size();
 typesets[dim].insert( id );
 // assert(typesets[dim].size() - oldsize == 1);
 }
 for( int ii = 0; ii < 9; ++ii )
 {
 char tmp[64];
 sprintf( tmp, "%s (%lu) ", typenames[ii], (unsigned long)typesets[ii].size() );
 str->print( VB, tmp, typesets[ii] );
 }
 str->printf( VB, "Total: %lu\n", (unsigned long)sets.size() );
}
 
#define debug_barrier() debug_barrier_line( __LINE__ )
 
void WriteHDF5Parallel::debug_barrier_line( int lineno )
{
 const unsigned threshold = 2;
 static unsigned long count = 0;
 if( dbgOut.get_verbosity() >= threshold && myPcomm )
 {
 dbgOut.printf( threshold, "*********** Debug Barrier %lu (@%d)***********\n", ++count, lineno );
 MPI_Barrier( myPcomm->proc_config().proc_comm() );
 }
}
 
WriterIface* WriteHDF5Parallel::factory( Interface* iface )
{
 return new WriteHDF5Parallel( iface );
}
 
WriteHDF5Parallel::WriteHDF5Parallel( Interface* iface )
 : WriteHDF5( iface ), myPcomm( NULL ), pcommAllocated( false ), hslabOp( H5S_SELECT_OR )
{
}
 
WriteHDF5Parallel::~WriteHDF5Parallel()
{
 if( pcommAllocated && myPcomm ) delete myPcomm;
}
 
// The parent WriteHDF5 class has ExportSet structs that are
// populated with the entities to be written, grouped by type
// (and for elements, connectivity length). This function:
// o determines which entities are to be written by a remote processor
// o removes those entities from the ExportSet structs in WriteMesh
// o passes them back in a Range
ErrorCode WriteHDF5Parallel::gather_interface_meshes( Range& nonowned )
{
 ErrorCode result;
 
 // START_SERIAL;
 dbgOut.print( 3, "Pre-interface mesh:\n" );
 dbgOut.print( 3, nodeSet.range );
 for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
 dbgOut.print( 3, eiter->range );
 dbgOut.print( 3, setSet.range );
 
 // Move handles of non-owned entities from lists of entities
 // that this processor will write to the 'nonowned' list.
 
 nonowned.clear();
 result = myPcomm->filter_pstatus( nodeSet.range, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonowned );
 if( MB_SUCCESS != result ) return error( result );
 nodeSet.range = subtract( nodeSet.range, nonowned );
 
 for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
 {
 Range tmpset;
 result = myPcomm->filter_pstatus( eiter->range, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &tmpset );
 if( MB_SUCCESS != result ) return error( result );
 eiter->range = subtract( eiter->range, tmpset );
 nonowned.merge( tmpset );
 }
 
 dbgOut.print( 3, "Post-interface mesh:\n" );
 dbgOut.print( 3, nodeSet.range );
 for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
 dbgOut.print( 3, eiter->range );
 dbgOut.print( 3, setSet.range );
 
 // END_SERIAL;
 
 return MB_SUCCESS;
}
 
ErrorCode WriteHDF5Parallel::parallel_create_file( const char* filename,
 bool overwrite,
 const std::vector< std::string >& qa_records,
 const FileOptions& opts,
 const Tag* user_tag_list,
 int user_tag_count,
 int dimension,
 double* times )
{
 ErrorCode rval;
 mhdf_Status status;
 
 int pcomm_no = 0;
 opts.get_int_option( "PARALLEL_COMM", pcomm_no );
 
 myPcomm = ParallelComm::get_pcomm( iFace, pcomm_no );
 if( 0 == myPcomm )
 {
 myPcomm = new ParallelComm( iFace, MPI_COMM_WORLD );
 pcommAllocated = true;
 }
 
 MPI_Info info = MPI_INFO_NULL;
 std::string cb_size;
 rval = opts.get_str_option( "CB_BUFFER_SIZE", cb_size );
 if( MB_SUCCESS == rval )
 {
 MPI_Info_create( &info );
 MPI_Info_set( info, const_cast< char* >( "cb_buffer_size" ), const_cast< char* >( cb_size.c_str() ) );
 }
 
 dbgOut.set_rank( myPcomm->proc_config().proc_rank() );
 dbgOut.limit_output_to_first_N_procs( 32 );
 Range nonlocal;
 debug_barrier();
 dbgOut.tprint( 1, "Gathering interface meshes\n" );
 rval = gather_interface_meshes( nonlocal );
 if( MB_SUCCESS != rval ) return error( rval );
 
 /**************** Get tag names for sets likely to be shared ***********/
 // debug_barrier();
 // dbgOut.tprint(1, "Getting shared entity sets\n");
 // rval = get_sharedset_tags();
 // if (MB_SUCCESS != rval) return error(rval);
 
 /**************** Create actual file and write meta info ***************/
 
 debug_barrier();
 if( myPcomm->proc_config().proc_rank() == 0 )
 {
 dbgOut.tprintf( 1, "Creating file: %s\n", filename );
 
 // Create the file
 const char* type_names[MBMAXTYPE];
 memset( type_names, 0, MBMAXTYPE * sizeof( char* ) );
 for( EntityType i = MBEDGE; i < MBENTITYSET; ++i )
 type_names[i] = CN::EntityTypeName( i );
 
 dbgOut.tprint( 1, "call mhdf_createFile\n" );
 filePtr = mhdf_createFile( filename, overwrite, type_names, MBMAXTYPE, id_type, &status );
 if( !filePtr )
 {
 MB_SET_ERR( MB_FAILURE, mhdf_message( &status ) );
 }
 
 dbgOut.tprint( 1, "call write_qa\n" );
 rval = write_qa( qa_records );
 if( MB_SUCCESS != rval ) return error( rval );
 }
 
 /**************** Create node coordinate table ***************/
 CpuTimer timer;
 debug_barrier();
 dbgOut.tprint( 1, "creating node table\n" );
 topState.start( "creating node table" );
 rval = create_node_table( dimension );
 topState.end( rval );
 if( MB_SUCCESS != rval ) return error( rval );
 if( times ) times[CREATE_NODE_TIME] = timer.time_elapsed();
 
 /**************** Create element tables ***************/
 
 debug_barrier();
 dbgOut.tprint( 1, "negotiating element types\n" );
 topState.start( "negotiating element types" );
 rval = negotiate_type_list();
 topState.end( rval );
 if( MB_SUCCESS != rval ) return error( rval );
 if( times ) times[NEGOTIATE_TYPES_TIME] = timer.time_elapsed();
 dbgOut.tprint( 1, "creating element table\n" );
 topState.start( "creating element tables" );
 rval = create_element_tables();
 topState.end( rval );
 if( MB_SUCCESS != rval ) return error( rval );
 if( times ) times[CREATE_ELEM_TIME] = timer.time_elapsed();
 
 /*************** Exchange file IDs *****************/
 
 debug_barrier();
 dbgOut.tprint( 1, "communicating file ids\n" );
 topState.start( "communicating file ids" );
 rval = exchange_file_ids( nonlocal );
 topState.end( rval );
 if( MB_SUCCESS != rval ) return error( rval );
 if( times ) times[FILEID_EXCHANGE_TIME] = timer.time_elapsed();
 
 /**************** Create meshset tables *********************/
 
 debug_barrier();
 dbgOut.tprint( 1, "creating meshset table\n" );
 topState.start( "creating meshset tables" );
 rval = create_meshset_tables( times );
 topState.end( rval );
 if( MB_SUCCESS != rval ) return error( rval );
 if( times ) times[CREATE_SET_TIME] = timer.time_elapsed();
 
 /**************** Create adjacency tables *********************/
 
 debug_barrier();
 dbgOut.tprint( 1, "creating adjacency table\n" );
 topState.start( "creating adjacency tables" );
 rval = create_adjacency_tables();
 topState.end( rval );
 if( MB_SUCCESS != rval ) return error( rval );
 if( times ) times[CREATE_ADJ_TIME] = timer.time_elapsed();
 
 /**************** Create tag data *********************/
 
 debug_barrier();
 dbgOut.tprint( 1, "creating tag tables\n" );
 topState.start( "creating tag tables" );
 rval = gather_tags( user_tag_list, user_tag_count );
 if( MB_SUCCESS != rval ) return error( rval );
 rval = create_tag_tables();
 topState.end( rval );
 if( MB_SUCCESS != rval ) return error( rval );
 if( times ) times[CREATE_TAG_TIME] = timer.time_elapsed();
 
 /************** Close serial file and reopen parallel *****************/
 
 if( 0 == myPcomm->proc_config().proc_rank() ) mhdf_closeFile( filePtr, &status );
 
 MPI_Barrier( myPcomm->proc_config().proc_comm() );
 dbgOut.tprint( 1, "(re)opening file in parallel mode\n" );
 unsigned long junk;
 hid_t hdf_opt = H5Pcreate( H5P_FILE_ACCESS );
 H5Pset_fapl_mpio( hdf_opt, myPcomm->proc_config().proc_comm(), info );
 filePtr = mhdf_openFileWithOpt( filename, 1, &junk, id_type, hdf_opt, &status );
 H5Pclose( hdf_opt );
 if( !filePtr )
 {
 MB_SET_ERR( MB_FAILURE, mhdf_message( &status ) );
 }
 
 if( collectiveIO )
 {
 dbgOut.print( 1, "USING COLLECTIVE IO\n" );
 writeProp = H5Pcreate( H5P_DATASET_XFER );
 H5Pset_dxpl_mpio( writeProp, H5FD_MPIO_COLLECTIVE );
 }
 
 /* Test if we can use H5S_APPEND when selecting hyperslabs */
 if( MB_SUCCESS != opts.get_null_option( "HYPERSLAB_OR" ) &&
 ( MB_SUCCESS == opts.get_null_option( "HYPERSLAB_APPEND" ) || HDF5_can_append_hyperslabs() ) )
 {
 dbgOut.print( 1, "HDF5 library supports H5Sselect_hyperlsab with H5S_SELECT_APPEND\n" );
 hslabOp = H5S_SELECT_APPEND;
 }
 
 dbgOut.tprint( 1, "Exiting parallel_create_file\n" );
 return MB_SUCCESS;
}
 
class TagNameCompare
{
 Interface* iFace;
 std::string name1, name2;
 
 public:
 TagNameCompare( Interface* iface ) : iFace( iface ) {}
 bool operator()( const WriteHDF5::TagDesc& t1, const WriteHDF5::TagDesc& t2 );
};
 
bool TagNameCompare::operator()( const WriteHDF5::TagDesc& t1, const WriteHDF5::TagDesc& t2 )
{
 iFace->tag_get_name( t1.tag_id, name1 );
 iFace->tag_get_name( t2.tag_id, name2 );
 return name1 < name2;
}
 
struct serial_tag_data
{
 TagType storage;
 DataType type;
 int size;
 int name_len;
 int def_val_len;
 char name[sizeof( unsigned long )];
 
 static size_t pad( size_t len )
 {
 if( len % sizeof( unsigned long ) )
 return len + sizeof( unsigned long ) - len % sizeof( unsigned long );
 else
 return len;
 }
 
 static size_t def_val_bytes( int def_val_len, DataType type )
 {
 switch( type )
 {
 case MB_TYPE_BIT:
 return def_val_len ? 1 : 0;
 case MB_TYPE_OPAQUE:
 return def_val_len;
 case MB_TYPE_INTEGER:
 return def_val_len * sizeof( int );
 case MB_TYPE_DOUBLE:
 return def_val_len * sizeof( double );
 case MB_TYPE_HANDLE:
 return def_val_len * sizeof( EntityHandle );
 }
 return 0;
 }
 
 static size_t len( int name_len, int def_val_len, DataType type )
 {
 return sizeof( serial_tag_data ) + pad( name_len + def_val_bytes( def_val_len, type ) ) -
 sizeof( unsigned long );
 }
 size_t len() const
 {
 return len( name_len, def_val_len, type );
 }
 void* default_value()
 {
 return def_val_len ? name + name_len : 0;
 }
 const void* default_value() const
 {
 return const_cast< serial_tag_data* >( this )->default_value();
 }
 void set_default_value( const void* val )
 {
 memcpy( default_value(), val, def_val_bytes( def_val_len, type ) );
 }
};
 
ErrorCode WriteHDF5Parallel::append_serial_tag_data( std::vector< unsigned char >& buffer,
 const WriteHDF5::TagDesc& tag )
{
 ErrorCode rval;
 
 std::string name;
 rval = iFace->tag_get_name( tag.tag_id, name );
 if( MB_SUCCESS != rval ) return error( rval );
 
 // Get name length, including space for null char
 size_t name_len = name.size() + 1;
 if( name_len == 1 ) return MB_SUCCESS; // Skip tags with no name
 
 DataType data_type;
 rval = iFace->tag_get_data_type( tag.tag_id, data_type );
 if( MB_SUCCESS != rval ) return error( rval );
 
 // Get default value
 int def_val_len;
 const void* def_val;
 if( MB_SUCCESS != iFace->tag_get_default_value( tag.tag_id, def_val, def_val_len ) )
 {
 def_val_len = 0;
 def_val = 0;
 }
 
 // Allocate struct within buffer
 size_t init_size = buffer.size();
 buffer.resize( init_size + serial_tag_data::len( name_len, def_val_len, data_type ) );
 serial_tag_data* ptr = reinterpret_cast< serial_tag_data* >( &buffer[init_size] );
 
 // Populate struct
 rval = iFace->tag_get_type( tag.tag_id, ptr->storage );
 if( MB_SUCCESS != rval ) return error( rval );
 ptr->type = data_type;
 rval = iFace->tag_get_length( tag.tag_id, ptr->size );
 if( MB_VARIABLE_DATA_LENGTH == rval )
 ptr->size = MB_VARIABLE_LENGTH;
 else if( MB_SUCCESS != rval )
 return error( rval );
 ptr->name_len = name_len;
 Range range;
 memset( ptr->name, 0, ptr->name_len );
 memcpy( ptr->name, name.data(), name.size() );
 ptr->def_val_len = def_val_len;
 ptr->set_default_value( def_val );
 
 return MB_SUCCESS;
}
 
ErrorCode WriteHDF5Parallel::check_serial_tag_data( const std::vector< unsigned char >& buffer,
 std::vector< TagDesc* >* missing,
 std::vector< TagDesc* >* newlist )
{
 ErrorCode rval;
 
 // Use 'write_sparse' field as a 'visited' mark
 std::list< TagDesc >::iterator tag_iter;
 if( missing )
 for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
 tag_iter->write_sparse = true;
 
 // Use a set as a temporary for what will ultimately go in
 // newlist because we need to pass back newlist in the order
 // of the tagList member.
 std::set< TagDesc* > newset;
 
 // Iterate over data from, updating the local list of tags.
 // Be careful to keep tagList sorted such that in the end all
 // procs have the same list in the same order.
 std::vector< unsigned char >::const_iterator diter = buffer.begin();
 tag_iter = tagList.begin();
 while( diter < buffer.end() )
 {
 // Get struct from buffer
 const serial_tag_data* ptr = reinterpret_cast< const serial_tag_data* >( &*diter );
 
 // Find local struct for tag
 std::string name( ptr->name );
 std::string n;
 iFace->tag_get_name( tag_iter->tag_id, n ); // Second time we've called, so shouldn't fail
 if( n > name )
 {
 tag_iter = tagList.begin(); // New proc, start search from beginning
 }
 iFace->tag_get_name( tag_iter->tag_id, n );
 while( n < name )
 {
 ++tag_iter;
 if( tag_iter == tagList.end() ) break;
 iFace->tag_get_name( tag_iter->tag_id, n );
 }
 if( tag_iter == tagList.end() || n != name )
 { // New tag
 TagDesc newtag;
 
 if( ptr->size == MB_VARIABLE_LENGTH )
 rval = iFace->tag_get_handle( name.c_str(), ptr->def_val_len, ptr->type, newtag.tag_id,
 MB_TAG_VARLEN | MB_TAG_CREAT | ptr->storage, ptr->default_value() );
 else
 rval = iFace->tag_get_handle( name.c_str(), ptr->size, ptr->type, newtag.tag_id,
 MB_TAG_CREAT | ptr->storage, ptr->default_value() );
 if( MB_SUCCESS != rval ) return error( rval );
 
 newtag.sparse_offset = 0;
 newtag.var_data_offset = 0;
 newtag.write_sparse = false;
 newtag.max_num_ents = 0;
 newtag.max_num_vals = 0;
 
 tag_iter = tagList.insert( tag_iter, newtag );
 if( newlist ) newset.insert( &*tag_iter );
 }
 else
 { // Check that tag is as expected
 DataType type;
 iFace->tag_get_data_type( tag_iter->tag_id, type );
 if( type != ptr->type )
 {
 MB_SET_ERR( MB_FAILURE, "Processes have inconsistent data type for tag \"" << name << "\"" );
 }
 int size;
 iFace->tag_get_length( tag_iter->tag_id, size );
 if( size != ptr->size )
 {
 MB_SET_ERR( MB_FAILURE, "Processes have inconsistent size for tag \"" << name << "\"" );
 }
 tag_iter->write_sparse = false;
 }
 
 // Step to next variable-length struct.
 diter += ptr->len();
 }
 
 // Now pass back any local tags that weren't in the buffer
 if( missing )
 {
 for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
 {
 if( tag_iter->write_sparse )
 {
 tag_iter->write_sparse = false;
 missing->push_back( &*tag_iter );
 }
 }
 }
 
 // Be careful to populate newlist in the same, sorted, order as tagList
 if( newlist )
 {
 for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
 if( newset.find( &*tag_iter ) != newset.end() ) newlist->push_back( &*tag_iter );
 }
 
 return MB_SUCCESS;
}
 
static void set_bit( int position, unsigned char* bytes )
{
 int byte = position / 8;
 int bit = position % 8;
 bytes[byte] |= ( ( (unsigned char)1 ) << bit );
}
 
static bool get_bit( int position, const unsigned char* bytes )
{
 int byte = position / 8;
 int bit = position % 8;
 return 0 != ( bytes[byte] & ( ( (unsigned char)1 ) << bit ) );
}
 
ErrorCode WriteHDF5Parallel::create_tag_tables()
{
 std::list< TagDesc >::iterator tag_iter;
 ErrorCode rval;
 int err;
 const int rank = myPcomm->proc_config().proc_rank();
 const MPI_Comm comm = myPcomm->proc_config().proc_comm();
 
 subState.start( "negotiating tag list" );
 
 dbgOut.tprint( 1, "communicating tag metadata\n" );
 
 dbgOut.printf( 2, "Exchanging tag data for %d tags.\n", (int)tagList.size() );
 
 // Sort tagList contents in alphabetical order by tag name
 tagList.sort( TagNameCompare( iFace ) );
 
 // Negotiate total list of tags to write
 
 // Build concatenated list of all tag data
 std::vector< unsigned char > tag_buffer;
 for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
 {
 rval = append_serial_tag_data( tag_buffer, *tag_iter );
 CHECK_MB( rval );
 }
 
 // Broadcast list from root to all other procs
 unsigned long size = tag_buffer.size();
 err = MPI_Bcast( &size, 1, MPI_UNSIGNED_LONG, 0, comm );
 CHECK_MPI( err );
 tag_buffer.resize( size );
 err = MPI_Bcast( &tag_buffer[0], size, MPI_UNSIGNED_CHAR, 0, comm );
 CHECK_MPI( err );
 
 // Update local tag list
 std::vector< TagDesc* > missing;
 rval = check_serial_tag_data( tag_buffer, &missing, 0 );
 CHECK_MB( rval );
 
 // Check if we're done (0->done, 1->more, 2+->error)
 int code, lcode = ( MB_SUCCESS != rval ) ? rval + 2 : missing.empty() ? 0 : 1;
 err = MPI_Allreduce( &lcode, &code, 1, MPI_INT, MPI_MAX, comm );
 CHECK_MPI( err );
 if( code > 1 )
 {
 MB_SET_ERR_CONT( "Inconsistent tag definitions between procs" );
 return error( (ErrorCode)( code - 2 ) );
 }
 
 // If not done...
 if( code )
 {
 dbgOut.print( 1, "Not all procs had same tag definitions, negotiating...\n" );
 
 // Get tags defined on this proc but not on root proc
 tag_buffer.clear();
 for( size_t i = 0; i < missing.size(); ++i )
 {
 rval = append_serial_tag_data( tag_buffer, *missing[i] );
 CHECK_MB( rval );
 }
 
 // Gather extra tag definitions on root processor
 std::vector< int > junk; // don't care how many from each proc
 assert( rank || tag_buffer.empty() ); // must be empty on root
 err = my_Gatherv( &tag_buffer[0], tag_buffer.size(), MPI_UNSIGNED_CHAR, tag_buffer, junk, 0, comm );
 CHECK_MPI( err );
 
 // Process serialized tag descriptions on root, and
 rval = MB_SUCCESS;
 if( 0 == rank )
 {
 // Process serialized tag descriptions on root, and
 std::vector< TagDesc* > newlist;
 rval = check_serial_tag_data( tag_buffer, 0, &newlist );
 tag_buffer.clear();
 // re-serialize a unique list of new tag definitions
 for( size_t i = 0; MB_SUCCESS == rval && i != newlist.size(); ++i )
 {
 rval = append_serial_tag_data( tag_buffer, *newlist[i] );
 CHECK_MB( rval );
 }
 }
 
 // Broadcast any new tag definitions from root to other procs
 long this_size = tag_buffer.size();
 if( MB_SUCCESS != rval ) this_size = -rval;
 err = MPI_Bcast( &this_size, 1, MPI_LONG, 0, comm );
 CHECK_MPI( err );
 if( this_size < 0 )
 {
 MB_SET_ERR_CONT( "Inconsistent tag definitions between procs" );
 return error( (ErrorCode)-this_size );
 }
 tag_buffer.resize( this_size );
 err = MPI_Bcast( &tag_buffer[0], this_size, MPI_UNSIGNED_CHAR, 0, comm );
 CHECK_MPI( err );
 
 // Process new tag definitions
 rval = check_serial_tag_data( tag_buffer, 0, 0 );
 CHECK_MB( rval );
 }
 
 subState.end();
 subState.start( "negotiate which element/tag combinations are dense" );
 
 // Figure out for which tag/element combinations we can
 // write dense tag data.
 
 // Construct a table of bits,
 // where each row of the table corresponds to a tag
 // and each column to an element group.
 
 // Two extra, because first is nodes and last is sets.
 // (n+7)/8 is ceil(n/8)
 const int bytes_per_tag = ( exportList.size() + 9 ) / 8;
 std::vector< unsigned char > data( bytes_per_tag * tagList.size(), 0 );
 std::vector< unsigned char > recv( data.size(), 0 );
 unsigned char* iter = &data[0];
 if( writeTagDense && !data.empty() )
 {
 for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, iter += bytes_per_tag )
 {
 
 Range tagged;
 rval = get_sparse_tagged_entities( *tag_iter, tagged );
 CHECK_MB( rval );
 
 int s;
 if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) ) continue;
 
 std::string n;
 iFace->tag_get_name( tag_iter->tag_id,
 n ); // Second time we've called, so shouldn't fail
 
 // Check if we want to write this tag in dense format even if not
 // all of the entities have a tag value. The criterion of this
 // is that the tag be dense, have a default value, and have at
 // least 2/3 of the entities tagged.
 bool prefer_dense = false;
 TagType type;
 rval = iFace->tag_get_type( tag_iter->tag_id, type );
 CHECK_MB( rval );
 if( MB_TAG_DENSE == type )
 {
 const void* defval = 0;
 rval = iFace->tag_get_default_value( tag_iter->tag_id, defval, s );
 if( MB_SUCCESS == rval ) prefer_dense = true;
 }
 
 int i = 0;
 if( check_dense_format_tag( nodeSet, tagged, prefer_dense ) )
 {
 set_bit( i, iter );
 dbgOut.printf( 2, "Can write dense data for \"%s\"/Nodes\n", n.c_str() );
 }
 std::list< ExportSet >::const_iterator ex_iter = exportList.begin();
 for( ++i; ex_iter != exportList.end(); ++i, ++ex_iter )
 {
 // when writing in parallel, on some partitions, some of these element ranges might
 // be empty so do not turn this tag as sparse, just because of that, leave it dense,
 // if we prefer dense
 if( ( prefer_dense && ex_iter->range.empty() ) ||
 check_dense_format_tag( *ex_iter, tagged, prefer_dense ) )
 {
 set_bit( i, iter );
 dbgOut.printf( 2, "Can write dense data for \"%s\"/%s\n", n.c_str(), ex_iter->name() );
 }
 }
 if( check_dense_format_tag( setSet, tagged, prefer_dense ) )
 {
 set_bit( i, iter );
 dbgOut.printf( 2, "Can write dense data for \"%s\"/Sets\n", n.c_str() );
 }
 }
 
 // Do bit-wise AND of list over all processors (only write dense format
 // if all proccesses want dense format for this group of entities).
 err = MPI_Allreduce( &data[0], &recv[0], data.size(), MPI_UNSIGNED_CHAR, MPI_BAND,
 myPcomm->proc_config().proc_comm() );
 CHECK_MPI( err );
 } // if (writeTagDense)
 
 // Store initial counts for sparse-formatted tag data.
 // The total number of values to send and receive will be the number of
 // tags plus the number of var-len tags because we need to negotiate
 // offsets into two different tables for the var-len tags.
 std::vector< long > counts;
 
 // Record dense tag/element combinations
 iter = &recv[0];
 const unsigned char* iter2 = &data[0];
 for( tag_iter = tagList.begin(); tag_iter != tagList.end();
 ++tag_iter, iter += bytes_per_tag, iter2 += bytes_per_tag )
 {
 
 Range tagged;
 rval = get_sparse_tagged_entities( *tag_iter, tagged );
 CHECK_MB( rval );
 
 std::string n;
 iFace->tag_get_name( tag_iter->tag_id, n ); // Second time we've called, so shouldn't fail
 
 int i = 0;
 if( get_bit( i, iter ) )
 {
 assert( get_bit( i, iter2 ) );
 tag_iter->dense_list.push_back( nodeSet );
 tagged -= nodeSet.range;
 dbgOut.printf( 2, "Will write dense data for \"%s\"/Nodes\n", n.c_str() );
 }
 std::list< ExportSet >::const_iterator ex_iter = exportList.begin();
 for( ++i; ex_iter != exportList.end(); ++i, ++ex_iter )
 {
 if( get_bit( i, iter ) )
 {
 assert( get_bit( i, iter2 ) );
 tag_iter->dense_list.push_back( *ex_iter );
 dbgOut.printf( 2, "WIll write dense data for \"%s\"/%s\n", n.c_str(), ex_iter->name() );
 tagged -= ex_iter->range;
 }
 }
 if( get_bit( i, iter ) )
 {
 assert( get_bit( i, iter2 ) );
 tag_iter->dense_list.push_back( setSet );
 dbgOut.printf( 2, "Will write dense data for \"%s\"/Sets\n", n.c_str() );
 tagged -= setSet.range;
 }
 
 counts.push_back( tagged.size() );
 
 int s;
 if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
 {
 unsigned long data_len;
 rval = get_tag_data_length( *tag_iter, tagged, data_len );
 CHECK_MB( rval );
 counts.push_back( data_len );
 }
 }
 
 subState.end();
 subState.start( "Negotiate offsets for sparse tag info" );
 
 std::vector< long > offsets( counts.size() ), maxima( counts.size() ), totals( counts.size() );
 rval = create_dataset( counts.size(), &counts[0], &offsets[0], &maxima[0], &totals[0] );
 CHECK_MB( rval );
 
 // Copy values into local structs and if root then create tables
 size_t idx = 0;
 for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, ++idx )
 {
 assert( idx < counts.size() );
 tag_iter->sparse_offset = offsets[idx];
 tag_iter->max_num_ents = maxima[idx];
 tag_iter->write_sparse = ( 0 != totals[idx] );
 int s;
 if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
 {
 ++idx;
 assert( idx < counts.size() );
 tag_iter->var_data_offset = offsets[idx];
 tag_iter->max_num_vals = maxima[idx];
 }
 else
 {
 tag_iter->var_data_offset = 0;
 tag_iter->max_num_vals = 0;
 }
 }
 
 subState.end();
 
 // Create tag tables on root process
 if( 0 == myPcomm->proc_config().proc_rank() )
 {
 size_t iidx = 0;
 for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, ++iidx )
 {
 assert( iidx < totals.size() );
 unsigned long num_ents = totals[iidx];
 unsigned long num_val = 0;
 int s;
 if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
 {
 ++iidx;
 assert( iidx < totals.size() );
 num_val = totals[iidx];
 }
 dbgOut.printf( 2, "Writing tag description for tag 0x%lx with %lu values\n",
 (unsigned long)tag_iter->tag_id, num_val ? num_val : num_ents );
 
 rval = create_tag( *tag_iter, num_ents, num_val );
 if( MB_SUCCESS != rval ) return error( rval );
 }
 }
 
 if( dbgOut.get_verbosity() > 1 )
 {
 dbgOut.printf( 2, "Tags: %12s %8s %8s %8s %8s %8s\n", "Name", "Count", "Offset", "Var Off", "Max Ent",
 "Handle" );
 
 for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
 {
 std::string name;
 iFace->tag_get_name( tag_iter->tag_id, name );
 size_t this_size;
 get_num_sparse_tagged_entities( *tag_iter, this_size );
 dbgOut.printf( 2, "%18s %8lu %8lu %8lu %8lu 0x%7lx\n", name.c_str(), (unsigned long)this_size,
 (unsigned long)tag_iter->sparse_offset, (unsigned long)tag_iter->var_data_offset,
 (unsigned long)tag_iter->max_num_ents, (unsigned long)tag_iter->tag_id );
 }
 }
 
 return MB_SUCCESS;
}
 
struct DatasetVals
{
 long start_id;
 long max_count;
 long total;
};
STATIC_ASSERT( ( sizeof( DatasetVals ) == 3 * sizeof( long ) ) );
 
ErrorCode WriteHDF5Parallel::create_dataset( int num_datasets,
 const long* num_owned,
 long* offsets_out,
 long* max_proc_entities,
 long* total_entities,
 const DataSetCreator& creator,
 ExportSet* groups[],
 wid_t* first_ids_out )
{
 int result;
 ErrorCode rval;
 const unsigned rank = myPcomm->proc_config().proc_rank();
 const unsigned nproc = myPcomm->proc_config().proc_size();
 const MPI_Comm comm = myPcomm->proc_config().proc_comm();
 
 // Gather entity counts for each processor on root
 std::vector< long > counts( rank ? 0 : nproc * num_datasets );
 (void)VALGRIND_CHECK_MEM_IS_DEFINED( &num_owned, sizeof( long ) );
 result = MPI_Gather( const_cast< long* >( num_owned ), num_datasets, MPI_LONG, &counts[0], num_datasets, MPI_LONG,
 0, comm );
 CHECK_MPI( result );
 
 // Create node data in file
 DatasetVals zero_val = { 0, 0, 0 };
 std::vector< DatasetVals > cumulative( num_datasets, zero_val );
 if( rank == 0 )
 {
 for( unsigned i = 0; i < nproc; i++ )
 {
 const long* proc_data = &counts[i * num_datasets];
 for( int index = 0; index < num_datasets; ++index )
 {
 cumulative[index].total += proc_data[index];
 if( proc_data[index] > cumulative[index].max_count ) cumulative[index].max_count = proc_data[index];
 }
 }
 
 for( int index = 0; index < num_datasets; ++index )
 {
 if( cumulative[index].total )
 {
 rval = creator( this, cumulative[index].total, groups ? groups[index] : 0, cumulative[index].start_id );
 CHECK_MB( rval );
 }
 else
 {
 cumulative[index].start_id = -1;
 }
 }
 }
 
 // Send id offset to every proc
 result = MPI_Bcast( (void*)&cumulative[0], 3 * num_datasets, MPI_LONG, 0, comm );
 CHECK_MPI( result );
 for( int index = 0; index < num_datasets; ++index )
 {
 if( first_ids_out ) first_ids_out[index] = (wid_t)cumulative[index].start_id;
 max_proc_entities[index] = cumulative[index].max_count;
 total_entities[index] = cumulative[index].total;
 }
 
 // Convert array of per-process counts to per-process offsets
 if( rank == 0 )
 {
 // Initialize prev_size with data sizes for root process
 std::vector< long > prev_size( counts.begin(), counts.begin() + num_datasets );
 // Root process gets offset zero
 std::fill( counts.begin(), counts.begin() + num_datasets, 0L );
 // For each proc other than this one (root)
 for( unsigned i = 1; i < nproc; ++i )
 {
 // Get pointer to offsets for previous process in list
 long* prev_data = &counts[( i - 1 ) * num_datasets];
 // Get pointer to offsets for this process in list
 long* proc_data = &counts[i * num_datasets];
 // For each data set
 for( int j = 0; j < num_datasets; ++j )
 {
 // Get size of data in dataset from process i
 long mysize = proc_data[j];
 // Offset for process i is offset of previous process plus
 // number of values previous process will write
 proc_data[j] = prev_data[j] + prev_size[j];
 // Store my size, as it is no longer available in 'counts'
 prev_size[j] = mysize;
 }
 }
 }
 
 // Send each proc it's offset in the table
 if( rank == 0 )
 {
 (void)VALGRIND_CHECK_MEM_IS_DEFINED( &counts[0], num_datasets * nproc * sizeof( long ) );
 }
 result = MPI_Scatter( &counts[0], num_datasets, MPI_LONG, offsets_out, num_datasets, MPI_LONG, 0, comm );
 CHECK_MPI( result );
 
 return MB_SUCCESS;
}
 
ErrorCode WriteHDF5Parallel::create_node_table( int dimension )
{
 nodeSet.num_nodes = dimension; // Put it here so NodeSetCreator can access it
 struct NodeSetCreator : public DataSetCreator
 {
 ErrorCode operator()( WriteHDF5* file, long count, const ExportSet* group, long& start_id ) const
 {
 mhdf_Status status;
 hid_t handle = mhdf_createNodeCoords( file->file_ptr(), group->num_nodes, count, &start_id, &status );
 CHECK_HDFN( status );
 mhdf_closeData( file->file_ptr(), handle, &status );
 CHECK_HDFN( status );
 return MB_SUCCESS;
 }
 };
 
 const long count = nodeSet.range.size();
 ExportSet* array[] = { &nodeSet };
 ErrorCode rval = create_dataset( 1, &count, &nodeSet.offset, &nodeSet.max_num_ents, &nodeSet.total_num_ents,
 NodeSetCreator(), array, &nodeSet.first_id );
 CHECK_MB( rval );
 return assign_ids( nodeSet.range, nodeSet.first_id + nodeSet.offset );
}
 
struct elemtype
{
 int mbtype;
 int numnode;
 
 elemtype( int vals[2] ) : mbtype( vals[0] ), numnode( vals[1] ) {}
 elemtype( int t, int n ) : mbtype( t ), numnode( n ) {}
 
 bool operator==( const elemtype& other ) const
 {
 return mbtype == other.mbtype && ( mbtype == MBENTITYSET || numnode == other.numnode );
 }
 bool operator<( const elemtype& other ) const
 {
 if( mbtype > other.mbtype ) return false;
 
 return mbtype < other.mbtype || ( mbtype != MBENTITYSET && numnode < other.numnode );
 }
 bool operator!=( const elemtype& other ) const
 {
 return !this->operator==( other );
 }
};
 
ErrorCode WriteHDF5Parallel::negotiate_type_list()
{
 int result;
 const MPI_Comm comm = myPcomm->proc_config().proc_comm();
 
 exportList.sort();
 int num_types = exportList.size();
 
 // Get list of types on this processor
 typedef std::vector< std::pair< int, int > > typelist;
 typelist my_types( num_types );
 (void)VALGRIND_MAKE_VEC_UNDEFINED( my_types );
 typelist::iterator viter = my_types.begin();
 for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
 {
 viter->first = eiter->type;
 viter->second = eiter->num_nodes;
 ++viter;
 }
 
 dbgOut.print( 2, "Local Element Types:\n" );
 for( viter = my_types.begin(); viter != my_types.end(); ++viter )
 {
 int type = viter->first;
 int count = viter->second;
 dbgOut.printf( 2, " %s : %d\n", CN::EntityTypeName( (EntityType)type ), count );
 }
 
 // Broadcast number of types from root to all nodes
 int num_types0 = num_types;
 result = MPI_Bcast( &num_types0, 1, MPI_INT, 0, comm );
 CHECK_MPI( result );
 // Broadcast type list from root to all nodes
 typelist root_types( num_types0 );
 if( 0 == myPcomm->proc_config().proc_rank() ) root_types = my_types;
 result = MPI_Bcast( (void*)&root_types[0], 2 * num_types0, MPI_INT, 0, comm );
 CHECK_MPI( result );
 
 // Build local list of any types that root did not know about
 typelist non_root_types;
 viter = root_types.begin();
 for( typelist::iterator iter = my_types.begin(); iter != my_types.end(); ++iter )
 {
 if( viter == root_types.end() || *viter != *iter )
 non_root_types.push_back( *iter );
 else
 ++viter;
 }
 
 // Determine if any process had types not defined on the root
 int non_root_count = non_root_types.size();
 int not_done;
 result = MPI_Allreduce( &non_root_count, &not_done, 1, MPI_INT, MPI_LOR, comm );
 CHECK_MPI( result );
 if( not_done )
 {
 // Get number of types each processor has that root does not
 std::vector< int > counts( myPcomm->proc_config().proc_size() );
 int two_count = 2 * non_root_count;
 result = MPI_Gather( &two_count, 1, MPI_INT, &counts[0], 1, MPI_INT, 0, comm );
 CHECK_MPI( result );
 
 // Get list of types from each processor
 std::vector< int > displs( myPcomm->proc_config().proc_size() + 1 );
 (void)VALGRIND_MAKE_VEC_UNDEFINED( displs );
 displs[0] = 0;
 for( unsigned long i = 1; i <= myPcomm->proc_config().proc_size(); ++i )
 displs[i] = displs[i - 1] + counts[i - 1];
 int total = displs[myPcomm->proc_config().proc_size()];
 typelist alltypes( total / 2 );
 (void)VALGRIND_MAKE_VEC_UNDEFINED( alltypes );
 (void)VALGRIND_CHECK_MEM_IS_DEFINED( &non_root_types[0], non_root_types.size() * sizeof( int ) );
 result = MPI_Gatherv( (void*)&non_root_types[0], 2 * non_root_count, MPI_INT, (int*)&alltypes[0], &counts[0],
 &displs[0], MPI_INT, 0, comm );
 CHECK_MPI( result );
 
 // Merge type lists.
 // Prefer O(n) insertions with O(ln n) search time because
 // we expect data from a potentially large number of processes,
 // but with a small total number of element types.
 if( 0 == myPcomm->proc_config().proc_rank() )
 {
 for( viter = alltypes.begin(); viter != alltypes.end(); ++viter )
 {
 typelist::iterator titer = std::lower_bound( my_types.begin(), my_types.end(), *viter );
 if( titer == my_types.end() || *titer != *viter ) my_types.insert( titer, *viter );
 }
 
 dbgOut.print( 2, "Global Element Types:\n" );
 for( viter = my_types.begin(); viter != my_types.end(); ++viter )
 dbgOut.printf( 2, " %s : %d\n", CN::EntityTypeName( (EntityType)viter->first ), viter->second );
 }
 
 // Send total number of types to each processor
 total = my_types.size();
 result = MPI_Bcast( &total, 1, MPI_INT, 0, comm );
 CHECK_MPI( result );
 
 // Send list of types to each processor
 my_types.resize( total );
 result = MPI_Bcast( (void*)&my_types[0], 2 * total, MPI_INT, 0, comm );
 CHECK_MPI( result );
 }
 else
 {
 // Special case: if root had types but some subset of procs did not
 // have those types, but there are no types that the root doesn't
 // know about then we still need to update processes that are missing
 // types.
 my_types.swap( root_types );
 }
 
 // Insert missing types into exportList, with an empty
 // range of entities to export.
 std::list< ExportSet >::iterator ex_iter = exportList.begin();
 for( viter = my_types.begin(); viter != my_types.end(); ++viter )
 {
 while( ex_iter != exportList.end() && *ex_iter < *viter )
 ++ex_iter;
 
 if( ex_iter == exportList.end() || !( *ex_iter == *viter ) )
 {
 ExportSet insert;
 insert.type = (EntityType)viter->first;
 insert.num_nodes = viter->second;
 insert.first_id = 0;
 insert.offset = 0;
 insert.adj_offset = 0;
 ex_iter = exportList.insert( ex_iter, insert );
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode WriteHDF5Parallel::create_element_tables()
{
 struct ElemSetCreator : public DataSetCreator
 {
 ErrorCode operator()( WriteHDF5* file, long size, const ExportSet* ex, long& start_id ) const
 {
 return file->create_elem_table( *ex, size, start_id );
 }
 };
 
 const int numtypes = exportList.size();
 std::vector< ExportSet* > groups( numtypes );
 std::vector< long > counts( numtypes ), offsets( numtypes ), max_ents( numtypes ), total_ents( numtypes );
 std::vector< wid_t > start_ids( numtypes );
 
 size_t idx = 0;
 std::list< ExportSet >::iterator ex_iter;
 for( ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter, ++idx )
 {
 groups[idx] = &*ex_iter;
 counts[idx] = ex_iter->range.size();
 }
 ErrorCode rval = create_dataset( numtypes, &counts[0], &offsets[0], &max_ents[0], &total_ents[0], ElemSetCreator(),
 &groups[0], &start_ids[0] );
 CHECK_MB( rval );
 
 for( idx = 0, ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter, ++idx )
 {
 ex_iter->first_id = start_ids[idx];
 ex_iter->offset = offsets[idx];
 ex_iter->max_num_ents = max_ents[idx];
 ex_iter->total_num_ents = total_ents[idx];
 rval = assign_ids( ex_iter->range, ex_iter->first_id + ex_iter->offset );
 CHECK_MB( rval );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode WriteHDF5Parallel::create_adjacency_tables()
{
 struct AdjSetCreator : public DataSetCreator
 {
 ErrorCode operator()( WriteHDF5* file, long size, const ExportSet* ex, long& start_id ) const
 {
 mhdf_Status status;
 hid_t handle = mhdf_createAdjacency( file->file_ptr(), ex->name(), size, &status );
 CHECK_HDFN( status );
 mhdf_closeData( file->file_ptr(), handle, &status );
 CHECK_HDFN( status );
 start_id = -1;
 return MB_SUCCESS;
 }
 };
 
 std::vector< ExportSet* > groups;
#ifdef WRITE_NODE_ADJACENCIES
 groups.push_back( &nodeSet );
#endif
 for( std::list< ExportSet >::iterator ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter )
 groups.push_back( &*ex_iter );
 
 ErrorCode rval;
 const int numtypes = groups.size();
 std::vector< long > counts( numtypes );
 std::vector< long > offsets( numtypes );
 std::vector< long > max_ents( numtypes );
 std::vector< long > totals( numtypes );
 for( int i = 0; i < numtypes; ++i )
 {
 wid_t count;
 rval = count_adjacencies( groups[i]->range, count );
 CHECK_MB( rval );
 counts[i] = count;
 }
 
 rval = create_dataset( numtypes, &counts[0], &offsets[0], &max_ents[0], &totals[0], AdjSetCreator(), &groups[0] );
 CHECK_MB( rval );
 
 // Cppcheck warning (false positive): variable groups is assigned a value that is never used
 for( int i = 0; i < numtypes; ++i )
 {
 groups[i]->max_num_adjs = max_ents[i];
 groups[i]->adj_offset = offsets[i];
 }
 return MB_SUCCESS;
}
 
const unsigned SSVB = 3;
 
void WriteHDF5Parallel::print_set_sharing_data( const Range& range, const char* label, Tag idt )
{
 dbgOut.printf( SSVB, "set\tid\towner\t%-*s\tfid\tshared\n", (int)( sizeof( EntityHandle ) * 2 ), "handle" );
 for( Range::iterator it = range.begin(); it != range.end(); ++it )
 {
 int id;
 iFace->tag_get_data( idt, &*it, 1, &id );
 EntityHandle handle = 0;
 unsigned owner = 0;
 wid_t file_id = 0;
 myPcomm->get_entityset_owner( *it, owner, &handle );
 if( !idMap.find( *it, file_id ) ) file_id = 0;
 dbgOut.printf( SSVB, "%s\t%d\t%u\t%lx\t%lu\t", label, id, owner, (unsigned long)handle,
 (unsigned long)file_id );
 std::vector< unsigned > procs;
 myPcomm->get_entityset_procs( *it, procs );
 if( procs.empty() )
 dbgOut.print( SSVB, "<none>\n" );
 else
 {
 for( unsigned i = 0; i < procs.size() - 1; ++i )
 dbgOut.printf( SSVB, "%u,", procs[i] );
 dbgOut.printf( SSVB, "%u\n", procs.back() );
 }
 }
}
 
void WriteHDF5Parallel::print_shared_sets()
{
 const char* tag_names[][2] = { { MATERIAL_SET_TAG_NAME, "block" },
 { DIRICHLET_SET_TAG_NAME, "nodeset" },
 { NEUMANN_SET_TAG_NAME, "sideset" },
 { 0, 0 } };
 
 for( int i = 0; tag_names[i][0]; ++i )
 {
 Tag tag;
 if( MB_SUCCESS != iFace->tag_get_handle( tag_names[i][0], 1, MB_TYPE_INTEGER, tag ) ) continue;
 
 Range tagged;
 iFace->get_entities_by_type_and_tag( 0, MBENTITYSET, &tag, 0, 1, tagged );
 print_set_sharing_data( tagged, tag_names[i][1], tag );
 }
 
 Tag geom, id;
 if( MB_SUCCESS != iFace->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom ) ) return;
 id = iFace->globalId_tag();
 
 const char* geom_names[] = { "vertex", "curve", "surface", "volume" };
 for( int d = 0; d <= 3; ++d )
 {
 Range tagged;
 const void* vals[] = { &d };
 iFace->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom, vals, 1, tagged );
 print_set_sharing_data( tagged, geom_names[d], id );
 }
}
 
ErrorCode WriteHDF5Parallel::communicate_shared_set_ids( const Range& owned, const Range& remote )
{
 ErrorCode rval;
 int mperr;
 const int TAG = 0xD0E;
 // const unsigned rank = myPcomm->proc_config().proc_rank();
 const MPI_Comm comm = myPcomm->proc_config().proc_comm();
 
 dbgOut.tprint( 1, "COMMUNICATING SHARED SET IDS\n" );
 dbgOut.print( 6, "Owned, shared sets: ", owned );
 
 // Post receive buffers for all procs for which we share sets
 
 std::vector< unsigned > procs;
 rval = myPcomm->get_entityset_owners( procs );
 CHECK_MB( rval );
 std::vector< unsigned >::iterator it = std::find( procs.begin(), procs.end(), myPcomm->proc_config().proc_rank() );
 if( it != procs.end() ) procs.erase( it );
 
 std::vector< MPI_Request > recv_req( procs.size(), MPI_REQUEST_NULL );
 std::vector< std::vector< unsigned long > > recv_buf( procs.size() );
 
 size_t recv_count = 0;
 for( size_t i = 0; i < procs.size(); ++i )
 {
 Range tmp;
 rval = myPcomm->get_owned_sets( procs[i], tmp );
 CHECK_MB( rval );
 size_t count =
 intersect( tmp, remote ).size(); // Necessary because we might not be writing all of the database
 if( count )
 {
 dbgOut.printf( 6, "Sets owned by proc %u (remote handles): ", procs[i] );
 if( dbgOut.get_verbosity() >= 6 )
 {
 Range remote_handles;
 tmp = intersect( tmp, remote );
 for( Range::iterator j = tmp.begin(); j != tmp.end(); ++j )
 {
 unsigned r;
 EntityHandle h;
 myPcomm->get_entityset_owner( *j, r, &h );
 assert( r == procs[i] );
 remote_handles.insert( h );
 }
 dbgOut.print( 6, remote_handles );
 }
 recv_count++;
 recv_buf[i].resize( 2 * count + 1 );
 dbgOut.printf( 5, "Posting receive buffer of size %lu for proc %u (%lu of %lu owned sets)\n",
 (unsigned long)recv_buf[i].size(), procs[i], count, tmp.size() );
 mperr =
 MPI_Irecv( &recv_buf[i][0], recv_buf[i].size(), MPI_UNSIGNED_LONG, procs[i], TAG, comm, &recv_req[i] );
 CHECK_MPI( mperr );
 }
 }
 
 // Send set ids to all procs with which we share them
 
 // First build per-process lists of sets for which we need to send data
 std::map< unsigned, Range > send_sets;
 std::vector< unsigned > set_procs;
 for( Range::reverse_iterator i = owned.rbegin(); i != owned.rend(); ++i )
 {
 set_procs.clear();
 rval = myPcomm->get_entityset_procs( *i, set_procs );
 CHECK_MB( rval );
 for( size_t j = 0; j < set_procs.size(); ++j )
 if( set_procs[j] != myPcomm->proc_config().proc_rank() ) send_sets[set_procs[j]].insert( *i );
 }
 assert( send_sets.find( myPcomm->proc_config().proc_rank() ) == send_sets.end() );
 
 // Now send the data
 std::vector< std::vector< unsigned long > > send_buf( send_sets.size() );
 std::vector< MPI_Request > send_req( send_sets.size() );
 std::map< unsigned, Range >::iterator si = send_sets.begin();
 for( size_t i = 0; si != send_sets.end(); ++si, ++i )
 {
 dbgOut.printf( 6, "Sending data for shared sets to proc %u: ", si->first );
 dbgOut.print( 6, si->second );
 
 send_buf[i].reserve( 2 * si->second.size() + 1 );
 send_buf[i].push_back( si->second.size() );
 for( Range::iterator j = si->second.begin(); j != si->second.end(); ++j )
 {
 send_buf[i].push_back( *j );
 send_buf[i].push_back( idMap.find( *j ) );
 }
 dbgOut.printf( 5, "Sending buffer of size %lu to proc %u (%lu of %lu owned sets)\n",
 (unsigned long)send_buf[i].size(), si->first, si->second.size(), owned.size() );
 mperr = MPI_Isend( &send_buf[i][0], send_buf[i].size(), MPI_UNSIGNED_LONG, si->first, TAG, comm, &send_req[i] );
 }
 
 // Process received data
 MPI_Status status;
 int idx;
 while( recv_count-- )
 {
 mperr = MPI_Waitany( recv_req.size(), &recv_req[0], &idx, &status );
 CHECK_MPI( mperr );
 
 assert( (unsigned)status.MPI_SOURCE == procs[idx] );
 assert( 2 * recv_buf[idx].front() + 1 == recv_buf[idx].size() );
 const size_t n = std::min< size_t >( recv_buf[idx].front(), ( recv_buf[idx].size() - 1 ) / 2 );
 dbgOut.printf( 5, "Received buffer of size %lu from proc %d\n", (unsigned long)( 2 * n + 1 ),
 (int)status.MPI_SOURCE );
 
 for( size_t i = 0; i < n; ++i )
 {
 EntityHandle handle = 0;
 rval = myPcomm->get_entityset_local_handle( procs[idx], recv_buf[idx][2 * i + 1], handle );
 CHECK_MB( rval );
 assert( handle != 0 );
 if( !idMap.insert( handle, recv_buf[idx][2 * i + 2], 1 ).second )
 error( MB_FAILURE ); // Conflicting IDs??????
 }
 
 recv_req[idx] = MPI_REQUEST_NULL;
 }
 assert( MPI_SUCCESS == MPI_Waitany( recv_req.size(), &recv_req[0], &idx, &status ) &&
 MPI_UNDEFINED == idx ); // Check that we got them all
 
 // Wait for all sends to complete before we release send
 // buffers (implicitly releases when we return from this function)
 
 std::vector< MPI_Status > stats( send_req.size() );
 mperr = MPI_Waitall( send_req.size(), &send_req[0], &stats[0] );
 CHECK_MPI( mperr );
 
 if( dbgOut.get_verbosity() >= SSVB ) print_shared_sets();
 
 return MB_SUCCESS;
}
 
// void get_global_ids(Interface* iFace, const unsigned long* ptr,
// size_t len, unsigned flags,
// std::vector<int>& ids)
//{
// Tag idtag;
// iFace->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, idtag);
// for (size_t i = 0; i < len; ++i) {
// if (flags & MESHSET_ORDERED) {
// int tmp;
// iFace->tag_get_data(idtag, ptr + i, 1, &tmp);
// ids.push_back(tmp);
// continue;
// }
//
// EntityHandle s = ptr[i];
// EntityHandle e = ptr[++i];
// for (; s <= e; ++s) {
// int tmp;
// iFace->tag_get_data(idtag, &s, 1, &tmp);
// ids.push_back(tmp);
// }
// }
//}
 
ErrorCode WriteHDF5Parallel::pack_set( Range::const_iterator it, unsigned long* buffer, size_t buffer_size )
{
 ErrorCode rval;
 const EntityHandle* ptr;
 int len;
 unsigned char flags;
 std::vector< wid_t > tmp;
 size_t newlen;
 
 // Buffer must always contain at least flags and desired sizes
 assert( buffer_size >= 4 );
 buffer_size -= 4;
 
 Range::const_iterator nd = it;
 ++nd;
 rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::CONTENTS, &len, &flags );
 CHECK_MB( rval );
 
 // Tag mattag;
 // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
 // int block;
 // if (MB_SUCCESS != iFace->tag_get_data(mattag, &*it, 1, &block))
 // block = 0;
 //
 // if (block) {
 // std::vector<int> ids;
 // get_global_ids(iFace, ptr, len, flags, ids);
 //}
 
 if( len && !( flags & MESHSET_ORDERED ) )
 {
 tmp.clear();
 bool blocked = false;
 assert( ( 0 == len % 2 ) );
 rval = range_to_blocked_list( ptr, len / 2, tmp, blocked );
 CHECK_MB( rval );
 if( blocked ) flags |= mhdf_SET_RANGE_BIT;
 }
 else
 {
 tmp.resize( len );
 rval = vector_to_id_list( ptr, len, &tmp[0], newlen, true );
 CHECK_MB( rval );
 tmp.resize( newlen );
 }
 
 buffer[0] = flags;
 buffer[1] = tmp.size();
 if( tmp.size() <= buffer_size ) std::copy( tmp.begin(), tmp.end(), buffer + 4 );
 
 rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::CHILDREN, &len );
 CHECK_MB( rval );
 tmp.resize( len );
 rval = vector_to_id_list( ptr, len, &tmp[0], newlen, true );
 tmp.resize( newlen );
 buffer[2] = tmp.size();
 if( tmp.size() <= buffer_size - buffer[1] ) std::copy( tmp.begin(), tmp.end(), buffer + 4 + buffer[1] );
 
 rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::PARENTS, &len );
 CHECK_MB( rval );
 tmp.resize( len );
 rval = vector_to_id_list( ptr, len, &tmp[0], newlen, true );
 tmp.resize( newlen );
 buffer[3] = tmp.size();
 if( tmp.size() <= buffer_size - buffer[1] - buffer[2] )
 std::copy( tmp.begin(), tmp.end(), buffer + 4 + buffer[1] + buffer[2] );
 
 return MB_SUCCESS;
}
 
template < typename TYPE >
static void convert_to_ranged_ids( const TYPE* buffer, size_t len, std::vector< WriteHDF5::wid_t >& result )
{
 if( !len )
 {
 result.clear();
 return;
 }
 
 result.resize( len * 2 );
 Range tmp;
 for( size_t i = 0; i < len; i++ )
 tmp.insert( (EntityHandle)buffer[i] );
 result.resize( tmp.psize() * 2 );
 int j = 0;
 for( Range::const_pair_iterator pit = tmp.const_pair_begin(); pit != tmp.const_pair_end(); ++pit, j++ )
 {
 result[2 * j] = pit->first;
 result[2 * j + 1] = pit->second - pit->first + 1;
 }
}
 
static void merge_ranged_ids( const unsigned long* range_list, size_t len, std::vector< WriteHDF5::wid_t >& result )
{
 typedef WriteHDF5::wid_t wid_t;
 assert( 0 == len % 2 );
 assert( 0 == result.size() % 2 );
 STATIC_ASSERT( sizeof( std::pair< wid_t, wid_t > ) == 2 * sizeof( wid_t ) );
 
 result.insert( result.end(), range_list, range_list + len );
 size_t plen = result.size() / 2;
 Range tmp;
 for( size_t i = 0; i < plen; i++ )
 {
 EntityHandle starth = (EntityHandle)result[2 * i];
 EntityHandle endh = (EntityHandle)result[2 * i] + (wid_t)result[2 * i + 1] - 1; // id + count - 1
 tmp.insert( starth, endh );
 }
 // Now convert back to std::vector<WriteHDF5::wid_t>, compressed range format
 result.resize( tmp.psize() * 2 );
 int j = 0;
 for( Range::const_pair_iterator pit = tmp.const_pair_begin(); pit != tmp.const_pair_end(); ++pit, j++ )
 {
 result[2 * j] = pit->first;
 result[2 * j + 1] = pit->second - pit->first + 1;
 }
}
 
static void merge_vector_ids( const unsigned long* list, size_t len, std::vector< WriteHDF5::wid_t >& result )
{
 result.insert( result.end(), list, list + len );
}
 
ErrorCode WriteHDF5Parallel::unpack_set( EntityHandle set, const unsigned long* buffer, size_t buffer_size )
{
 // Use local variables for readability
 assert( buffer_size >= 4 );
 assert( buffer[1] + buffer[2] + buffer[3] <= buffer_size );
 const unsigned long flags = buffer[0];
 unsigned long num_content = buffer[1];
 const unsigned long num_child = buffer[2];
 const unsigned long num_parent = buffer[3];
 const unsigned long* contents = buffer + 4;
 const unsigned long* children = contents + num_content;
 const unsigned long* parents = children + num_child;
 
 SpecialSetData* data = find_set_data( set );
 assert( NULL != data );
 if( NULL == data ) return MB_FAILURE;
 
 // Tag mattag;
 // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
 // int block;
 // if (MB_SUCCESS != iFace->tag_get_data(mattag, &set, 1, &block))
 // block = 0;
 
 // If either the current data or the new data is in ranged format,
 // then change the other to ranged format if it isn't already
 // in both cases when they differ, the data will end up "compressed range"
 std::vector< wid_t > tmp;
 if( ( flags & mhdf_SET_RANGE_BIT ) != ( data->setFlags & mhdf_SET_RANGE_BIT ) )
 {
 if( flags & mhdf_SET_RANGE_BIT )
 {
 tmp = data->contentIds;
 convert_to_ranged_ids( &tmp[0], tmp.size(), data->contentIds );
 data->setFlags |= mhdf_SET_RANGE_BIT;
 }
 else
 {
 tmp.clear();
 convert_to_ranged_ids( contents, num_content, tmp );
 num_content = tmp.size();
 if( sizeof( wid_t ) < sizeof( long ) )
 {
 size_t old_size = tmp.size();
 tmp.resize( sizeof( long ) * old_size / sizeof( wid_t ) );
 unsigned long* array = reinterpret_cast< unsigned long* >( &tmp[0] );
 for( long i = ( (long)old_size ) - 1; i >= 0; --i )
 array[i] = tmp[i];
 contents = array;
 }
 else if( sizeof( wid_t ) > sizeof( long ) )
 {
 unsigned long* array = reinterpret_cast< unsigned long* >( &tmp[0] );
 std::copy( tmp.begin(), tmp.end(), array );
 }
 contents = reinterpret_cast< unsigned long* >( &tmp[0] );
 }
 }
 
 if( data->setFlags & mhdf_SET_RANGE_BIT )
 merge_ranged_ids( contents, num_content, data->contentIds );
 else
 merge_vector_ids( contents, num_content, data->contentIds );
 
 merge_vector_ids( children, num_child, data->childIds );
 merge_vector_ids( parents, num_parent, data->parentIds );
 return MB_SUCCESS;
}
 
ErrorCode WriteHDF5Parallel::communicate_shared_set_data( const Range& owned, const Range& remote )
{
 ErrorCode rval;
 int mperr;
 const unsigned rank = myPcomm->proc_config().proc_rank();
 const MPI_Comm comm = myPcomm->proc_config().proc_comm();
 
 dbgOut.tprintf( 1, "COMMUNICATING SHARED SET DATA (%lu owned & %lu remote)\n", (unsigned long)owned.size(),
 (unsigned long)remote.size() );
 
 // Calculate the total number of messages to be in transit (send and receive)
 size_t nummess = 0;
 std::vector< unsigned > procs;
 ;
 Range shared( owned );
 shared.merge( remote );
 for( Range::iterator i = shared.begin(); i != shared.end(); ++i )
 {
 procs.clear();
 rval = myPcomm->get_entityset_procs( *i, procs );
 CHECK_MB( rval );
 nummess += procs.size();
 }
 
 // Choose a receive buffer size. We need 4*sizeof(long) minimum,
 // but that is almost useless so use 16*sizeof(long) as the minimum
 // instead. Choose an upper limit such that we don't exceed 32 MB
 // of allocated memory (unless we absolutely must to meet the minimum.)
 // Also, don't initially choose buffers larger than 128*sizeof(long).
 const size_t MAX_BUFFER_MEM = 32 * 1024 * 1024 / sizeof( long );
 // const size_t INIT_BUFFER_SIZE = 128;
 const size_t INIT_BUFFER_SIZE = 1024;
 const size_t MIN_BUFFER_SIZE = 16;
 size_t init_buff_size = INIT_BUFFER_SIZE;
 if( init_buff_size * nummess > MAX_BUFFER_MEM ) init_buff_size = MAX_BUFFER_MEM / nummess;
 if( init_buff_size < MIN_BUFFER_SIZE ) init_buff_size = MIN_BUFFER_SIZE;
 
 dbgOut.printf( 2, "Using buffer size of %lu for an expected message count of %lu\n", (unsigned long)init_buff_size,
 (unsigned long)nummess );
 
 // Count number of recvs
 size_t numrecv = 0;
 for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
 {
 procs.clear();
 rval = myPcomm->get_entityset_procs( *i, procs );
 CHECK_MB( rval );
 numrecv += procs.size();
 if( std::find( procs.begin(), procs.end(), rank ) != procs.end() ) --numrecv;
 }
 
 // Post receive buffers for all owned sets for all sharing procs
 std::vector< MPI_Request > recv_req( numrecv, MPI_REQUEST_NULL );
 std::vector< MPI_Request > lrecv_req( numrecv, MPI_REQUEST_NULL );
 
 std::vector< std::vector< unsigned long > > recv_buf( numrecv, std::vector< unsigned long >( init_buff_size ) );
 int idx = 0;
 for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
 {
 procs.clear();
 rval = myPcomm->get_entityset_procs( *i, procs );
 CHECK_MB( rval );
 for( size_t j = 0; j < procs.size(); ++j )
 {
 if( procs[j] == rank ) continue;
 int tag = ID_FROM_HANDLE( *i );
 if( *i != CREATE_HANDLE( MBENTITYSET, tag ) )
 {
#ifndef NDEBUG
 abort();
#endif
 CHECK_MB( MB_FAILURE );
 }
 dbgOut.printf( 5, "Posting buffer to receive set %d from proc %u\n", tag, procs[j] );
 mperr =
 MPI_Irecv( &recv_buf[idx][0], init_buff_size, MPI_UNSIGNED_LONG, procs[j], tag, comm, &recv_req[idx] );
 CHECK_MPI( mperr );
 ++idx;
 }
 }
 assert( (size_t)idx == numrecv );
 
 // Now send set data for all remote sets that I know about
 std::vector< MPI_Request > send_req( remote.size() );
 std::vector< std::vector< unsigned long > > send_buf( remote.size() );
 idx = 0;
 for( Range::iterator i = remote.begin(); i != remote.end(); ++i, ++idx )
 {
 send_buf[idx].resize( init_buff_size );
 rval = pack_set( i, &send_buf[idx][0], init_buff_size );
 CHECK_MB( rval );
 EntityHandle remote_handle;
 unsigned owner;
 rval = myPcomm->get_entityset_owner( *i, owner, &remote_handle );
 CHECK_MB( rval );
 
 int tag = ID_FROM_HANDLE( remote_handle );
 assert( remote_handle == CREATE_HANDLE( MBENTITYSET, tag ) );
 dbgOut.printf( 5, "Sending %lu values for set %d to proc %u\n",
 send_buf[idx][1] + send_buf[idx][2] + send_buf[idx][3] + 4, tag, owner );
 mperr = MPI_Isend( &send_buf[idx][0], init_buff_size, MPI_UNSIGNED_LONG, owner, tag, comm, &send_req[idx] );
 CHECK_MPI( mperr );
 }
 
 // Tag mattag;
 // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
 
 // Now initialize local data for managing contents of owned, shared sets
 assert( specialSets.empty() );
 specialSets.clear();
 specialSets.reserve( owned.size() );
 for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
 {
 // int block;
 // if (MB_SUCCESS != iFace->tag_get_data(mattag, &*i, 1, &block))
 // block = 0;
 // std::vector<int> ids;
 
 SpecialSetData data;
 data.setHandle = *i;
 rval = iFace->get_meshset_options( *i, data.setFlags );
 CHECK_MB( rval );
 specialSets.push_back( data );
 std::vector< EntityHandle > list;
 if( data.setFlags & MESHSET_ORDERED )
 {
 list.clear();
 rval = iFace->get_entities_by_handle( *i, list );
 CHECK_MB( rval );
 rval = vector_to_id_list( list, specialSets.back().contentIds, true );
 CHECK_MB( rval );
 // if (block)
 // get_global_ids(iFace, &list[0], list.size(), MESHSET_ORDERED, ids);
 }
 else
 {
 Range range;
 rval = iFace->get_entities_by_handle( *i, range );
 CHECK_MB( rval );
 bool ranged;
 rval = range_to_blocked_list( range, specialSets.back().contentIds, ranged );
 if( ranged ) specialSets.back().setFlags |= mhdf_SET_RANGE_BIT;
 // if (block) {
 // std::vector<EntityHandle> tmp;
 // for (Range::const_pair_iterator pi = range.const_pair_begin(); pi !=
 // range.const_pair_end(); ++pi) {
 // tmp.push_back(pi->first);
 // tmp.push_back(pi->second);
 // }
 // get_global_ids(iFace, &tmp[0], tmp.size(), ranged ? 0 : MESHSET_ORDERED, ids);
 //}
 }
 
 list.clear();
 rval = iFace->get_parent_meshsets( *i, list );
 CHECK_MB( rval );
 rval = vector_to_id_list( list, specialSets.back().parentIds, true );
 CHECK_MB( rval );
 rval = iFace->get_child_meshsets( *i, list );
 CHECK_MB( rval );
 rval = vector_to_id_list( list, specialSets.back().childIds, true );
 CHECK_MB( rval );
 }
 
 // Process received buffers, repost larger buffers where necessary
 size_t remaining = numrecv;
 numrecv = 0;
 while( remaining-- )
 {
 std::vector< unsigned long > dead;
 MPI_Status status;
 mperr = MPI_Waitany( recv_req.size(), &recv_req[0], &idx, &status );
 CHECK_MPI( mperr );
 EntityHandle handle = CREATE_HANDLE( MBENTITYSET, status.MPI_TAG );
 std::vector< unsigned long >& buff = recv_buf[idx];
 size_t size = buff[1] + buff[2] + buff[3] + 4;
 dbgOut.printf( 5, "Received %lu values for set %d from proc %d\n", (unsigned long)size, status.MPI_TAG,
 status.MPI_SOURCE );
 if( size <= init_buff_size )
 {
 rval = unpack_set( handle, &buff[0], init_buff_size );
 CHECK_MB( rval );
 dead.swap( buff ); // Release memory
 }
 else
 {
 // Data was too big for init_buff_size
 // repost with larger buffer
 buff.resize( size );
 dbgOut.printf( 5, "Re-Posting buffer to receive set %d from proc %d with size %lu\n", status.MPI_TAG,
 status.MPI_SOURCE, (unsigned long)size );
 mperr = MPI_Irecv( &buff[0], size, MPI_UNSIGNED_LONG, status.MPI_SOURCE, status.MPI_TAG, comm,
 &lrecv_req[idx] );
 CHECK_MPI( mperr );
 ++numrecv;
 }
 recv_req[idx] = MPI_REQUEST_NULL;
 }
 
 // Wait for sends to complete
 MPI_Waitall( send_req.size(), &send_req[0], MPI_STATUSES_IGNORE );
 
 // Re-send sets that didn't fit initial buffer size
 idx = 0;
 for( Range::iterator i = remote.begin(); i != remote.end(); ++i, ++idx )
 {
 std::vector< unsigned long >& buff = send_buf[idx];
 size_t size = buff[1] + buff[2] + buff[3] + 4;
 if( size <= init_buff_size ) continue;
 
 buff.resize( size );
 rval = pack_set( i, &buff[0], size );
 CHECK_MB( rval );
 EntityHandle remote_handle;
 unsigned owner;
 rval = myPcomm->get_entityset_owner( *i, owner, &remote_handle );
 CHECK_MB( rval );
 
 int tag = ID_FROM_HANDLE( remote_handle );
 assert( remote_handle == CREATE_HANDLE( MBENTITYSET, tag ) );
 dbgOut.printf( 5, "Sending %lu values for set %d to proc %u\n", (unsigned long)size, tag, owner );
 mperr = MPI_Isend( &buff[0], size, MPI_UNSIGNED_LONG, owner, tag, comm, &send_req[idx] );
 CHECK_MPI( mperr );
 }
 
 // Process received buffers
 remaining = numrecv;
 while( remaining-- )
 {
 std::vector< unsigned long > dead;
 MPI_Status status;
 mperr = MPI_Waitany( lrecv_req.size(), &lrecv_req[0], &idx, &status );
 CHECK_MPI( mperr );
 EntityHandle handle = CREATE_HANDLE( MBENTITYSET, status.MPI_TAG );
 std::vector< unsigned long >& buff = recv_buf[idx];
 dbgOut.printf( 5, "Received %lu values for set %d from proc %d\n", 4 + buff[1] + buff[2] + buff[3],
 status.MPI_TAG, status.MPI_SOURCE );
 rval = unpack_set( handle, &buff[0], buff.size() );
 CHECK_MB( rval );
 dead.swap( buff ); // Release memory
 
 lrecv_req[idx] = MPI_REQUEST_NULL;
 }
 
 // Wait for sends to complete
 MPI_Waitall( send_req.size(), &send_req[0], MPI_STATUSES_IGNORE );
 
 return MB_SUCCESS;
}
 
ErrorCode WriteHDF5Parallel::create_meshset_tables( double* times )
{
 Range::const_iterator riter;
 const unsigned rank = myPcomm->proc_config().proc_rank();
 
 START_SERIAL;
 print_type_sets( iFace, &dbgOut, setSet.range );
 END_SERIAL;
 CpuTimer timer;
 
 // Remove remote sets from setSets
 Range shared, owned, remote;
 ErrorCode rval = myPcomm->get_shared_sets( shared );
 CHECK_MB( rval );
 shared = intersect( shared, setSet.range );
 rval = myPcomm->get_owned_sets( rank, owned );
 CHECK_MB( rval );
 owned = intersect( owned, setSet.range );
 remote = subtract( shared, owned );
 setSet.range = subtract( setSet.range, remote );
 
 // Create set meta table
 struct SetDescCreator : public DataSetCreator
 {
 ErrorCode operator()( WriteHDF5* writer, long size, const ExportSet*, long& start_id ) const
 {
 return writer->create_set_meta( size, start_id );
 }
 };
 long count = setSet.range.size();
 rval = create_dataset( 1, &count, &setSet.offset, &setSet.max_num_ents, &setSet.total_num_ents, SetDescCreator(),
 NULL, &setSet.first_id );
 CHECK_MB( rval );
 writeSets = setSet.max_num_ents > 0;
 
 rval = assign_ids( setSet.range, setSet.first_id + setSet.offset );
 CHECK_MB( rval );
 if( times ) times[SET_OFFSET_TIME] = timer.time_elapsed();
 
 // Exchange file IDS for sets between all procs
 rval = communicate_shared_set_ids( owned, remote );
 CHECK_MB( rval );
 if( times ) times[SHARED_SET_IDS] = timer.time_elapsed();
 
 // Communicate remote set contents, children, etc.
 rval = communicate_shared_set_data( owned, remote );
 CHECK_MB( rval );
 if( times ) times[SHARED_SET_CONTENTS] = timer.time_elapsed();
 
 // Communicate counts for owned sets
 long data_counts[3]; // { #contents, #children, #parents }
 rval = count_set_size( setSet.range, data_counts[0], data_counts[1], data_counts[2] );
 CHECK_MB( rval );
 if( times ) times[SET_OFFSET_TIME] += timer.time_elapsed();
 
 long offsets[3], max_counts[3], totals[3];
 rval = create_dataset( 3, data_counts, offsets, max_counts, totals );
 CHECK_MB( rval );
 
 // Create the datasets
 if( 0 == myPcomm->proc_config().proc_rank() )
 {
 rval = create_set_tables( totals[0], totals[1], totals[2] );
 CHECK_MB( rval );
 }
 
 // Store communicated global data
 setContentsOffset = offsets[0];
 setChildrenOffset = offsets[1];
 setParentsOffset = offsets[2];
 maxNumSetContents = max_counts[0];
 maxNumSetChildren = max_counts[1];
 maxNumSetParents = max_counts[2];
 writeSetContents = totals[0] > 0;
 writeSetChildren = totals[1] > 0;
 writeSetParents = totals[2] > 0;
 
 dbgOut.printf( 2, "set contents: %ld local, %ld global, offset = %ld\n", data_counts[0], totals[0], offsets[0] );
 dbgOut.printf( 2, "set children: %ld local, %ld global, offset = %ld\n", data_counts[1], totals[1], offsets[1] );
 dbgOut.printf( 2, "set parents: %ld local, %ld global, offset = %ld\n", data_counts[2], totals[2], offsets[2] );
 
 return MB_SUCCESS;
}
 
void WriteHDF5Parallel::remove_remote_entities( EntityHandle relative, Range& range )
{
 Range result;
 result.merge( intersect( range, nodeSet.range ) );
 result.merge( intersect( range, setSet.range ) );
 for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
 result.merge( intersect( range, eiter->range ) );
 
 // result.merge(intersect(range, myParallelSets));
 Range sets;
 int junk;
 sets.merge( Range::lower_bound( range.begin(), range.end(), CREATE_HANDLE( MBENTITYSET, 0, junk ) ), range.end() );
 remove_remote_sets( relative, sets );
 result.merge( sets );
 range.swap( result );
}
 
void WriteHDF5Parallel::remove_remote_sets( EntityHandle /* relative */, Range& range )
{
 Range result( intersect( range, setSet.range ) );
 // Store the non-intersecting entities separately if needed
 // Range remaining(subtract(range, result));
 range.swap( result );
}
 
void WriteHDF5Parallel::remove_remote_entities( EntityHandle relative, std::vector< EntityHandle >& vect )
{
 Range intrsct;
 for( std::vector< EntityHandle >::const_iterator iter = vect.begin(); iter != vect.end(); ++iter )
 intrsct.insert( *iter );
 remove_remote_entities( relative, intrsct );
 
 unsigned int read, write;
 for( read = write = 0; read < vect.size(); ++read )
 {
 if( intrsct.find( vect[read] ) != intrsct.end() )
 {
 if( read != write ) vect[write] = vect[read];
 ++write;
 }
 }
 if( write != vect.size() ) vect.resize( write );
}
 
void WriteHDF5Parallel::remove_remote_sets( EntityHandle relative, std::vector< EntityHandle >& vect )
{
 Range intrsct;
 for( std::vector< EntityHandle >::const_iterator iter = vect.begin(); iter != vect.end(); ++iter )
 intrsct.insert( *iter );
 remove_remote_sets( relative, intrsct );
 
 unsigned int read, write;
 for( read = write = 0; read < vect.size(); ++read )
 {
 if( intrsct.find( vect[read] ) != intrsct.end() )
 {
 if( read != write ) vect[write] = vect[read];
 ++write;
 }
 }
 if( write != vect.size() ) vect.resize( write );
}
 
ErrorCode WriteHDF5Parallel::exchange_file_ids( const Range& nonlocal )
{
 ErrorCode rval;
 
 // For each entity owned on the interface, write its file id to
 // a tag. The sets of entities to be written should already contain
 // only owned entities so by intersecting with them we not only
 // filter by entities to be written, but also restrict to entities
 // owned by the proc
 
 // Get list of interface entities
 Range imesh, tmp;
 for( std::list< ExportSet >::reverse_iterator i = exportList.rbegin(); i != exportList.rend(); ++i )
 {
 tmp.clear();
 rval = myPcomm->filter_pstatus( i->range, PSTATUS_SHARED, PSTATUS_AND, -1, &tmp );
 if( MB_SUCCESS != rval ) return error( rval );
 imesh.merge( tmp );
 }
 tmp.clear();
 rval = myPcomm->filter_pstatus( nodeSet.range, PSTATUS_SHARED, PSTATUS_AND, -1, &tmp );
 if( MB_SUCCESS != rval ) return error( rval );
 imesh.merge( tmp );
 
 // Create tag to store file IDs
 EntityHandle default_val = 0;
 Tag file_id_tag = 0;
 rval = iFace->tag_get_handle( "__hdf5_ll_fileid", 1, MB_TYPE_HANDLE, file_id_tag, MB_TAG_DENSE | MB_TAG_CREAT,
 &default_val );
 if( MB_SUCCESS != rval ) return error( rval );
 
 // Copy file IDs into tag
 std::vector< EntityHandle > file_id_vect( imesh.size() );
 Range::const_iterator i;
 std::vector< EntityHandle >::iterator j = file_id_vect.begin();
 for( i = imesh.begin(); i != imesh.end(); ++i, ++j )
 {
 *j = idMap.find( *i );
 if( !*j )
 {
 iFace->tag_delete( file_id_tag );
 return error( MB_FAILURE );
 }
 }
 rval = iFace->tag_set_data( file_id_tag, imesh, &file_id_vect[0] );
 if( MB_SUCCESS != rval )
 {
 iFace->tag_delete( file_id_tag );
 return error( rval );
 }
 
 // Do communication
 rval = myPcomm->exchange_tags( file_id_tag, imesh );
 if( MB_SUCCESS != rval )
 {
 iFace->tag_delete( file_id_tag );
 return error( rval );
 }
 
 // Copy file IDs from tag into idMap for remote entities
 file_id_vect.resize( nonlocal.size() );
 rval = iFace->tag_get_data( file_id_tag, nonlocal, &file_id_vect[0] );
 if( MB_SUCCESS != rval )
 {
 iFace->tag_delete( file_id_tag );
 return error( rval );
 }
 
 j = file_id_vect.begin();
 for( i = nonlocal.begin(); i != nonlocal.end(); ++i, ++j )
 {
 if( *j == 0 )
 {
 int owner = -1;
 myPcomm->get_owner( *i, owner );
 const char* name = CN::EntityTypeName( TYPE_FROM_HANDLE( *i ) );
 int id = ID_FROM_HANDLE( *i );
 MB_SET_ERR_CONT( "Process " << myPcomm->proc_config().proc_rank()
 << " did not receive valid id handle for shared " << name << " " << id
 << " owned by process " << owner );
 dbgOut.printf( 1,
 "Did not receive valid remote id for "
 "shared %s %d owned by process %d",
 name, id, owner );
 iFace->tag_delete( file_id_tag );
 return error( MB_FAILURE );
 }
 else
 {
 if( !idMap.insert( *i, *j, 1 ).second )
 {
 iFace->tag_delete( file_id_tag );
 return error( MB_FAILURE );
 }
 }
 }
 
#ifndef NDEBUG
 // Check that writer is correct with regards to which entities
 // that it owns by verifying that the file ids that we thought
 // we were sending where not received instead
 file_id_vect.resize( imesh.size() );
 rval = iFace->tag_get_data( file_id_tag, imesh, &file_id_vect[0] );
 if( MB_SUCCESS != rval )
 {
 iFace->tag_delete( file_id_tag );
 return error( rval );
 }
 int invalid_count = 0;
 j = file_id_vect.begin();
 for( i = imesh.begin(); i != imesh.end(); ++i, ++j )
 {
 EntityHandle h = idMap.find( *i );
 if( *j != h )
 {
 ++invalid_count;
 dbgOut.printf( 1, "Conflicting ownership for %s %ld\n", CN::EntityTypeName( TYPE_FROM_HANDLE( *i ) ),
 (long)ID_FROM_HANDLE( *i ) );
 }
 }
 if( invalid_count )
 {
 iFace->tag_delete( file_id_tag );
 MB_SET_ERR( MB_FAILURE, invalid_count << " entities with conflicting ownership found by process "
 << myPcomm->proc_config().proc_rank()
 << ". This will result in duplicate entities written to file" );
 }
#endif
 
 return iFace->tag_delete( file_id_tag );
}
 
void WriteHDF5Parallel::print_times( const double* times ) const
{
 if( !myPcomm )
 {
 WriteHDF5::print_times( times );
 }
 else
 {
 double recv[NUM_TIMES];
 MPI_Reduce( (void*)times, recv, NUM_TIMES, MPI_DOUBLE, MPI_MAX, 0, myPcomm->proc_config().proc_comm() );
 if( 0 == myPcomm->proc_config().proc_rank() ) WriteHDF5::print_times( recv );
 }
}
 
} // namespace moab