ReadParallel.cpp

Go to the documentation of this file.

#include "ReadParallel.hpp"
#include "moab/Core.hpp"
#include "moab/ProcConfig.hpp"
#include "moab/FileOptions.hpp"
#include "moab/Error.hpp"
#include "moab/ReaderWriterSet.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
 
#include <iostream>
#include <iomanip>
#include <iterator>
#include <sstream>
#include <algorithm>
#include <cassert>
 
namespace moab
{
 
const bool debug = false;
 
const char* ReadParallel::ParallelActionsNames[] = { "PARALLEL READ",
                                                     "PARALLEL READ PART",
                                                     "PARALLEL BROADCAST",
                                                     "PARALLEL DELETE NONLOCAL",
                                                     "PARALLEL CHECK_GIDS_SERIAL",
                                                     "PARALLEL GET_FILESET_ENTS",
                                                     "PARALLEL RESOLVE_SHARED_ENTS",
                                                     "PARALLEL EXCHANGE_GHOSTS",
                                                     "PARALLEL RESOLVE_SHARED_SETS",
                                                     "PARALLEL_AUGMENT_SETS_WITH_GHOSTS",
                                                     "PARALLEL PRINT_PARALLEL",
                                                     "PARALLEL_CREATE_TRIVIAL_PARTITION",
                                                     "PARALLEL_CORRECT_THIN_GHOST_LAYERS" };
 
const char* ReadParallel::parallelOptsNames[] = { "NONE", "BCAST", "BCAST_DELETE", "READ_DELETE", "READ_PART", "", 0 };
 
ReadParallel::ReadParallel( Interface* impl, ParallelComm* pc )
    : mbImpl( impl ), myPcomm( pc ), myDebug( "ReadPara", std::cerr )
{
    if( !myPcomm )
    {
        myPcomm = ParallelComm::get_pcomm( mbImpl, 0 );
        if( NULL == myPcomm ) myPcomm = new ParallelComm( mbImpl, MPI_COMM_WORLD );
    }
    myDebug.set_rank( myPcomm->proc_config().proc_rank() );
    if( debug )  // For backwards compatibility, enable all debug output if constant is true
        myDebug.set_verbosity( 10 );
 
    impl->query_interface( mError );
}
 
ErrorCode ReadParallel::load_file( const char** file_names,
                                   const int num_files,
                                   const EntityHandle* file_set,
                                   const FileOptions& opts,
                                   const ReaderIface::SubsetList* subset_list,
                                   const Tag* file_id_tag )
{
    int tmpval;
    if( MB_SUCCESS == opts.get_int_option( "DEBUG_PIO", 1, tmpval ) )
    {
        myDebug.set_verbosity( tmpval );
        myPcomm->set_debug_verbosity( tmpval );
    }
    myDebug.tprint( 1, "Setting up...\n" );
 
    // Get parallel settings
    int parallel_mode;
    ErrorCode result = opts.match_option( "PARALLEL", parallelOptsNames, parallel_mode );
    if( MB_FAILURE == result )
    {
        MB_SET_ERR( MB_FAILURE, "Unexpected value for 'PARALLEL' option" );
    }
    else if( MB_ENTITY_NOT_FOUND == result )
    {
        parallel_mode = 0;
    }
 
    // Get partition setting
    bool distrib;
    std::string partition_tag_name;
    result = opts.get_option( "PARTITION", partition_tag_name );
    if( MB_ENTITY_NOT_FOUND == result )
    {
        distrib            = false;
        partition_tag_name = "";
    }
    else
    {
        distrib = true;
        if( partition_tag_name.empty() ) partition_tag_name = PARALLEL_PARTITION_TAG_NAME;
 
        // Also get deprecated PARTITION_DISTRIBUTE option
        // so that higher-level code doesn't return an error
        // due to an unrecognized option
        opts.get_null_option( "PARTITION_DISTRIBUTE" );
    }
 
    // Get partition tag value(s), if any, and whether they're to be
    // distributed or assigned
    std::vector< int > partition_tag_vals;
    opts.get_ints_option( "PARTITION_VAL", partition_tag_vals );
 
    // see if partition tag name is "TRIVIAL", if the tag exists
    bool create_trivial_partition = false;
    if( partition_tag_name == std::string( "TRIVIAL" ) )
    {
        Tag ttag;  // see if the trivial tag exists
        result = mbImpl->tag_get_handle( partition_tag_name.c_str(), ttag );
        if( MB_TAG_NOT_FOUND == result ) create_trivial_partition = true;
    }
    // See if we need to report times
    bool cputime = false;
    result       = opts.get_null_option( "CPUTIME" );
    if( MB_SUCCESS == result ) cputime = true;
 
    // See if we need to report times
    bool print_parallel = false;
    result              = opts.get_null_option( "PRINT_PARALLEL" );
    if( MB_SUCCESS == result ) print_parallel = true;
 
    // Get ghosting options
    std::string ghost_str;
    int bridge_dim, ghost_dim = -1, num_layers, addl_ents = 0;
    result = opts.get_str_option( "PARALLEL_GHOSTS", ghost_str );
    if( MB_TYPE_OUT_OF_RANGE == result )
    {
        ghost_dim  = 3;
        bridge_dim = 0;
        num_layers = 1;
    }
    else if( MB_SUCCESS == result )
    {
        int num_fields = sscanf( ghost_str.c_str(), "%d.%d.%d.%d", &ghost_dim, &bridge_dim, &num_layers, &addl_ents );
        if( 3 > num_fields )
        {
            MB_SET_ERR( MB_FAILURE, "Didn't read 3 fields from PARALLEL_GHOSTS string" );
        }
    }
 
    // Get resolve_shared_ents option
    std::string shared_str;
    int resolve_dim = -2, shared_dim = -1;
    result = opts.get_str_option( "PARALLEL_RESOLVE_SHARED_ENTS", shared_str );
    if( MB_TYPE_OUT_OF_RANGE == result )
    {
        resolve_dim = -1;
        shared_dim  = -1;
    }
    else if( MB_SUCCESS == result )
    {
        int num_fields = sscanf( shared_str.c_str(), "%d.%d", &resolve_dim, &shared_dim );
        if( 2 != num_fields )
        {
            MB_SET_ERR( MB_FAILURE, "Didn't read 2 fields from PARALLEL_RESOLVE_SHARED_ENTS string" );
        }
    }
 
    // Get skip augmenting with ghosts option
    bool skip_augment = false;
    result            = opts.get_null_option( "SKIP_AUGMENT_WITH_GHOSTS" );
    if( MB_SUCCESS == result ) skip_augment = true;
 
    bool correct_thin_ghosts = false;
    result                   = opts.get_null_option( "PARALLEL_THIN_GHOST_LAYER" );
    if( MB_SUCCESS == result ) correct_thin_ghosts = true;
 
    // Get MPI IO processor rank
    int reader_rank;
    result = opts.get_int_option( "MPI_IO_RANK", reader_rank );
    if( MB_ENTITY_NOT_FOUND == result )
        reader_rank = 0;
    else if( MB_SUCCESS != result )
    {
        MB_SET_ERR( MB_FAILURE, "Unexpected value for 'MPI_IO_RANK' option" );
    }
 
    // Now that we've parsed all the parallel options, make an instruction queue
    std::vector< int > pa_vec;
    bool is_reader = ( reader_rank == (int)myPcomm->proc_config().proc_rank() );
 
    bool partition_by_rank = false;
    if( MB_SUCCESS == opts.get_null_option( "PARTITION_BY_RANK" ) )
    {
        partition_by_rank = true;
        if( !partition_tag_vals.empty() )
        {
            MB_SET_ERR( MB_FAILURE, "Cannot specify both PARTITION_VALS and PARTITION_BY_RANK" );
        }
    }
 
    double factor_seq;
    if( MB_SUCCESS == opts.get_real_option( "PARALLEL_SEQUENCE_FACTOR", factor_seq ) )
    {
        if( factor_seq < 1. ) MB_SET_ERR( MB_FAILURE, "cannot have sequence factor less than 1." );
        mbImpl->set_sequence_multiplier( factor_seq );
    }
    switch( parallel_mode )
    {
        case POPT_BCAST:
            myDebug.print( 1, "Read mode is BCAST\n" );
            if( is_reader )
            {
                pa_vec.push_back( PA_READ );
                pa_vec.push_back( PA_CHECK_GIDS_SERIAL );
                pa_vec.push_back( PA_GET_FILESET_ENTS );
            }
            pa_vec.push_back( PA_BROADCAST );
            if( !is_reader ) pa_vec.push_back( PA_GET_FILESET_ENTS );
            break;
 
        case POPT_BCAST_DELETE:
            myDebug.print( 1, "Read mode is BCAST_DELETE\n" );
            if( is_reader )
            {
                pa_vec.push_back( PA_READ );
                pa_vec.push_back( PA_CHECK_GIDS_SERIAL );
                pa_vec.push_back( PA_GET_FILESET_ENTS );
                if( create_trivial_partition ) pa_vec.push_back( PA_CREATE_TRIVIAL_PARTITION );
            }
            pa_vec.push_back( PA_BROADCAST );
            if( !is_reader ) pa_vec.push_back( PA_GET_FILESET_ENTS );
            pa_vec.push_back( PA_DELETE_NONLOCAL );
            break;
 
        case POPT_DEFAULT:
        case POPT_READ_DELETE:
            myDebug.print( 1, "Read mode is READ_DELETE\n" );
            pa_vec.push_back( PA_READ );
            pa_vec.push_back( PA_CHECK_GIDS_SERIAL );
            pa_vec.push_back( PA_GET_FILESET_ENTS );
            pa_vec.push_back( PA_DELETE_NONLOCAL );
            break;
 
        case POPT_READ_PART:
            myDebug.print( 1, "Read mode is READ_PART\n" );
            pa_vec.push_back( PA_READ_PART );
            break;
 
        default:
            MB_SET_ERR( MB_FAILURE, "Unexpected parallel read mode" );
    }
 
    if( -2 != resolve_dim ) pa_vec.push_back( PA_RESOLVE_SHARED_ENTS );
 
    if( -1 != ghost_dim ) pa_vec.push_back( PA_EXCHANGE_GHOSTS );
 
    if( -2 != resolve_dim )
    {
        pa_vec.push_back( PA_RESOLVE_SHARED_SETS );
        if( -1 != ghost_dim && !skip_augment ) pa_vec.push_back( PA_AUGMENT_SETS_WITH_GHOSTS );
        if( -1 != ghost_dim && correct_thin_ghosts ) pa_vec.push_back( PA_CORRECT_THIN_GHOSTS );
    }
 
    if( print_parallel ) pa_vec.push_back( PA_PRINT_PARALLEL );
 
    result = load_file( file_names, num_files, file_set, parallel_mode, partition_tag_name, partition_tag_vals, distrib,
                        partition_by_rank, pa_vec, opts, subset_list, file_id_tag, reader_rank, cputime, resolve_dim,
                        shared_dim, ghost_dim, bridge_dim, num_layers, addl_ents );MB_CHK_ERR( result );
 
    if( parallel_mode == POPT_BCAST_DELETE && !is_reader ) opts.mark_all_seen();
 
    return MB_SUCCESS;
}
 
ErrorCode ReadParallel::load_file( const char** file_names,
                                   const int num_files,
                                   const EntityHandle* file_set_ptr,
                                   int /*parallel_mode*/,
                                   std::string& partition_tag_name,
                                   std::vector< int >& partition_tag_vals,
                                   bool distrib,
                                   bool partition_by_rank,
                                   std::vector< int >& pa_vec,
                                   const FileOptions& opts,
                                   const ReaderIface::SubsetList* subset_list,
                                   const Tag* file_id_tag,
                                   const int reader_rank,
                                   const bool cputime,
                                   const int resolve_dim,
                                   const int shared_dim,
                                   const int ghost_dim,
                                   const int bridge_dim,
                                   const int num_layers,
                                   const int addl_ents )
{
    ErrorCode result = MB_SUCCESS;
    if( myPcomm == NULL ) myPcomm = new ParallelComm( mbImpl, MPI_COMM_WORLD );
 
    Range entities;
    Tag file_set_tag = 0;
    int other_sets   = 0;
    ReaderWriterSet::iterator iter;
    Range other_file_sets, file_sets;
    Core* impl = dynamic_cast< Core* >( mbImpl );
 
    std::vector< double > act_times( pa_vec.size() + 1 );
    std::vector< int >::iterator vit;
    int i, j;
    act_times[0] = MPI_Wtime();
 
    // Make a new set for the parallel read
    EntityHandle file_set;
    if( !file_set_ptr || !( *file_set_ptr ) )
    {
        result = mbImpl->create_meshset( MESHSET_SET, file_set );MB_CHK_SET_ERR( result, "Trouble creating file set" );
    }
    else
        file_set = *file_set_ptr;
 
    bool i_read     = false;
    Tag id_tag      = 0;
    bool use_id_tag = false;
    Range ents;
 
    for( i = 1, vit = pa_vec.begin(); vit != pa_vec.end(); ++vit, i++ )
    {
        ErrorCode tmp_result = MB_SUCCESS;
        switch( *vit )
        {
                //==================
            case PA_READ:
                i_read = true;
 
                for( j = 0; j < num_files; j++ )
                {
                    myDebug.tprintf( 1, "Reading file: \"%s\"\n", file_names[j] );
 
                    EntityHandle new_file_set;
                    result = mbImpl->create_meshset( MESHSET_SET, new_file_set );MB_CHK_SET_ERR( result, "Trouble creating file set" );
                    tmp_result = impl->serial_load_file( file_names[j], &new_file_set, opts, subset_list, file_id_tag );
                    if( MB_SUCCESS != tmp_result ) break;
 
                    // Put the contents of each file set for the reader into the
                    // file set for the parallel read
                    assert( 0 != new_file_set );
                    Range all_ents;
                    tmp_result = mbImpl->get_entities_by_handle( new_file_set, all_ents );
                    if( MB_SUCCESS != tmp_result ) break;
                    all_ents.insert( new_file_set );
                    tmp_result = mbImpl->add_entities( file_set, all_ents );
                    if( MB_SUCCESS != tmp_result ) break;
                }
                if( MB_SUCCESS != tmp_result ) break;
 
                // Mark the file set for this parallel reader
                tmp_result = mbImpl->tag_get_handle( "__file_set", 1, MB_TYPE_INTEGER, file_set_tag,
                                                     MB_TAG_SPARSE | MB_TAG_CREAT );
                if( MB_SUCCESS != tmp_result ) break;
 
                tmp_result = mbImpl->tag_set_data( file_set_tag, &file_set, 1, &other_sets );
                break;
 
                //==================
            case PA_READ_PART: {
                myDebug.tprintf( 1, "Reading file: \"%s\"\n", file_names[0] );
 
                i_read = true;
                if( num_files != 1 )
                {
                    MB_SET_ERR( MB_NOT_IMPLEMENTED, "Multiple file read not supported for READ_PART" );
                }
 
                // If we're going to resolve shared entities, then we need
                // to ask the file reader to populate a tag with unique ids
                // (typically file ids/indices/whatever.)
                if( std::find( pa_vec.begin(), pa_vec.end(), PA_RESOLVE_SHARED_ENTS ) != pa_vec.end() )
                {
                    use_id_tag = true;
                    if( !file_id_tag )
                    {
                        // This tag is really used for resolving shared entities with crystal router
                        // In the end, this is an identifier that gets converted to long
                        // In hdf5 file reader, we also convert from hdf5 file id type to long
                        tmp_result = mbImpl->tag_get_handle( "", sizeof( long ), MB_TYPE_OPAQUE, id_tag,
                                                             MB_TAG_DENSE | MB_TAG_CREAT );
                        if( MB_SUCCESS != tmp_result ) break;
                        file_id_tag = &id_tag;
                    }
                }
 
                ReaderIface::IDTag parts = { partition_tag_name.c_str(), 0, 0 };
                ReaderIface::SubsetList sl;
                sl.num_parts = 0;
                int rank     = myPcomm->rank();
                if( partition_by_rank )
                {
                    assert( partition_tag_vals.empty() );
                    parts.tag_values     = &rank;
                    parts.num_tag_values = 1;
                }
                else
                {
                    sl.num_parts   = myPcomm->size();
                    sl.part_number = myPcomm->rank();
                    if( !partition_tag_vals.empty() )
                    {
                        parts.tag_values     = &partition_tag_vals[0];
                        parts.num_tag_values = partition_tag_vals.size();
                    }
                }
                std::vector< ReaderIface::IDTag > subset;
                if( subset_list )
                {
                    std::vector< ReaderIface::IDTag > tmplist( subset_list->tag_list,
                                                               subset_list->tag_list + subset_list->tag_list_length );
                    tmplist.push_back( parts );
                    subset.swap( tmplist );
                    sl.tag_list        = &subset[0];
                    sl.tag_list_length = subset.size();
                }
                else
                {
                    sl.tag_list        = &parts;
                    sl.tag_list_length = 1;
                }
                tmp_result = impl->serial_load_file( *file_names, &file_set, opts, &sl, file_id_tag );
                if( MB_SUCCESS != tmp_result ) break;
 
                if( !partition_tag_name.empty() )
                {
                    Tag part_tag;
                    tmp_result = impl->tag_get_handle( partition_tag_name.c_str(), 1, MB_TYPE_INTEGER, part_tag );
                    if( MB_SUCCESS != tmp_result ) break;
 
                    tmp_result = impl->get_entities_by_type_and_tag( file_set, MBENTITYSET, &part_tag, 0, 1,
                                                                     myPcomm->partition_sets() );
                }
            }
            break;
 
                //==================
            case PA_GET_FILESET_ENTS:
                myDebug.tprint( 1, "Getting fileset entities.\n" );
 
                // Get entities in the file set, and add actual file set to it;
                // mark the file set to make sure any receiving procs know which it is
                tmp_result = mbImpl->get_entities_by_handle( file_set, entities );
                if( MB_SUCCESS != tmp_result )
                {
                    entities.clear();
                    break;
                }
 
                // Add actual file set to entities too
                entities.insert( file_set );
                break;
                //==================
            case PA_CREATE_TRIVIAL_PARTITION: {
                myDebug.tprint( 1, "create trivial partition, for higher dim entities.\n" );
                // get high dim entities (2 or 3)
                Range hi_dim_ents = entities.subset_by_dimension( 3 );
                if( hi_dim_ents.empty() ) hi_dim_ents = entities.subset_by_dimension( 2 );
                if( hi_dim_ents.empty() ) hi_dim_ents = entities.subset_by_dimension( 1 );
                if( hi_dim_ents.empty() ) MB_SET_ERR( MB_FAILURE, "there are no elements of dim 1-3" );
 
                size_t num_hi_ents     = hi_dim_ents.size();
                unsigned int num_parts = myPcomm->size();
 
                // create first the trivial partition tag
                int dum_id = -1;
                Tag ttag;  // trivial tag
                tmp_result = mbImpl->tag_get_handle( partition_tag_name.c_str(), 1, MB_TYPE_INTEGER, ttag,
                                                     MB_TAG_CREAT | MB_TAG_SPARSE, &dum_id );MB_CHK_SET_ERR( tmp_result, "Can't create trivial partition tag" );
 
                // Compute the number of high dim entities on each part
                size_t nPartEnts = num_hi_ents / num_parts;
 
                // Number of extra entities after equal split over parts
                int iextra          = num_hi_ents % num_parts;
                Range::iterator itr = hi_dim_ents.begin();
                for( int k = 0; k < (int)num_parts; k++ )
                {
                    // create a mesh set, insert a subrange of entities
                    EntityHandle part_set;
                    tmp_result = mbImpl->create_meshset( MESHSET_SET, part_set );MB_CHK_SET_ERR( tmp_result, "Can't create part set" );
                    entities.insert( part_set );
 
                    tmp_result = mbImpl->tag_set_data( ttag, &part_set, 1, &k );MB_CHK_SET_ERR( tmp_result, "Can't set trivial partition tag" );
                    Range subrange;
                    size_t num_ents_in_part = nPartEnts;
                    if( i < iextra ) num_ents_in_part++;
                    for( size_t i1 = 0; i1 < num_ents_in_part; i1++, itr++ )
                        subrange.insert( *itr );
                    tmp_result = mbImpl->add_entities( part_set, subrange );MB_CHK_SET_ERR( tmp_result, "Can't add entities to trivial part " << k );
                    myDebug.tprintf( 1, "create trivial part %d with %lu entities \n", k, num_ents_in_part );
                    tmp_result = mbImpl->add_entities( file_set, &part_set, 1 );MB_CHK_SET_ERR( tmp_result, "Can't add trivial part to file set " << k );
                }
            }
 
            break;
                //==================
            case PA_BROADCAST:
                // Do the actual broadcast; if single-processor, ignore error
                myDebug.tprint( 1, "Broadcasting mesh.\n" );
 
                if( myPcomm->proc_config().proc_size() > 1 )
                {
                    tmp_result = myPcomm->broadcast_entities( reader_rank, entities );
                    if( MB_SUCCESS != tmp_result ) break;
                }
 
                if( debug )
                {
                    std::cerr << "Bcast done; entities:" << std::endl;
                    mbImpl->list_entities( 0, 0 );
                }
 
                // Add the received entities to this fileset if I wasn't the reader
                if( !i_read && MB_SUCCESS == tmp_result ) tmp_result = mbImpl->add_entities( file_set, entities );
 
                break;
 
                //==================
            case PA_DELETE_NONLOCAL:
                myDebug.tprint( 1, "Deleting nonlocal entities.\n" );
 
                tmp_result = delete_nonlocal_entities( partition_tag_name, partition_tag_vals, distrib, file_set );
                if( debug )
                {
                    std::cerr << "Delete nonlocal done; entities:" << std::endl;
                    mbImpl->list_entities( 0, 0 );
                }
 
                if( MB_SUCCESS == tmp_result ) tmp_result = create_partition_sets( partition_tag_name, file_set );
 
                break;
 
                //==================
            case PA_CHECK_GIDS_SERIAL:
                myDebug.tprint( 1, "Checking global IDs.\n" );
 
                tmp_result = myPcomm->check_global_ids( file_set, 0, 1, true, false );
                break;
 
                //==================
            case PA_RESOLVE_SHARED_ENTS:
                myDebug.tprint( 1, "Resolving shared entities.\n" );
 
                if( 1 == myPcomm->size() )
                    tmp_result = MB_SUCCESS;
                else
                    tmp_result =
                        myPcomm->resolve_shared_ents( file_set, resolve_dim, shared_dim, use_id_tag ? file_id_tag : 0 );
                if( MB_SUCCESS != tmp_result ) break;
 
#ifndef NDEBUG
                // check number of owned vertices through pcomm's public interface
                tmp_result = mbImpl->get_entities_by_type( 0, MBVERTEX, ents );
                if( MB_SUCCESS == tmp_result )
                    tmp_result = myPcomm->filter_pstatus( ents, PSTATUS_NOT_OWNED, PSTATUS_NOT );
                if( MB_SUCCESS == tmp_result )
                    myDebug.tprintf( 1, "Proc %u reports %lu owned vertices.\n", myPcomm->proc_config().proc_rank(),
                                     ents.size() );
#endif
                break;
 
                //==================
            case PA_EXCHANGE_GHOSTS:
                myDebug.tprint( 1, "Exchanging ghost entities.\n" );
 
                tmp_result = myPcomm->exchange_ghost_cells( ghost_dim, bridge_dim, num_layers, addl_ents, true, true,
                                                            &file_set );
                break;
 
                //==================
            case PA_RESOLVE_SHARED_SETS:
                myDebug.tprint( 1, "Resolving shared sets.\n" );
 
                if( 1 == myPcomm->size() )
                    tmp_result = MB_SUCCESS;
                else
                    tmp_result = myPcomm->resolve_shared_sets( file_set, use_id_tag ? file_id_tag : 0 );
                break;
 
                //==================
            case PA_AUGMENT_SETS_WITH_GHOSTS:
                myDebug.tprint( 1, "Augmenting sets with ghost entities.\n" );
 
                if( 1 == myPcomm->size() )
                    tmp_result = MB_SUCCESS;
                else
                    tmp_result = myPcomm->augment_default_sets_with_ghosts( file_set );
                break;
                //==================
            case PA_CORRECT_THIN_GHOSTS:
                myDebug.tprint( 1, "correcting thin ghost layers.\n" );
                if( 2 >= myPcomm->size() )  // it is a problem only for multi-shared entities
                    tmp_result = MB_SUCCESS;
                else
                    tmp_result = myPcomm->correct_thin_ghost_layers();
                break;
            case PA_PRINT_PARALLEL:
                myDebug.tprint( 1, "Printing parallel information.\n" );
 
                tmp_result = myPcomm->list_entities( 0, -1 );
                break;
 
                //==================
            default:
                MB_SET_ERR( MB_FAILURE, "Unexpected parallel action" );
        }  // switch (*vit)
 
        if( MB_SUCCESS != tmp_result )
        {
            MB_SET_ERR( MB_FAILURE, "Failed in step " << ParallelActionsNames[*vit] );
        }
 
        if( cputime ) act_times[i] = MPI_Wtime();
    }  // for (i = 1, vit = pa_vec.begin(); vit != pa_vec.end(); ++vit, i++)
 
    if( use_id_tag )
    {
        result = mbImpl->tag_delete( id_tag );MB_CHK_SET_ERR( result, "Trouble deleting id tag" );
    }
 
    if( cputime )
    {
        for( i = pa_vec.size(); i > 0; i-- )
            act_times[i] -= act_times[i - 1];
 
        // Replace initial time with overall time
        act_times[0] = MPI_Wtime() - act_times[0];
        // Get the maximum over all procs
        if( 0 != myPcomm->proc_config().proc_rank() )
        {
            MPI_Reduce( &act_times[0], 0, pa_vec.size() + 1, MPI_DOUBLE, MPI_MAX, 0,
                        myPcomm->proc_config().proc_comm() );
        }
        else
        {
#if( MPI_VERSION >= 2 )
            MPI_Reduce( MPI_IN_PLACE, &act_times[0], pa_vec.size() + 1, MPI_DOUBLE, MPI_MAX, 0,
                        myPcomm->proc_config().proc_comm() );
#else
            // Note, extra comm-size allocation is required
            std::vector< double > act_times_tmp( pa_vec.size() + 1 );
            MPI_Reduce( &act_times[0], &act_times_tmp[0], pa_vec.size() + 1, MPI_DOUBLE, MPI_MAX, 0,
                        myPcomm->proc_config().proc_comm() );
            act_times = act_times_tmp;  // extra copy here too
#endif
            std::cout << "Parallel Read times: " << std::endl;
            for( i = 1, vit = pa_vec.begin(); vit != pa_vec.end(); ++vit, i++ )
                std::cout << "  " << act_times[i] << " " << ParallelActionsNames[*vit] << std::endl;
            std::cout << "  " << act_times[0] << " PARALLEL TOTAL" << std::endl;
        }
    }
 
    return MB_SUCCESS;
}
 
ErrorCode ReadParallel::delete_nonlocal_entities( std::string& ptag_name,
                                                  std::vector< int >& ptag_vals,
                                                  bool distribute,
                                                  EntityHandle file_set )
{
    Range partition_sets;
 
    Tag ptag;
    ErrorCode result = mbImpl->tag_get_handle( ptag_name.c_str(), 1, MB_TYPE_INTEGER, ptag );MB_CHK_SET_ERR( result, "Failed getting tag handle in delete_nonlocal_entities" );
 
    result = mbImpl->get_entities_by_type_and_tag( file_set, MBENTITYSET, &ptag, NULL, 1, myPcomm->partition_sets() );MB_CHK_SET_ERR( result, "Failed to get sets with partition-type tag" );
 
    int proc_sz = myPcomm->proc_config().proc_size();
    int proc_rk = myPcomm->proc_config().proc_rank();
 
    if( !ptag_vals.empty() )
    {
        // Values input, get sets with those values
        Range tmp_sets;
        std::vector< int > tag_vals( myPcomm->partition_sets().size() );
        result = mbImpl->tag_get_data( ptag, myPcomm->partition_sets(), &tag_vals[0] );MB_CHK_SET_ERR( result, "Failed to get tag data for partition vals tag" );
        for( std::vector< int >::iterator pit = tag_vals.begin(); pit != tag_vals.end(); ++pit )
        {
            std::vector< int >::iterator pit2 = std::find( ptag_vals.begin(), ptag_vals.end(), *pit );
            if( pit2 != ptag_vals.end() ) tmp_sets.insert( myPcomm->partition_sets()[pit - tag_vals.begin()] );
        }
 
        myPcomm->partition_sets().swap( tmp_sets );
    }
 
    if( distribute )
    {
        // For now, require that number of partition sets be greater
        // than number of procs
        if( myPcomm->partition_sets().size() < (unsigned int)proc_sz )
        {
            MB_SET_ERR( MB_FAILURE, "Too few parts; P = " << proc_rk << ", tag = " << ptag
                                                          << ", # sets = " << myPcomm->partition_sets().size() );
        }
 
        Range tmp_sets;
        // Distribute the partition sets
        unsigned int num_sets     = myPcomm->partition_sets().size() / proc_sz;
        unsigned int num_leftover = myPcomm->partition_sets().size() % proc_sz;
        int begin_set             = 0;
        if( proc_rk < (int)num_leftover )
        {
            num_sets++;
            begin_set = num_sets * proc_rk;
        }
        else
            begin_set = proc_rk * num_sets + num_leftover;
 
        for( unsigned int i = 0; i < num_sets; i++ )
            tmp_sets.insert( myPcomm->partition_sets()[begin_set + i] );
 
        myPcomm->partition_sets().swap( tmp_sets );
    }
 
    myDebug.print( 1, "My partition sets: ", myPcomm->partition_sets() );
 
    result = delete_nonlocal_entities( file_set );MB_CHK_ERR( result );
 
    return MB_SUCCESS;
}
 
ErrorCode ReadParallel::create_partition_sets( std::string& ptag_name, EntityHandle file_set )
{
    int proc_rk      = myPcomm->proc_config().proc_rank();
    ErrorCode result = MB_SUCCESS;
 
    Tag ptag;
 
    // Tag the partition sets with a standard tag name
    if( ptag_name.empty() ) ptag_name = PARALLEL_PARTITION_TAG_NAME;
    bool tag_created = false;
    result = mbImpl->tag_get_handle( ptag_name.c_str(), 1, MB_TYPE_INTEGER, ptag, MB_TAG_SPARSE | MB_TAG_CREAT, 0,
                                     &tag_created );MB_CHK_SET_ERR( result, "Trouble getting PARALLEL_PARTITION tag" );
 
    if( !tag_created )
    {
        // This tag already exists; better check to see that tagged sets
        // agree with this partition
        Range tagged_sets;
        int* proc_rk_ptr = &proc_rk;
        result = mbImpl->get_entities_by_type_and_tag( file_set, MBENTITYSET, &ptag, (const void* const*)&proc_rk_ptr,
                                                       1, tagged_sets );MB_CHK_SET_ERR( result, "Trouble getting tagged sets" );
        if( !tagged_sets.empty() && tagged_sets != myPcomm->partition_sets() )
        {
            result = mbImpl->tag_delete_data( ptag, tagged_sets );MB_CHK_SET_ERR( result, "Trouble deleting data of PARALLEL_PARTITION tag" );
        }
        else if( tagged_sets == myPcomm->partition_sets() )
            return MB_SUCCESS;
    }
 
    // If we get here, we need to assign the tag
    std::vector< int > values( myPcomm->partition_sets().size() );
    for( unsigned int i = 0; i < myPcomm->partition_sets().size(); i++ )
        values[i] = proc_rk;
    result = mbImpl->tag_set_data( ptag, myPcomm->partition_sets(), &values[0] );MB_CHK_SET_ERR( result, "Trouble setting data to PARALLEL_PARTITION tag" );
 
    return MB_SUCCESS;
}
 
ErrorCode ReadParallel::delete_nonlocal_entities( EntityHandle file_set )
{
    // Get partition entities and ents related to/used by those
    // get ents in the partition
    ReadUtilIface* read_iface;
    mbImpl->query_interface( read_iface );
    Range partition_ents, all_sets;
 
    myDebug.tprint( 2, "Gathering related entities.\n" );
 
    ErrorCode result = read_iface->gather_related_ents( myPcomm->partition_sets(), partition_ents, &file_set );MB_CHK_SET_ERR( result, "Failure gathering related entities" );
 
    // Get pre-existing entities
    Range file_ents;
    result = mbImpl->get_entities_by_handle( file_set, file_ents );MB_CHK_SET_ERR( result, "Couldn't get pre-existing entities" );
 
    if( 0 == myPcomm->proc_config().proc_rank() )
    {
        myDebug.print( 2, "File entities: ", file_ents );
    }
 
    // Get deletable entities by subtracting partition ents from file ents
    Range deletable_ents = subtract( file_ents, partition_ents );
 
    // Cache deletable vs. keepable sets
    Range deletable_sets = deletable_ents.subset_by_type( MBENTITYSET );
    Range keepable_sets  = subtract( file_ents.subset_by_type( MBENTITYSET ), deletable_sets );
 
    myDebug.tprint( 2, "Removing deletable entities from keepable sets.\n" );
 
    // Remove deletable ents from all keepable sets
    for( Range::iterator rit = keepable_sets.begin(); rit != keepable_sets.end(); ++rit )
    {
        result = mbImpl->remove_entities( *rit, deletable_ents );MB_CHK_SET_ERR( result, "Failure removing deletable entities" );
    }
    result = mbImpl->remove_entities( file_set, deletable_ents );MB_CHK_SET_ERR( result, "Failure removing deletable entities" );
 
    myDebug.tprint( 2, "Deleting deletable entities.\n" );
 
    if( 0 == myPcomm->proc_config().proc_rank() )
    {
        myDebug.print( 2, "Deletable sets: ", deletable_sets );
    }
 
    // Delete sets, then ents
    if( !deletable_sets.empty() )
    {
        result = mbImpl->delete_entities( deletable_sets );MB_CHK_SET_ERR( result, "Failure deleting sets in delete_nonlocal_entities" );
    }
 
    deletable_ents -= deletable_sets;
 
    if( 0 == myPcomm->proc_config().proc_rank() )
    {
        myDebug.print( 2, "Deletable entities: ", deletable_ents );
    }
 
    if( !deletable_ents.empty() )
    {
        result = mbImpl->delete_entities( deletable_ents );MB_CHK_SET_ERR( result, "Failure deleting entities in delete_nonlocal_entities" );
    }
 
    return MB_SUCCESS;
}
 
}  // namespace moab