ReadHDF5Dataset.cpp

Go to the documentation of this file.

/** \file   ReadHDF5Dataset.cpp
 *  \author Jason Kraftcheck
 *  \date   2010-07-09
 */

#include <cassert>
#include <cstring>
#include <cstdlib>
#include <iostream>
#include "moab/MOABConfig.h"
#include "ReadHDF5Dataset.hpp"

#include "moab_mpe.h"

#include <H5Dpublic.h>
#include <H5Tpublic.h>
#include <H5Ppublic.h>
#include <H5Spublic.h>
#ifdef MOAB_HAVE_HDF5_PARALLEL
#include <H5FDmpi.h>
#include <H5FDmpio.h>
#endif

#define HDF5_16API ( H5_VERS_MAJOR < 2 && H5_VERS_MINOR < 8 )

namespace moab
{

// Selection of hyperslabs appears to be superlinear.  Don't try to select
// more than a few thousand at a time or things start to get real slow.
const size_t DEFAULT_HYPERSLAB_SELECTION_LIMIT  = 200;
size_t ReadHDF5Dataset::hyperslabSelectionLimit = DEFAULT_HYPERSLAB_SELECTION_LIMIT;
void ReadHDF5Dataset::default_hyperslab_selection_limit()
{
    hyperslabSelectionLimit = DEFAULT_HYPERSLAB_SELECTION_LIMIT;
}

H5S_seloper_t ReadHDF5Dataset::hyperslabSelectOp = H5S_SELECT_OR;

#ifdef MOAB_HAVE_LIBMPE
static std::pair< int, int > allocate_mpe_state( const char* name, const char* color )
{
    std::pair< int, int > result;
    result.first  = MPE_Log_get_event_number();
    result.second = MPE_Log_get_event_number();
    MPE_Describe_state( result.first, result.second, name, color );
    return result;
}
#else
static std::pair< int, int > allocate_mpe_state( const char*, const char* )
{
    return std::pair< int, int >();
}
#endif

bool ReadHDF5Dataset::haveMPEEvents = false;
std::pair< int, int > ReadHDF5Dataset::mpeReadEvent;
std::pair< int, int > ReadHDF5Dataset::mpeReduceEvent;

ReadHDF5Dataset::ReadHDF5Dataset( const char* debug_desc, bool parallel, const Comm* communicator )
    : closeDataSet( false ), dataSet( -1 ), dataSpace( -1 ), dataType( -1 ), fileType( -1 ), ioProp( H5P_DEFAULT ),
      dataSpaceRank( 0 ), rowsInTable( 0 ), doConversion( false ), nativeParallel( parallel ), readCount( 0 ),
      bufferSize( 0 ), mpiComm( communicator ), mpeDesc( debug_desc )
{
    if( !haveMPEEvents )
    {
        haveMPEEvents  = true;
        mpeReadEvent   = allocate_mpe_state( "ReadHDF5Dataset::read", "yellow" );
        mpeReduceEvent = allocate_mpe_state( "ReadHDF5Dataset::all_reduce", "yellow" );
    }

#ifndef MOAB_HAVE_HDF5_PARALLEL
    if( nativeParallel ) throw Exception( __LINE__ );
#else
    if( nativeParallel && !mpiComm ) throw Exception( __LINE__ );

    if( mpiComm )
    {
        ioProp = H5Pcreate( H5P_DATASET_XFER );
        H5Pset_dxpl_mpio( ioProp, H5FD_MPIO_COLLECTIVE );
    }
#endif
}

ReadHDF5Dataset::ReadHDF5Dataset( const char* debug_desc,
                                  hid_t data_set_handle,
                                  bool parallel,
                                  const Comm* communicator,
                                  bool close_data_set )
    : closeDataSet( close_data_set ), dataSet( data_set_handle ), dataSpace( -1 ), dataType( -1 ), fileType( -1 ),
      ioProp( H5P_DEFAULT ), dataSpaceRank( 0 ), rowsInTable( 0 ), doConversion( false ), nativeParallel( parallel ),
      readCount( 0 ), bufferSize( 0 ), mpiComm( communicator ), mpeDesc( debug_desc )
{
    if( !haveMPEEvents )
    {
        haveMPEEvents  = true;
        mpeReadEvent   = allocate_mpe_state( "ReadHDF5Dataset::read", "yellow" );
        mpeReduceEvent = allocate_mpe_state( "ReadHDF5Dataset::all_reduce", "yellow" );
    }

    init( data_set_handle, close_data_set );

#ifndef MOAB_HAVE_HDF5_PARALLEL
    if( nativeParallel ) throw Exception( __LINE__ );
#else
    if( nativeParallel && !mpiComm ) throw Exception( __LINE__ );

    if( mpiComm )
    {
        ioProp = H5Pcreate( H5P_DATASET_XFER );
        H5Pset_dxpl_mpio( ioProp, H5FD_MPIO_COLLECTIVE );
    }
#endif
}

void ReadHDF5Dataset::init( hid_t data_set_handle, bool close_data_set )
{
    closeDataSet = close_data_set;
    dataSet      = data_set_handle;

    fileType = H5Dget_type( data_set_handle );
    if( fileType < 0 ) throw Exception( __LINE__ );

    dataSpace = H5Dget_space( dataSet );
    if( dataSpace < 0 ) throw Exception( __LINE__ );

    dataSpaceRank = H5Sget_simple_extent_dims( dataSpace, dataSetCount, dataSetOffset );
    if( dataSpaceRank < 0 ) throw Exception( __LINE__ );
    rowsInTable = dataSetCount[0];

    for( int i = 0; i < dataSpaceRank; ++i )
        dataSetOffset[i] = 0;

    currOffset = rangeEnd = internalRange.end();
}

unsigned ReadHDF5Dataset::columns() const
{
    if( dataSpaceRank == 1 )
        return 1;
    else if( dataSpaceRank == 2 )
        return dataSetCount[1];

    throw Exception( __LINE__ );
}

void ReadHDF5Dataset::set_column( unsigned column )
{
    if( dataSpaceRank != 2 || column >= dataSetCount[1] ) throw Exception( __LINE__ );
    dataSetCount[1]  = 1;
    dataSetOffset[1] = column;
}

Range::const_iterator ReadHDF5Dataset::next_end( Range::const_iterator iter )
{
    size_t slabs_remaining = hyperslabSelectionLimit;
    size_t avail           = bufferSize;
    while( iter != rangeEnd && slabs_remaining )
    {
        size_t count = *( iter.end_of_block() ) - *iter + 1;
        if( count >= avail )
        {
            iter += avail;
            break;
        }

        avail -= count;
        iter += count;
        --slabs_remaining;
    }
    return iter;
}

void ReadHDF5Dataset::set_file_ids( const Range& file_ids, EntityHandle start_id, hsize_t row_count, hid_t data_type )
{
    startID    = start_id;
    currOffset = file_ids.begin();
    rangeEnd   = file_ids.end();
    readCount  = 0;
    bufferSize = row_count;

    // if a) user specified buffer size and b) we're doing a true
    // parallel partial read and c) we're doing collective I/O, then
    // we need to know the maximum number of reads that will be done.
#ifdef MOAB_HAVE_HDF5_PARALLEL
    if( nativeParallel )
    {
        Range::const_iterator iter = currOffset;
        while( iter != rangeEnd )
        {
            ++readCount;
            iter = next_end( iter );
        }

        MPE_Log_event( mpeReduceEvent.first, (int)readCount, mpeDesc.c_str() );
        unsigned long recv = readCount, send = readCount;
        MPI_Allreduce( &send, &recv, 1, MPI_UNSIGNED_LONG, MPI_MAX, *mpiComm );
        readCount = recv;
        MPE_Log_event( mpeReduceEvent.second, (int)readCount, mpeDesc.c_str() );
    }
#endif

    dataType     = data_type;
    htri_t equal = H5Tequal( fileType, dataType );
    if( equal < 0 ) throw Exception( __LINE__ );
    doConversion = !equal;

    // We always read in the format of the file to avoid stupind HDF5
    // library behavior when reading in parallel.  We call H5Tconvert
    // ourselves to do the data conversion.  If the type we're reading
    // from the file is larger than the type we want in memory, then
    // we need to reduce num_rows so that we can read the larger type
    // from the file into the passed buffer mean to accomodate num_rows
    // of values of the smaller in-memory type.
    if( doConversion )
    {
        size_t mem_size, file_size;
        mem_size  = H5Tget_size( dataType );
        file_size = H5Tget_size( fileType );
        if( file_size > mem_size ) bufferSize = bufferSize * mem_size / file_size;
    }
}

void ReadHDF5Dataset::set_all_file_ids( hsize_t row_count, hid_t data_type )
{
    internalRange.clear();
    internalRange.insert( (EntityHandle)1, (EntityHandle)( rowsInTable ) );
    set_file_ids( internalRange, 1, row_count, data_type );
}

ReadHDF5Dataset::~ReadHDF5Dataset()
{
    if( fileType >= 0 ) H5Tclose( fileType );
    if( dataSpace >= 0 ) H5Sclose( dataSpace );
    if( closeDataSet && dataSet >= 0 ) H5Dclose( dataSet );
    dataSpace = dataSet = -1;
    if( ioProp != H5P_DEFAULT ) H5Pclose( ioProp );
}

void ReadHDF5Dataset::read( void* buffer, size_t& rows_read )
{
    herr_t err;
    rows_read = 0;

    MPE_Log_event( mpeReadEvent.first, (int)readCount, mpeDesc.c_str() );
    if( currOffset != rangeEnd )
    {

        // Build H5S hyperslab selection describing the portions of the
        // data set to read
        H5S_seloper_t sop       = H5S_SELECT_SET;
        Range::iterator new_end = next_end( currOffset );
        while( currOffset != new_end )
        {
            size_t count = *( currOffset.end_of_block() ) - *currOffset + 1;
            if( new_end != rangeEnd && *currOffset + count > *new_end )
            {
                count = *new_end - *currOffset;
            }
            rows_read += count;

            dataSetOffset[0] = *currOffset - startID;
            dataSetCount[0]  = count;
            err              = H5Sselect_hyperslab( dataSpace, sop, dataSetOffset, NULL, dataSetCount, 0 );
            if( err < 0 ) throw Exception( __LINE__ );
            sop = hyperslabSelectOp;  // subsequent calls to select_hyperslab append

            currOffset += count;
        }

        // Create a data space describing the memory in which to read the data
        dataSetCount[0] = rows_read;
        hid_t mem_id    = H5Screate_simple( dataSpaceRank, dataSetCount, NULL );
        if( mem_id < 0 ) throw Exception( __LINE__ );

        // Do the actual read
        err = H5Dread( dataSet, fileType, mem_id, dataSpace, ioProp, buffer );
        H5Sclose( mem_id );
        if( err < 0 ) throw Exception( __LINE__ );

        if( readCount ) --readCount;

        if( doConversion )
        {
            err = H5Tconvert( fileType, dataType, rows_read * columns(), buffer, 0, H5P_DEFAULT );
            if( err < 0 ) throw Exception( __LINE__ );
        }
    }
    else if( readCount )
    {
        null_read();
        --readCount;
    }
    MPE_Log_event( mpeReadEvent.second, (int)readCount, mpeDesc.c_str() );
}

void ReadHDF5Dataset::null_read()
{
    herr_t err;
    err = H5Sselect_none( dataSpace );
    if( err < 0 ) throw Exception( __LINE__ );

    //#if HDF5_16API
    hsize_t one  = 1;
    hid_t mem_id = H5Screate_simple( 1, &one, NULL );
    if( mem_id < 0 ) throw Exception( __LINE__ );
    err = H5Sselect_none( mem_id );
    if( err < 0 )
    {
        H5Sclose( mem_id );
        throw Exception( __LINE__ );
    }
    //#else
    //  hid_t mem_id = H5Screate(H5S_NULL);
    //  if (mem_id < 0)
    //    throw Exception(__LINE__);
    //#endif

    err = H5Dread( dataSet, fileType, mem_id, dataSpace, ioProp, 0 );
    H5Sclose( mem_id );
    if( err < 0 ) throw Exception( __LINE__ );
}

}  // namespace moab