MOAB: Mesh Oriented datABase  (version 5.3.1)
WriteHDF5Parallel.cpp
Go to the documentation of this file.
00001 #undef DEBUG
00002 #undef TIME_DEBUG
00003 
00004 #include <cstdarg>
00005 #include <ctime>
00006 #include <cstdlib>
00007 
00008 #include <cstring>
00009 #include <cassert>
00010 
00011 #include <vector>
00012 #include <set>
00013 #include <map>
00014 #include <utility>
00015 #include <iostream>
00016 #include <sstream>
00017 #include <string>
00018 
00019 #include "moab/Interface.hpp"
00020 #include "Internals.hpp"
00021 #include "MBTagConventions.hpp"
00022 #include "MBParallelConventions.h"
00023 #include "moab/ParallelComm.hpp"
00024 #include "moab/CN.hpp"
00025 #include "moab/Range.hpp"
00026 #include "moab/CpuTimer.hpp"
00027 
00028 #include "WriteHDF5Parallel.hpp"
00029 
00030 #ifndef MOAB_HAVE_HDF5
00031 #error Attempt to compile WriteHDF5Parallel with HDF5 support disabled
00032 #endif
00033 
00034 #include <H5Tpublic.h>
00035 #include <H5Ppublic.h>
00036 #include <H5FDmpi.h>
00037 #include <H5FDmpio.h>
00038 
00039 #include "mhdf.h"
00040 
00041 #include "IODebugTrack.hpp"
00042 #include "moab/FileOptions.hpp"
00043 
00044 namespace
00045 {
00046 template < bool Condition >
00047 struct STATIC_ASSERTION;
00048 template <>
00049 struct STATIC_ASSERTION< true >
00050 {
00051 };
00052 }  // namespace
00053 
00054 #define PP_CAT_( a, b ) a##b
00055 #define PP_CAT( a, b )  PP_CAT_( a, b )
00056 #define STATIC_ASSERT( Condition )                                                      \
00057     enum                                                                                \
00058     {                                                                                   \
00059         PP_CAT( dummy, __LINE__ ) = sizeof( ::STATIC_ASSERTION< (bool)( Condition ) > ) \
00060     }
00061 
00062 namespace moab
00063 {
00064 
00065 #ifndef _WIN32  // problematic for windows
00066 // Need an MPI type that we can put handles in
00067 STATIC_ASSERT( sizeof( unsigned long ) >= sizeof( EntityHandle ) );
00068 
00069 // Need an MPI type that we can put file IDs in
00070 STATIC_ASSERT( sizeof( unsigned long ) >= sizeof( WriteHDF5::wid_t ) );
00071 #endif
00072 
00073 // This function doesn't do anything useful. It's just a nice
00074 // place to set a break point to determine why the reader fails.
00075 static inline ErrorCode error( ErrorCode rval )
00076 {
00077     return rval;
00078 }
00079 
00080 const char* mpi_err_str( int errorcode )
00081 {
00082     static char buffer[2048];
00083     int len = sizeof( buffer );
00084     MPI_Error_string( errorcode, buffer, &len );
00085     buffer[std::min( (size_t)len, sizeof( buffer ) - 1 )] = '\0';
00086     return buffer;
00087 }
00088 
00089 #define MPI_FAILURE_MSG( A ) \
00090     "MPI Failure at " __FILE__ ":%d : (Code %d) %s\n", __LINE__, (int)( A ), mpi_err_str( ( A ) )
00091 
00092 #define CHECK_MPI( A )                                                                               \
00093     do                                                                                               \
00094     {                                                                                                \
00095         if( MPI_SUCCESS != ( A ) )                                                                   \
00096         {                                                                                            \
00097             MB_SET_ERR_CONT( "MPI Failure : (Code " << (int)( A ) << ") " << mpi_err_str( ( A ) ) ); \
00098             dbgOut.printf( 1, MPI_FAILURE_MSG( ( A ) ) );                                            \
00099             return error( MB_FAILURE );                                                              \
00100         }                                                                                            \
00101     } while( false )
00102 
00103 #define MB_FAILURE_MSG( A ) "MOAB_Failure at " __FILE__ ":%d : %s (%d)\n", __LINE__, ErrorCodeStr[( A )], (int)( A )
00104 
00105 #define CHECK_MB( A )                                                    \
00106     do                                                                   \
00107     {                                                                    \
00108         if( MB_SUCCESS != ( A ) )                                        \
00109         {                                                                \
00110             MB_SET_ERR_CONT( "MOAB Failure : " << ErrorCodeStr[( A )] ); \
00111             dbgOut.printf( 1, MB_FAILURE_MSG( ( A ) ) );                 \
00112             return error( A );                                           \
00113         }                                                                \
00114     } while( false )
00115 
00116 #define HDF_FAILURE_MSG( A ) "MHDF Failure at " __FILE__ ":%d : %s\n", __LINE__, mhdf_message( &( A ) )
00117 
00118 #define CHECK_HDF( A )                                                      \
00119     do                                                                      \
00120     {                                                                       \
00121         if( mhdf_isError( &( A ) ) )                                        \
00122         {                                                                   \
00123             MB_SET_ERR_CONT( "MHDF Failure : " << mhdf_message( &( A ) ) ); \
00124             dbgOut.printf( 1, HDF_FAILURE_MSG( ( A ) ) );                   \
00125             return error( MB_FAILURE );                                     \
00126         }                                                                   \
00127     } while( false )
00128 
00129 #define CHECK_HDFN( A )                                                     \
00130     do                                                                      \
00131     {                                                                       \
00132         if( mhdf_isError( &( A ) ) )                                        \
00133         {                                                                   \
00134             MB_SET_ERR_CONT( "MHDF Failure : " << mhdf_message( &( A ) ) ); \
00135             return error( MB_FAILURE );                                     \
00136         }                                                                   \
00137     } while( false )
00138 
00139 #ifdef VALGRIND
00140 #include <valgrind/memcheck.h>
00141 
00142 template < typename T >
00143 inline void VALGRIND_MAKE_VEC_UNDEFINED( std::vector< T >& v )
00144 {
00145     if( v.size() ) {}
00146     (void)VALGRIND_MAKE_MEM_UNDEFINED( &v[0], v.size() * sizeof( T ) );
00147 }
00148 
00149 #else
00150 #ifndef VALGRIND_CHECK_MEM_IS_DEFINED
00151 #define VALGRIND_CHECK_MEM_IS_DEFINED( a, b ) ( (void)0 )
00152 #endif
00153 #ifndef VALGRIND_CHECK_MEM_IS_ADDRESSABLE
00154 #define VALGRIND_CHECK_MEM_IS_ADDRESSABLE( a, b ) ( (void)0 )
00155 #endif
00156 #ifndef VALGRIND_MAKE_MEM_UNDEFINED
00157 #define VALGRIND_MAKE_MEM_UNDEFINED( a, b ) ( (void)0 )
00158 #endif
00159 
00160 template < typename T >
00161 inline void VALGRIND_MAKE_VEC_UNDEFINED( std::vector< T >& )
00162 {
00163     /* Nothing to do */
00164 }
00165 
00166 #endif
00167 
00168 #ifndef NDEBUG
00169 #define START_SERIAL                                                      \
00170     for( unsigned _x = 0; _x < myPcomm->proc_config().proc_size(); ++_x ) \
00171     {                                                                     \
00172         MPI_Barrier( myPcomm->proc_config().proc_comm() );                \
00173         if( _x != myPcomm->proc_config().proc_rank() ) continue
00174 #define END_SERIAL \
00175     }              \
00176     MPI_Barrier( myPcomm->proc_config().proc_comm() )
00177 #else
00178 #define START_SERIAL
00179 #define END_SERIAL
00180 #endif
00181 
00182 static int my_Gatherv( void* sendbuf, int sendcount, MPI_Datatype sendtype, std::vector< unsigned char >& recvbuf,
00183                        std::vector< int >& recvcounts, int root, MPI_Comm comm )
00184 {
00185     int nproc, rank, bytes, err;
00186     MPI_Comm_size( comm, &nproc );
00187     MPI_Comm_rank( comm, &rank );
00188     MPI_Type_size( sendtype, &bytes );
00189 
00190     recvcounts.resize( rank == root ? nproc : 0 );
00191     err = MPI_Gather( &sendcount, 1, MPI_INT, &recvcounts[0], 1, MPI_INT, root, comm );
00192     if( MPI_SUCCESS != err ) return err;
00193 
00194     std::vector< int > disp( recvcounts.size() );
00195     if( root == rank )
00196     {
00197         disp[0] = 0;
00198         for( int i = 1; i < nproc; ++i )
00199             disp[i] = disp[i - 1] + recvcounts[i - 1];
00200         recvbuf.resize( bytes * ( disp.back() + recvcounts.back() ) );
00201     }
00202 
00203     return MPI_Gatherv( sendbuf, sendcount, sendtype, &recvbuf[0], &recvcounts[0], &disp[0], sendtype, root, comm );
00204 }
00205 
00206 static void print_type_sets( Interface* iFace, DebugOutput* str, Range& sets )
00207 {
00208     const unsigned VB = 2;
00209     if( str->get_verbosity() < VB ) return;
00210 
00211     Tag gid, did, bid, sid, nid;
00212     gid = iFace->globalId_tag();
00213     iFace->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, did );
00214     iFace->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, bid );
00215     iFace->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, nid );
00216     iFace->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, sid );
00217     Range typesets[10];
00218     const char* typenames[] = {
00219         "Block ", "Sideset ", "NodeSet", "Vertex", "Curve", "Surface", "Volume", "Body", "Other"
00220     };
00221     for( Range::iterator riter = sets.begin(); riter != sets.end(); ++riter )
00222     {
00223         unsigned dim, id;  //, oldsize;
00224         if( MB_SUCCESS == iFace->tag_get_data( bid, &*riter, 1, &id ) )
00225             dim = 0;
00226         else if( MB_SUCCESS == iFace->tag_get_data( sid, &*riter, 1, &id ) )
00227             dim = 1;
00228         else if( MB_SUCCESS == iFace->tag_get_data( nid, &*riter, 1, &id ) )
00229             dim = 2;
00230         else if( MB_SUCCESS == iFace->tag_get_data( did, &*riter, 1, &dim ) )
00231         {
00232             id = 0;
00233             iFace->tag_get_data( gid, &*riter, 1, &id );
00234             dim += 3;
00235         }
00236         else
00237         {
00238             id  = *riter;
00239             dim = 9;
00240         }
00241 
00242         // oldsize = typesets[dim].size();
00243         typesets[dim].insert( id );
00244         // assert(typesets[dim].size() - oldsize == 1);
00245     }
00246     for( int ii = 0; ii < 9; ++ii )
00247     {
00248         char tmp[64];
00249         sprintf( tmp, "%s (%lu) ", typenames[ii], (unsigned long)typesets[ii].size() );
00250         str->print( VB, tmp, typesets[ii] );
00251     }
00252     str->printf( VB, "Total: %lu\n", (unsigned long)sets.size() );
00253 }
00254 
00255 #define debug_barrier() debug_barrier_line( __LINE__ )
00256 
00257 void WriteHDF5Parallel::debug_barrier_line( int lineno )
00258 {
00259     const unsigned threshold   = 2;
00260     static unsigned long count = 0;
00261     if( dbgOut.get_verbosity() >= threshold && myPcomm )
00262     {
00263         dbgOut.printf( threshold, "*********** Debug Barrier %lu (@%d)***********\n", ++count, lineno );
00264         MPI_Barrier( myPcomm->proc_config().proc_comm() );
00265     }
00266 }
00267 
00268 WriterIface* WriteHDF5Parallel::factory( Interface* iface )
00269 {
00270     return new WriteHDF5Parallel( iface );
00271 }
00272 
00273 WriteHDF5Parallel::WriteHDF5Parallel( Interface* iface )
00274     : WriteHDF5( iface ), myPcomm( NULL ), pcommAllocated( false ), hslabOp( H5S_SELECT_OR )
00275 {
00276 }
00277 
00278 WriteHDF5Parallel::~WriteHDF5Parallel()
00279 {
00280     if( pcommAllocated && myPcomm ) delete myPcomm;
00281 }
00282 
00283 // The parent WriteHDF5 class has ExportSet structs that are
00284 // populated with the entities to be written, grouped by type
00285 // (and for elements, connectivity length).  This function:
00286 //  o determines which entities are to be written by a remote processor
00287 //  o removes those entities from the ExportSet structs in WriteMesh
00288 //  o passes them back in a Range
00289 ErrorCode WriteHDF5Parallel::gather_interface_meshes( Range& nonowned )
00290 {
00291     ErrorCode result;
00292 
00293     // START_SERIAL;
00294     dbgOut.print( 3, "Pre-interface mesh:\n" );
00295     dbgOut.print( 3, nodeSet.range );
00296     for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
00297         dbgOut.print( 3, eiter->range );
00298     dbgOut.print( 3, setSet.range );
00299 
00300     // Move handles of non-owned entities from lists of entities
00301     // that this processor will write to the 'nonowned' list.
00302 
00303     nonowned.clear();
00304     result = myPcomm->filter_pstatus( nodeSet.range, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonowned );
00305     if( MB_SUCCESS != result ) return error( result );
00306     nodeSet.range = subtract( nodeSet.range, nonowned );
00307 
00308     for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
00309     {
00310         Range tmpset;
00311         result = myPcomm->filter_pstatus( eiter->range, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &tmpset );
00312         if( MB_SUCCESS != result ) return error( result );
00313         eiter->range = subtract( eiter->range, tmpset );
00314         nonowned.merge( tmpset );
00315     }
00316 
00317     dbgOut.print( 3, "Post-interface mesh:\n" );
00318     dbgOut.print( 3, nodeSet.range );
00319     for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
00320         dbgOut.print( 3, eiter->range );
00321     dbgOut.print( 3, setSet.range );
00322 
00323     // END_SERIAL;
00324 
00325     return MB_SUCCESS;
00326 }
00327 
00328 ErrorCode WriteHDF5Parallel::parallel_create_file( const char* filename, bool overwrite,
00329                                                    const std::vector< std::string >& qa_records,
00330                                                    const FileOptions& opts, const Tag* user_tag_list,
00331                                                    int user_tag_count, int dimension, double* times )
00332 {
00333     ErrorCode rval;
00334     mhdf_Status status;
00335 
00336     int pcomm_no = 0;
00337     opts.get_int_option( "PARALLEL_COMM", pcomm_no );
00338 
00339     myPcomm = ParallelComm::get_pcomm( iFace, pcomm_no );
00340     if( 0 == myPcomm )
00341     {
00342         myPcomm        = new ParallelComm( iFace, MPI_COMM_WORLD );
00343         pcommAllocated = true;
00344     }
00345 
00346     MPI_Info info = MPI_INFO_NULL;
00347     std::string cb_size;
00348     rval = opts.get_str_option( "CB_BUFFER_SIZE", cb_size );
00349     if( MB_SUCCESS == rval )
00350     {
00351         MPI_Info_create( &info );
00352         MPI_Info_set( info, const_cast< char* >( "cb_buffer_size" ), const_cast< char* >( cb_size.c_str() ) );
00353     }
00354 
00355     dbgOut.set_rank( myPcomm->proc_config().proc_rank() );
00356     dbgOut.limit_output_to_first_N_procs( 32 );
00357     Range nonlocal;
00358     debug_barrier();
00359     dbgOut.tprint( 1, "Gathering interface meshes\n" );
00360     rval = gather_interface_meshes( nonlocal );
00361     if( MB_SUCCESS != rval ) return error( rval );
00362 
00363     /**************** Get tag names for sets likely to be shared ***********/
00364     // debug_barrier();
00365     // dbgOut.tprint(1, "Getting shared entity sets\n");
00366     // rval = get_sharedset_tags();
00367     // if (MB_SUCCESS != rval) return error(rval);
00368 
00369     /**************** Create actual file and write meta info ***************/
00370 
00371     debug_barrier();
00372     if( myPcomm->proc_config().proc_rank() == 0 )
00373     {
00374         dbgOut.tprintf( 1, "Creating file: %s\n", filename );
00375 
00376         // Create the file
00377         const char* type_names[MBMAXTYPE];
00378         memset( type_names, 0, MBMAXTYPE * sizeof( char* ) );
00379         for( EntityType i = MBEDGE; i < MBENTITYSET; ++i )
00380             type_names[i] = CN::EntityTypeName( i );
00381 
00382         dbgOut.tprint( 1, "call mhdf_createFile\n" );
00383         filePtr = mhdf_createFile( filename, overwrite, type_names, MBMAXTYPE, id_type, &status );
00384         if( !filePtr ) { MB_SET_ERR( MB_FAILURE, mhdf_message( &status ) ); }
00385 
00386         dbgOut.tprint( 1, "call write_qa\n" );
00387         rval = write_qa( qa_records );
00388         if( MB_SUCCESS != rval ) return error( rval );
00389     }
00390 
00391     /**************** Create node coordinate table ***************/
00392     CpuTimer timer;
00393     debug_barrier();
00394     dbgOut.tprint( 1, "creating node table\n" );
00395     topState.start( "creating node table" );
00396     rval = create_node_table( dimension );
00397     topState.end( rval );
00398     if( MB_SUCCESS != rval ) return error( rval );
00399     if( times ) times[CREATE_NODE_TIME] = timer.time_elapsed();
00400 
00401     /**************** Create element tables ***************/
00402 
00403     debug_barrier();
00404     dbgOut.tprint( 1, "negotiating element types\n" );
00405     topState.start( "negotiating element types" );
00406     rval = negotiate_type_list();
00407     topState.end( rval );
00408     if( MB_SUCCESS != rval ) return error( rval );
00409     if( times ) times[NEGOTIATE_TYPES_TIME] = timer.time_elapsed();
00410     dbgOut.tprint( 1, "creating element table\n" );
00411     topState.start( "creating element tables" );
00412     rval = create_element_tables();
00413     topState.end( rval );
00414     if( MB_SUCCESS != rval ) return error( rval );
00415     if( times ) times[CREATE_ELEM_TIME] = timer.time_elapsed();
00416 
00417     /*************** Exchange file IDs *****************/
00418 
00419     debug_barrier();
00420     dbgOut.tprint( 1, "communicating file ids\n" );
00421     topState.start( "communicating file ids" );
00422     rval = exchange_file_ids( nonlocal );
00423     topState.end( rval );
00424     if( MB_SUCCESS != rval ) return error( rval );
00425     if( times ) times[FILEID_EXCHANGE_TIME] = timer.time_elapsed();
00426 
00427     /**************** Create meshset tables *********************/
00428 
00429     debug_barrier();
00430     dbgOut.tprint( 1, "creating meshset table\n" );
00431     topState.start( "creating meshset tables" );
00432     rval = create_meshset_tables( times );
00433     topState.end( rval );
00434     if( MB_SUCCESS != rval ) return error( rval );
00435     if( times ) times[CREATE_SET_TIME] = timer.time_elapsed();
00436 
00437     /**************** Create adjacency tables *********************/
00438 
00439     debug_barrier();
00440     dbgOut.tprint( 1, "creating adjacency table\n" );
00441     topState.start( "creating adjacency tables" );
00442     rval = create_adjacency_tables();
00443     topState.end( rval );
00444     if( MB_SUCCESS != rval ) return error( rval );
00445     if( times ) times[CREATE_ADJ_TIME] = timer.time_elapsed();
00446 
00447     /**************** Create tag data *********************/
00448 
00449     debug_barrier();
00450     dbgOut.tprint( 1, "creating tag tables\n" );
00451     topState.start( "creating tag tables" );
00452     rval = gather_tags( user_tag_list, user_tag_count );
00453     if( MB_SUCCESS != rval ) return error( rval );
00454     rval = create_tag_tables();
00455     topState.end( rval );
00456     if( MB_SUCCESS != rval ) return error( rval );
00457     if( times ) times[CREATE_TAG_TIME] = timer.time_elapsed();
00458 
00459     /************** Close serial file and reopen parallel *****************/
00460 
00461     if( 0 == myPcomm->proc_config().proc_rank() ) mhdf_closeFile( filePtr, &status );
00462 
00463     MPI_Barrier( myPcomm->proc_config().proc_comm() );
00464     dbgOut.tprint( 1, "(re)opening file in parallel mode\n" );
00465     unsigned long junk;
00466     hid_t hdf_opt = H5Pcreate( H5P_FILE_ACCESS );
00467     H5Pset_fapl_mpio( hdf_opt, myPcomm->proc_config().proc_comm(), info );
00468     filePtr = mhdf_openFileWithOpt( filename, 1, &junk, id_type, hdf_opt, &status );
00469     H5Pclose( hdf_opt );
00470     if( !filePtr ) { MB_SET_ERR( MB_FAILURE, mhdf_message( &status ) ); }
00471 
00472     if( collectiveIO )
00473     {
00474         dbgOut.print( 1, "USING COLLECTIVE IO\n" );
00475         writeProp = H5Pcreate( H5P_DATASET_XFER );
00476         H5Pset_dxpl_mpio( writeProp, H5FD_MPIO_COLLECTIVE );
00477     }
00478 
00479     /* Test if we can use H5S_APPEND when selecting hyperslabs */
00480     if( MB_SUCCESS != opts.get_null_option( "HYPERSLAB_OR" ) &&
00481         ( MB_SUCCESS == opts.get_null_option( "HYPERSLAB_APPEND" ) || HDF5_can_append_hyperslabs() ) )
00482     {
00483         dbgOut.print( 1, "HDF5 library supports H5Sselect_hyperlsab with H5S_SELECT_APPEND\n" );
00484         hslabOp = H5S_SELECT_APPEND;
00485     }
00486 
00487     dbgOut.tprint( 1, "Exiting parallel_create_file\n" );
00488     return MB_SUCCESS;
00489 }
00490 
00491 class TagNameCompare
00492 {
00493     Interface* iFace;
00494     std::string name1, name2;
00495 
00496   public:
00497     TagNameCompare( Interface* iface ) : iFace( iface ) {}
00498     bool operator()( const WriteHDF5::TagDesc& t1, const WriteHDF5::TagDesc& t2 );
00499 };
00500 
00501 bool TagNameCompare::operator()( const WriteHDF5::TagDesc& t1, const WriteHDF5::TagDesc& t2 )
00502 {
00503     iFace->tag_get_name( t1.tag_id, name1 );
00504     iFace->tag_get_name( t2.tag_id, name2 );
00505     return name1 < name2;
00506 }
00507 
00508 struct serial_tag_data
00509 {
00510     TagType storage;
00511     DataType type;
00512     int size;
00513     int name_len;
00514     int def_val_len;
00515     char name[sizeof( unsigned long )];
00516 
00517     static size_t pad( size_t len )
00518     {
00519         if( len % sizeof( unsigned long ) )
00520             return len + sizeof( unsigned long ) - len % sizeof( unsigned long );
00521         else
00522             return len;
00523     }
00524 
00525     static size_t def_val_bytes( int def_val_len, DataType type )
00526     {
00527         switch( type )
00528         {
00529             case MB_TYPE_BIT:
00530                 return def_val_len ? 1 : 0;
00531             case MB_TYPE_OPAQUE:
00532                 return def_val_len;
00533             case MB_TYPE_INTEGER:
00534                 return def_val_len * sizeof( int );
00535             case MB_TYPE_DOUBLE:
00536                 return def_val_len * sizeof( double );
00537             case MB_TYPE_HANDLE:
00538                 return def_val_len * sizeof( EntityHandle );
00539         }
00540         return 0;
00541     }
00542 
00543     static size_t len( int name_len, int def_val_len, DataType type )
00544     {
00545         return sizeof( serial_tag_data ) + pad( name_len + def_val_bytes( def_val_len, type ) ) -
00546                sizeof( unsigned long );
00547     }
00548     size_t len() const
00549     {
00550         return len( name_len, def_val_len, type );
00551     }
00552     void* default_value()
00553     {
00554         return def_val_len ? name + name_len : 0;
00555     }
00556     const void* default_value() const
00557     {
00558         return const_cast< serial_tag_data* >( this )->default_value();
00559     }
00560     void set_default_value( const void* val )
00561     {
00562         memcpy( default_value(), val, def_val_bytes( def_val_len, type ) );
00563     }
00564 };
00565 
00566 ErrorCode WriteHDF5Parallel::append_serial_tag_data( std::vector< unsigned char >& buffer,
00567                                                      const WriteHDF5::TagDesc& tag )
00568 {
00569     ErrorCode rval;
00570 
00571     std::string name;
00572     rval = iFace->tag_get_name( tag.tag_id, name );
00573     if( MB_SUCCESS != rval ) return error( rval );
00574 
00575     // Get name length, including space for null char
00576     size_t name_len = name.size() + 1;
00577     if( name_len == 1 ) return MB_SUCCESS;  // Skip tags with no name
00578 
00579     DataType data_type;
00580     rval = iFace->tag_get_data_type( tag.tag_id, data_type );
00581     if( MB_SUCCESS != rval ) return error( rval );
00582 
00583     // Get default value
00584     int def_val_len;
00585     const void* def_val;
00586     if( MB_SUCCESS != iFace->tag_get_default_value( tag.tag_id, def_val, def_val_len ) )
00587     {
00588         def_val_len = 0;
00589         def_val     = 0;
00590     }
00591 
00592     // Allocate struct within buffer
00593     size_t init_size = buffer.size();
00594     buffer.resize( init_size + serial_tag_data::len( name_len, def_val_len, data_type ) );
00595     serial_tag_data* ptr = reinterpret_cast< serial_tag_data* >( &buffer[init_size] );
00596 
00597     // Populate struct
00598     rval = iFace->tag_get_type( tag.tag_id, ptr->storage );
00599     if( MB_SUCCESS != rval ) return error( rval );
00600     ptr->type = data_type;
00601     rval      = iFace->tag_get_length( tag.tag_id, ptr->size );
00602     if( MB_VARIABLE_DATA_LENGTH == rval )
00603         ptr->size = MB_VARIABLE_LENGTH;
00604     else if( MB_SUCCESS != rval )
00605         return error( rval );
00606     ptr->name_len = name_len;
00607     Range range;
00608     memset( ptr->name, 0, ptr->name_len );
00609     memcpy( ptr->name, name.data(), name.size() );
00610     ptr->def_val_len = def_val_len;
00611     ptr->set_default_value( def_val );
00612 
00613     return MB_SUCCESS;
00614 }
00615 
00616 ErrorCode WriteHDF5Parallel::check_serial_tag_data( const std::vector< unsigned char >& buffer,
00617                                                     std::vector< TagDesc* >* missing, std::vector< TagDesc* >* newlist )
00618 {
00619     ErrorCode rval;
00620 
00621     // Use 'write_sparse' field as a 'visited' mark
00622     std::list< TagDesc >::iterator tag_iter;
00623     if( missing )
00624         for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
00625             tag_iter->write_sparse = true;
00626 
00627     // Use a set as a temporary for what will ultimately go in
00628     // newlist because we need to pass back newlist in the order
00629     // of the tagList member.
00630     std::set< TagDesc* > newset;
00631 
00632     // Iterate over data from, updating the local list of tags.
00633     // Be careful to keep tagList sorted such that in the end all
00634     // procs have the same list in the same order.
00635     std::vector< unsigned char >::const_iterator diter = buffer.begin();
00636     tag_iter                                           = tagList.begin();
00637     while( diter < buffer.end() )
00638     {
00639         // Get struct from buffer
00640         const serial_tag_data* ptr = reinterpret_cast< const serial_tag_data* >( &*diter );
00641 
00642         // Find local struct for tag
00643         std::string name( ptr->name );
00644         std::string n;
00645         iFace->tag_get_name( tag_iter->tag_id, n );  // Second time we've called, so shouldn't fail
00646         if( n > name )
00647         {
00648             tag_iter = tagList.begin();  // New proc, start search from beginning
00649         }
00650         iFace->tag_get_name( tag_iter->tag_id, n );
00651         while( n < name )
00652         {
00653             ++tag_iter;
00654             if( tag_iter == tagList.end() ) break;
00655             iFace->tag_get_name( tag_iter->tag_id, n );
00656         }
00657         if( tag_iter == tagList.end() || n != name )
00658         {  // New tag
00659             TagDesc newtag;
00660 
00661             if( ptr->size == MB_VARIABLE_LENGTH )
00662                 rval = iFace->tag_get_handle( name.c_str(), ptr->def_val_len, ptr->type, newtag.tag_id,
00663                                               MB_TAG_VARLEN | MB_TAG_CREAT | ptr->storage, ptr->default_value() );
00664             else
00665                 rval = iFace->tag_get_handle( name.c_str(), ptr->size, ptr->type, newtag.tag_id,
00666                                               MB_TAG_CREAT | ptr->storage, ptr->default_value() );
00667             if( MB_SUCCESS != rval ) return error( rval );
00668 
00669             newtag.sparse_offset   = 0;
00670             newtag.var_data_offset = 0;
00671             newtag.write_sparse    = false;
00672             newtag.max_num_ents    = 0;
00673             newtag.max_num_vals    = 0;
00674 
00675             tag_iter = tagList.insert( tag_iter, newtag );
00676             if( newlist ) newset.insert( &*tag_iter );
00677         }
00678         else
00679         {  // Check that tag is as expected
00680             DataType type;
00681             iFace->tag_get_data_type( tag_iter->tag_id, type );
00682             if( type != ptr->type )
00683             {
00684                 MB_SET_ERR( MB_FAILURE, "Processes have inconsistent data type for tag \"" << name << "\"" );
00685             }
00686             int size;
00687             iFace->tag_get_length( tag_iter->tag_id, size );
00688             if( size != ptr->size )
00689             {
00690                 MB_SET_ERR( MB_FAILURE, "Processes have inconsistent size for tag \"" << name << "\"" );
00691             }
00692             tag_iter->write_sparse = false;
00693         }
00694 
00695         // Step to next variable-length struct.
00696         diter += ptr->len();
00697     }
00698 
00699     // Now pass back any local tags that weren't in the buffer
00700     if( missing )
00701     {
00702         for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
00703         {
00704             if( tag_iter->write_sparse )
00705             {
00706                 tag_iter->write_sparse = false;
00707                 missing->push_back( &*tag_iter );
00708             }
00709         }
00710     }
00711 
00712     // Be careful to populate newlist in the same, sorted, order as tagList
00713     if( newlist )
00714     {
00715         for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
00716             if( newset.find( &*tag_iter ) != newset.end() ) newlist->push_back( &*tag_iter );
00717     }
00718 
00719     return MB_SUCCESS;
00720 }
00721 
00722 static void set_bit( int position, unsigned char* bytes )
00723 {
00724     int byte = position / 8;
00725     int bit  = position % 8;
00726     bytes[byte] |= ( ( (unsigned char)1 ) << bit );
00727 }
00728 
00729 static bool get_bit( int position, const unsigned char* bytes )
00730 {
00731     int byte = position / 8;
00732     int bit  = position % 8;
00733     return 0 != ( bytes[byte] & ( ( (unsigned char)1 ) << bit ) );
00734 }
00735 
00736 ErrorCode WriteHDF5Parallel::create_tag_tables()
00737 {
00738     std::list< TagDesc >::iterator tag_iter;
00739     ErrorCode rval;
00740     int err;
00741     const int rank      = myPcomm->proc_config().proc_rank();
00742     const MPI_Comm comm = myPcomm->proc_config().proc_comm();
00743 
00744     subState.start( "negotiating tag list" );
00745 
00746     dbgOut.tprint( 1, "communicating tag metadata\n" );
00747 
00748     dbgOut.printf( 2, "Exchanging tag data for %d tags.\n", (int)tagList.size() );
00749 
00750     // Sort tagList contents in alphabetical order by tag name
00751     tagList.sort( TagNameCompare( iFace ) );
00752 
00753     // Negotiate total list of tags to write
00754 
00755     // Build concatenated list of all tag data
00756     std::vector< unsigned char > tag_buffer;
00757     for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
00758     {
00759         rval = append_serial_tag_data( tag_buffer, *tag_iter );
00760         CHECK_MB( rval );
00761     }
00762 
00763     // Broadcast list from root to all other procs
00764     unsigned long size = tag_buffer.size();
00765     err                = MPI_Bcast( &size, 1, MPI_UNSIGNED_LONG, 0, comm );
00766     CHECK_MPI( err );
00767     tag_buffer.resize( size );
00768     err = MPI_Bcast( &tag_buffer[0], size, MPI_UNSIGNED_CHAR, 0, comm );
00769     CHECK_MPI( err );
00770 
00771     // Update local tag list
00772     std::vector< TagDesc* > missing;
00773     rval = check_serial_tag_data( tag_buffer, &missing, 0 );
00774     CHECK_MB( rval );
00775 
00776     // Check if we're done (0->done, 1->more, 2+->error)
00777     int code, lcode = ( MB_SUCCESS != rval ) ? rval + 2 : missing.empty() ? 0 : 1;
00778     err = MPI_Allreduce( &lcode, &code, 1, MPI_INT, MPI_MAX, comm );
00779     CHECK_MPI( err );
00780     if( code > 1 )
00781     {
00782         MB_SET_ERR_CONT( "Inconsistent tag definitions between procs" );
00783         return error( ( ErrorCode )( code - 2 ) );
00784     }
00785 
00786     // If not done...
00787     if( code )
00788     {
00789         dbgOut.print( 1, "Not all procs had same tag definitions, negotiating...\n" );
00790 
00791         // Get tags defined on this proc but not on root proc
00792         tag_buffer.clear();
00793         for( size_t i = 0; i < missing.size(); ++i )
00794         {
00795             rval = append_serial_tag_data( tag_buffer, *missing[i] );
00796             CHECK_MB( rval );
00797         }
00798 
00799         // Gather extra tag definitions on root processor
00800         std::vector< int > junk;               // don't care how many from each proc
00801         assert( rank || tag_buffer.empty() );  // must be empty on root
00802         err = my_Gatherv( &tag_buffer[0], tag_buffer.size(), MPI_UNSIGNED_CHAR, tag_buffer, junk, 0, comm );
00803         CHECK_MPI( err );
00804 
00805         // Process serialized tag descriptions on root, and
00806         rval = MB_SUCCESS;
00807         if( 0 == rank )
00808         {
00809             // Process serialized tag descriptions on root, and
00810             std::vector< TagDesc* > newlist;
00811             rval = check_serial_tag_data( tag_buffer, 0, &newlist );
00812             tag_buffer.clear();
00813             // re-serialize a unique list of new tag definitions
00814             for( size_t i = 0; MB_SUCCESS == rval && i != newlist.size(); ++i )
00815             {
00816                 rval = append_serial_tag_data( tag_buffer, *newlist[i] );
00817                 CHECK_MB( rval );
00818             }
00819         }
00820 
00821         // Broadcast any new tag definitions from root to other procs
00822         long this_size = tag_buffer.size();
00823         if( MB_SUCCESS != rval ) this_size = -rval;
00824         err = MPI_Bcast( &this_size, 1, MPI_LONG, 0, comm );
00825         CHECK_MPI( err );
00826         if( this_size < 0 )
00827         {
00828             MB_SET_ERR_CONT( "Inconsistent tag definitions between procs" );
00829             return error( (ErrorCode)-this_size );
00830         }
00831         tag_buffer.resize( this_size );
00832         err = MPI_Bcast( &tag_buffer[0], this_size, MPI_UNSIGNED_CHAR, 0, comm );
00833         CHECK_MPI( err );
00834 
00835         // Process new tag definitions
00836         rval = check_serial_tag_data( tag_buffer, 0, 0 );
00837         CHECK_MB( rval );
00838     }
00839 
00840     subState.end();
00841     subState.start( "negotiate which element/tag combinations are dense" );
00842 
00843     // Figure out for which tag/element combinations we can
00844     // write dense tag data.
00845 
00846     // Construct a table of bits,
00847     // where each row of the table corresponds to a tag
00848     // and each column to an element group.
00849 
00850     // Two extra, because first is nodes and last is sets.
00851     // (n+7)/8 is ceil(n/8)
00852     const int bytes_per_tag = ( exportList.size() + 9 ) / 8;
00853     std::vector< unsigned char > data( bytes_per_tag * tagList.size(), 0 );
00854     std::vector< unsigned char > recv( data.size(), 0 );
00855     unsigned char* iter = &data[0];
00856     if( writeTagDense && !data.empty() )
00857     {
00858         for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, iter += bytes_per_tag )
00859         {
00860 
00861             Range tagged;
00862             rval = get_sparse_tagged_entities( *tag_iter, tagged );
00863             CHECK_MB( rval );
00864 
00865             int s;
00866             if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) ) continue;
00867 
00868             std::string n;
00869             iFace->tag_get_name( tag_iter->tag_id,
00870                                  n );  // Second time we've called, so shouldn't fail
00871 
00872             // Check if we want to write this tag in dense format even if not
00873             // all of the entities have a tag value.  The criterion of this
00874             // is that the tag be dense, have a default value, and have at
00875             // least 2/3 of the entities tagged.
00876             bool prefer_dense = false;
00877             TagType type;
00878             rval = iFace->tag_get_type( tag_iter->tag_id, type );
00879             CHECK_MB( rval );
00880             if( MB_TAG_DENSE == type )
00881             {
00882                 const void* defval = 0;
00883                 rval               = iFace->tag_get_default_value( tag_iter->tag_id, defval, s );
00884                 if( MB_SUCCESS == rval ) prefer_dense = true;
00885             }
00886 
00887             int i = 0;
00888             if( check_dense_format_tag( nodeSet, tagged, prefer_dense ) )
00889             {
00890                 set_bit( i, iter );
00891                 dbgOut.printf( 2, "Can write dense data for \"%s\"/Nodes\n", n.c_str() );
00892             }
00893             std::list< ExportSet >::const_iterator ex_iter = exportList.begin();
00894             for( ++i; ex_iter != exportList.end(); ++i, ++ex_iter )
00895             {
00896                 // when writing in parallel, on some partitions, some of these element ranges might
00897                 // be empty so do not turn this tag as sparse, just because of that, leave it dense,
00898                 // if we prefer dense
00899                 if( ( prefer_dense && ex_iter->range.empty() ) ||
00900                     check_dense_format_tag( *ex_iter, tagged, prefer_dense ) )
00901                 {
00902                     set_bit( i, iter );
00903                     dbgOut.printf( 2, "Can write dense data for \"%s\"/%s\n", n.c_str(), ex_iter->name() );
00904                 }
00905             }
00906             if( check_dense_format_tag( setSet, tagged, prefer_dense ) )
00907             {
00908                 set_bit( i, iter );
00909                 dbgOut.printf( 2, "Can write dense data for \"%s\"/Sets\n", n.c_str() );
00910             }
00911         }
00912 
00913         // Do bit-wise AND of list over all processors (only write dense format
00914         // if all proccesses want dense format for this group of entities).
00915         err = MPI_Allreduce( &data[0], &recv[0], data.size(), MPI_UNSIGNED_CHAR, MPI_BAND,
00916                              myPcomm->proc_config().proc_comm() );
00917         CHECK_MPI( err );
00918     }  // if (writeTagDense)
00919 
00920     // Store initial counts for sparse-formatted tag data.
00921     // The total number of values to send and receive will be the number of
00922     // tags plus the number of var-len tags because we need to negotiate
00923     // offsets into two different tables for the var-len tags.
00924     std::vector< long > counts;
00925 
00926     // Record dense tag/element combinations
00927     iter                       = &recv[0];
00928     const unsigned char* iter2 = &data[0];
00929     for( tag_iter = tagList.begin(); tag_iter != tagList.end();
00930          ++tag_iter, iter += bytes_per_tag, iter2 += bytes_per_tag )
00931     {
00932 
00933         Range tagged;
00934         rval = get_sparse_tagged_entities( *tag_iter, tagged );
00935         CHECK_MB( rval );
00936 
00937         std::string n;
00938         iFace->tag_get_name( tag_iter->tag_id, n );  // Second time we've called, so shouldn't fail
00939 
00940         int i = 0;
00941         if( get_bit( i, iter ) )
00942         {
00943             assert( get_bit( i, iter2 ) );
00944             tag_iter->dense_list.push_back( nodeSet );
00945             tagged -= nodeSet.range;
00946             dbgOut.printf( 2, "Will write dense data for \"%s\"/Nodes\n", n.c_str() );
00947         }
00948         std::list< ExportSet >::const_iterator ex_iter = exportList.begin();
00949         for( ++i; ex_iter != exportList.end(); ++i, ++ex_iter )
00950         {
00951             if( get_bit( i, iter ) )
00952             {
00953                 assert( get_bit( i, iter2 ) );
00954                 tag_iter->dense_list.push_back( *ex_iter );
00955                 dbgOut.printf( 2, "WIll write dense data for \"%s\"/%s\n", n.c_str(), ex_iter->name() );
00956                 tagged -= ex_iter->range;
00957             }
00958         }
00959         if( get_bit( i, iter ) )
00960         {
00961             assert( get_bit( i, iter2 ) );
00962             tag_iter->dense_list.push_back( setSet );
00963             dbgOut.printf( 2, "Will write dense data for \"%s\"/Sets\n", n.c_str() );
00964             tagged -= setSet.range;
00965         }
00966 
00967         counts.push_back( tagged.size() );
00968 
00969         int s;
00970         if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
00971         {
00972             unsigned long data_len;
00973             rval = get_tag_data_length( *tag_iter, tagged, data_len );
00974             CHECK_MB( rval );
00975             counts.push_back( data_len );
00976         }
00977     }
00978 
00979     subState.end();
00980     subState.start( "Negotiate offsets for sparse tag info" );
00981 
00982     std::vector< long > offsets( counts.size() ), maxima( counts.size() ), totals( counts.size() );
00983     rval = create_dataset( counts.size(), &counts[0], &offsets[0], &maxima[0], &totals[0] );
00984     CHECK_MB( rval );
00985 
00986     // Copy values into local structs and if root then create tables
00987     size_t idx = 0;
00988     for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, ++idx )
00989     {
00990         assert( idx < counts.size() );
00991         tag_iter->sparse_offset = offsets[idx];
00992         tag_iter->max_num_ents  = maxima[idx];
00993         tag_iter->write_sparse  = ( 0 != totals[idx] );
00994         int s;
00995         if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
00996         {
00997             ++idx;
00998             assert( idx < counts.size() );
00999             tag_iter->var_data_offset = offsets[idx];
01000             tag_iter->max_num_vals    = maxima[idx];
01001         }
01002         else
01003         {
01004             tag_iter->var_data_offset = 0;
01005             tag_iter->max_num_vals    = 0;
01006         }
01007     }
01008 
01009     subState.end();
01010 
01011     // Create tag tables on root process
01012     if( 0 == myPcomm->proc_config().proc_rank() )
01013     {
01014         size_t iidx = 0;
01015         for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter, ++iidx )
01016         {
01017             assert( iidx < totals.size() );
01018             unsigned long num_ents = totals[iidx];
01019             unsigned long num_val  = 0;
01020             int s;
01021             if( MB_VARIABLE_DATA_LENGTH == iFace->tag_get_length( tag_iter->tag_id, s ) )
01022             {
01023                 ++iidx;
01024                 assert( iidx < totals.size() );
01025                 num_val = totals[iidx];
01026             }
01027             dbgOut.printf( 2, "Writing tag description for tag 0x%lx with %lu values\n",
01028                            (unsigned long)tag_iter->tag_id, num_val ? num_val : num_ents );
01029 
01030             rval = create_tag( *tag_iter, num_ents, num_val );
01031             if( MB_SUCCESS != rval ) return error( rval );
01032         }
01033     }
01034 
01035     if( dbgOut.get_verbosity() > 1 )
01036     {
01037         dbgOut.printf( 2, "Tags: %12s %8s %8s %8s %8s %8s\n", "Name", "Count", "Offset", "Var Off", "Max Ent",
01038                        "Handle" );
01039 
01040         for( tag_iter = tagList.begin(); tag_iter != tagList.end(); ++tag_iter )
01041         {
01042             std::string name;
01043             iFace->tag_get_name( tag_iter->tag_id, name );
01044             size_t this_size;
01045             get_num_sparse_tagged_entities( *tag_iter, this_size );
01046             dbgOut.printf( 2, "%18s %8lu %8lu %8lu %8lu 0x%7lx\n", name.c_str(), (unsigned long)this_size,
01047                            (unsigned long)tag_iter->sparse_offset, (unsigned long)tag_iter->var_data_offset,
01048                            (unsigned long)tag_iter->max_num_ents, (unsigned long)tag_iter->tag_id );
01049         }
01050     }
01051 
01052     return MB_SUCCESS;
01053 }
01054 
01055 struct DatasetVals
01056 {
01057     long start_id;
01058     long max_count;
01059     long total;
01060 };
01061 STATIC_ASSERT( ( sizeof( DatasetVals ) == 3 * sizeof( long ) ) );
01062 
01063 ErrorCode WriteHDF5Parallel::create_dataset( int num_datasets, const long* num_owned, long* offsets_out,
01064                                              long* max_proc_entities, long* total_entities,
01065                                              const DataSetCreator& creator, ExportSet* groups[], wid_t* first_ids_out )
01066 {
01067     int result;
01068     ErrorCode rval;
01069     const unsigned rank  = myPcomm->proc_config().proc_rank();
01070     const unsigned nproc = myPcomm->proc_config().proc_size();
01071     const MPI_Comm comm  = myPcomm->proc_config().proc_comm();
01072 
01073     // Gather entity counts for each processor on root
01074     std::vector< long > counts( rank ? 0 : nproc * num_datasets );
01075     (void)VALGRIND_CHECK_MEM_IS_DEFINED( &num_owned, sizeof( long ) );
01076     result = MPI_Gather( const_cast< long* >( num_owned ), num_datasets, MPI_LONG, &counts[0], num_datasets, MPI_LONG,
01077                          0, comm );
01078     CHECK_MPI( result );
01079 
01080     // Create node data in file
01081     DatasetVals zero_val = { 0, 0, 0 };
01082     std::vector< DatasetVals > cumulative( num_datasets, zero_val );
01083     if( rank == 0 )
01084     {
01085         for( unsigned i = 0; i < nproc; i++ )
01086         {
01087             const long* proc_data = &counts[i * num_datasets];
01088             for( int index = 0; index < num_datasets; ++index )
01089             {
01090                 cumulative[index].total += proc_data[index];
01091                 if( proc_data[index] > cumulative[index].max_count ) cumulative[index].max_count = proc_data[index];
01092             }
01093         }
01094 
01095         for( int index = 0; index < num_datasets; ++index )
01096         {
01097             if( cumulative[index].total )
01098             {
01099                 rval = creator( this, cumulative[index].total, groups ? groups[index] : 0, cumulative[index].start_id );
01100                 CHECK_MB( rval );
01101             }
01102             else
01103             {
01104                 cumulative[index].start_id = -1;
01105             }
01106         }
01107     }
01108 
01109     // Send id offset to every proc
01110     result = MPI_Bcast( (void*)&cumulative[0], 3 * num_datasets, MPI_LONG, 0, comm );
01111     CHECK_MPI( result );
01112     for( int index = 0; index < num_datasets; ++index )
01113     {
01114         if( first_ids_out ) first_ids_out[index] = (wid_t)cumulative[index].start_id;
01115         max_proc_entities[index] = cumulative[index].max_count;
01116         total_entities[index]    = cumulative[index].total;
01117     }
01118 
01119     // Convert array of per-process counts to per-process offsets
01120     if( rank == 0 )
01121     {
01122         // Initialize prev_size with data sizes for root process
01123         std::vector< long > prev_size( counts.begin(), counts.begin() + num_datasets );
01124         // Root process gets offset zero
01125         std::fill( counts.begin(), counts.begin() + num_datasets, 0L );
01126         // For each proc other than this one (root)
01127         for( unsigned i = 1; i < nproc; ++i )
01128         {
01129             // Get pointer to offsets for previous process in list
01130             long* prev_data = &counts[( i - 1 ) * num_datasets];
01131             // Get pointer to offsets for this process in list
01132             long* proc_data = &counts[i * num_datasets];
01133             // For each data set
01134             for( int j = 0; j < num_datasets; ++j )
01135             {
01136                 // Get size of data in dataset from process i
01137                 long mysize = proc_data[j];
01138                 // Offset for process i is offset of previous process plus
01139                 // number of values previous process will write
01140                 proc_data[j] = prev_data[j] + prev_size[j];
01141                 // Store my size, as it is no longer available in 'counts'
01142                 prev_size[j] = mysize;
01143             }
01144         }
01145     }
01146 
01147     // Send each proc it's offset in the table
01148     if( rank == 0 ) { (void)VALGRIND_CHECK_MEM_IS_DEFINED( &counts[0], num_datasets * nproc * sizeof( long ) ); }
01149     result = MPI_Scatter( &counts[0], num_datasets, MPI_LONG, offsets_out, num_datasets, MPI_LONG, 0, comm );
01150     CHECK_MPI( result );
01151 
01152     return MB_SUCCESS;
01153 }
01154 
01155 ErrorCode WriteHDF5Parallel::create_node_table( int dimension )
01156 {
01157     nodeSet.num_nodes = dimension;  // Put it here so NodeSetCreator can access it
01158     struct NodeSetCreator : public DataSetCreator
01159     {
01160         ErrorCode operator()( WriteHDF5* file, long count, const ExportSet* group, long& start_id ) const
01161         {
01162             mhdf_Status status;
01163             hid_t handle = mhdf_createNodeCoords( file->file_ptr(), group->num_nodes, count, &start_id, &status );
01164             CHECK_HDFN( status );
01165             mhdf_closeData( file->file_ptr(), handle, &status );
01166             CHECK_HDFN( status );
01167             return MB_SUCCESS;
01168         }
01169     };
01170 
01171     const long count   = nodeSet.range.size();
01172     ExportSet* array[] = { &nodeSet };
01173     ErrorCode rval     = create_dataset( 1, &count, &nodeSet.offset, &nodeSet.max_num_ents, &nodeSet.total_num_ents,
01174                                      NodeSetCreator(), array, &nodeSet.first_id );
01175     CHECK_MB( rval );
01176     return assign_ids( nodeSet.range, nodeSet.first_id + nodeSet.offset );
01177 }
01178 
01179 struct elemtype
01180 {
01181     int mbtype;
01182     int numnode;
01183 
01184     elemtype( int vals[2] ) : mbtype( vals[0] ), numnode( vals[1] ) {}
01185     elemtype( int t, int n ) : mbtype( t ), numnode( n ) {}
01186 
01187     bool operator==( const elemtype& other ) const
01188     {
01189         return mbtype == other.mbtype && ( mbtype == MBENTITYSET || numnode == other.numnode );
01190     }
01191     bool operator<( const elemtype& other ) const
01192     {
01193         if( mbtype > other.mbtype ) return false;
01194 
01195         return mbtype < other.mbtype || ( mbtype != MBENTITYSET && numnode < other.numnode );
01196     }
01197     bool operator!=( const elemtype& other ) const
01198     {
01199         return !this->operator==( other );
01200     }
01201 };
01202 
01203 ErrorCode WriteHDF5Parallel::negotiate_type_list()
01204 {
01205     int result;
01206     const MPI_Comm comm = myPcomm->proc_config().proc_comm();
01207 
01208     exportList.sort();
01209     int num_types = exportList.size();
01210 
01211     // Get list of types on this processor
01212     typedef std::vector< std::pair< int, int > > typelist;
01213     typelist my_types( num_types );
01214     (void)VALGRIND_MAKE_VEC_UNDEFINED( my_types );
01215     typelist::iterator viter = my_types.begin();
01216     for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
01217     {
01218         viter->first  = eiter->type;
01219         viter->second = eiter->num_nodes;
01220         ++viter;
01221     }
01222 
01223     dbgOut.print( 2, "Local Element Types:\n" );
01224     for( viter = my_types.begin(); viter != my_types.end(); ++viter )
01225     {
01226         int type  = viter->first;
01227         int count = viter->second;
01228         dbgOut.printf( 2, "  %s : %d\n", CN::EntityTypeName( (EntityType)type ), count );
01229     }
01230 
01231     // Broadcast number of types from root to all nodes
01232     int num_types0 = num_types;
01233     result         = MPI_Bcast( &num_types0, 1, MPI_INT, 0, comm );
01234     CHECK_MPI( result );
01235     // Broadcast type list from root to all nodes
01236     typelist root_types( num_types0 );
01237     if( 0 == myPcomm->proc_config().proc_rank() ) root_types = my_types;
01238     result = MPI_Bcast( (void*)&root_types[0], 2 * num_types0, MPI_INT, 0, comm );
01239     CHECK_MPI( result );
01240 
01241     // Build local list of any types that root did not know about
01242     typelist non_root_types;
01243     viter = root_types.begin();
01244     for( typelist::iterator iter = my_types.begin(); iter != my_types.end(); ++iter )
01245     {
01246         if( viter == root_types.end() || *viter != *iter )
01247             non_root_types.push_back( *iter );
01248         else
01249             ++viter;
01250     }
01251 
01252     // Determine if any process had types not defined on the root
01253     int non_root_count = non_root_types.size();
01254     int not_done;
01255     result = MPI_Allreduce( &non_root_count, &not_done, 1, MPI_INT, MPI_LOR, comm );
01256     CHECK_MPI( result );
01257     if( not_done )
01258     {
01259         // Get number of types each processor has that root does not
01260         std::vector< int > counts( myPcomm->proc_config().proc_size() );
01261         int two_count = 2 * non_root_count;
01262         result        = MPI_Gather( &two_count, 1, MPI_INT, &counts[0], 1, MPI_INT, 0, comm );
01263         CHECK_MPI( result );
01264 
01265         // Get list of types from each processor
01266         std::vector< int > displs( myPcomm->proc_config().proc_size() + 1 );
01267         (void)VALGRIND_MAKE_VEC_UNDEFINED( displs );
01268         displs[0] = 0;
01269         for( unsigned long i = 1; i <= myPcomm->proc_config().proc_size(); ++i )
01270             displs[i] = displs[i - 1] + counts[i - 1];
01271         int total = displs[myPcomm->proc_config().proc_size()];
01272         typelist alltypes( total / 2 );
01273         (void)VALGRIND_MAKE_VEC_UNDEFINED( alltypes );
01274         (void)VALGRIND_CHECK_MEM_IS_DEFINED( &non_root_types[0], non_root_types.size() * sizeof( int ) );
01275         result = MPI_Gatherv( (void*)&non_root_types[0], 2 * non_root_count, MPI_INT, (int*)&alltypes[0], &counts[0],
01276                               &displs[0], MPI_INT, 0, comm );
01277         CHECK_MPI( result );
01278 
01279         // Merge type lists.
01280         // Prefer O(n) insertions with O(ln n) search time because
01281         // we expect data from a potentially large number of processes,
01282         // but with a small total number of element types.
01283         if( 0 == myPcomm->proc_config().proc_rank() )
01284         {
01285             for( viter = alltypes.begin(); viter != alltypes.end(); ++viter )
01286             {
01287                 typelist::iterator titer = std::lower_bound( my_types.begin(), my_types.end(), *viter );
01288                 if( titer == my_types.end() || *titer != *viter ) my_types.insert( titer, *viter );
01289             }
01290 
01291             dbgOut.print( 2, "Global Element Types:\n" );
01292             for( viter = my_types.begin(); viter != my_types.end(); ++viter )
01293                 dbgOut.printf( 2, "  %s : %d\n", CN::EntityTypeName( (EntityType)viter->first ), viter->second );
01294         }
01295 
01296         // Send total number of types to each processor
01297         total  = my_types.size();
01298         result = MPI_Bcast( &total, 1, MPI_INT, 0, comm );
01299         CHECK_MPI( result );
01300 
01301         // Send list of types to each processor
01302         my_types.resize( total );
01303         result = MPI_Bcast( (void*)&my_types[0], 2 * total, MPI_INT, 0, comm );
01304         CHECK_MPI( result );
01305     }
01306     else
01307     {
01308         // Special case: if root had types but some subset of procs did not
01309         // have those types, but there are no types that the root doesn't
01310         // know about then we still need to update processes that are missing
01311         // types.
01312         my_types.swap( root_types );
01313     }
01314 
01315     // Insert missing types into exportList, with an empty
01316     // range of entities to export.
01317     std::list< ExportSet >::iterator ex_iter = exportList.begin();
01318     for( viter = my_types.begin(); viter != my_types.end(); ++viter )
01319     {
01320         while( ex_iter != exportList.end() && *ex_iter < *viter )
01321             ++ex_iter;
01322 
01323         if( ex_iter == exportList.end() || !( *ex_iter == *viter ) )
01324         {
01325             ExportSet insert;
01326             insert.type       = (EntityType)viter->first;
01327             insert.num_nodes  = viter->second;
01328             insert.first_id   = 0;
01329             insert.offset     = 0;
01330             insert.adj_offset = 0;
01331             ex_iter           = exportList.insert( ex_iter, insert );
01332         }
01333     }
01334 
01335     return MB_SUCCESS;
01336 }
01337 
01338 ErrorCode WriteHDF5Parallel::create_element_tables()
01339 {
01340     struct ElemSetCreator : public DataSetCreator
01341     {
01342         ErrorCode operator()( WriteHDF5* file, long size, const ExportSet* ex, long& start_id ) const
01343         {
01344             return file->create_elem_table( *ex, size, start_id );
01345         }
01346     };
01347 
01348     const int numtypes = exportList.size();
01349     std::vector< ExportSet* > groups( numtypes );
01350     std::vector< long > counts( numtypes ), offsets( numtypes ), max_ents( numtypes ), total_ents( numtypes );
01351     std::vector< wid_t > start_ids( numtypes );
01352 
01353     size_t idx = 0;
01354     std::list< ExportSet >::iterator ex_iter;
01355     for( ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter, ++idx )
01356     {
01357         groups[idx] = &*ex_iter;
01358         counts[idx] = ex_iter->range.size();
01359     }
01360     ErrorCode rval = create_dataset( numtypes, &counts[0], &offsets[0], &max_ents[0], &total_ents[0], ElemSetCreator(),
01361                                      &groups[0], &start_ids[0] );
01362     CHECK_MB( rval );
01363 
01364     for( idx = 0, ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter, ++idx )
01365     {
01366         ex_iter->first_id       = start_ids[idx];
01367         ex_iter->offset         = offsets[idx];
01368         ex_iter->max_num_ents   = max_ents[idx];
01369         ex_iter->total_num_ents = total_ents[idx];
01370         rval                    = assign_ids( ex_iter->range, ex_iter->first_id + ex_iter->offset );
01371         CHECK_MB( rval );
01372     }
01373 
01374     return MB_SUCCESS;
01375 }
01376 
01377 ErrorCode WriteHDF5Parallel::create_adjacency_tables()
01378 {
01379     struct AdjSetCreator : public DataSetCreator
01380     {
01381         ErrorCode operator()( WriteHDF5* file, long size, const ExportSet* ex, long& start_id ) const
01382         {
01383             mhdf_Status status;
01384             hid_t handle = mhdf_createAdjacency( file->file_ptr(), ex->name(), size, &status );
01385             CHECK_HDFN( status );
01386             mhdf_closeData( file->file_ptr(), handle, &status );
01387             CHECK_HDFN( status );
01388             start_id = -1;
01389             return MB_SUCCESS;
01390         }
01391     };
01392 
01393     std::vector< ExportSet* > groups;
01394 #ifdef WRITE_NODE_ADJACENCIES
01395     groups.push_back( &nodeSet );
01396 #endif
01397     for( std::list< ExportSet >::iterator ex_iter = exportList.begin(); ex_iter != exportList.end(); ++ex_iter )
01398         groups.push_back( &*ex_iter );
01399 
01400     ErrorCode rval;
01401     const int numtypes = groups.size();
01402     std::vector< long > counts( numtypes );
01403     std::vector< long > offsets( numtypes );
01404     std::vector< long > max_ents( numtypes );
01405     std::vector< long > totals( numtypes );
01406     for( int i = 0; i < numtypes; ++i )
01407     {
01408         wid_t count;
01409         rval = count_adjacencies( groups[i]->range, count );
01410         CHECK_MB( rval );
01411         counts[i] = count;
01412     }
01413 
01414     rval = create_dataset( numtypes, &counts[0], &offsets[0], &max_ents[0], &totals[0], AdjSetCreator(), &groups[0] );
01415     CHECK_MB( rval );
01416 
01417     // Cppcheck warning (false positive): variable groups is assigned a value that is never used
01418     for( int i = 0; i < numtypes; ++i )
01419     {
01420         groups[i]->max_num_adjs = max_ents[i];
01421         groups[i]->adj_offset   = offsets[i];
01422     }
01423     return MB_SUCCESS;
01424 }
01425 
01426 const unsigned SSVB = 3;
01427 
01428 void WriteHDF5Parallel::print_set_sharing_data( const Range& range, const char* label, Tag idt )
01429 {
01430     dbgOut.printf( SSVB, "set\tid\towner\t%-*s\tfid\tshared\n", (int)( sizeof( EntityHandle ) * 2 ), "handle" );
01431     for( Range::iterator it = range.begin(); it != range.end(); ++it )
01432     {
01433         int id;
01434         iFace->tag_get_data( idt, &*it, 1, &id );
01435         EntityHandle handle = 0;
01436         unsigned owner      = 0;
01437         wid_t file_id       = 0;
01438         myPcomm->get_entityset_owner( *it, owner, &handle );
01439         if( !idMap.find( *it, file_id ) ) file_id = 0;
01440         dbgOut.printf( SSVB, "%s\t%d\t%u\t%lx\t%lu\t", label, id, owner, (unsigned long)handle,
01441                        (unsigned long)file_id );
01442         std::vector< unsigned > procs;
01443         myPcomm->get_entityset_procs( *it, procs );
01444         if( procs.empty() )
01445             dbgOut.print( SSVB, "<none>\n" );
01446         else
01447         {
01448             for( unsigned i = 0; i < procs.size() - 1; ++i )
01449                 dbgOut.printf( SSVB, "%u,", procs[i] );
01450             dbgOut.printf( SSVB, "%u\n", procs.back() );
01451         }
01452     }
01453 }
01454 
01455 void WriteHDF5Parallel::print_shared_sets()
01456 {
01457     const char* tag_names[][2] = { { MATERIAL_SET_TAG_NAME, "block" },
01458                                    { DIRICHLET_SET_TAG_NAME, "nodeset" },
01459                                    { NEUMANN_SET_TAG_NAME, "sideset" },
01460                                    { 0, 0 } };
01461 
01462     for( int i = 0; tag_names[i][0]; ++i )
01463     {
01464         Tag tag;
01465         if( MB_SUCCESS != iFace->tag_get_handle( tag_names[i][0], 1, MB_TYPE_INTEGER, tag ) ) continue;
01466 
01467         Range tagged;
01468         iFace->get_entities_by_type_and_tag( 0, MBENTITYSET, &tag, 0, 1, tagged );
01469         print_set_sharing_data( tagged, tag_names[i][1], tag );
01470     }
01471 
01472     Tag geom, id;
01473     if( MB_SUCCESS != iFace->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geom ) ) return;
01474     id = iFace->globalId_tag();
01475 
01476     const char* geom_names[] = { "vertex", "curve", "surface", "volume" };
01477     for( int d = 0; d <= 3; ++d )
01478     {
01479         Range tagged;
01480         const void* vals[] = { &d };
01481         iFace->get_entities_by_type_and_tag( 0, MBENTITYSET, &geom, vals, 1, tagged );
01482         print_set_sharing_data( tagged, geom_names[d], id );
01483     }
01484 }
01485 
01486 ErrorCode WriteHDF5Parallel::communicate_shared_set_ids( const Range& owned, const Range& remote )
01487 {
01488     ErrorCode rval;
01489     int mperr;
01490     const int TAG = 0xD0E;
01491     // const unsigned rank = myPcomm->proc_config().proc_rank();
01492     const MPI_Comm comm = myPcomm->proc_config().proc_comm();
01493 
01494     dbgOut.tprint( 1, "COMMUNICATING SHARED SET IDS\n" );
01495     dbgOut.print( 6, "Owned, shared sets: ", owned );
01496 
01497     // Post receive buffers for all procs for which we share sets
01498 
01499     std::vector< unsigned > procs;
01500     rval = myPcomm->get_entityset_owners( procs );
01501     CHECK_MB( rval );
01502     std::vector< unsigned >::iterator it = std::find( procs.begin(), procs.end(), myPcomm->proc_config().proc_rank() );
01503     if( it != procs.end() ) procs.erase( it );
01504 
01505     std::vector< MPI_Request > recv_req( procs.size(), MPI_REQUEST_NULL );
01506     std::vector< std::vector< unsigned long > > recv_buf( procs.size() );
01507 
01508     size_t recv_count = 0;
01509     for( size_t i = 0; i < procs.size(); ++i )
01510     {
01511         Range tmp;
01512         rval = myPcomm->get_owned_sets( procs[i], tmp );
01513         CHECK_MB( rval );
01514         size_t count =
01515             intersect( tmp, remote ).size();  // Necessary because we might not be writing all of the database
01516         if( count )
01517         {
01518             dbgOut.printf( 6, "Sets owned by proc %u (remote handles): ", procs[i] );
01519             if( dbgOut.get_verbosity() >= 6 )
01520             {
01521                 Range remote_handles;
01522                 tmp = intersect( tmp, remote );
01523                 for( Range::iterator j = tmp.begin(); j != tmp.end(); ++j )
01524                 {
01525                     unsigned r;
01526                     EntityHandle h;
01527                     myPcomm->get_entityset_owner( *j, r, &h );
01528                     assert( r == procs[i] );
01529                     remote_handles.insert( h );
01530                 }
01531                 dbgOut.print( 6, remote_handles );
01532             }
01533             recv_count++;
01534             recv_buf[i].resize( 2 * count + 1 );
01535             dbgOut.printf( 5, "Posting receive buffer of size %lu for proc %u (%lu of %lu owned sets)\n",
01536                            (unsigned long)recv_buf[i].size(), procs[i], count, tmp.size() );
01537             mperr =
01538                 MPI_Irecv( &recv_buf[i][0], recv_buf[i].size(), MPI_UNSIGNED_LONG, procs[i], TAG, comm, &recv_req[i] );
01539             CHECK_MPI( mperr );
01540         }
01541     }
01542 
01543     // Send set ids to all procs with which we share them
01544 
01545     // First build per-process lists of sets for which we need to send data
01546     std::map< unsigned, Range > send_sets;
01547     std::vector< unsigned > set_procs;
01548     for( Range::reverse_iterator i = owned.rbegin(); i != owned.rend(); ++i )
01549     {
01550         set_procs.clear();
01551         rval = myPcomm->get_entityset_procs( *i, set_procs );
01552         CHECK_MB( rval );
01553         for( size_t j = 0; j < set_procs.size(); ++j )
01554             if( set_procs[j] != myPcomm->proc_config().proc_rank() ) send_sets[set_procs[j]].insert( *i );
01555     }
01556     assert( send_sets.find( myPcomm->proc_config().proc_rank() ) == send_sets.end() );
01557 
01558     // Now send the data
01559     std::vector< std::vector< unsigned long > > send_buf( send_sets.size() );
01560     std::vector< MPI_Request > send_req( send_sets.size() );
01561     std::map< unsigned, Range >::iterator si = send_sets.begin();
01562     for( size_t i = 0; si != send_sets.end(); ++si, ++i )
01563     {
01564         dbgOut.printf( 6, "Sending data for shared sets to proc %u: ", si->first );
01565         dbgOut.print( 6, si->second );
01566 
01567         send_buf[i].reserve( 2 * si->second.size() + 1 );
01568         send_buf[i].push_back( si->second.size() );
01569         for( Range::iterator j = si->second.begin(); j != si->second.end(); ++j )
01570         {
01571             send_buf[i].push_back( *j );
01572             send_buf[i].push_back( idMap.find( *j ) );
01573         }
01574         dbgOut.printf( 5, "Sending buffer of size %lu to proc %u (%lu of %lu owned sets)\n",
01575                        (unsigned long)send_buf[i].size(), si->first, si->second.size(), owned.size() );
01576         mperr = MPI_Isend( &send_buf[i][0], send_buf[i].size(), MPI_UNSIGNED_LONG, si->first, TAG, comm, &send_req[i] );
01577     }
01578 
01579     // Process received data
01580     MPI_Status status;
01581     int idx;
01582     while( recv_count-- )
01583     {
01584         mperr = MPI_Waitany( recv_req.size(), &recv_req[0], &idx, &status );
01585         CHECK_MPI( mperr );
01586 
01587         assert( (unsigned)status.MPI_SOURCE == procs[idx] );
01588         assert( 2 * recv_buf[idx].front() + 1 == recv_buf[idx].size() );
01589         const size_t n = std::min< size_t >( recv_buf[idx].front(), ( recv_buf[idx].size() - 1 ) / 2 );
01590         dbgOut.printf( 5, "Received buffer of size %lu from proc %d\n", (unsigned long)( 2 * n + 1 ),
01591                        (int)status.MPI_SOURCE );
01592 
01593         for( size_t i = 0; i < n; ++i )
01594         {
01595             EntityHandle handle = 0;
01596             rval                = myPcomm->get_entityset_local_handle( procs[idx], recv_buf[idx][2 * i + 1], handle );
01597             CHECK_MB( rval );
01598             assert( handle != 0 );
01599             if( !idMap.insert( handle, recv_buf[idx][2 * i + 2], 1 ).second )
01600                 error( MB_FAILURE );  // Conflicting IDs??????
01601         }
01602 
01603         recv_req[idx] = MPI_REQUEST_NULL;
01604     }
01605     assert( MPI_SUCCESS == MPI_Waitany( recv_req.size(), &recv_req[0], &idx, &status ) &&
01606             MPI_UNDEFINED == idx );  // Check that we got them all
01607 
01608     // Wait for all sends to complete before we release send
01609     // buffers (implicitly releases when we return from this function)
01610 
01611     std::vector< MPI_Status > stats( send_req.size() );
01612     mperr = MPI_Waitall( send_req.size(), &send_req[0], &stats[0] );
01613     CHECK_MPI( mperr );
01614 
01615     if( dbgOut.get_verbosity() >= SSVB ) print_shared_sets();
01616 
01617     return MB_SUCCESS;
01618 }
01619 
01620 // void get_global_ids(Interface* iFace, const unsigned long* ptr,
01621 //                    size_t len, unsigned flags,
01622 //                    std::vector<int>& ids)
01623 //{
01624 //  Tag idtag;
01625 //  iFace->tag_get_handle(GLOBAL_ID_TAG_NAME, 1, MB_TYPE_INTEGER, idtag);
01626 //  for (size_t i = 0; i < len; ++i) {
01627 //    if (flags & MESHSET_ORDERED) {
01628 //      int tmp;
01629 //      iFace->tag_get_data(idtag, ptr + i, 1, &tmp);
01630 //      ids.push_back(tmp);
01631 //      continue;
01632 //    }
01633 //
01634 //    EntityHandle s = ptr[i];
01635 //    EntityHandle e = ptr[++i];
01636 //    for (; s <= e; ++s) {
01637 //      int tmp;
01638 //      iFace->tag_get_data(idtag, &s, 1, &tmp);
01639 //      ids.push_back(tmp);
01640 //    }
01641 //  }
01642 //}
01643 
01644 ErrorCode WriteHDF5Parallel::pack_set( Range::const_iterator it, unsigned long* buffer, size_t buffer_size )
01645 {
01646     ErrorCode rval;
01647     const EntityHandle* ptr;
01648     int len;
01649     unsigned char flags;
01650     std::vector< wid_t > tmp;
01651     size_t newlen;
01652 
01653     // Buffer must always contain at least flags and desired sizes
01654     assert( buffer_size >= 4 );
01655     buffer_size -= 4;
01656 
01657     Range::const_iterator nd = it;
01658     ++nd;
01659     rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::CONTENTS, &len, &flags );
01660     CHECK_MB( rval );
01661 
01662     // Tag mattag;
01663     // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
01664     // int block;
01665     // if (MB_SUCCESS != iFace->tag_get_data(mattag, &*it, 1, &block))
01666     //  block = 0;
01667     //
01668     // if (block) {
01669     //  std::vector<int> ids;
01670     //  get_global_ids(iFace, ptr, len, flags, ids);
01671     //}
01672 
01673     if( len && !( flags & MESHSET_ORDERED ) )
01674     {
01675         tmp.clear();
01676         bool blocked = false;
01677         assert( ( 0 == len % 2 ) );
01678         rval = range_to_blocked_list( ptr, len / 2, tmp, blocked );
01679         CHECK_MB( rval );
01680         if( blocked ) flags |= mhdf_SET_RANGE_BIT;
01681     }
01682     else
01683     {
01684         tmp.resize( len );
01685         rval = vector_to_id_list( ptr, len, &tmp[0], newlen, true );
01686         CHECK_MB( rval );
01687         tmp.resize( newlen );
01688     }
01689 
01690     buffer[0] = flags;
01691     buffer[1] = tmp.size();
01692     if( tmp.size() <= buffer_size ) std::copy( tmp.begin(), tmp.end(), buffer + 4 );
01693 
01694     rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::CHILDREN, &len );
01695     CHECK_MB( rval );
01696     tmp.resize( len );
01697     rval = vector_to_id_list( ptr, len, &tmp[0], newlen, true );
01698     tmp.resize( newlen );
01699     buffer[2] = tmp.size();
01700     if( tmp.size() <= buffer_size - buffer[1] ) std::copy( tmp.begin(), tmp.end(), buffer + 4 + buffer[1] );
01701 
01702     rval = writeUtil->get_entity_list_pointers( it, nd, &ptr, WriteUtilIface::PARENTS, &len );
01703     CHECK_MB( rval );
01704     tmp.resize( len );
01705     rval = vector_to_id_list( ptr, len, &tmp[0], newlen, true );
01706     tmp.resize( newlen );
01707     buffer[3] = tmp.size();
01708     if( tmp.size() <= buffer_size - buffer[1] - buffer[2] )
01709         std::copy( tmp.begin(), tmp.end(), buffer + 4 + buffer[1] + buffer[2] );
01710 
01711     return MB_SUCCESS;
01712 }
01713 
01714 template < typename TYPE >
01715 static void convert_to_ranged_ids( const TYPE* buffer, size_t len, std::vector< WriteHDF5::wid_t >& result )
01716 {
01717     if( !len )
01718     {
01719         result.clear();
01720         return;
01721     }
01722 
01723     result.resize( len * 2 );
01724     Range tmp;
01725     for( size_t i = 0; i < len; i++ )
01726         tmp.insert( (EntityHandle)buffer[i] );
01727     result.resize( tmp.psize() * 2 );
01728     int j = 0;
01729     for( Range::const_pair_iterator pit = tmp.const_pair_begin(); pit != tmp.const_pair_end(); ++pit, j++ )
01730     {
01731         result[2 * j]     = pit->first;
01732         result[2 * j + 1] = pit->second - pit->first + 1;
01733     }
01734 }
01735 
01736 static void merge_ranged_ids( const unsigned long* range_list, size_t len, std::vector< WriteHDF5::wid_t >& result )
01737 {
01738     typedef WriteHDF5::wid_t wid_t;
01739     assert( 0 == len % 2 );
01740     assert( 0 == result.size() % 2 );
01741     STATIC_ASSERT( sizeof( std::pair< wid_t, wid_t > ) == 2 * sizeof( wid_t ) );
01742 
01743     result.insert( result.end(), range_list, range_list + len );
01744     size_t plen = result.size() / 2;
01745     Range tmp;
01746     for( size_t i = 0; i < plen; i++ )
01747     {
01748         EntityHandle starth = (EntityHandle)result[2 * i];
01749         EntityHandle endh   = (EntityHandle)result[2 * i] + (wid_t)result[2 * i + 1] - 1;  // id + count - 1
01750         tmp.insert( starth, endh );
01751     }
01752     // Now convert back to std::vector<WriteHDF5::wid_t>, compressed range format
01753     result.resize( tmp.psize() * 2 );
01754     int j = 0;
01755     for( Range::const_pair_iterator pit = tmp.const_pair_begin(); pit != tmp.const_pair_end(); ++pit, j++ )
01756     {
01757         result[2 * j]     = pit->first;
01758         result[2 * j + 1] = pit->second - pit->first + 1;
01759     }
01760 }
01761 
01762 static void merge_vector_ids( const unsigned long* list, size_t len, std::vector< WriteHDF5::wid_t >& result )
01763 {
01764     result.insert( result.end(), list, list + len );
01765 }
01766 
01767 ErrorCode WriteHDF5Parallel::unpack_set( EntityHandle set, const unsigned long* buffer, size_t buffer_size )
01768 {
01769     // Use local variables for readability
01770     assert( buffer_size >= 4 );
01771     assert( buffer[1] + buffer[2] + buffer[3] <= buffer_size );
01772     const unsigned long flags      = buffer[0];
01773     unsigned long num_content      = buffer[1];
01774     const unsigned long num_child  = buffer[2];
01775     const unsigned long num_parent = buffer[3];
01776     const unsigned long* contents  = buffer + 4;
01777     const unsigned long* children  = contents + num_content;
01778     const unsigned long* parents   = children + num_child;
01779 
01780     SpecialSetData* data = find_set_data( set );
01781     assert( NULL != data );
01782     if( NULL == data ) return MB_FAILURE;
01783 
01784     // Tag mattag;
01785     // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
01786     // int block;
01787     // if (MB_SUCCESS != iFace->tag_get_data(mattag, &set, 1, &block))
01788     //  block = 0;
01789 
01790     // If either the current data or the new data is in ranged format,
01791     // then change the other to ranged format if it isn't already
01792     // in both cases when they differ, the data will end up "compressed range"
01793     std::vector< wid_t > tmp;
01794     if( ( flags & mhdf_SET_RANGE_BIT ) != ( data->setFlags & mhdf_SET_RANGE_BIT ) )
01795     {
01796         if( flags & mhdf_SET_RANGE_BIT )
01797         {
01798             tmp = data->contentIds;
01799             convert_to_ranged_ids( &tmp[0], tmp.size(), data->contentIds );
01800             data->setFlags |= mhdf_SET_RANGE_BIT;
01801         }
01802         else
01803         {
01804             tmp.clear();
01805             convert_to_ranged_ids( contents, num_content, tmp );
01806             num_content = tmp.size();
01807             if( sizeof( wid_t ) < sizeof( long ) )
01808             {
01809                 size_t old_size = tmp.size();
01810                 tmp.resize( sizeof( long ) * old_size / sizeof( wid_t ) );
01811                 unsigned long* array = reinterpret_cast< unsigned long* >( &tmp[0] );
01812                 for( long i = ( (long)old_size ) - 1; i >= 0; --i )
01813                     array[i] = tmp[i];
01814                 contents = array;
01815             }
01816             else if( sizeof( wid_t ) > sizeof( long ) )
01817             {
01818                 unsigned long* array = reinterpret_cast< unsigned long* >( &tmp[0] );
01819                 std::copy( tmp.begin(), tmp.end(), array );
01820             }
01821             contents = reinterpret_cast< unsigned long* >( &tmp[0] );
01822         }
01823     }
01824 
01825     if( data->setFlags & mhdf_SET_RANGE_BIT )
01826         merge_ranged_ids( contents, num_content, data->contentIds );
01827     else
01828         merge_vector_ids( contents, num_content, data->contentIds );
01829 
01830     merge_vector_ids( children, num_child, data->childIds );
01831     merge_vector_ids( parents, num_parent, data->parentIds );
01832     return MB_SUCCESS;
01833 }
01834 
01835 ErrorCode WriteHDF5Parallel::communicate_shared_set_data( const Range& owned, const Range& remote )
01836 {
01837     ErrorCode rval;
01838     int mperr;
01839     const unsigned rank = myPcomm->proc_config().proc_rank();
01840     const MPI_Comm comm = myPcomm->proc_config().proc_comm();
01841 
01842     dbgOut.tprintf( 1, "COMMUNICATING SHARED SET DATA (%lu owned & %lu remote)\n", (unsigned long)owned.size(),
01843                     (unsigned long)remote.size() );
01844 
01845     // Calculate the total number of messages to be in transit (send and receive)
01846     size_t nummess = 0;
01847     std::vector< unsigned > procs;
01848     ;
01849     Range shared( owned );
01850     shared.merge( remote );
01851     for( Range::iterator i = shared.begin(); i != shared.end(); ++i )
01852     {
01853         procs.clear();
01854         rval = myPcomm->get_entityset_procs( *i, procs );
01855         CHECK_MB( rval );
01856         nummess += procs.size();
01857     }
01858 
01859     // Choose a receive buffer size. We need 4*sizeof(long) minimum,
01860     // but that is almost useless so use 16*sizeof(long) as the minimum
01861     // instead. Choose an upper limit such that we don't exceed 32 MB
01862     // of allocated memory (unless we absolutely must to meet the minimum.)
01863     // Also, don't initially choose buffers larger than 128*sizeof(long).
01864     const size_t MAX_BUFFER_MEM = 32 * 1024 * 1024 / sizeof( long );
01865     // const size_t INIT_BUFFER_SIZE = 128;
01866     const size_t INIT_BUFFER_SIZE = 1024;
01867     const size_t MIN_BUFFER_SIZE  = 16;
01868     size_t init_buff_size         = INIT_BUFFER_SIZE;
01869     if( init_buff_size * nummess > MAX_BUFFER_MEM ) init_buff_size = MAX_BUFFER_MEM / nummess;
01870     if( init_buff_size < MIN_BUFFER_SIZE ) init_buff_size = MIN_BUFFER_SIZE;
01871 
01872     dbgOut.printf( 2, "Using buffer size of %lu for an expected message count of %lu\n", (unsigned long)init_buff_size,
01873                    (unsigned long)nummess );
01874 
01875     // Count number of recvs
01876     size_t numrecv = 0;
01877     for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
01878     {
01879         procs.clear();
01880         rval = myPcomm->get_entityset_procs( *i, procs );
01881         CHECK_MB( rval );
01882         numrecv += procs.size();
01883         if( std::find( procs.begin(), procs.end(), rank ) != procs.end() ) --numrecv;
01884     }
01885 
01886     // Post receive buffers for all owned sets for all sharing procs
01887     std::vector< MPI_Request > recv_req( numrecv, MPI_REQUEST_NULL );
01888     std::vector< MPI_Request > lrecv_req( numrecv, MPI_REQUEST_NULL );
01889 
01890     std::vector< std::vector< unsigned long > > recv_buf( numrecv, std::vector< unsigned long >( init_buff_size ) );
01891     int idx = 0;
01892     for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
01893     {
01894         procs.clear();
01895         rval = myPcomm->get_entityset_procs( *i, procs );
01896         CHECK_MB( rval );
01897         for( size_t j = 0; j < procs.size(); ++j )
01898         {
01899             if( procs[j] == rank ) continue;
01900             int tag = ID_FROM_HANDLE( *i );
01901             if( *i != CREATE_HANDLE( MBENTITYSET, tag ) )
01902             {
01903 #ifndef NDEBUG
01904                 abort();
01905 #endif
01906                 CHECK_MB( MB_FAILURE );
01907             }
01908             dbgOut.printf( 5, "Posting buffer to receive set %d from proc %u\n", tag, procs[j] );
01909             mperr =
01910                 MPI_Irecv( &recv_buf[idx][0], init_buff_size, MPI_UNSIGNED_LONG, procs[j], tag, comm, &recv_req[idx] );
01911             CHECK_MPI( mperr );
01912             ++idx;
01913         }
01914     }
01915     assert( (size_t)idx == numrecv );
01916 
01917     // Now send set data for all remote sets that I know about
01918     std::vector< MPI_Request > send_req( remote.size() );
01919     std::vector< std::vector< unsigned long > > send_buf( remote.size() );
01920     idx = 0;
01921     for( Range::iterator i = remote.begin(); i != remote.end(); ++i, ++idx )
01922     {
01923         send_buf[idx].resize( init_buff_size );
01924         rval = pack_set( i, &send_buf[idx][0], init_buff_size );
01925         CHECK_MB( rval );
01926         EntityHandle remote_handle;
01927         unsigned owner;
01928         rval = myPcomm->get_entityset_owner( *i, owner, &remote_handle );
01929         CHECK_MB( rval );
01930 
01931         int tag = ID_FROM_HANDLE( remote_handle );
01932         assert( remote_handle == CREATE_HANDLE( MBENTITYSET, tag ) );
01933         dbgOut.printf( 5, "Sending %lu values for set %d to proc %u\n",
01934                        send_buf[idx][1] + send_buf[idx][2] + send_buf[idx][3] + 4, tag, owner );
01935         mperr = MPI_Isend( &send_buf[idx][0], init_buff_size, MPI_UNSIGNED_LONG, owner, tag, comm, &send_req[idx] );
01936         CHECK_MPI( mperr );
01937     }
01938 
01939     // Tag mattag;
01940     // iFace->tag_get_handle(MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mattag);
01941 
01942     // Now initialize local data for managing contents of owned, shared sets
01943     assert( specialSets.empty() );
01944     specialSets.clear();
01945     specialSets.reserve( owned.size() );
01946     for( Range::iterator i = owned.begin(); i != owned.end(); ++i )
01947     {
01948         // int block;
01949         // if (MB_SUCCESS != iFace->tag_get_data(mattag, &*i, 1, &block))
01950         //  block = 0;
01951         // std::vector<int> ids;
01952 
01953         SpecialSetData data;
01954         data.setHandle = *i;
01955         rval           = iFace->get_meshset_options( *i, data.setFlags );
01956         CHECK_MB( rval );
01957         specialSets.push_back( data );
01958         std::vector< EntityHandle > list;
01959         if( data.setFlags & MESHSET_ORDERED )
01960         {
01961             list.clear();
01962             rval = iFace->get_entities_by_handle( *i, list );
01963             CHECK_MB( rval );
01964             rval = vector_to_id_list( list, specialSets.back().contentIds, true );
01965             CHECK_MB( rval );
01966             // if (block)
01967             //  get_global_ids(iFace, &list[0], list.size(), MESHSET_ORDERED, ids);
01968         }
01969         else
01970         {
01971             Range range;
01972             rval = iFace->get_entities_by_handle( *i, range );
01973             CHECK_MB( rval );
01974             bool ranged;
01975             rval = range_to_blocked_list( range, specialSets.back().contentIds, ranged );
01976             if( ranged ) specialSets.back().setFlags |= mhdf_SET_RANGE_BIT;
01977             // if (block) {
01978             //  std::vector<EntityHandle> tmp;
01979             //  for (Range::const_pair_iterator pi = range.const_pair_begin(); pi !=
01980             //  range.const_pair_end(); ++pi) {
01981             //    tmp.push_back(pi->first);
01982             //    tmp.push_back(pi->second);
01983             //  }
01984             //  get_global_ids(iFace, &tmp[0], tmp.size(), ranged ? 0 : MESHSET_ORDERED, ids);
01985             //}
01986         }
01987 
01988         list.clear();
01989         rval = iFace->get_parent_meshsets( *i, list );
01990         CHECK_MB( rval );
01991         rval = vector_to_id_list( list, specialSets.back().parentIds, true );
01992         CHECK_MB( rval );
01993         rval = iFace->get_child_meshsets( *i, list );
01994         CHECK_MB( rval );
01995         rval = vector_to_id_list( list, specialSets.back().childIds, true );
01996         CHECK_MB( rval );
01997     }
01998 
01999     // Process received buffers, repost larger buffers where necessary
02000     size_t remaining = numrecv;
02001     numrecv          = 0;
02002     while( remaining-- )
02003     {
02004         std::vector< unsigned long > dead;
02005         MPI_Status status;
02006         mperr = MPI_Waitany( recv_req.size(), &recv_req[0], &idx, &status );
02007         CHECK_MPI( mperr );
02008         EntityHandle handle                = CREATE_HANDLE( MBENTITYSET, status.MPI_TAG );
02009         std::vector< unsigned long >& buff = recv_buf[idx];
02010         size_t size                        = buff[1] + buff[2] + buff[3] + 4;
02011         dbgOut.printf( 5, "Received %lu values for set %d from proc %d\n", (unsigned long)size, status.MPI_TAG,
02012                        status.MPI_SOURCE );
02013         if( size <= init_buff_size )
02014         {
02015             rval = unpack_set( handle, &buff[0], init_buff_size );
02016             CHECK_MB( rval );
02017             dead.swap( buff );  // Release memory
02018         }
02019         else
02020         {
02021             // Data was too big for init_buff_size
02022             // repost with larger buffer
02023             buff.resize( size );
02024             dbgOut.printf( 5, "Re-Posting buffer to receive set %d from proc %d with size %lu\n", status.MPI_TAG,
02025                            status.MPI_SOURCE, (unsigned long)size );
02026             mperr = MPI_Irecv( &buff[0], size, MPI_UNSIGNED_LONG, status.MPI_SOURCE, status.MPI_TAG, comm,
02027                                &lrecv_req[idx] );
02028             CHECK_MPI( mperr );
02029             ++numrecv;
02030         }
02031         recv_req[idx] = MPI_REQUEST_NULL;
02032     }
02033 
02034     // Wait for sends to complete
02035     MPI_Waitall( send_req.size(), &send_req[0], MPI_STATUSES_IGNORE );
02036 
02037     // Re-send sets that didn't fit initial buffer size
02038     idx = 0;
02039     for( Range::iterator i = remote.begin(); i != remote.end(); ++i, ++idx )
02040     {
02041         std::vector< unsigned long >& buff = send_buf[idx];
02042         size_t size                        = buff[1] + buff[2] + buff[3] + 4;
02043         if( size <= init_buff_size ) continue;
02044 
02045         buff.resize( size );
02046         rval = pack_set( i, &buff[0], size );
02047         CHECK_MB( rval );
02048         EntityHandle remote_handle;
02049         unsigned owner;
02050         rval = myPcomm->get_entityset_owner( *i, owner, &remote_handle );
02051         CHECK_MB( rval );
02052 
02053         int tag = ID_FROM_HANDLE( remote_handle );
02054         assert( remote_handle == CREATE_HANDLE( MBENTITYSET, tag ) );
02055         dbgOut.printf( 5, "Sending %lu values for set %d to proc %u\n", (unsigned long)size, tag, owner );
02056         mperr = MPI_Isend( &buff[0], size, MPI_UNSIGNED_LONG, owner, tag, comm, &send_req[idx] );
02057         CHECK_MPI( mperr );
02058     }
02059 
02060     // Process received buffers
02061     remaining = numrecv;
02062     while( remaining-- )
02063     {
02064         std::vector< unsigned long > dead;
02065         MPI_Status status;
02066         mperr = MPI_Waitany( lrecv_req.size(), &lrecv_req[0], &idx, &status );
02067         CHECK_MPI( mperr );
02068         EntityHandle handle                = CREATE_HANDLE( MBENTITYSET, status.MPI_TAG );
02069         std::vector< unsigned long >& buff = recv_buf[idx];
02070         dbgOut.printf( 5, "Received %lu values for set %d from proc %d\n", 4 + buff[1] + buff[2] + buff[3],
02071                        status.MPI_TAG, status.MPI_SOURCE );
02072         rval = unpack_set( handle, &buff[0], buff.size() );
02073         CHECK_MB( rval );
02074         dead.swap( buff );  // Release memory
02075 
02076         lrecv_req[idx] = MPI_REQUEST_NULL;
02077     }
02078 
02079     // Wait for sends to complete
02080     MPI_Waitall( send_req.size(), &send_req[0], MPI_STATUSES_IGNORE );
02081 
02082     return MB_SUCCESS;
02083 }
02084 
02085 ErrorCode WriteHDF5Parallel::create_meshset_tables( double* times )
02086 {
02087     Range::const_iterator riter;
02088     const unsigned rank = myPcomm->proc_config().proc_rank();
02089 
02090     START_SERIAL;
02091     print_type_sets( iFace, &dbgOut, setSet.range );
02092     END_SERIAL;
02093     CpuTimer timer;
02094 
02095     // Remove remote sets from setSets
02096     Range shared, owned, remote;
02097     ErrorCode rval = myPcomm->get_shared_sets( shared );
02098     CHECK_MB( rval );
02099     shared = intersect( shared, setSet.range );
02100     rval   = myPcomm->get_owned_sets( rank, owned );
02101     CHECK_MB( rval );
02102     owned        = intersect( owned, setSet.range );
02103     remote       = subtract( shared, owned );
02104     setSet.range = subtract( setSet.range, remote );
02105 
02106     // Create set meta table
02107     struct SetDescCreator : public DataSetCreator
02108     {
02109         ErrorCode operator()( WriteHDF5* writer, long size, const ExportSet*, long& start_id ) const
02110         {
02111             return writer->create_set_meta( size, start_id );
02112         }
02113     };
02114     long count = setSet.range.size();
02115     rval = create_dataset( 1, &count, &setSet.offset, &setSet.max_num_ents, &setSet.total_num_ents, SetDescCreator(),
02116                            NULL, &setSet.first_id );
02117     CHECK_MB( rval );
02118     writeSets = setSet.max_num_ents > 0;
02119 
02120     rval = assign_ids( setSet.range, setSet.first_id + setSet.offset );
02121     CHECK_MB( rval );
02122     if( times ) times[SET_OFFSET_TIME] = timer.time_elapsed();
02123 
02124     // Exchange file IDS for sets between all procs
02125     rval = communicate_shared_set_ids( owned, remote );
02126     CHECK_MB( rval );
02127     if( times ) times[SHARED_SET_IDS] = timer.time_elapsed();
02128 
02129     // Communicate remote set contents, children, etc.
02130     rval = communicate_shared_set_data( owned, remote );
02131     CHECK_MB( rval );
02132     if( times ) times[SHARED_SET_CONTENTS] = timer.time_elapsed();
02133 
02134     // Communicate counts for owned sets
02135     long data_counts[3];  // { #contents, #children, #parents }
02136     rval = count_set_size( setSet.range, data_counts[0], data_counts[1], data_counts[2] );
02137     CHECK_MB( rval );
02138     if( times ) times[SET_OFFSET_TIME] += timer.time_elapsed();
02139 
02140     long offsets[3], max_counts[3], totals[3];
02141     rval = create_dataset( 3, data_counts, offsets, max_counts, totals );
02142     CHECK_MB( rval );
02143 
02144     // Create the datasets
02145     if( 0 == myPcomm->proc_config().proc_rank() )
02146     {
02147         rval = create_set_tables( totals[0], totals[1], totals[2] );
02148         CHECK_MB( rval );
02149     }
02150 
02151     // Store communicated global data
02152     setContentsOffset = offsets[0];
02153     setChildrenOffset = offsets[1];
02154     setParentsOffset  = offsets[2];
02155     maxNumSetContents = max_counts[0];
02156     maxNumSetChildren = max_counts[1];
02157     maxNumSetParents  = max_counts[2];
02158     writeSetContents  = totals[0] > 0;
02159     writeSetChildren  = totals[1] > 0;
02160     writeSetParents   = totals[2] > 0;
02161 
02162     dbgOut.printf( 2, "set contents: %ld local, %ld global, offset = %ld\n", data_counts[0], totals[0], offsets[0] );
02163     dbgOut.printf( 2, "set children: %ld local, %ld global, offset = %ld\n", data_counts[1], totals[1], offsets[1] );
02164     dbgOut.printf( 2, "set parents: %ld local, %ld global, offset = %ld\n", data_counts[2], totals[2], offsets[2] );
02165 
02166     return MB_SUCCESS;
02167 }
02168 
02169 void WriteHDF5Parallel::remove_remote_entities( EntityHandle relative, Range& range )
02170 {
02171     Range result;
02172     result.merge( intersect( range, nodeSet.range ) );
02173     result.merge( intersect( range, setSet.range ) );
02174     for( std::list< ExportSet >::iterator eiter = exportList.begin(); eiter != exportList.end(); ++eiter )
02175         result.merge( intersect( range, eiter->range ) );
02176 
02177     // result.merge(intersect(range, myParallelSets));
02178     Range sets;
02179     int junk;
02180     sets.merge( Range::lower_bound( range.begin(), range.end(), CREATE_HANDLE( MBENTITYSET, 0, junk ) ), range.end() );
02181     remove_remote_sets( relative, sets );
02182     result.merge( sets );
02183     range.swap( result );
02184 }
02185 
02186 void WriteHDF5Parallel::remove_remote_sets( EntityHandle /* relative */, Range& range )
02187 {
02188     Range result( intersect( range, setSet.range ) );
02189     // Store the non-intersecting entities separately if needed
02190     // Range remaining(subtract(range, result));
02191     range.swap( result );
02192 }
02193 
02194 void WriteHDF5Parallel::remove_remote_entities( EntityHandle relative, std::vector< EntityHandle >& vect )
02195 {
02196     Range intrsct;
02197     for( std::vector< EntityHandle >::const_iterator iter = vect.begin(); iter != vect.end(); ++iter )
02198         intrsct.insert( *iter );
02199     remove_remote_entities( relative, intrsct );
02200 
02201     unsigned int read, write;
02202     for( read = write = 0; read < vect.size(); ++read )
02203     {
02204         if( intrsct.find( vect[read] ) != intrsct.end() )
02205         {
02206             if( read != write ) vect[write] = vect[read];
02207             ++write;
02208         }
02209     }
02210     if( write != vect.size() ) vect.resize( write );
02211 }
02212 
02213 void WriteHDF5Parallel::remove_remote_sets( EntityHandle relative, std::vector< EntityHandle >& vect )
02214 {
02215     Range intrsct;
02216     for( std::vector< EntityHandle >::const_iterator iter = vect.begin(); iter != vect.end(); ++iter )
02217         intrsct.insert( *iter );
02218     remove_remote_sets( relative, intrsct );
02219 
02220     unsigned int read, write;
02221     for( read = write = 0; read < vect.size(); ++read )
02222     {
02223         if( intrsct.find( vect[read] ) != intrsct.end() )
02224         {
02225             if( read != write ) vect[write] = vect[read];
02226             ++write;
02227         }
02228     }
02229     if( write != vect.size() ) vect.resize( write );
02230 }
02231 
02232 ErrorCode WriteHDF5Parallel::exchange_file_ids( const Range& nonlocal )
02233 {
02234     ErrorCode rval;
02235 
02236     // For each entity owned on the interface, write its file id to
02237     // a tag. The sets of entities to be written should already contain
02238     // only owned entities so by intersecting with them we not only
02239     // filter by entities to be written, but also restrict to entities
02240     // owned by the proc
02241 
02242     // Get list of interface entities
02243     Range imesh, tmp;
02244     for( std::list< ExportSet >::reverse_iterator i = exportList.rbegin(); i != exportList.rend(); ++i )
02245     {
02246         tmp.clear();
02247         rval = myPcomm->filter_pstatus( i->range, PSTATUS_SHARED, PSTATUS_AND, -1, &tmp );
02248         if( MB_SUCCESS != rval ) return error( rval );
02249         imesh.merge( tmp );
02250     }
02251     tmp.clear();
02252     rval = myPcomm->filter_pstatus( nodeSet.range, PSTATUS_SHARED, PSTATUS_AND, -1, &tmp );
02253     if( MB_SUCCESS != rval ) return error( rval );
02254     imesh.merge( tmp );
02255 
02256     // Create tag to store file IDs
02257     EntityHandle default_val = 0;
02258     Tag file_id_tag          = 0;
02259     rval = iFace->tag_get_handle( "__hdf5_ll_fileid", 1, MB_TYPE_HANDLE, file_id_tag, MB_TAG_DENSE | MB_TAG_CREAT,
02260                                   &default_val );
02261     if( MB_SUCCESS != rval ) return error( rval );
02262 
02263     // Copy file IDs into tag
02264     std::vector< EntityHandle > file_id_vect( imesh.size() );
02265     Range::const_iterator i;
02266     std::vector< EntityHandle >::iterator j = file_id_vect.begin();
02267     for( i = imesh.begin(); i != imesh.end(); ++i, ++j )
02268     {
02269         *j = idMap.find( *i );
02270         if( !*j )
02271         {
02272             iFace->tag_delete( file_id_tag );
02273             return error( MB_FAILURE );
02274         }
02275     }
02276     rval = iFace->tag_set_data( file_id_tag, imesh, &file_id_vect[0] );
02277     if( MB_SUCCESS != rval )
02278     {
02279         iFace->tag_delete( file_id_tag );
02280         return error( rval );
02281     }
02282 
02283     // Do communication
02284     rval = myPcomm->exchange_tags( file_id_tag, imesh );
02285     if( MB_SUCCESS != rval )
02286     {
02287         iFace->tag_delete( file_id_tag );
02288         return error( rval );
02289     }
02290 
02291     // Copy file IDs from tag into idMap for remote entities
02292     file_id_vect.resize( nonlocal.size() );
02293     rval = iFace->tag_get_data( file_id_tag, nonlocal, &file_id_vect[0] );
02294     if( MB_SUCCESS != rval )
02295     {
02296         iFace->tag_delete( file_id_tag );
02297         return error( rval );
02298     }
02299 
02300     j = file_id_vect.begin();
02301     for( i = nonlocal.begin(); i != nonlocal.end(); ++i, ++j )
02302     {
02303         if( *j == 0 )
02304         {
02305             int owner = -1;
02306             myPcomm->get_owner( *i, owner );
02307             const char* name = CN::EntityTypeName( TYPE_FROM_HANDLE( *i ) );
02308             int id           = ID_FROM_HANDLE( *i );
02309             MB_SET_ERR_CONT( "Process " << myPcomm->proc_config().proc_rank()
02310                                         << " did not receive valid id handle for shared " << name << " " << id
02311                                         << " owned by process " << owner );
02312             dbgOut.printf( 1,
02313                            "Did not receive valid remote id for "
02314                            "shared %s %d owned by process %d",
02315                            name, id, owner );
02316             iFace->tag_delete( file_id_tag );
02317             return error( MB_FAILURE );
02318         }
02319         else
02320         {
02321             if( !idMap.insert( *i, *j, 1 ).second )
02322             {
02323                 iFace->tag_delete( file_id_tag );
02324                 return error( MB_FAILURE );
02325             }
02326         }
02327     }
02328 
02329 #ifndef NDEBUG
02330     // Check that writer is correct with regards to which entities
02331     // that it owns by verifying that the file ids that we thought
02332     // we were sending where not received instead
02333     file_id_vect.resize( imesh.size() );
02334     rval = iFace->tag_get_data( file_id_tag, imesh, &file_id_vect[0] );
02335     if( MB_SUCCESS != rval )
02336     {
02337         iFace->tag_delete( file_id_tag );
02338         return error( rval );
02339     }
02340     int invalid_count = 0;
02341     j                 = file_id_vect.begin();
02342     for( i = imesh.begin(); i != imesh.end(); ++i, ++j )
02343     {
02344         EntityHandle h = idMap.find( *i );
02345         if( *j != h )
02346         {
02347             ++invalid_count;
02348             dbgOut.printf( 1, "Conflicting ownership for %s %ld\n", CN::EntityTypeName( TYPE_FROM_HANDLE( *i ) ),
02349                            (long)ID_FROM_HANDLE( *i ) );
02350         }
02351     }
02352     if( invalid_count )
02353     {
02354         iFace->tag_delete( file_id_tag );
02355         MB_SET_ERR( MB_FAILURE, invalid_count << " entities with conflicting ownership found by process "
02356                                               << myPcomm->proc_config().proc_rank()
02357                                               << ". This will result in duplicate entities written to file" );
02358     }
02359 #endif
02360 
02361     return iFace->tag_delete( file_id_tag );
02362 }
02363 
02364 void WriteHDF5Parallel::print_times( const double* times ) const
02365 {
02366     if( !myPcomm ) { WriteHDF5::print_times( times ); }
02367     else
02368     {
02369         double recv[NUM_TIMES];
02370         MPI_Reduce( (void*)times, recv, NUM_TIMES, MPI_DOUBLE, MPI_MAX, 0, myPcomm->proc_config().proc_comm() );
02371         if( 0 == myPcomm->proc_config().proc_rank() ) WriteHDF5::print_times( recv );
02372     }
02373 }
02374 
02375 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines