From 96ef15c5c6cdf2649e14bc742b521cc825351df9 Mon Sep 17 00:00:00 2001 From: Mark Berrill Date: Wed, 2 Jun 2021 15:36:44 -0400 Subject: [PATCH] Creating HDF5 writer using Xmdf for visualization --- IO/HDF5Writer.cpp | 283 +++++++++++++++++ IO/HDF5_IO.cpp | 585 +++++++++++++++++++++++++++++++++++ IO/HDF5_IO.h | 169 +++++++++++ IO/HDF5_IO.hpp | 348 +++++++++++++++++++++ IO/Mesh.cpp | 6 +- IO/Mesh.h | 2 +- IO/MeshDatabase.cpp | 2 +- IO/Reader.cpp | 142 +++++++-- IO/SiloWriter.cpp | 251 +++++++++++++++ IO/Writer.cpp | 297 +++--------------- IO/Writer.h | 11 +- IO/Xdmf.cpp | 620 ++++++++++++++++++++++++++++++++++++++ IO/Xdmf.h | 130 ++++++++ IO/netcdf.cpp | 2 +- IO/netcdf.h | 4 +- IO/silo.cpp | 6 +- IO/silo.hpp | 2 +- StackTrace/StackTrace.cpp | 2 +- common/Array.cpp | 143 +++------ common/Array.h | 232 ++++++++++---- common/Array.hpp | 600 +++++++++++++++++++----------------- common/ArraySize.h | 233 +++++++++++--- common/FunctionTable.cpp | 147 +++++++++ common/FunctionTable.h | 144 +++++++-- common/FunctionTable.hpp | 239 ++++++++++----- tests/TestWriter.cpp | 24 +- 26 files changed, 3751 insertions(+), 873 deletions(-) create mode 100644 IO/HDF5Writer.cpp create mode 100644 IO/HDF5_IO.cpp create mode 100644 IO/HDF5_IO.h create mode 100644 IO/HDF5_IO.hpp create mode 100644 IO/SiloWriter.cpp create mode 100644 IO/Xdmf.cpp create mode 100644 IO/Xdmf.h create mode 100644 common/FunctionTable.cpp diff --git a/IO/HDF5Writer.cpp b/IO/HDF5Writer.cpp new file mode 100644 index 00000000..3fca33bd --- /dev/null +++ b/IO/HDF5Writer.cpp @@ -0,0 +1,283 @@ +#include "IO/HDF5_IO.h" +#include "IO/IOHelpers.h" +#include "IO/MeshDatabase.h" +#include "IO/Writer.h" +#include "IO/Xdmf.h" +#include "common/MPI.h" +#include "common/Utilities.h" + +#include +#include +#include +#include +#include + + +#ifdef USE_HDF5 + + +std::string to_string( const ArraySize &s ) +{ + std::string out = "[" + std::to_string( s[0] ); + for ( size_t i = 1; i < s.ndim(); i++ ) + out += "," + to_string( s[i] ); + out += "]"; + return out; +} + + +Xdmf::Center getXdmfType( IO::VariableType type ) +{ + if ( type == IO::VariableType::NodeVariable ) { + return Xdmf::Center::Node; + } else if ( type == IO::VariableType::VolumeVariable ) { + return Xdmf::Center::Cell; + } else { + ERROR( "Variable type not supported" ); + } + return Xdmf::Center::Null; +} + + +// Write a PointList mesh (and variables) to a file +template +static void writeCoordinates( hid_t fid, const std::vector &points ) +{ + std::vector x( points.size() ), y( points.size() ), z( points.size() ); + for ( size_t i = 0; i < x.size(); i++ ) { + x[i] = points[i].x; + y[i] = points[i].y; + z[i] = points[i].z; + } + IO::HDF5::writeHDF5( fid, "x", x ); + IO::HDF5::writeHDF5( fid, "y", y ); + IO::HDF5::writeHDF5( fid, "z", z ); +} +static void writeHDF5PointList( hid_t fid, const std::string &filename, + const IO::MeshDataStruct &meshData, IO::MeshDatabase database, Xdmf &xmf ) +{ + auto meshname = database.domains[0].name; + const auto &mesh = dynamic_cast( *meshData.mesh ); + auto gid = IO::HDF5::createGroup( fid, meshname ); + if ( meshData.precision == IO::DataType::Double ) { + writeCoordinates( gid, mesh.getPoints() ); + } else if ( meshData.precision == IO::DataType::Float ) { + writeCoordinates( gid, mesh.getPoints() ); + } else { + ERROR( "Unsupported format" ); + } + auto path = filename + ":/" + meshname + "/"; + auto domain = Xdmf::createPointMesh( + meshname, 3, mesh.getPoints().size(), path + "x", path + "y", path + "z" ); + for ( size_t i = 0; i < meshData.vars.size(); i++ ) { + auto &var = *meshData.vars[i]; + auto data = var.data; + auto rankType = Xdmf::RankType::Null; + if ( data.ndim() == 1 ) { + rankType = Xdmf::RankType::Scalar; + } else if ( data.ndim() == 2 && data.size( 1 ) == 3 ) { + // Vector data, need to permute for visit + rankType = Xdmf::RankType::Vector; + data = data.permute( { 1, 0 } ); + } else { + ERROR( "Unable to determine variable rank: " + to_string( var.data.size() ) ); + } + if ( var.precision == IO::DataType::Double ) { + IO::HDF5::writeHDF5( gid, var.name, data ); + } else if ( var.precision == IO::DataType::Float ) { + IO::HDF5::writeHDF5( gid, var.name, data.cloneTo() ); + } else if ( var.precision == IO::DataType::Int ) { + IO::HDF5::writeHDF5( gid, var.name, data.cloneTo() ); + } else { + ERROR( "Unsupported format" ); + } + domain.addVariable( + meshname, var.name, data.size(), rankType, Xdmf::Center::Node, path + var.name ); + } + xmf.addMesh( meshData.meshName, domain ); +} +// Write a TriMesh mesh (and variables) to a file +static void writeHDF5TriMesh2( hid_t fid, const std::string &filename, + const IO::MeshDataStruct &meshData, const IO::TriMesh &mesh, IO::MeshDatabase database, + Xdmf &xmf ) +{ + auto meshname = database.domains[0].name; + auto gid = IO::HDF5::createGroup( fid, meshname ); + auto path = filename + ":/" + meshname + "/"; + // Write the verticies + if ( meshData.precision == IO::DataType::Double ) { + writeCoordinates( gid, mesh.vertices->getPoints() ); + } else if ( meshData.precision == IO::DataType::Float ) { + writeCoordinates( gid, mesh.vertices->getPoints() ); + } else { + ERROR( "Unsupported format" ); + } + // Write the connectivity + Array tri( 3, mesh.A.size() ); + for ( size_t i = 0; i < mesh.A.size(); i++ ) { + tri( 0, i ) = mesh.A[i]; + tri( 1, i ) = mesh.B[i]; + tri( 2, i ) = mesh.C[i]; + } + IO::HDF5::writeHDF5( gid, "tri", tri ); + auto domain = + Xdmf::createUnstructuredMesh( meshname, 3, Xdmf::TopologyType::Triangle, tri.size( 1 ), + path + "tri", mesh.vertices->getPoints().size(), path + "x", path + "y", path + "z" ); + // Write the variables + for ( size_t i = 0; i < meshData.vars.size(); i++ ) { + auto &var = *meshData.vars[i]; + auto data = var.data; + auto rankType = Xdmf::RankType::Null; + if ( data.ndim() == 1 ) { + rankType = Xdmf::RankType::Scalar; + } else if ( data.ndim() == 2 && data.size( 1 ) == 3 ) { + // Vector data, need to permute for visit + rankType = Xdmf::RankType::Vector; + data = data.permute( { 1, 0 } ); + } else { + ERROR( "Unable to determine variable rank: " + to_string( var.data.size() ) ); + } + if ( var.precision == IO::DataType::Double ) { + IO::HDF5::writeHDF5( gid, var.name, data ); + } else if ( var.precision == IO::DataType::Float ) { + IO::HDF5::writeHDF5( gid, var.name, data.cloneTo() ); + } else if ( var.precision == IO::DataType::Int ) { + IO::HDF5::writeHDF5( gid, var.name, data.cloneTo() ); + } else { + ERROR( "Unsupported format" ); + } + domain.addVariable( + meshname, var.name, data.size(), rankType, getXdmfType( var.type ), path + var.name ); + } + xmf.addMesh( meshData.meshName, domain ); +} +static void writeHDF5TriMesh( hid_t fid, const std::string &filename, + const IO::MeshDataStruct &meshData, IO::MeshDatabase database, Xdmf &xmf ) +{ + const IO::TriMesh &mesh = dynamic_cast( *meshData.mesh ); + writeHDF5TriMesh2( fid, filename, meshData, mesh, database, xmf ); +} +static void writeHDF5TriList( hid_t fid, const std::string &filename, + const IO::MeshDataStruct &meshData, IO::MeshDatabase database, Xdmf &xmf ) +{ + auto mesh = getTriMesh( meshData.mesh ); + writeHDF5TriMesh2( fid, filename, meshData, *mesh, database, xmf ); +} +// Write a DomainMesh mesh (and variables) to a file +static void writeHDF5DomainMesh( hid_t fid, const std::string &filename, + const IO::MeshDataStruct &meshData, IO::MeshDatabase database, Xdmf &xmf ) +{ + auto &mesh = dynamic_cast( *meshData.mesh ); + auto meshname = database.domains[0].name; + auto gid = IO::HDF5::createGroup( fid, meshname ); + auto path = filename + ":/" + meshname + "/"; + // Write the mesh + RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); + std::vector range = { info.ix * mesh.Lx / info.nx, ( info.ix + 1 ) * mesh.Lx / info.nx, + info.jy * mesh.Ly / info.ny, ( info.jy + 1 ) * mesh.Ly / info.ny, + info.kz * mesh.Lz / info.nz, ( info.kz + 1 ) * mesh.Lz / info.nz }; + std::vector N = { mesh.nx, mesh.ny, mesh.nz }; + std::vector rankinfo = { mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz }; + IO::HDF5::writeHDF5( gid, "range", range ); + IO::HDF5::writeHDF5( gid, "N", N ); + IO::HDF5::writeHDF5( gid, "rankinfo", rankinfo ); + // xmf.addUniformMesh( meshname, range, ArraySize( N[0], N[1], N[2] ) ); + // Write a curvilinear mesh due to bug with vector data on nodes loading into visit + Array x( N[0] + 1, N[1] + 1, N[2] + 1 ); + Array y( N[0] + 1, N[1] + 1, N[2] + 1 ); + Array z( N[0] + 1, N[1] + 1, N[2] + 1 ); + double dx = ( range[1] - range[0] ) / N[0]; + double dy = ( range[3] - range[2] ) / N[1]; + double dz = ( range[5] - range[4] ) / N[2]; + for ( int k = 0; k <= N[2]; k++ ) { + for ( int j = 0; j <= N[1]; j++ ) { + for ( int i = 0; i <= N[0]; i++ ) { + x( i, j, k ) = range[0] + dx * i; + y( i, j, k ) = range[2] + dy * j; + z( i, j, k ) = range[4] + dz * k; + } + } + } + IO::HDF5::writeHDF5( gid, "x", x ); + IO::HDF5::writeHDF5( gid, "y", y ); + IO::HDF5::writeHDF5( gid, "z", z ); + auto domain = Xdmf::createCurvilinearMesh( + meshname, ArraySize( N[0], N[1], N[2] ), path + "x", path + "y", path + "z" ); + // Write the variables + for ( size_t i = 0; i < meshData.vars.size(); i++ ) { + auto &var = *meshData.vars[i]; + auto data = var.data; + auto rankType = Xdmf::RankType::Null; + if ( data.ndim() == 3 ) { + rankType = Xdmf::RankType::Scalar; + } else if ( data.ndim() == 4 && data.size( 3 ) == 3 ) { + // Vector data, need to permute for visit + rankType = Xdmf::RankType::Vector; + data = data.permute( { 3, 0, 1, 2 } ); + } else { + ERROR( "Unable to determine variable rank: " + to_string( var.data.size() ) ); + } + if ( var.precision == IO::DataType::Double ) { + IO::HDF5::writeHDF5( gid, var.name, data ); + } else if ( var.precision == IO::DataType::Float ) { + IO::HDF5::writeHDF5( gid, var.name, data.cloneTo() ); + } else if ( var.precision == IO::DataType::Int ) { + IO::HDF5::writeHDF5( gid, var.name, data.cloneTo() ); + } else { + ERROR( "Unsupported format" ); + } + domain.addVariable( + meshname, var.name, data.size(), rankType, getXdmfType( var.type ), path + var.name ); + } + IO::HDF5::closeGroup( gid ); + xmf.addMesh( meshData.meshName, domain ); +} +// Write a mesh (and variables) to a file +static IO::MeshDatabase write_domain_hdf5( hid_t fid, const std::string &filename, + const IO::MeshDataStruct &mesh, IO::FileFormat format, int rank, Xdmf &xmf ) +{ + // Create the MeshDatabase + auto database = getDatabase( filename, mesh, format, rank ); + if ( database.meshClass == "PointList" ) { + writeHDF5PointList( fid, filename, mesh, database, xmf ); + } else if ( database.meshClass == "TriMesh" ) { + writeHDF5TriMesh( fid, filename, mesh, database, xmf ); + } else if ( database.meshClass == "TriList" ) { + writeHDF5TriList( fid, filename, mesh, database, xmf ); + } else if ( database.meshClass == "DomainMesh" ) { + writeHDF5DomainMesh( fid, filename, mesh, database, xmf ); + } else { + ERROR( "Unknown mesh class" ); + } + return database; +} +// Write the mesh data to hdf5 +std::vector writeMeshesHDF5( const std::vector &meshData, + const std::string &path, IO::FileFormat format, int rank, Xdmf &xmf ) +{ + + std::vector meshes_written; + char filename[100], fullpath[200]; + sprintf( filename, "%05i.h5", rank ); + sprintf( fullpath, "%s/%s", path.c_str(), filename ); + auto fid = IO::HDF5::openHDF5( fullpath, "w", IO::HDF5::Compression::GZIP ); + for ( size_t i = 0; i < meshData.size(); i++ ) { + meshes_written.push_back( + write_domain_hdf5( fid, filename, meshData[i], format, rank, xmf ) ); + } + IO::HDF5::closeHDF5( fid ); + return meshes_written; +} + + +#else + + +std::vector writeMeshesHDF5( + const std::vector &, const std::string &, IO::FileFormat, int ); +{ + return std::vector(); +} + + +#endif diff --git a/IO/HDF5_IO.cpp b/IO/HDF5_IO.cpp new file mode 100644 index 00000000..c659f845 --- /dev/null +++ b/IO/HDF5_IO.cpp @@ -0,0 +1,585 @@ +#include "IO/HDF5_IO.h" +#include "IO/HDF5_IO.hpp" +#include "common/Array.h" +#include "common/Utilities.h" + +#include +#include +#include +#include + + +namespace IO { +namespace HDF5 { + + +#ifdef USE_HDF5 // USE HDF5 + + +/************************************************************************ + * HDF5 helper routines * + ************************************************************************/ +inline const void *H5Ptr( const void *x ) { return x == nullptr ? ( (void *) 1 ) : x; } +bool H5Gexists( hid_t fid, const std::string &name ) +{ + H5E_auto2_t func; + void *client; + H5Eget_auto2( H5E_DEFAULT, &func, &client ); + H5Eset_auto2( H5E_DEFAULT, nullptr, nullptr ); + int status = H5Gget_objinfo( fid, name.data(), 0, nullptr ); + H5Eset_auto2( H5E_DEFAULT, func, client ); + return status == 0; +} +bool H5Dexists( hid_t fid, const std::string &name ) +{ + H5E_auto2_t func; + void *client; + H5Eget_auto2( H5E_DEFAULT, &func, &client ); + H5Eset_auto2( H5E_DEFAULT, nullptr, nullptr ); + hid_t dataset = H5Dopen2( fid, name.data(), H5P_DEFAULT ); + H5Eset_auto2( H5E_DEFAULT, func, client ); + bool exists = dataset > 0; + // if ( exists ) + // H5Dclose( dataset ); + return exists; +} +hid_t createGroup( hid_t fid, const std::string &name ) +{ + return H5Gcreate2( fid, name.data(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT ); +} +hid_t openGroup( hid_t fid, const std::string &name ) +{ + INSIST( H5Gexists( fid, name ), "Group " + name + " does not exist" ); + return H5Gopen2( fid, name.data(), H5P_DEFAULT ); +} +void closeGroup( hid_t gid ) { H5Gclose( gid ); } + + +/************************************************************************ + * Complex struct that is compatible with HDF5 * + ************************************************************************/ +typedef struct { + double re; + double im; +} complex_t; +inline void convert( size_t N, const std::complex *x, complex_t *y ) +{ + for ( size_t i = 0; i < N; i++ ) { + y[i].re = x[i].real(); + y[i].im = x[i].imag(); + } +} +inline void convert( size_t N, const complex_t *x, std::complex *y ) +{ + for ( size_t i = 0; i < N; i++ ) { + y[i] = std::complex( x[i].re, x[i].im ); + } +} + + +/************************************************************************ + * Get the HDF5 data type * + ************************************************************************/ +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_UCHAR ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_CHAR ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_UINT8 ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_INT8 ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_UINT16 ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_INT16 ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_INT ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_UINT ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_LONG ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_ULONG ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_FLOAT ); +} +template<> +hid_t getHDF5datatype() +{ + return H5Tcopy( H5T_NATIVE_DOUBLE ); +} +template<> +hid_t getHDF5datatype>() +{ + hid_t datatype = H5Tcreate( H5T_COMPOUND, sizeof( complex_t ) ); + H5Tinsert( datatype, "real", HOFFSET( complex_t, re ), H5T_NATIVE_DOUBLE ); + H5Tinsert( datatype, "imag", HOFFSET( complex_t, im ), H5T_NATIVE_DOUBLE ); + return datatype; +} +template<> +hid_t getHDF5datatype() +{ + hid_t datatype = H5Tcopy( H5T_C_S1 ); + H5Tset_size( datatype, H5T_VARIABLE ); + return datatype; +} + + +/************************************************************************ + * Read/write Array types * + ************************************************************************/ +template<> +void readHDF5>( hid_t fid, const std::string &name, Array &data ) +{ + if ( !H5Dexists( fid, name ) ) { + // Dataset does not exist + data.resize( 0 ); + return; + } + hid_t dataset = H5Dopen2( fid, name.data(), H5P_DEFAULT ); + hid_t datatype = H5Dget_type( dataset ); + hid_t dataspace = H5Dget_space( dataset ); + hsize_t dims0[10]; + int ndim = H5Sget_simple_extent_dims( dataspace, dims0, nullptr ); + auto dims = convertSize( ndim, dims0 ); + data.resize( dims ); + hid_t datatype2 = getHDF5datatype(); + if ( data.empty() ) { + // The data is empty + } else if ( H5Tequal( datatype, datatype2 ) ) { + // The type of Array and the data in HDF5 match + auto **tmp = new char *[data.length() * sizeof( char * )]; + memset( tmp, 0, data.length() * sizeof( char * ) ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, tmp ); + for ( size_t i = 0; i < data.length(); i++ ) + data( i ) = std::string( tmp[i] ); + H5Dvlen_reclaim( datatype, dataspace, H5P_DEFAULT, tmp ); + delete[] tmp; + } else { + ERROR( "Unknown format for std::string" ); + } + H5Dclose( dataset ); + H5Tclose( datatype ); + H5Tclose( datatype2 ); + H5Sclose( dataspace ); +} +template<> +void readHDF5>>( + hid_t fid, const std::string &name, Array> &data ) +{ + if ( !H5Dexists( fid, name ) ) { + // Dataset does not exist + data.resize( 0 ); + return; + } + hid_t dataset = H5Dopen2( fid, name.data(), H5P_DEFAULT ); + hid_t datatype = H5Dget_type( dataset ); + hid_t dataspace = H5Dget_space( dataset ); + hsize_t dims0[10]; + int ndim = H5Sget_simple_extent_dims( dataspace, dims0, nullptr ); + auto dims = convertSize( ndim, dims0 ); + data.resize( dims ); + hid_t datatype2 = getHDF5datatype>(); + if ( data.empty() ) { + // The data is empty + } else if ( H5Tequal( datatype, datatype2 ) ) { + // The type of Array and the data in HDF5 match + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data.data() ); + } else { + ERROR( "We need to convert data formats" ); + } + H5Dclose( dataset ); + H5Tclose( datatype ); + H5Tclose( datatype2 ); + H5Sclose( dataspace ); +} +// clang-format off +#define readWriteHDF5Array( TYPE ) \ + template<> \ + void writeHDF5>( hid_t fid, const std::string &name, const Array &data ) \ + { \ + writeHDF5ArrayDefault( fid, name, data ); \ + } \ + template<> \ + void readHDF5>( hid_t fid, const std::string &name, Array &data ) \ + { \ + readHDF5ArrayDefault( fid, name, data ); \ + } +readWriteHDF5Array( bool ) +readWriteHDF5Array( char ) +readWriteHDF5Array( int8_t ) +readWriteHDF5Array( int16_t ) +readWriteHDF5Array( int32_t ) +readWriteHDF5Array( int64_t ) +readWriteHDF5Array( uint8_t ) +readWriteHDF5Array( uint16_t ) +readWriteHDF5Array( uint32_t ) +readWriteHDF5Array( uint64_t ) +readWriteHDF5Array( float ) +readWriteHDF5Array( double ) + // clang-format on + + + /************************************************************************ + * Read/write scalar types * + ************************************************************************/ + template<> + void readHDF5( hid_t fid, const std::string &name, std::string &data ) +{ + hid_t dataset = H5Dopen2( fid, name.data(), H5P_DEFAULT ); + hid_t datatype = H5Dget_type( dataset ); + hid_t datatype0 = getHDF5datatype(); + if ( H5Tequal( datatype, datatype0 ) ) { + hid_t dataspace = H5Dget_space( dataset ); + char *tmp[1] = { nullptr }; + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, tmp ); + data = std::string( tmp[0] ); + H5Dvlen_reclaim( datatype, dataspace, H5P_DEFAULT, tmp ); + H5Sclose( dataspace ); + } else { + Array tmp; + readHDF5( fid, name, tmp ); + data = std::string( tmp.data(), tmp.length() ); + } + H5Dclose( dataset ); + H5Tclose( datatype ); + H5Tclose( datatype0 ); +} +template<> +void writeHDF5( hid_t fid, const std::string &name, const std::string &data ) +{ + Array tmp; + tmp.viewRaw( { data.length() }, (char *) data.data() ); + writeHDF5( fid, name, tmp ); +} +// clang-format off +#define readWriteHDF5Scalar( TYPE ) \ + template<> \ + void writeHDF5( hid_t fid, const std::string &name, const TYPE &data ) \ + { \ + hid_t dataspace = H5Screate( H5S_SCALAR ); \ + hid_t datatype = getHDF5datatype(); \ + hid_t dataset = H5Dcreate2( \ + fid, name.data(), datatype, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT ); \ + H5Dwrite( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, H5Ptr( &data ) ); \ + H5Dclose( dataset ); \ + H5Tclose( datatype ); \ + H5Sclose( dataspace ); \ + } \ + template<> \ + void readHDF5( hid_t fid, const std::string &name, TYPE &data ) \ + { \ + Array tmp; \ + readHDF5( fid, name, tmp ); \ + INSIST( tmp.ndim() == 1 && tmp.length() == 1, "Error loading " + std::string( name ) ); \ + data = tmp( 0 ); \ + } +readWriteHDF5Scalar( bool ) +readWriteHDF5Scalar( char ) +readWriteHDF5Scalar( int8_t ) +readWriteHDF5Scalar( int16_t ) +readWriteHDF5Scalar( int32_t ) +readWriteHDF5Scalar( int64_t ) +readWriteHDF5Scalar( uint8_t ) +readWriteHDF5Scalar( uint16_t ) +readWriteHDF5Scalar( uint32_t ) +readWriteHDF5Scalar( uint64_t ) +readWriteHDF5Scalar( float ) +readWriteHDF5Scalar( double ) +readWriteHDF5Scalar( std::complex ) + // clang-format on + + + /****************************************************************** + * Create custom error handler * + ******************************************************************/ + herr_t hdf5_error_handler( hid_t err_stack, void * ) +{ + FILE *fid = tmpfile(); + H5Eprint2( err_stack, fid ); + H5Eclear2( err_stack ); + rewind( fid ); + char msg[1024]; + size_t N = fread( msg, 1, sizeof( msg ) - 1, fid ); + fclose( fid ); + msg[N] = 0; + std::string msg2 = "Error calling HDF5 routine:\n"; + ERROR( msg2 + msg ); + return 0; +} +bool set_hdf5_error_handler() +{ + hid_t error_stack = 0; + H5E_auto2_t fun = hdf5_error_handler; + H5Eset_auto2( error_stack, fun, nullptr ); + return true; +} +bool global_is_hdf5_error_handler_set = set_hdf5_error_handler(); + + +/****************************************************************** + * Open/close HDF5 files * + ******************************************************************/ +hid_t openHDF5( const std::string &filename, const char *mode, Compression compress ) +{ + // Set cache size to 3MBs and instruct the cache to discard the fully read chunk + auto pid = H5P_DEFAULT; + /*auto pid = H5Pcreate( H5P_FILE_ACCESS ); + int nelemts; + size_t nslots, nbytes; + double w0; + H5Pget_cache(pid,& nelemts,& nslots,& nbytes,& w0); + H5Pset_cache(pid, nelemts, 1999, 3*1024*1024, 1.0); */ + // Open the file + hid_t fid = 0; + if ( strcmp( mode, "r" ) == 0 ) { + fid = H5Fopen( filename.data(), H5F_ACC_RDONLY, pid ); + } else if ( strcmp( mode, "w" ) == 0 ) { + fid = H5Fcreate( filename.data(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT ); + } else if ( strcmp( mode, "rw" ) == 0 ) { + fid = H5Fopen( filename.data(), H5F_ACC_RDWR, H5P_DEFAULT ); + } else { + ERROR( "Invalid mode for opening HDF5 file" ); + } + if ( strcmp( mode, "w" ) == 0 ) { + if ( compress == Compression::None ) { + writeHDF5( fid, "DefaultCompression", 0 ); + } else if ( compress == Compression::GZIP ) { + writeHDF5( fid, "DefaultCompression", 1 ); + } else if ( compress == Compression::SZIP ) { + writeHDF5( fid, "DefaultCompression", 2 ); + } else { + ERROR( "Internal error" ); + } + } + // H5Pclose( pid ); + return fid; +} +void closeHDF5( hid_t fid ) +{ + // Try to close any remaining objects (needed to ensure we can reopen the data if desired) + hid_t file[1000], set[1000], group[1000], type[1000], attr[1000]; + size_t N_file = H5Fget_obj_ids( fid, H5F_OBJ_FILE, 1000, file ); + size_t N_set = H5Fget_obj_ids( fid, H5F_OBJ_DATASET, 1000, set ); + size_t N_group = H5Fget_obj_ids( fid, H5F_OBJ_GROUP, 1000, group ); + size_t N_type = H5Fget_obj_ids( fid, H5F_OBJ_DATATYPE, 1000, type ); + size_t N_attr = H5Fget_obj_ids( fid, H5F_OBJ_ATTR, 1000, attr ); + for ( size_t i = 0; i < N_file; i++ ) { + if ( file[i] != fid ) + H5Fclose( file[i] ); + } + for ( size_t i = 0; i < N_set; i++ ) + H5Dclose( set[i] ); + for ( size_t i = 0; i < N_group; i++ ) + H5Gclose( group[i] ); + for ( size_t i = 0; i < N_type; i++ ) + H5Tclose( type[i] ); + for ( size_t i = 0; i < N_attr; i++ ) + H5Aclose( attr[i] ); + // Flush the data (needed to ensure we can reopen the data if desired) + unsigned intent; + H5Fget_intent( fid, &intent ); + if ( intent == H5F_ACC_RDWR || intent == H5F_ACC_TRUNC ) + H5Fflush( fid, H5F_SCOPE_GLOBAL ); + // Close the file + H5Fclose( fid ); +} + + +/************************************************************************ + * Check if we support compression * + ************************************************************************/ +Compression defaultCompression( hid_t fid ) +{ + hid_t root = H5Gopen2( fid, "/", H5P_DEFAULT ); + if ( !H5Dexists( root, "DefaultCompression" ) ) + return Compression::None; + int tmp; + readHDF5( root, "DefaultCompression", tmp ); + Compression compress = Compression::None; + if ( tmp == 0 ) { + compress = Compression::None; + } else if ( tmp == 1 ) { + compress = Compression::GZIP; + } else if ( tmp == 2 ) { + compress = Compression::SZIP; + } else { + ERROR( "Internal error" ); + } + return compress; +} + + +/************************************************************************ + * Create a default chunk size * + ************************************************************************/ +hid_t createChunk( const std::vector &dims, Compression compress ) +{ + if ( compress == Compression::None || dims.empty() ) + return H5P_DEFAULT; + hsize_t length = 1; + for ( auto d : dims ) + length *= d; + if ( length < 512 ) + return H5P_DEFAULT; + hid_t plist = H5Pcreate( H5P_DATASET_CREATE ); + auto status = H5Pset_chunk( plist, dims.size(), dims.data() ); + ASSERT( status == 0 ); + if ( compress == Compression::GZIP ) { + status = H5Pset_deflate( plist, 7 ); + ASSERT( status == 0 ); + } else if ( compress == Compression::SZIP ) { + status = H5Pset_szip( plist, H5_SZIP_NN_OPTION_MASK, 16 ); + ASSERT( status == 0 ); + } + return plist; +} + + +/************************************************************************ + * Write Array * + ************************************************************************/ +template<> +void writeHDF5>>( + hid_t fid, const std::string &name, const Array> &data ) +{ + hid_t datatype = getHDF5datatype>(); + // Copy the data + size_t N = data.length(); + auto *y = new complex_t[N]; + convert( N, data.data(), y ); + // Save the array + auto dim = arraySize( data ); + hid_t dataspace = H5Screate_simple( dim.size(), dim.data(), nullptr ); + hid_t dataset = + H5Dcreate2( fid, name.data(), datatype, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT ); + H5Dwrite( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, H5Ptr( y ) ); + H5Dclose( dataset ); + H5Tclose( datatype ); + H5Sclose( dataspace ); + delete[] y; +} +template<> +void writeHDF5>( + hid_t fid, const std::string &name, const Array &data ) +{ + auto dim = arraySize( data ); + hid_t dataspace = H5Screate_simple( dim.size(), dim.data(), nullptr ); + auto **tmp = new char *[data.length() + 1]; + memset( tmp, 0, ( data.length() + 1 ) * sizeof( char * ) ); + for ( size_t i = 0; i < data.length(); i++ ) { + tmp[i] = const_cast( data( i ).data() ); + } + hid_t datatype = getHDF5datatype(); + hid_t props = H5Pcreate( H5P_DATASET_CREATE ); + hid_t dataset = H5Dcreate1( fid, name.data(), datatype, dataspace, props ); + H5Dwrite( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, tmp ); + H5Pclose( props ); + H5Dclose( dataset ); + H5Tclose( datatype ); + H5Sclose( dataspace ); + delete[] tmp; +} + + +/************************************************************************ + * Specializations for std::vector * + ************************************************************************/ +template<> +void readHDF5>( hid_t fid, const std::string &name, std::vector &data ) +{ + Array tmp; + readHDF5( fid, name, tmp ); + data.resize( tmp.length() ); + for ( size_t i = 0; i < data.size(); i++ ) + data[i] = tmp( i ); +} +template<> +void writeHDF5>( hid_t fid, const std::string &name, const std::vector &x ) +{ + Array y( x.size() ); + for ( size_t i = 0; i < x.size(); i++ ) + y( i ) = x[i]; + writeHDF5( fid, name, y ); +} + + +/************************************************************************ + * Explicit instantiations for std::vector * + ***********************************************************************/ +// clang-format off +#define INSTANTIATE_STD_VECTOR( TYPE ) \ + template<> void readHDF5>( hid_t fid, const std::string &name, std::vector &x ) \ + { \ + Array y; \ + readHDF5( fid, name, y ); \ + x.resize( y.length() ); \ + for ( size_t i = 0; i < x.size(); i++ ) \ + x[i] = y( i ); \ + } \ + template<> void writeHDF5>( hid_t fid, const std::string &name, const std::vector &x ) \ + { \ + Array y; \ + y.viewRaw( { x.size() }, const_cast( x.data() ) ); \ + writeHDF5( fid, name, y ); \ + } +INSTANTIATE_STD_VECTOR( char ) +INSTANTIATE_STD_VECTOR( unsigned char ) +INSTANTIATE_STD_VECTOR( int ) +INSTANTIATE_STD_VECTOR( unsigned int ) +INSTANTIATE_STD_VECTOR( int16_t ) +INSTANTIATE_STD_VECTOR( uint16_t ) +INSTANTIATE_STD_VECTOR( int64_t ) +INSTANTIATE_STD_VECTOR( uint64_t ) +INSTANTIATE_STD_VECTOR( float ) +INSTANTIATE_STD_VECTOR( double ) +INSTANTIATE_STD_VECTOR( std::string ) +// clang-format on + + +#else // No HDF5 +// Dummy implimentations for no HDF5 +hid_t openHDF5( const std::string &, const char *, Compression ) { return 0; } +void closeHDF5( hid_t ) {} +bool H5Gexists( hid_t, const std::string & ) { return false; } +bool H5Dexists( hid_t, const std::string & ) { return false; } +hid_t createGroup( hid_t, const std::string & ) { return 0; } +hid_t openGroup( hid_t, const std::string & ) { return 0; } +void closeGroup( hid_t ) {} +#endif + + +} // namespace HDF5 +} // namespace IO diff --git a/IO/HDF5_IO.h b/IO/HDF5_IO.h new file mode 100644 index 00000000..7521a0e0 --- /dev/null +++ b/IO/HDF5_IO.h @@ -0,0 +1,169 @@ +// This file contains helper functions and interfaces for reading/writing HDF5 +#ifndef included_HDF5_h +#define included_HDF5_h + +#include "common/ArraySize.h" + +#include +#include + + +// Include the headers and define some basic types +#ifdef USE_HDF5 +// Using HDF5 +#include "hdf5.h" +#else +// Not using HDF5 +typedef int hid_t; +typedef size_t hsize_t; +#endif + + +namespace IO { +namespace HDF5 { + + +enum class Compression : uint8_t { None, GZIP, SZIP }; + + +/** + * \brief Open an HDF5 file + * \details This function opens and HDF5 file for reading/writing. + * Once complete, we must close the file using closeHDF5 + * @param[in] filename File to open + * @param[in] mode C string containing a file access mode. It can be: + * "r" read: Open file for input operations. The file must exist. + * "w" write: Create an empty file for output operations. + * If a file with the same name already exists, its contents + * are discarded and the file is treated as a new empty file. + * "rw" read+write: Open file for reading and writing. The file must exist. + * @param[in] compress Default compression + * @return Return a handle to the file. + */ +hid_t openHDF5( + const std::string &filename, const char *mode, Compression compress = Compression::None ); + + +/** + * \brief Open an HDF5 file + * \details This function opens and HDF5 file for reading/writing + * @param[in] fid File to open + */ +void closeHDF5( hid_t fid ); + + +/** + * \brief Retrun the the default compression + * \details This function returns the default compression used when the file was created + * @param[in] fid File/Group id + */ +Compression defaultCompression( hid_t fid ); + + +/** + * \brief Open an HDF5 file + * \details This function create a chunk for HDF5 + * @param[in] dims Chunk size + * @param[in] compress Compression to use + * @return Return a handle to the file. + */ +hid_t createChunk( const std::vector &dims, Compression compress ); + + +/** + * \brief Write a structure to HDF5 + * \details This function writes a C++ class/struct to HDF5. + * This is a templated function and users can impliment their own data + * types by creating explicit instantiations for a given type. + * There is no default instantiation except when compiled without HDF5 which is a no-op. + * @param[in] fid File or group to write to + * @param[in] name The name of the variable + * @param[in] data The structure to write + */ +template +void writeHDF5( hid_t fid, const std::string &name, const T &data ); + + +/** + * \brief Read a structure from HDF5 + * \details This function reads a C++ class/struct from HDF5. + * This is a templated function and users can impliment their own data + * types by creating explicit instantiations for a given type. + * There is no default instantiation except when compiled without HDF5 which is a no-op. + * @param[in] fid File or group to read from + * @param[in] name The name of the variable + * @param[out] data The structure to read + */ +template +void readHDF5( hid_t fid, const std::string &name, T &data ); + + +/** + * \brief Check if group exists + * \details This function checks if an HDF5 group exists in the file + * @param[in] fid ID of group or database to read + * @param[in] name The name of the group + */ +bool H5Gexists( hid_t fid, const std::string &name ); + + +/** + * \brief Check if dataset exists + * \details This function checks if an HDF5 dataset exists in the file + * @param[in] fid File to open + * @param[in] name The name of the dataset + */ +bool H5Dexists( hid_t fid, const std::string &name ); + + +/** + * \brief Create a group + * \details This function creates a new HDF5 group + * @param[in] fid File or group to write to + * @param[in] name The name of the group + */ +hid_t createGroup( hid_t fid, const std::string &name ); + + +/** + * \brief Open a group + * \details This function opens an HDF5 group + * @param[in] fid File or group to write to + * @param[in] name The name of the group + */ +hid_t openGroup( hid_t fid, const std::string &name ); + + +/** + * \brief Close a group + * \details This function closes an HDF5 group + * @param[in] fid Group to close + */ +void closeGroup( hid_t fid ); + + +/** + * \brief Get HDF5 data type + * \details This function returns the id of the data type + */ +template +hid_t getHDF5datatype(); + + +// Default no-op implimentations for use without HDF5 +// clang-format off +#ifndef USE_HDF5 +template void readHDF5( hid_t, const std::string&, T& ) {} +template void writeHDF5( hid_t, const std::string&, const T& ) {} +template void readHDF5Array( hid_t, const std::string&, Array& ) {} +template void writeHDF5Array( hid_t, const std::string&, const Array& ) {} +template hid_t getHDF5datatype() { return 0; } +#endif +// clang-format on + + +} // namespace HDF5 +} // namespace IO + + +#endif diff --git a/IO/HDF5_IO.hpp b/IO/HDF5_IO.hpp new file mode 100644 index 00000000..c9e01d90 --- /dev/null +++ b/IO/HDF5_IO.hpp @@ -0,0 +1,348 @@ +// This file contains helper functions and interfaces for reading/writing HDF5 +#ifndef included_HDF5_hpp +#define included_HDF5_hpp +#ifdef USE_HDF5 + +#include "IO/HDF5_IO.h" +#include "common/Array.h" +#include "common/Array.hpp" +#include "common/Utilities.h" + +#include +#include +#include +#include +#include + + +namespace IO { +namespace HDF5 { + + +/******************************************************** + * External instantiations (scalar) * + ********************************************************/ +// clang-format off +template<> void writeHDF5( hid_t, const std::string &, const char & ); +template<> void readHDF5( hid_t, const std::string &, char & ); +template<> void writeHDF5( hid_t, const std::string &, const bool & ); +template<> void readHDF5( hid_t, const std::string &, bool & ); +template<> void writeHDF5( hid_t, const std::string &, const int & ); +template<> void readHDF5( hid_t, const std::string &, int & ); +template<> void writeHDF5( hid_t, const std::string &, const long & ); +template<> void readHDF5( hid_t, const std::string &, long & ); +template<> void writeHDF5( hid_t, const std::string &, const float & ); +template<> void readHDF5( hid_t, const std::string &, float & ); +template<> void writeHDF5( hid_t, const std::string &, const double & ); +template<> void readHDF5( hid_t, const std::string &, double & ); +template<> void writeHDF5( hid_t, const std::string &, const unsigned char & ); +template<> void readHDF5( hid_t, const std::string &, unsigned char & ); +template<> void writeHDF5( hid_t, const std::string &, const unsigned int & ); +template<> void readHDF5( hid_t, const std::string &, unsigned int & ); +template<> void writeHDF5( hid_t, const std::string &, const unsigned long & ); +template<> void readHDF5( hid_t, const std::string &, unsigned long & ); +template<> void writeHDF5( hid_t, const std::string &, const std::string & ); +template<> void readHDF5( hid_t, const std::string &, std::string & ); +template<> void writeHDF5>( hid_t, const std::string &, const std::complex & ); +template<> void readHDF5>( hid_t, const std::string &, std::complex & ); +template<> void writeHDF5>( hid_t, const std::string &, const std::complex & ); +template<> void readHDF5>( hid_t, const std::string &, std::complex & ); +// clang-format on + + +/******************************************************** + * External instantiations (Array) * + ********************************************************/ +// clang-format off +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>( hid_t, const std::string &, const Array & ); +template<> void readHDF5>( hid_t, const std::string &, Array & ); +template<> void writeHDF5>>( hid_t, const std::string &, const Array> & ); +template<> void readHDF5>>( hid_t, const std::string &, Array> & ); +template<> void writeHDF5>>( hid_t, const std::string &, const Array> & ); +template<> void readHDF5>>( hid_t, const std::string &, Array> & ); +// clang-format on + + +/****************************************************************** + * Default implimentation * + ******************************************************************/ +/*template +void writeHDF5( hid_t fid, const std::string &name, const TYPE &x ) +{ + NULL_USE( fid ); + if constexpr ( is_shared_ptr::value ) { + // We are dealing with a std::shared_ptr + writeHDF5( fid, name, *x ); + } else if constexpr ( is_vector::value ) { + // We are dealing with a std::vector + typedef decltype( *x.begin() ) TYPE2; + typedef typename std::remove_reference::type TYPE3; + typedef typename std::remove_cv::type TYPE4; + Array y; + y.viewRaw( { x.size() }, const_cast( x.data() ) ); + writeHDF5( fid, name, y ); + } else if constexpr ( std::is_array::value ) { + // We are dealing with a std::array + typedef decltype( *x.begin() ) TYPE2; + typedef typename std::remove_reference::type TYPE3; + typedef typename std::remove_cv::type TYPE4; + Array y; + y.viewRaw( { x.size() }, const_cast( x.data() ) ); + writeHDF5( fid, name, y ); + } else if constexpr ( is_Array::value ) { + // We are dealing with an Array + std::string typeName = Utilities::demangle( typeid( TYPE ).name() ); + throw std::logic_error( "Unsupported type writeHDF5>" ); + } else if constexpr ( std::is_same::value ) { + // We are dealing with a std::string (should be handled through specialization) + throw std::logic_error( "Internal error" ); + } else if constexpr ( std::is_same::value || + std::is_same::value || + std::is_same::value ) { + // We are dealing with a string or char array + writeHDF5( fid, name, std::string( x ) ); + } else if constexpr ( has_size::value ) { + // We are dealing with a container + typedef decltype( *x.begin() ) TYPE2; + typedef typename std::remove_reference::type TYPE3; + typedef typename std::remove_cv::type TYPE4; + std::vector x2( x.begin(), x.end() ); + writeHDF5>( fid, name, x2 ); + } else { + throw std::logic_error( "Unsupported type" ); + } +} +template +void readHDF5( hid_t fid, const std::string &name, TYPE &x ) +{ + NULL_USE( fid ); + if constexpr ( is_shared_ptr::value ) { + // We are dealing with a std::shared_ptr + readHDF5( fid, name, *x ); + } else if constexpr ( is_vector::value ) { + // We are dealing with a std::vector + typedef typename std::remove_reference::type TYPE2; + Array y; + readHDF5( fid, name, y ); + x.resize( y.length() ); + // Swap the elements in the arrays to use the move operator + for ( size_t i = 0; i < x.size(); i++ ) + std::swap( x[i], y( i ) ); + } else if constexpr ( std::is_array::value ) { + // We are dealing with a std::array + typedef typename std::remove_reference::type TYPE2; + Array y; + readHDF5( fid, name, y ); + ASSERT( y.length() == x.size() ); + // Swap the elements in the arrays to use the move operator + for ( size_t i = 0; i < x.size(); i++ ) + std::swap( x[i], y( i ) ); + } else if constexpr ( is_Array::value ) { + // We are dealing with an Array + std::string typeName = Utilities::demangle( typeid( TYPE ).name() ); + throw std::logic_error( "Unsupported type readHDF5>" ); + } else if constexpr ( std::is_same::value ) { + // We are dealing with a std::string (should be handled through specialization) + throw std::logic_error( "Internal error" ); + } else if constexpr ( std::is_same::value || + std::is_same::value || + std::is_same::value ) { + // We are dealing with a string or char array + throw std::logic_error( + "Reading data into a string, char*, const char* is not supported" ); + } else if constexpr ( has_size::value ) { + // We are dealing with a container + typedef typename std::remove_reference::type TYPE2; + Array y; + readHDF5( fid, name, y ); + if ( x.size() == y.length() ) { + auto it = x.begin(); + for ( size_t i = 0; i < y.length(); i++, ++it ) + *it = y( i ); + } else { + throw std::logic_error( "Reading data into an arbitrary container is not finished" ); + } + } else { + throw std::logic_error( "Unsupported type" ); + } +}*/ + + +/************************************************************************ + * Helper function to get the size of an Array * + * Note that HDF5 uses C ordered arrays so we need to flip the dimensions* + ************************************************************************/ +template +inline std::vector arraySize( const Array &x ) +{ + int N = x.ndim(); + auto s1 = x.size(); + std::vector s2( std::max( N, 1 ), 0 ); + for ( int i = 0; i < N; i++ ) + s2[N - i - 1] = static_cast( s1[i] ); + return s2; +} +inline std::vector convertSize( int N, const hsize_t *dims ) +{ + if ( N == 0 ) + return std::vector( 1, 1 ); + std::vector size( N, 0 ); + for ( int i = 0; i < N; i++ ) + size[N - i - 1] = static_cast( dims[i] ); + return size; +} + + +/************************************************************************ + * readAndConvertHDF5Data * + ************************************************************************/ +template +typename std::enable_if::value || std::is_floating_point::value, void>::type +readAndConvertHDF5Data( hid_t dataset, hid_t datatype, Array &data ) +{ + if ( H5Tequal( datatype, H5T_NATIVE_CHAR ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_UCHAR ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_INT8 ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_UINT8 ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_INT ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_UINT ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_LONG ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_ULONG ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_FLOAT ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else if ( H5Tequal( datatype, H5T_NATIVE_DOUBLE ) ) { + Array data2( data.size() ); + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data2.data() ); + data.copy( data2 ); + } else { + ERROR( "We need to convert unknown data format" ); + } +} +template +typename std::enable_if::value && !std::is_floating_point::value, + void>::type +readAndConvertHDF5Data( hid_t, hid_t, Array & ) +{ + ERROR( "Unable to convert data" ); +} + + +/************************************************************************ + * Default writeHDF5Array * + ************************************************************************/ +template +void writeHDF5ArrayDefault( hid_t fid, const std::string &name, const Array &data ) +{ + size_t N_bytes = data.length() * sizeof( T ); + auto dim = arraySize( data ); + hid_t plist = H5P_DEFAULT; + if ( N_bytes < 0x7500 ) { + // Use compact storage (limited to < 30K) + plist = H5Pcreate( H5P_DATASET_CREATE ); + auto status = H5Pset_layout( plist, H5D_COMPACT ); + ASSERT( status == 0 ); + } else if ( std::is_same::value || std::is_same::value ) { + // Use compression if availible + plist = createChunk( dim, defaultCompression( fid ) ); + } + hid_t dataspace = H5Screate_simple( dim.size(), dim.data(), NULL ); + hid_t datatype = getHDF5datatype(); + hid_t dataset = + H5Dcreate2( fid, name.data(), datatype, dataspace, H5P_DEFAULT, plist, H5P_DEFAULT ); + const void *ptr = data.data() == NULL ? ( (void *) 1 ) : data.data(); + H5Dwrite( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, ptr ); + H5Dclose( dataset ); + H5Tclose( datatype ); + H5Sclose( dataspace ); + if ( plist != H5P_DEFAULT ) + H5Pclose( plist ); +} + + +/************************************************************************ + * Default readHDF5Array * + ************************************************************************/ +template +void readHDF5ArrayDefault( hid_t fid, const std::string &name, Array &data ) +{ + if ( !H5Dexists( fid, name ) ) { + // Dataset does not exist + data.resize( 0 ); + return; + } + hid_t dataset = H5Dopen2( fid, name.data(), H5P_DEFAULT ); + hid_t datatype = H5Dget_type( dataset ); + hid_t dataspace = H5Dget_space( dataset ); + hsize_t dims0[10]; + int ndim = H5Sget_simple_extent_dims( dataspace, dims0, NULL ); + auto dims = convertSize( ndim, dims0 ); + data.resize( dims ); + hid_t datatype2 = getHDF5datatype(); + if ( data.empty() ) { + // The data is empty + } else if ( H5Tequal( datatype, datatype2 ) ) { + // The type of Array and the data in HDF5 match + H5Dread( dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data.data() ); + } else { + // Try to convert the data + readAndConvertHDF5Data( dataset, datatype, data ); + } + H5Dclose( dataset ); + H5Tclose( datatype ); + H5Tclose( datatype2 ); + H5Sclose( dataspace ); +} + + +} // namespace HDF5 +} // namespace IO + + +#endif +#endif diff --git a/IO/Mesh.cpp b/IO/Mesh.cpp index 9966bf52..5e5eb96a 100644 --- a/IO/Mesh.cpp +++ b/IO/Mesh.cpp @@ -133,7 +133,7 @@ size_t PointList::numberPointsVar( VariableType type ) const } std::pair PointList::pack( int level ) const { - std::pair data_out( 0, NULL ); + std::pair data_out( 0, nullptr ); if ( level == 0 ) { data_out.first = ( 2 + 3 * points.size() ) * sizeof( double ); double *data_ptr = new double[2 + 3 * points.size()]; @@ -596,6 +596,8 @@ std::string getString( FileFormat type ) return "new(single)"; else if ( type == FileFormat::SILO ) return "silo"; + else if ( type == FileFormat::HDF5 ) + return "hdf5"; else ERROR( "Invalid type" ); return ""; @@ -611,6 +613,8 @@ FileFormat getFileFormat( const std::string &type_in ) return FileFormat::NEW_SINGLE; else if ( type == "silo" || type == "4" ) return FileFormat::SILO; + else if ( type == "hdf5" || type == "5" ) + return FileFormat::HDF5; else ERROR( "Invalid type: " + type ); return FileFormat::SILO; diff --git a/IO/Mesh.h b/IO/Mesh.h index a420f95d..9e5f32e6 100644 --- a/IO/Mesh.h +++ b/IO/Mesh.h @@ -24,7 +24,7 @@ enum class VariableType { }; enum class DataType { Double, Float, Int, Null }; enum class MeshType { PointMesh, SurfaceMesh, VolumeMesh, Unknown }; -enum class FileFormat { OLD, NEW, NEW_SINGLE, SILO }; +enum class FileFormat { OLD, NEW, NEW_SINGLE, SILO, HDF5 }; //! Convert enums to/from strings (more future-proof than static_cast) diff --git a/IO/MeshDatabase.cpp b/IO/MeshDatabase.cpp index 63702c7b..975c9e27 100644 --- a/IO/MeshDatabase.cpp +++ b/IO/MeshDatabase.cpp @@ -366,7 +366,7 @@ std::vector read( const std::string &filename ) PROFILE_START( "read" ); FILE *fid = fopen( filename.c_str(), "rb" ); if ( fid == NULL ) - ERROR( "Error opening file" ); + ERROR( "Error opening file: " + filename ); char *line = new char[10000]; while ( std::fgets( line, 1000, fid ) != NULL ) { if ( line[0] < 32 ) { diff --git a/IO/Reader.cpp b/IO/Reader.cpp index e63b8dd3..014e1e7b 100644 --- a/IO/Reader.cpp +++ b/IO/Reader.cpp @@ -1,13 +1,10 @@ #include "IO/Reader.h" +#include "IO/HDF5_IO.h" #include "IO/IOHelpers.h" #include "IO/Mesh.h" #include "IO/MeshDatabase.h" -#include "common/Utilities.h" - -#ifdef USE_SILO #include "IO/silo.h" -#endif - +#include "common/Utilities.h" #include #include @@ -62,14 +59,16 @@ std::vector IO::readTimesteps( const std::string &path, const std:: filename += "summary.LBM"; } else if ( format == "silo" ) { filename += "LBM.visit"; + } else if ( format == "hdf5" ) { + filename += "LBM.visit"; } else if ( format == "auto" ) { - bool test_old = fileExists( path + "/summary.LBM" ); - bool test_silo = fileExists( path + "/LBM.visit" ); - if ( test_old && test_silo ) { + bool test_old = fileExists( path + "/summary.LBM" ); + bool test_new = fileExists( path + "/LBM.visit" ); + if ( test_old && test_new ) { ERROR( "Unable to determine format (both summary.LBM and LBM.visit exist)" ); } else if ( test_old ) { filename += "summary.LBM"; - } else if ( test_silo ) { + } else if ( test_new ) { filename += "LBM.visit"; } else { ERROR( "Unable to determine format (neither summary.LBM or LBM.visit exist)" ); @@ -88,6 +87,9 @@ std::vector IO::readTimesteps( const std::string &path, const std:: std::string line( buf ); line.resize( line.size() - 1 ); auto pos = line.find( "summary.silo" ); + if ( pos != std::string::npos ) + line.resize( pos ); + pos = line.find( "summary.xmf" ); if ( pos != std::string::npos ) line.resize( pos ); if ( line.empty() ) @@ -170,8 +172,8 @@ std::shared_ptr IO::getMesh( const std::string &path, const std::strin if ( count % 3 != 0 ) ERROR( "Error reading file" ); if ( meshDatabase.type == IO::MeshType::PointMesh ) { - size_t N = count / 3; - std::shared_ptr pointlist( new PointList( N ) ); + size_t N = count / 3; + auto pointlist = std::make_shared( N ); std::vector &P = pointlist->points; for ( size_t i = 0; i < N; i++ ) { P[i].x = data[3 * i + 0]; @@ -182,8 +184,8 @@ std::shared_ptr IO::getMesh( const std::string &path, const std::strin } else if ( meshDatabase.type == IO::MeshType::SurfaceMesh ) { if ( count % 9 != 0 ) ERROR( "Error reading file (2)" ); - size_t N_tri = count / 9; - std::shared_ptr trilist( new TriList( N_tri ) ); + size_t N_tri = count / 9; + auto trilist = std::make_shared( N_tri ); std::vector &A = trilist->A; std::vector &B = trilist->B; std::vector &C = trilist->C; @@ -201,7 +203,7 @@ std::shared_ptr IO::getMesh( const std::string &path, const std::strin mesh = trilist; } else if ( meshDatabase.type == IO::MeshType::VolumeMesh ) { // this was never supported in the old format - mesh = std::shared_ptr( new DomainMesh() ); + mesh = std::make_shared(); } else { ERROR( "Unknown mesh type" ); } @@ -222,13 +224,13 @@ std::shared_ptr IO::getMesh( const std::string &path, const std::strin fclose( fid ); ASSERT( count == bytes ); if ( meshDatabase.meshClass == "PointList" ) { - mesh.reset( new IO::PointList() ); + mesh = std::make_shared(); } else if ( meshDatabase.meshClass == "TriMesh" ) { - mesh.reset( new IO::TriMesh() ); + mesh = std::make_shared(); } else if ( meshDatabase.meshClass == "TriList" ) { - mesh.reset( new IO::TriList() ); + mesh = std::make_shared(); } else if ( meshDatabase.meshClass == "DomainMesh" ) { - mesh.reset( new IO::DomainMesh() ); + mesh = std::make_shared(); } else { ERROR( "Unknown mesh class" ); } @@ -243,7 +245,7 @@ std::shared_ptr IO::getMesh( const std::string &path, const std::strin if ( meshDatabase.meshClass == "PointList" ) { Array coords = silo::readPointMesh( fid, database.name ); ASSERT( coords.size( 1 ) == 3 ); - std::shared_ptr mesh2( new IO::PointList( coords.size( 0 ) ) ); + auto mesh2 = std::make_shared( coords.size( 0 ) ); for ( size_t i = 0; i < coords.size( 1 ); i++ ) { mesh2->points[i].x = coords( i, 0 ); mesh2->points[i].y = coords( i, 1 ); @@ -257,7 +259,7 @@ std::shared_ptr IO::getMesh( const std::string &path, const std::strin ASSERT( tri.size( 1 ) == 3 && coords.size( 1 ) == 3 ); int N_tri = tri.size( 0 ); int N_point = coords.size( 0 ); - std::shared_ptr mesh2( new IO::TriMesh( N_tri, N_point ) ); + auto mesh2 = std::make_shared( N_tri, N_point ); for ( int i = 0; i < N_point; i++ ) { mesh2->vertices->points[i].x = coords( i, 0 ); mesh2->vertices->points[i].y = coords( i, 1 ); @@ -280,14 +282,81 @@ std::shared_ptr IO::getMesh( const std::string &path, const std::strin silo::readUniformMesh( fid, database.name, range, N ); auto rankinfo = silo::read( fid, database.name + "_rankinfo" ); RankInfoStruct rank_data( rankinfo[0], rankinfo[1], rankinfo[2], rankinfo[3] ); - mesh.reset( new IO::DomainMesh( rank_data, N[0], N[1], N[2], range[1] - range[0], - range[3] - range[2], range[5] - range[4] ) ); + mesh = std::make_shared( rank_data, N[0], N[1], N[2], + range[1] - range[0], range[3] - range[2], range[5] - range[4] ); } else { ERROR( "Unknown mesh class" ); } silo::close( fid ); #else ERROR( "Build without silo support" ); +#endif + } else if ( meshDatabase.format == FileFormat::HDF5 ) { + // Reading an hdf5 file +#ifdef USE_HDF5 + auto &database = meshDatabase.domains[domain]; + auto filename = path + "/" + timestep + "/" + database.file; + auto fid = IO::HDF5::openHDF5( filename, "r" ); + auto gid = IO::HDF5::openGroup( fid, database.name ); + if ( meshDatabase.meshClass == "PointList" ) { + std::vector x, y, z; + IO::HDF5::readHDF5( gid, "x", x ); + IO::HDF5::readHDF5( gid, "y", y ); + IO::HDF5::readHDF5( gid, "z", z ); + ASSERT( y.size() == x.size() && z.size() == x.size() ); + auto mesh2 = std::make_shared( x.size() ); + for ( size_t i = 0; i < x.size(); i++ ) { + mesh2->points[i].x = x[i]; + mesh2->points[i].y = y[i]; + mesh2->points[i].z = z[i]; + } + mesh = mesh2; + } else if ( meshDatabase.meshClass == "TriMesh" || meshDatabase.meshClass == "TriList" ) { + // Read the points + std::vector x, y, z; + IO::HDF5::readHDF5( gid, "x", x ); + IO::HDF5::readHDF5( gid, "y", y ); + IO::HDF5::readHDF5( gid, "z", z ); + // Read the triangles + Array tri; + IO::HDF5::readHDF5( gid, "tri", tri ); + ASSERT( tri.size( 0 ) == 3 ); + size_t N_tri = tri.size( 1 ); + size_t N_point = x.size(); + auto mesh2 = std::make_shared( N_tri, N_point ); + for ( size_t i = 0; i < N_point; i++ ) { + mesh2->vertices->points[i].x = x[i]; + mesh2->vertices->points[i].y = y[i]; + mesh2->vertices->points[i].z = z[i]; + } + for ( size_t i = 0; i < N_tri; i++ ) { + mesh2->A[i] = tri( 0, i ); + mesh2->B[i] = tri( 1, i ); + mesh2->C[i] = tri( 2, i ); + } + if ( meshDatabase.meshClass == "TriMesh" ) { + mesh = mesh2; + } else if ( meshDatabase.meshClass == "TriList" ) { + auto trilist = IO::getTriList( std::dynamic_pointer_cast( mesh2 ) ); + mesh = trilist; + } + } else if ( meshDatabase.meshClass == "DomainMesh" ) { + std::vector range; + std::vector N; + std::vector rankinfo; + IO::HDF5::readHDF5( gid, "range", range ); + IO::HDF5::readHDF5( gid, "N", N ); + IO::HDF5::readHDF5( gid, "rankinfo", rankinfo ); + RankInfoStruct rank_data( rankinfo[0], rankinfo[1], rankinfo[2], rankinfo[3] ); + mesh = std::make_shared( rank_data, N[0], N[1], N[2], + range[1] - range[0], range[3] - range[2], range[5] - range[4] ); + } else { + ERROR( "Unknown mesh class" ); + } + IO::HDF5::closeGroup( gid ); + IO::HDF5::closeHDF5( fid ); +#else + ERROR( "Build without hdf5 support" ); #endif } else { ERROR( "Unknown format" ); @@ -322,7 +391,7 @@ std::shared_ptr IO::getVariable( const std::string &path, const st size_t N = atol( values[2].c_str() ); size_t bytes = atol( values[3].c_str() ); std::string precision = values[4]; - var = std::shared_ptr( new IO::Variable() ); + var = std::make_shared(); var->dim = dim; var->type = getVariableType( type ); var->name = variable; @@ -341,7 +410,7 @@ std::shared_ptr IO::getVariable( const std::string &path, const st auto variableDatabase = meshDatabase.getVariableDatabase( variable ); std::string filename = path + "/" + timestep + "/" + database.file; auto fid = silo::open( filename, silo::READ ); - var.reset( new Variable( variableDatabase.dim, variableDatabase.type, variable ) ); + var = std::make_shared( variableDatabase.dim, variableDatabase.type, variable ); if ( meshDatabase.meshClass == "PointList" ) { var->data = silo::readPointMeshVariable( fid, variable ); } else if ( meshDatabase.meshClass == "TriMesh" || meshDatabase.meshClass == "TriList" ) { @@ -355,7 +424,30 @@ std::shared_ptr IO::getVariable( const std::string &path, const st #else ERROR( "Build without silo support" ); #endif - + } else if ( meshDatabase.format == FileFormat::HDF5 ) { + // Reading an hdf5 file +#ifdef USE_HDF5 + auto &database = meshDatabase.domains[domain]; + auto varDatabase = meshDatabase.getVariableDatabase( variable ); + auto filename = path + "/" + timestep + "/" + database.file; + var = std::make_shared( varDatabase.dim, varDatabase.type, variable ); + auto fid = IO::HDF5::openHDF5( filename, "r" ); + auto gid = IO::HDF5::openGroup( fid, database.name ); + IO::HDF5::readHDF5( gid, var->name, var->data ); + IO::HDF5::closeHDF5( fid ); + if ( meshDatabase.meshClass == "PointList" || meshDatabase.meshClass == "TriMesh" || + meshDatabase.meshClass == "TriList" ) { + if ( var->data.ndim() == 2 && var->data.size( 0 ) == 3 ) + var->data = var->data.permute( { 1, 0 } ); + } else if ( meshDatabase.meshClass == "DomainMesh" ) { + if ( var->data.ndim() == 4 && var->data.size( 0 ) == 3 ) + var->data = var->data.permute( { 1, 2, 3, 0 } ); + } else { + ERROR( "Unknown mesh class" ); + } +#else + ERROR( "Build without silo support" ); +#endif } else { ERROR( "Unknown format" ); } diff --git a/IO/SiloWriter.cpp b/IO/SiloWriter.cpp new file mode 100644 index 00000000..c2e8df26 --- /dev/null +++ b/IO/SiloWriter.cpp @@ -0,0 +1,251 @@ +#include "IO/HDF5_IO.h" +#include "IO/IOHelpers.h" +#include "IO/MeshDatabase.h" +#include "IO/Writer.h" +#include "IO/silo.h" +#include "common/MPI.h" +#include "common/Utilities.h" + +#include +#include +#include +#include +#include + + +#ifdef USE_SILO + + +// Write a PointList mesh (and variables) to a file +template +static void writeSiloPointMesh( + DBfile *fid, const IO::PointList &mesh, const std::string &meshname ) +{ + const auto &points = mesh.getPoints(); + std::vector x( points.size() ), y( points.size() ), z( points.size() ); + for ( size_t i = 0; i < x.size(); i++ ) { + x[i] = points[i].x; + y[i] = points[i].y; + z[i] = points[i].z; + } + const TYPE *coords[] = { x.data(), y.data(), z.data() }; + IO::silo::writePointMesh( fid, meshname, 3, points.size(), coords ); +} +static void writeSiloPointList( + DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database ) +{ + const IO::PointList &mesh = dynamic_cast( *meshData.mesh ); + const std::string meshname = database.domains[0].name; + if ( meshData.precision == IO::DataType::Double ) { + writeSiloPointMesh( fid, mesh, meshname ); + } else if ( meshData.precision == IO::DataType::Float ) { + writeSiloPointMesh( fid, mesh, meshname ); + } else { + ERROR( "Unsupported format" ); + } + const auto &points = mesh.getPoints(); + std::vector x( points.size() ), y( points.size() ), z( points.size() ); + for ( size_t i = 0; i < x.size(); i++ ) { + x[i] = points[i].x; + y[i] = points[i].y; + z[i] = points[i].z; + } + const double *coords[] = { x.data(), y.data(), z.data() }; + IO::silo::writePointMesh( fid, meshname, 3, points.size(), coords ); + for ( size_t i = 0; i < meshData.vars.size(); i++ ) { + const IO::Variable &var = *meshData.vars[i]; + if ( var.precision == IO::DataType::Double ) { + IO::silo::writePointMeshVariable( fid, meshname, var.name, var.data ); + } else if ( var.precision == IO::DataType::Float ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + IO::silo::writePointMeshVariable( fid, meshname, var.name, data2 ); + } else if ( var.precision == IO::DataType::Int ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + IO::silo::writePointMeshVariable( fid, meshname, var.name, data2 ); + } else { + ERROR( "Unsupported format" ); + } + } +} +// Write a TriMesh mesh (and variables) to a file +template +static void writeSiloTriMesh( DBfile *fid, const IO::TriMesh &mesh, const std::string &meshname ) +{ + const auto &points = mesh.vertices->getPoints(); + std::vector x( points.size() ), y( points.size() ), z( points.size() ); + for ( size_t i = 0; i < x.size(); i++ ) { + x[i] = points[i].x; + y[i] = points[i].y; + z[i] = points[i].z; + } + const TYPE *coords[] = { x.data(), y.data(), z.data() }; + const int *tri[] = { mesh.A.data(), mesh.B.data(), mesh.C.data() }; + IO::silo::writeTriMesh( fid, meshname, 3, 2, points.size(), coords, mesh.A.size(), tri ); +} +static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct &meshData, + const IO::TriMesh &mesh, IO::MeshDatabase database ) +{ + const std::string meshname = database.domains[0].name; + if ( meshData.precision == IO::DataType::Double ) { + writeSiloTriMesh( fid, mesh, meshname ); + } else if ( meshData.precision == IO::DataType::Float ) { + writeSiloTriMesh( fid, mesh, meshname ); + } else { + ERROR( "Unsupported format" ); + } + for ( size_t i = 0; i < meshData.vars.size(); i++ ) { + const IO::Variable &var = *meshData.vars[i]; + if ( var.precision == IO::DataType::Double ) { + IO::silo::writeTriMeshVariable( fid, 3, meshname, var.name, var.data, var.type ); + } else if ( var.precision == IO::DataType::Float ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + IO::silo::writeTriMeshVariable( fid, 3, meshname, var.name, data2, var.type ); + } else if ( var.precision == IO::DataType::Int ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + IO::silo::writeTriMeshVariable( fid, 3, meshname, var.name, data2, var.type ); + } else { + ERROR( "Unsupported format" ); + } + } +} +static void writeSiloTriMesh( + DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database ) +{ + const IO::TriMesh &mesh = dynamic_cast( *meshData.mesh ); + writeSiloTriMesh2( fid, meshData, mesh, database ); +} +static void writeSiloTriList( + DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database ) +{ + auto mesh = getTriMesh( meshData.mesh ); + writeSiloTriMesh2( fid, meshData, *mesh, database ); +} +// Write a DomainMesh mesh (and variables) to a file +static void writeSiloDomainMesh( + DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database ) +{ + const IO::DomainMesh &mesh = dynamic_cast( *meshData.mesh ); + RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); + std::array range = { info.ix * mesh.Lx / info.nx, + ( info.ix + 1 ) * mesh.Lx / info.nx, info.jy * mesh.Ly / info.ny, + ( info.jy + 1 ) * mesh.Ly / info.ny, info.kz * mesh.Lz / info.nz, + ( info.kz + 1 ) * mesh.Lz / info.nz }; + std::array N = { mesh.nx, mesh.ny, mesh.nz }; + auto meshname = database.domains[0].name; + IO::silo::writeUniformMesh<3>( fid, meshname, range, N ); + IO::silo::write( + fid, meshname + "_rankinfo", { mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz } ); + for ( size_t i = 0; i < meshData.vars.size(); i++ ) { + const auto &var = *meshData.vars[i]; + if ( var.precision == IO::DataType::Double ) { + IO::silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, var.data, var.type ); + } else if ( var.precision == IO::DataType::Float ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + IO::silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, data2, var.type ); + } else if ( var.precision == IO::DataType::Int ) { + Array data2( var.data.size() ); + data2.copy( var.data ); + IO::silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, data2, var.type ); + } else { + ERROR( "Unsupported format" ); + } + } +} +// Write a mesh (and variables) to a file +static IO::MeshDatabase write_domain_silo( DBfile *fid, const std::string &filename, + const IO::MeshDataStruct &mesh, IO::FileFormat format, int rank ) +{ + // Create the MeshDatabase + auto database = getDatabase( filename, mesh, format, rank ); + if ( database.meshClass == "PointList" ) { + writeSiloPointList( fid, mesh, database ); + } else if ( database.meshClass == "TriMesh" ) { + writeSiloTriMesh( fid, mesh, database ); + } else if ( database.meshClass == "TriList" ) { + writeSiloTriList( fid, mesh, database ); + } else if ( database.meshClass == "DomainMesh" ) { + writeSiloDomainMesh( fid, mesh, database ); + } else { + ERROR( "Unknown mesh class" ); + } + return database; +} +// Write the summary file for silo +std::pair getSiloMeshType( const std::string &meshClass ) +{ + int meshType = 0; + int varType = 0; + if ( meshClass == "PointList" ) { + meshType = DB_POINTMESH; + varType = DB_POINTVAR; + } else if ( meshClass == "TriMesh" ) { + meshType = DB_UCDMESH; + varType = DB_UCDVAR; + } else if ( meshClass == "TriList" ) { + meshType = DB_UCDMESH; + varType = DB_UCDVAR; + } else if ( meshClass == "DomainMesh" ) { + meshType = DB_QUAD_RECT; + varType = DB_QUADVAR; + } else { + ERROR( "Unknown mesh class" ); + } + return std::make_pair( meshType, varType ); +} +void writeSiloSummary( + const std::vector &meshes_written, const std::string &filename ) +{ + auto fid = IO::silo::open( filename, IO::silo::CREATE ); + for ( const auto &data : meshes_written ) { + auto type = getSiloMeshType( data.meshClass ); + std::vector meshTypes( data.domains.size(), type.first ); + std::vector varTypes( data.domains.size(), type.second ); + std::vector meshNames; + for ( const auto &tmp : data.domains ) + meshNames.push_back( tmp.file + ":" + tmp.name ); + IO::silo::writeMultiMesh( fid, data.name, meshNames, meshTypes ); + for ( const auto &variable : data.variables ) { + std::vector varnames; + for ( const auto &tmp : data.domains ) + varnames.push_back( tmp.file + ":" + variable.name ); + IO::silo::writeMultiVar( fid, variable.name, varnames, varTypes ); + } + } + IO::silo::close( fid ); +} +// Write the mesh data to silo +std::vector writeMeshesSilo( const std::vector &meshData, + const std::string &path, IO::FileFormat format, int rank ) +{ + std::vector meshes_written; + char filename[100], fullpath[200]; + sprintf( filename, "%05i.silo", rank ); + sprintf( fullpath, "%s/%s", path.c_str(), filename ); + auto fid = IO::silo::open( fullpath, IO::silo::CREATE ); + for ( size_t i = 0; i < meshData.size(); i++ ) { + auto mesh = meshData[i].mesh; + meshes_written.push_back( write_domain_silo( fid, filename, meshData[i], format, rank ) ); + } + IO::silo::close( fid ); + return meshes_written; +} + + +#else + + +// Write the mesh data to silo +std::vector writeMeshesSilo( + const std::vector &, const std::string &, IO::FileFormat, int ) +{ + return std::vector(); +} +void writeSiloSummary( const std::vector &, const std::string & ); + + +#endif diff --git a/IO/Writer.cpp b/IO/Writer.cpp index 051db47d..0fb8d135 100644 --- a/IO/Writer.cpp +++ b/IO/Writer.cpp @@ -1,10 +1,13 @@ #include "IO/Writer.h" +#include "IO/HDF5_IO.h" #include "IO/IOHelpers.h" #include "IO/MeshDatabase.h" -#include "IO/silo.h" +#include "IO/Xdmf.h" #include "common/MPI.h" #include "common/Utilities.h" +#include "ProfilerApp.h" + #include #include #include @@ -12,7 +15,17 @@ #include -enum class Format { OLD, NEW, SILO, UNKNOWN }; +enum class Format { OLD, NEW, SILO, HDF5, UNKNOWN }; + + +/**************************************************** + * External declerations * + ****************************************************/ +std::vector writeMeshesSilo( + const std::vector &, const std::string &, IO::FileFormat, int ); +void writeSiloSummary( const std::vector &, const std::string & ); +std::vector writeMeshesHDF5( + const std::vector &, const std::string &, IO::FileFormat, int, Xdmf & ); /**************************************************** @@ -75,6 +88,8 @@ void IO::initialize( const std::string &path, const std::string &format, bool ap global_IO_format = Format::NEW; else if ( format == "silo" ) global_IO_format = Format::SILO; + else if ( format == "hdf5" ) + global_IO_format = Format::HDF5; else ERROR( "Unknown format" ); int rank = Utilities::MPI( MPI_COMM_WORLD ).getRank(); @@ -83,7 +98,7 @@ void IO::initialize( const std::string &path, const std::string &format, bool ap std::string filename; if ( global_IO_format == Format::OLD || global_IO_format == Format::NEW ) filename = global_IO_path + "/summary.LBM"; - else if ( global_IO_format == Format::SILO ) + else if ( global_IO_format == Format::SILO || global_IO_format == Format::HDF5 ) filename = global_IO_path + "/LBM.visit"; else ERROR( "Unknown format" ); @@ -116,10 +131,10 @@ static std::vector writeMeshesOrigFormat( domain.file = filename; domain.offset = 0; mesh_entry.domains.push_back( domain ); - if ( !meshData[i].vars.empty() ) { + static bool printVariableSupportMsg = true; + if ( !meshData[i].vars.empty() && printVariableSupportMsg ) { + printVariableSupportMsg = false; printf( "Warning: variables are not supported with this format (original)\n" ); - // for (size_t j=0; jname ); } const std::string meshClass = mesh->className(); if ( meshClass == "PointList" ) { @@ -170,7 +185,7 @@ static std::vector writeMeshesOrigFormat( // Create the database entry for the mesh data -static IO::MeshDatabase getDatabase( +IO::MeshDatabase IO::getDatabase( const std::string &filename, const IO::MeshDataStruct &mesh, IO::FileFormat format, int rank ) { char domainname[100]; @@ -212,6 +227,7 @@ static IO::MeshDatabase getDatabase( static IO::MeshDatabase write_domain( FILE *fid, const std::string &filename, const IO::MeshDataStruct &mesh, IO::FileFormat format, int rank ) { + ASSERT( !mesh.meshName.empty() ); const int level = 0; // Create the MeshDatabase IO::MeshDatabase database = getDatabase( filename, mesh, format, rank ); @@ -234,6 +250,8 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string &filename, size_t N = mesh.vars[i]->data.length(); size_t N_mesh = mesh.mesh->numberPointsVar( mesh.vars[i]->type ); ASSERT( N == dim * N_mesh ); + ASSERT( !type.empty() ); + ASSERT( !variable.name.empty() ); fprintf( fid, "Var: %s-%05i-%s: %i, %s, %lu, %lu, double\n", database.name.c_str(), rank, variable.name.c_str(), dim, type.data(), N_mesh, N * sizeof( double ) ); fwrite( mesh.vars[i]->data.data(), sizeof( double ), N, fid ); @@ -243,212 +261,6 @@ static IO::MeshDatabase write_domain( FILE *fid, const std::string &filename, } -#ifdef USE_SILO -// Write a PointList mesh (and variables) to a file -template -static void writeSiloPointMesh( - DBfile *fid, const IO::PointList &mesh, const std::string &meshname ) -{ - const auto &points = mesh.getPoints(); - std::vector x( points.size() ), y( points.size() ), z( points.size() ); - for ( size_t i = 0; i < x.size(); i++ ) { - x[i] = points[i].x; - y[i] = points[i].y; - z[i] = points[i].z; - } - const TYPE *coords[] = { x.data(), y.data(), z.data() }; - IO::silo::writePointMesh( fid, meshname, 3, points.size(), coords ); -} -static void writeSiloPointList( - DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database ) -{ - const IO::PointList &mesh = dynamic_cast( *meshData.mesh ); - const std::string meshname = database.domains[0].name; - if ( meshData.precision == IO::DataType::Double ) { - writeSiloPointMesh( fid, mesh, meshname ); - } else if ( meshData.precision == IO::DataType::Float ) { - writeSiloPointMesh( fid, mesh, meshname ); - } else { - ERROR( "Unsupported format" ); - } - const auto &points = mesh.getPoints(); - std::vector x( points.size() ), y( points.size() ), z( points.size() ); - for ( size_t i = 0; i < x.size(); i++ ) { - x[i] = points[i].x; - y[i] = points[i].y; - z[i] = points[i].z; - } - const double *coords[] = { x.data(), y.data(), z.data() }; - IO::silo::writePointMesh( fid, meshname, 3, points.size(), coords ); - for ( size_t i = 0; i < meshData.vars.size(); i++ ) { - const IO::Variable &var = *meshData.vars[i]; - if ( var.precision == IO::DataType::Double ) { - IO::silo::writePointMeshVariable( fid, meshname, var.name, var.data ); - } else if ( var.precision == IO::DataType::Float ) { - Array data2( var.data.size() ); - data2.copy( var.data ); - IO::silo::writePointMeshVariable( fid, meshname, var.name, data2 ); - } else if ( var.precision == IO::DataType::Int ) { - Array data2( var.data.size() ); - data2.copy( var.data ); - IO::silo::writePointMeshVariable( fid, meshname, var.name, data2 ); - } else { - ERROR( "Unsupported format" ); - } - } -} -// Write a TriMesh mesh (and variables) to a file -template -static void writeSiloTriMesh( DBfile *fid, const IO::TriMesh &mesh, const std::string &meshname ) -{ - const auto &points = mesh.vertices->getPoints(); - std::vector x( points.size() ), y( points.size() ), z( points.size() ); - for ( size_t i = 0; i < x.size(); i++ ) { - x[i] = points[i].x; - y[i] = points[i].y; - z[i] = points[i].z; - } - const TYPE *coords[] = { x.data(), y.data(), z.data() }; - const int *tri[] = { mesh.A.data(), mesh.B.data(), mesh.C.data() }; - IO::silo::writeTriMesh( fid, meshname, 3, 2, points.size(), coords, mesh.A.size(), tri ); -} -static void writeSiloTriMesh2( DBfile *fid, const IO::MeshDataStruct &meshData, - const IO::TriMesh &mesh, IO::MeshDatabase database ) -{ - const std::string meshname = database.domains[0].name; - if ( meshData.precision == IO::DataType::Double ) { - writeSiloTriMesh( fid, mesh, meshname ); - } else if ( meshData.precision == IO::DataType::Float ) { - writeSiloTriMesh( fid, mesh, meshname ); - } else { - ERROR( "Unsupported format" ); - } - for ( size_t i = 0; i < meshData.vars.size(); i++ ) { - const IO::Variable &var = *meshData.vars[i]; - if ( var.precision == IO::DataType::Double ) { - IO::silo::writeTriMeshVariable( fid, 3, meshname, var.name, var.data, var.type ); - } else if ( var.precision == IO::DataType::Float ) { - Array data2( var.data.size() ); - data2.copy( var.data ); - IO::silo::writeTriMeshVariable( fid, 3, meshname, var.name, data2, var.type ); - } else if ( var.precision == IO::DataType::Int ) { - Array data2( var.data.size() ); - data2.copy( var.data ); - IO::silo::writeTriMeshVariable( fid, 3, meshname, var.name, data2, var.type ); - } else { - ERROR( "Unsupported format" ); - } - } -} -static void writeSiloTriMesh( - DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database ) -{ - const IO::TriMesh &mesh = dynamic_cast( *meshData.mesh ); - writeSiloTriMesh2( fid, meshData, mesh, database ); -} -static void writeSiloTriList( - DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database ) -{ - auto mesh = getTriMesh( meshData.mesh ); - writeSiloTriMesh2( fid, meshData, *mesh, database ); -} -// Write a DomainMesh mesh (and variables) to a file -static void writeSiloDomainMesh( - DBfile *fid, const IO::MeshDataStruct &meshData, IO::MeshDatabase database ) -{ - const IO::DomainMesh &mesh = dynamic_cast( *meshData.mesh ); - RankInfoStruct info( mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz ); - std::array range = { info.ix * mesh.Lx / info.nx, - ( info.ix + 1 ) * mesh.Lx / info.nx, info.jy * mesh.Ly / info.ny, - ( info.jy + 1 ) * mesh.Ly / info.ny, info.kz * mesh.Lz / info.nz, - ( info.kz + 1 ) * mesh.Lz / info.nz }; - std::array N = { mesh.nx, mesh.ny, mesh.nz }; - auto meshname = database.domains[0].name; - IO::silo::writeUniformMesh<3>( fid, meshname, range, N ); - IO::silo::write( - fid, meshname + "_rankinfo", { mesh.rank, mesh.nprocx, mesh.nprocy, mesh.nprocz } ); - for ( size_t i = 0; i < meshData.vars.size(); i++ ) { - const auto &var = *meshData.vars[i]; - if ( var.precision == IO::DataType::Double ) { - IO::silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, var.data, var.type ); - } else if ( var.precision == IO::DataType::Float ) { - Array data2( var.data.size() ); - data2.copy( var.data ); - IO::silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, data2, var.type ); - } else if ( var.precision == IO::DataType::Int ) { - Array data2( var.data.size() ); - data2.copy( var.data ); - IO::silo::writeUniformMeshVariable<3>( fid, meshname, N, var.name, data2, var.type ); - } else { - ERROR( "Unsupported format" ); - } - } -} -// Write a mesh (and variables) to a file -static IO::MeshDatabase write_domain_silo( DBfile *fid, const std::string &filename, - const IO::MeshDataStruct &mesh, IO::FileFormat format, int rank ) -{ - // Create the MeshDatabase - auto database = getDatabase( filename, mesh, format, rank ); - if ( database.meshClass == "PointList" ) { - writeSiloPointList( fid, mesh, database ); - } else if ( database.meshClass == "TriMesh" ) { - writeSiloTriMesh( fid, mesh, database ); - } else if ( database.meshClass == "TriList" ) { - writeSiloTriList( fid, mesh, database ); - } else if ( database.meshClass == "DomainMesh" ) { - writeSiloDomainMesh( fid, mesh, database ); - } else { - ERROR( "Unknown mesh class" ); - } - return database; -} -// Write the summary file for silo -std::pair getSiloMeshType( const std::string &meshClass ) -{ - int meshType = 0; - int varType = 0; - if ( meshClass == "PointList" ) { - meshType = DB_POINTMESH; - varType = DB_POINTVAR; - } else if ( meshClass == "TriMesh" ) { - meshType = DB_UCDMESH; - varType = DB_UCDVAR; - } else if ( meshClass == "TriList" ) { - meshType = DB_UCDMESH; - varType = DB_UCDVAR; - } else if ( meshClass == "DomainMesh" ) { - meshType = DB_QUAD_RECT; - varType = DB_QUADVAR; - } else { - ERROR( "Unknown mesh class" ); - } - return std::make_pair( meshType, varType ); -} -void writeSiloSummary( - const std::vector &meshes_written, const std::string &filename ) -{ - auto fid = IO::silo::open( filename, IO::silo::CREATE ); - for ( const auto &data : meshes_written ) { - auto type = getSiloMeshType( data.meshClass ); - std::vector meshTypes( data.domains.size(), type.first ); - std::vector varTypes( data.domains.size(), type.second ); - std::vector meshNames; - for ( const auto &tmp : data.domains ) - meshNames.push_back( tmp.file + ":" + tmp.name ); - IO::silo::writeMultiMesh( fid, data.name, meshNames, meshTypes ); - for ( const auto &variable : data.variables ) { - std::vector varnames; - for ( const auto &tmp : data.domains ) - varnames.push_back( tmp.file + ":" + variable.name ); - IO::silo::writeMultiVar( fid, variable.name, varnames, varTypes ); - } - } - IO::silo::close( fid ); -} -#endif - - // Write the mesh data in the new format static std::vector writeMeshesNewFormat( const std::vector &meshData, const std::string &path, IO::FileFormat format, @@ -459,43 +271,14 @@ static std::vector writeMeshesNewFormat( sprintf( filename, "%05i", rank ); sprintf( fullpath, "%s/%s", path.c_str(), filename ); FILE *fid = fopen( fullpath, "wb" ); - for ( size_t i = 0; i < meshData.size(); i++ ) { - std::shared_ptr mesh = meshData[i].mesh; + ASSERT( fid != nullptr ); + for ( size_t i = 0; i < meshData.size(); i++ ) meshes_written.push_back( write_domain( fid, filename, meshData[i], format, rank ) ); - } fclose( fid ); return meshes_written; } -// Write the mesh data to silo -static std::vector writeMeshesSilo( - const std::vector &meshData, const std::string &path, IO::FileFormat format, - int rank ) -{ -#ifdef USE_SILO - std::vector meshes_written; - char filename[100], fullpath[200]; - sprintf( filename, "%05i.silo", rank ); - sprintf( fullpath, "%s/%s", path.c_str(), filename ); - auto fid = IO::silo::open( fullpath, IO::silo::CREATE ); - for ( size_t i = 0; i < meshData.size(); i++ ) { - auto mesh = meshData[i].mesh; - meshes_written.push_back( write_domain_silo( fid, filename, meshData[i], format, rank ) ); - } - IO::silo::close( fid ); - return meshes_written; -#else - NULL_USE( meshData ); - NULL_USE( path ); - NULL_USE( format ); - NULL_USE( rank ); - ERROR( "Application built without silo support" ); - return std::vector(); -#endif -} - - /**************************************************** * Write the mesh data * ****************************************************/ @@ -513,6 +296,7 @@ void IO::writeData( const std::string &subdir, const std::vector meshes_written; if ( global_IO_format == Format::OLD ) { // Write the original triangle format @@ -523,24 +307,28 @@ void IO::writeData( const std::string &subdir, const std::vector\n" ); + fprintf( xmf, "%s %s\n", indent.data(), location.data() ); + fprintf( xmf, "%s\n", indent.data() ); +} +template +static void addVariable( FILE *xmf, const std::string &indent, const std::string &name, + const std::string &type, const std::string ¢er, ArraySize size, + const std::string &location ) +{ + fprintf( xmf, "%s\n", indent.data(), + name.data(), type.data(), center.data() ); + addDataItem( xmf, indent + " ", size, location ); + fprintf( xmf, "%s\n", indent.data() ); +} + + +/**************************************************************** + * Enum functions * + ****************************************************************/ +/*template +static Xdmf::DataType getDataType() +{ + if ( std::is_same::value ) + return Xdmf::DataType::Char; + else if ( std::is_same::value ) + return Xdmf::DataType::Int32; + else if ( std::is_same::value ) + return Xdmf::DataType::Int64; + else if ( std::is_same::value ) + return Xdmf::DataType::Uint32; + else if ( std::is_same::value ) + return Xdmf::DataType::Uint64; + else if ( std::is_same::value ) + return Xdmf::DataType::Float; + else if ( std::is_same::value ) + return Xdmf::DataType::Double; + else + ERROR( "Invalid type" ); +}*/ +static const char *TopologyTypeNames[] = { "", "Polyvertex", "Polyline", "Polygon", "Triangle", + "Quadrilateral", "Tetrahedron", "Pyramid", "Wedge", "Hexahedron", "Edge_3", "Triangle_6", + "Quadrilateral_8", "Tetrahedron_10", "Pyramid_13", "Wedge_15", "Hexahedron_20", "Mixed", + "CurvilinearMesh2D", "CurvilinearMesh3D", "RectangularMesh2D", "RectangularMesh3D", + "UniformMesh2D", "UniformMesh3D" }; +static const uint8_t TopologyTypeDOFs[] = { 0, 1, 2, 0, 3, 4, 4, 5, 6, 8, 3, 6, 8, 10, 13, 15, 20, + 0, 0, 0, 0, 0, 0, 0 }; + + +/**************************************************************** + * Create a mesh * + ****************************************************************/ +Xdmf::MeshData Xdmf::createPointMesh( const std::string &name, uint8_t NDIM, size_t N, + const std::string &x, const std::string &y, const std::string &z ) +{ + return createUnstructuredMesh( name, NDIM, TopologyType::Polyvertex, N, "", N, x, y, z ); +} +Xdmf::MeshData Xdmf::createUniformMesh( + const std::string &name, const std::vector &range, ArraySize size ) +{ + ASSERT( range.size() == 2 * size.ndim() ); + MeshData data; + data.name = name; + data.size = size; + if ( size.ndim() == 2 ) + data.type = TopologyType::UniformMesh2D; + else if ( size.ndim() == 3 ) + data.type = TopologyType::UniformMesh3D; + else + ERROR( "# of dimensions != 2 or 3" ); + for ( int i = 0; i < 2 * size.ndim(); i++ ) + data.range[i] = range[i]; + return data; +} +Xdmf::MeshData Xdmf::createCurvilinearMesh( const std::string &name, ArraySize size, + const std::string &x, const std::string &y, const std::string &z ) +{ + MeshData data; + data.name = name; + if ( size.ndim() == 2 ) + data.type = TopologyType::CurvilinearMesh2D; + else if ( size.ndim() == 3 ) + data.type = TopologyType::CurvilinearMesh3D; + else + ERROR( "Invalid size for Curvilinear mesh" ); + data.size = size; + data.x = x; + data.y = y; + data.z = z; + return data; +} +Xdmf::MeshData Xdmf::createUnstructuredMesh( const std::string &name, uint8_t NDIM, + TopologyType type, size_t NumElements, const std::string &dofMap, size_t NumNodes, + const std::string &x, const std::string &y, const std::string &z ) +{ + ASSERT( type != TopologyType::Null ); + MeshData data; + data.name = name; + data.type = type; + data.size = { NDIM, NumElements, NumNodes }; + data.dofMap = dofMap; + data.x = x; + data.y = y; + data.z = z; + return data; +} + + +/**************************************************************** + * Add a variable * + ****************************************************************/ +void Xdmf::MeshData::addVariable( const std::string &meshName, const std::string &varName, + ArraySize varSize, RankType rank, Center center, const std::string &varData ) +{ + VarData var; + var.name = varName; + var.size = varSize; + var.data = varData; + var.rankType = rank; + var.center = center; + vars.push_back( std::move( var ) ); +} + + +/**************************************************************** + * Add a mesh domain * + ****************************************************************/ +void Xdmf::addMesh( const std::string &meshName, const MeshData &domain ) +{ + auto &domains = d_meshData[meshName]; + for ( const auto &domain2 : domains ) + ASSERT( domain2.name != domain.name ); + domains.push_back( domain ); +} + + +/**************************************************************** + * Write a variable * + ****************************************************************/ +static void writeVariable( FILE *fid, const Xdmf::VarData &var, const std::string &indent ) +{ + // Write the variable name + fprintf( fid, "%s\n" ); + } else if ( var.center == Xdmf::Center::Cell ) { + fprintf( fid, " Center=\"Cell\">\n" ); + } else if ( var.center == Xdmf::Center::Grid ) { + fprintf( fid, " Center=\"Grid\">\n" ); + } else if ( var.center == Xdmf::Center::Face ) { + fprintf( fid, " Center=\"Face\">\n" ); + } else if ( var.center == Xdmf::Center::Edge ) { + fprintf( fid, " Center=\"Edge\">\n" ); + } else if ( var.center == Xdmf::Center::Other ) { + fprintf( fid, " Center=\"Other\">\n" ); + } else { + ERROR( "Unknown center type" ); + } + // Write the data item + addDataItem( fid, indent + " ", var.size, var.data ); + // Finished + fprintf( fid, "%s\n", indent.data() ); +} + + +/**************************************************************** + * Write the mesh grid * + ****************************************************************/ +static void writeMeshGrid( FILE *fid, const Xdmf::MeshData &mesh, const std::string &indent ) +{ + const char *s = indent.data(); + double x0[3] = { mesh.range[0], mesh.range[2], mesh.range[4] }; + double dx[3] = { ( mesh.range[1] - mesh.range[0] ) / mesh.size[0], + ( mesh.range[3] - mesh.range[2] ) / mesh.size[1], + ( mesh.range[5] - mesh.range[4] ) / mesh.size[2] }; + switch ( mesh.type ) { + case Xdmf::TopologyType::UniformMesh2D: + // Write a uniform 2d mesh + fprintf( fid, "%s\n", s, mesh.name.data() ); + fprintf( fid, + "%s \n", s, + mesh.size[1] + 1, mesh.size[0] + 1 ); + fprintf( fid, "%s \n", s ); + fprintf( + fid, "%s \n", s ); + fprintf( fid, "%s %0.12e %0.12e\n", s, x0[0], x0[1] ); + fprintf( fid, "%s \n", s ); + fprintf( + fid, "%s \n", s ); + fprintf( fid, "%s %0.12e %0.12e\n", s, dx[0], dx[1] ); + fprintf( fid, "%s \n", s ); + fprintf( fid, "%s \n", s ); + break; + case Xdmf::TopologyType::UniformMesh3D: + // Write a uniform 3d mesh + fprintf( fid, "%s\n", s, mesh.name.data() ); + fprintf( fid, + "%s \n", s, + mesh.size[1] + 1, mesh.size[0] + 1 ); + fprintf( fid, "%s \n", s ); + fprintf( + fid, "%s \n", s ); + fprintf( fid, "%s %0.12e %0.12e %0.12e\n", s, x0[0], x0[1], x0[2] ); + fprintf( fid, "%s \n", s ); + fprintf( + fid, "%s \n", s ); + fprintf( fid, "%s %0.12e %0.12e %0.12e\n", s, dx[0], dx[1], dx[2] ); + fprintf( fid, "%s \n", s ); + fprintf( fid, "%s \n", s ); + break; + case Xdmf::TopologyType::CurvilinearMesh2D: + // Write a 2D curvillinear mesh + fprintf( fid, "%s\n", s, mesh.name.data() ); + fprintf( fid, "%s \n", s, + mesh.size[1] + 1, mesh.size[0] + 1 ); + fprintf( fid, "%s \n", s ); + addDataItem( fid, indent + " ", mesh.size + 1, mesh.x ); + addDataItem( fid, indent + " ", mesh.size + 1, mesh.y ); + fprintf( fid, "%s \n", s ); + break; + case Xdmf::TopologyType::CurvilinearMesh3D: + // Write a 3D curvillinear mesh + fprintf( fid, "%s\n", s, mesh.name.data() ); + fprintf( fid, "%s \n", + s, mesh.size[2] + 1, mesh.size[1] + 1, mesh.size[0] + 1 ); + fprintf( fid, "%s \n", s ); + addDataItem( fid, indent + " ", mesh.size + 1, mesh.x ); + addDataItem( fid, indent + " ", mesh.size + 1, mesh.y ); + addDataItem( fid, indent + " ", mesh.size + 1, mesh.z ); + fprintf( fid, "%s \n", s ); + break; + case Xdmf::TopologyType::Polyvertex: + case Xdmf::TopologyType::Polyline: + case Xdmf::TopologyType::Polygon: + case Xdmf::TopologyType::Triangle: + case Xdmf::TopologyType::Quadrilateral: + case Xdmf::TopologyType::Tetrahedron: + case Xdmf::TopologyType::Pyramid: + case Xdmf::TopologyType::Wedge: + case Xdmf::TopologyType::Hexahedron: + case Xdmf::TopologyType::Edge_3: + case Xdmf::TopologyType::Triangle_6: + case Xdmf::TopologyType::Quadrilateral_8: + case Xdmf::TopologyType::Tetrahedron_10: + case Xdmf::TopologyType::Pyramid_13: + case Xdmf::TopologyType::Wedge_15: + case Xdmf::TopologyType::Hexahedron_20: + // Write an unstructured mesh + { + int NDIM = mesh.size[0]; + size_t Nelem = mesh.size[1]; + size_t Nnode = mesh.size[2]; + uint8_t Ndofs = TopologyTypeDOFs[static_cast( mesh.type )]; + auto type = TopologyTypeNames[static_cast( mesh.type )]; + fprintf( fid, "%s\n", s, mesh.name.data() ); + fprintf( fid, "%s \n", Nelem ); + if ( !mesh.dofMap.empty() ) + addDataItem( fid, indent + " ", { Ndofs, Nelem }, mesh.dofMap ); + fprintf( fid, "%s \n", s ); + if ( NDIM == 2 ) { + if ( mesh.y.empty() ) { + fprintf( fid, "%s \n", s ); + addDataItem( fid, indent + " ", { 2, Nnode }, mesh.x ); + } else { + fprintf( fid, "%s \n", s ); + addDataItem( fid, indent + " ", Nnode, mesh.x ); + addDataItem( fid, indent + " ", Nnode, mesh.y ); + } + } else if ( NDIM == 3 ) { + if ( mesh.y.empty() ) { + fprintf( fid, "%s \n", s ); + addDataItem( fid, indent + " ", { 2, Nnode }, mesh.x ); + } else { + fprintf( fid, "%s \n", s ); + addDataItem( fid, indent + " ", Nnode, mesh.x ); + addDataItem( fid, indent + " ", Nnode, mesh.y ); + addDataItem( fid, indent + " ", Nnode, mesh.z ); + } + } else { + ERROR( "Dimensions other than 2 or 3 are not supported" ); + } + fprintf( fid, "%s \n", s ); + } + break; + default: { + auto msg = "Invalid mesh type: " + std::to_string( static_cast( mesh.type ) ) + " - " + + mesh.name; + ERROR( msg ); + } + } + // Write the variables + for ( const auto &var : mesh.vars ) + writeVariable( fid, var, indent + " " ); + fprintf( fid, "%s \n", s ); +} + + +/**************************************************************** + * Write the XDMF xml file * + ****************************************************************/ +void Xdmf::write( const std::string &filename ) const +{ + // Create XDMF file + auto fid = fopen( filename.data(), "w" ); + fprintf( fid, "\n" ); + fprintf( fid, "\n" ); + fprintf( fid, "\n" ); + fprintf( fid, "\n" ); + // Write an empty mesh to enable collections to work properly + fprintf( fid, " \n\n" ); + // Write each mesh + for ( const auto &data : d_meshData ) { + auto name = data.first; + auto domains = data.second; + if ( domains.empty() ) + continue; + if ( domains.size() == 1u && name == domains[0].name ) { + writeMeshGrid( fid, domains[0], " " ); + } else { + fprintf( fid, " \n", name.data() ); + for ( const auto &domain : domains ) + writeMeshGrid( fid, domain, " " ); + fprintf( fid, " \n\n" ); + } + } + fprintf( fid, "\n" ); + fprintf( fid, "\n" ); + fclose( fid ); +} + + +/**************************************************************** + * Pack/Unpack data * + ****************************************************************/ +template +typename std::enable_if::value, size_t>::type size( const T & ) +{ + return sizeof( T ); +} +template +typename std::enable_if::value, char *>::type pack( + char *ptr, const T &x ) +{ + memcpy( ptr, &x, sizeof( T ) ); + return ptr + sizeof( T ); +} +template +typename std::enable_if::value, char *>::type unpack( + char *ptr, T &x ) +{ + memcpy( &x, ptr, sizeof( T ) ); + return ptr + sizeof( T ); +} +static size_t size( const std::string &str ) { return sizeof( int ) + str.size(); } +static char *pack( char *ptr, const std::string &str ) +{ + int N = str.size(); + memcpy( ptr, &N, sizeof( int ) ); + ptr += sizeof( int ); + memcpy( ptr, str.data(), str.size() ); + ptr += str.size(); + return ptr; +} +static char *unpack( char *ptr, std::string &str ) +{ + int N = 0; + memcpy( &N, ptr, sizeof( int ) ); + ASSERT( N >= 0 && N < 1000 ); + ptr += sizeof( int ); + str = std::string( ptr, N ); + ptr += N; + return ptr; +} +static size_t size( const Xdmf::VarData &data ) +{ + size_t bytes = 0; + bytes += size( data.name ); + bytes += size( data.size ); + bytes += size( data.rankType ); + bytes += size( data.center ); + bytes += size( data.data ); + return bytes; +} +static char *pack( char *ptr, const Xdmf::VarData &data ) +{ + ptr = pack( ptr, data.name ); + ptr = pack( ptr, data.size ); + ptr = pack( ptr, data.rankType ); + ptr = pack( ptr, data.center ); + ptr = pack( ptr, data.data ); + return ptr; +} +static char *unpack( char *ptr, Xdmf::VarData &data ) +{ + int rankType = 0, center = 0; + ptr = unpack( ptr, data.name ); + ptr = unpack( ptr, data.size ); + ptr = unpack( ptr, rankType ); + ptr = unpack( ptr, center ); + ptr = unpack( ptr, data.data ); + data.rankType = static_cast( rankType ); + data.center = static_cast( center ); + return ptr; +} +static size_t size( const Xdmf::MeshData &data ) +{ + int N_vars = data.vars.size(); + size_t bytes = 0; + bytes += size( data.name ); + bytes += size( data.type ); + bytes += size( data.size ); + bytes += size( data.range ); + bytes += size( data.x ); + bytes += size( data.y ); + bytes += size( data.z ); + bytes += size( N_vars ); + for ( int i = 0; i < N_vars; i++ ) + bytes += size( data.vars[i] ); + return bytes; +} +static char *pack( char *ptr, const Xdmf::MeshData &data ) +{ + int N_vars = data.vars.size(); + ptr = pack( ptr, data.name ); + ptr = pack( ptr, data.type ); + ptr = pack( ptr, data.size ); + ptr = pack( ptr, data.range ); + ptr = pack( ptr, data.x ); + ptr = pack( ptr, data.y ); + ptr = pack( ptr, data.z ); + ptr = pack( ptr, N_vars ); + for ( int i = 0; i < N_vars; i++ ) + ptr = pack( ptr, data.vars[i] ); + return ptr; +} +static char *unpack( char *ptr, Xdmf::MeshData &data ) +{ + int N_vars = 0; + ptr = unpack( ptr, data.name ); + ptr = unpack( ptr, data.type ); + ptr = unpack( ptr, data.size ); + ptr = unpack( ptr, data.range ); + ptr = unpack( ptr, data.x ); + ptr = unpack( ptr, data.y ); + ptr = unpack( ptr, data.z ); + ptr = unpack( ptr, N_vars ); + data.vars.resize( N_vars ); + for ( int i = 0; i < N_vars; i++ ) + ptr = unpack( ptr, data.vars[i] ); + return ptr; +} +static size_t size( const std::vector &data ) +{ + size_t bytes = 0; + int N = data.size(); + bytes += size( N ); + for ( int i = 0; i < N; i++ ) + bytes += size( data[i] ); + return bytes; +} +static char *pack( char *ptr, const std::vector &data ) +{ + int N = data.size(); + ptr = pack( ptr, N ); + for ( int i = 0; i < N; i++ ) + ptr = pack( ptr, data[i] ); + return ptr; +} +static char *unpack( char *ptr, std::vector &data ) +{ + data.clear(); + int N = data.size(); + ptr = unpack( ptr, N ); + data.resize( N ); + for ( int i = 0; i < N; i++ ) + ptr = unpack( ptr, data[i] ); + return ptr; +} +static size_t size( const std::map> &data ) +{ + size_t bytes = 0; + int N_map = data.size(); + bytes += size( N_map ); + for ( const auto &tmp : data ) { + bytes += size( tmp.first ); + bytes += size( tmp.second ); + } + return bytes; +} +static char *pack( char *ptr, const std::map> &data ) +{ + int N_map = data.size(); + ptr = pack( ptr, N_map ); + for ( const auto &tmp : data ) { + ptr = pack( ptr, tmp.first ); + ptr = pack( ptr, tmp.second ); + } + return ptr; +} +static char *unpack( char *ptr, std::map> &data ) +{ + data.clear(); + int N_map = data.size(); + ptr = unpack( ptr, N_map ); + for ( int i = 0; i < N_map; i++ ) { + std::string name; + std::vector data2; + ptr = unpack( ptr, name ); + ptr = unpack( ptr, data2 ); + data[name] = std::move( data2 ); + } + return ptr; +} + + +/**************************************************************** + * Gather all data to rank 0 * + ****************************************************************/ +void Xdmf::gather( const Utilities::MPI &comm ) +{ + if ( comm.getRank() == 0 ) { + for ( int i = 1; i < comm.getSize(); i++ ) { + // Recieve the data + size_t N_meshes = 0, N_bytes = 0; + comm.recv( &N_meshes, 1, i, 717 ); + comm.recv( &N_bytes, 1, i, 718 ); + auto buf = new char[N_bytes]; + comm.recv( buf, N_bytes, i, 719 ); + // Unpack the data + std::map> data; + unpack( buf, data ); + delete[] buf; + // Add the meshes + for ( auto tmp : data ) { + const auto &name = tmp.first; + const auto &domains = tmp.second; + if ( domains.size() == 1u && domains[0].name == name ) { + // We are dealing with a single mesh + ASSERT( d_meshData.find( name ) == d_meshData.end() ); + d_meshData.insert( tmp ); + } else { + // Add the domains + auto &meshes = d_meshData[name]; + for ( auto domain : domains ) { + for ( const auto &tmp : meshes ) + ASSERT( tmp.name != domain.name ); + meshes.push_back( domain ); + } + } + } + } + } else { + // Send the number of meshes + size_t N_meshes = d_meshData.size(); + comm.send( &N_meshes, 1, 0, 717 ); + // Pack the send data + size_t N_bytes = size( d_meshData ); + comm.send( &N_bytes, 1, 0, 718 ); + auto buf = new char[N_bytes]; + pack( buf, d_meshData ); + // Send the data to rank 0 + comm.send( buf, N_bytes, 0, 719 ); + delete[] buf; + // Clear the internal data + d_meshData.clear(); + } +} diff --git a/IO/Xdmf.h b/IO/Xdmf.h new file mode 100644 index 00000000..c7d39801 --- /dev/null +++ b/IO/Xdmf.h @@ -0,0 +1,130 @@ +#include "IO/HDF5_IO.h" +#include "common/Array.h" +#include "common/MPI.h" +#include "common/Utilities.h" + +#include +#include +#include +#include + + +// Helper class to write/read XDMF files +class Xdmf +{ +public: + enum class TopologyType { + Null = 0, + Polyvertex, + Polyline, + Polygon, + Triangle, + Quadrilateral, + Tetrahedron, + Pyramid, + Wedge, + Hexahedron, + Edge_3, + Triangle_6, + Quadrilateral_8, + Tetrahedron_10, + Pyramid_13, + Wedge_15, + Hexahedron_20, + Mixed, + CurvilinearMesh2D, + CurvilinearMesh3D, + RectangularMesh2D, + RectangularMesh3D, + UniformMesh2D, + UniformMesh3D, + }; + enum class DataType { Null = 0, Char, Int32, Int64, Uint32, Uint64, Float, Double }; + enum class RankType { Null = 0, Scalar, Vector, Tensor, Tensor6, Matrix, GlobalID }; + enum class Center { Null = 0, Node, Edge, Face, Cell, Grid, Other }; + + struct VarData { + std::string name; // Variable name + ArraySize size; // Size of variable + RankType rankType; // Rank order of data + Center center; // Variable centering + std::string data; // Variable data + }; + + struct MeshData { + std::string name; // Name of mesh + TopologyType type; // Type of mesh + ArraySize size; // Size of mesh (meaning depends on mesh type) + double range[6]; // Range of the mesh (only used for UniformMesh2D/UniformMesh3D) + std::string x; // x coordinates (or xy/xyz coordinates) + std::string y; // y coordinates + std::string z; // z coordinates + std::string dofMap; // mesh connectivity + std::vector vars; // Variables + MeshData() : type( TopologyType::Null ), range{ 0 } {} + //! Add a variable + void addVariable( const std::string &meshName, const std::string &varName, + ArraySize varSize, RankType rank, Center center, const std::string &varData ); + }; + + +public: + //! Add a Point mesh + static MeshData createPointMesh( const std::string &name, uint8_t NDIM, size_t N, + const std::string &x, const std::string &y = "", const std::string &z = "" ); + + /*! + * @brief Add a uniform mesh + * @details This function adds a uniform rectangular mesh + * @param[in] name The name of the mesh + * @param[in] range The range of the mesh [ x_min, x_max, y_min, y_max, z_min, z_max ] + * @param[in] size The number of cells in the mesh + */ + static MeshData createUniformMesh( + const std::string &name, const std::vector &range, ArraySize size ); + + /*! + * @brief Add a Curvilinear mesh + * @details This function adds a curvilinear mesh + * @param[in] name The name of the mesh + * @param[in] size The number of cells in the mesh + * @param[in] x The x coordinates or the xy/xyz coordinates + * @param[in] y The y coordinates (may be null) + * @param[in] z The z coordinates (may be null) + */ + static MeshData createCurvilinearMesh( const std::string &name, ArraySize size, + const std::string &x, const std::string &y, const std::string &z = "" ); + + /*! + * @brief Add an unstructured mesh + * @details This function adds an unstructerd mesh to the class to write. + * The mesh may be one of several unsupported unstructured mesh types. + * This function does not support mixed elements. + * @param[in] name The name of the mesh + * @param[in] NDIM The number of physical dimensions + * @param[in] type The element type + * @param[in] NumElements The number of elements + * @param[in] dofMap The connectivity information (type x NumElements) + * @param[in] x The x coordinates or the xy/xyz coordinates + * @param[in] y The y coordinates (may be null) + * @param[in] z The z coordinates (may be null) + */ + static MeshData createUnstructuredMesh( const std::string &name, uint8_t NDIM, + TopologyType type, size_t NumElements, const std::string &dofMap, size_t NumNodes, + const std::string &x, const std::string &y = "", const std::string &z = "" ); + + +public: + //! Add a sub-domain + void addMesh( const std::string &meshName, const MeshData &domain ); + + //! Gather all data to rank 0 + void gather( const Utilities::MPI &comm ); + + //! Write the xml file + void write( const std::string &filename ) const; + + +private: + std::map> d_meshData; +}; diff --git a/IO/netcdf.cpp b/IO/netcdf.cpp index 06f41dba..12347709 100644 --- a/IO/netcdf.cpp +++ b/IO/netcdf.cpp @@ -496,7 +496,7 @@ template void write( int fid, const std::string &var, const std::vector< const Array &data, const RankInfoStruct &info ); -}; // namespace netcdf +} // namespace netcdf #else diff --git a/IO/netcdf.h b/IO/netcdf.h index eb77784d..5fcc5d82 100644 --- a/IO/netcdf.h +++ b/IO/netcdf.h @@ -138,5 +138,7 @@ void write( int fid, const std::string &var, const std::vector &dimids, const Array &data, const RankInfoStruct &rank_info ); -}; // namespace netcdf +} // namespace netcdf + + #endif diff --git a/IO/silo.cpp b/IO/silo.cpp index a9041e44..1a264eb6 100644 --- a/IO/silo.cpp +++ b/IO/silo.cpp @@ -10,7 +10,8 @@ #include -namespace IO::silo { +namespace IO { +namespace silo { /**************************************************** @@ -99,7 +100,8 @@ void writeMultiVar( DBfile *fid, const std::string &varname, } -}; // namespace IO::silo +} // namespace silo +} // namespace IO #else diff --git a/IO/silo.hpp b/IO/silo.hpp index 72ff7041..d86ad17f 100644 --- a/IO/silo.hpp +++ b/IO/silo.hpp @@ -245,11 +245,11 @@ Array readUniformMeshVariable( DBfile *fid, const std::string &varname ) copyData( data2, type, var->vals[i] ); memcpy( &data( 0, i ), data2.data(), var->nels * sizeof( TYPE ) ); } - DBFreeQuadvar( var ); std::vector dims( var->ndims + 1, var->nvals ); for ( int d = 0; d < var->ndims; d++ ) dims[d] = var->dims[d]; data.reshape( dims ); + DBFreeQuadvar( var ); return data; } diff --git a/StackTrace/StackTrace.cpp b/StackTrace/StackTrace.cpp index 55a24352..b288ed75 100644 --- a/StackTrace/StackTrace.cpp +++ b/StackTrace/StackTrace.cpp @@ -856,7 +856,7 @@ static void getFileAndLineObject( staticVectorfilename, info[i]->filenamePath ); diff --git a/common/Array.cpp b/common/Array.cpp index ecb3470a..5db23505 100644 --- a/common/Array.cpp +++ b/common/Array.cpp @@ -1,117 +1,74 @@ +// clang-format off #include "common/Array.h" #include "common/Array.hpp" +#include "common/Utilities.h" #include -/******************************************************** - * ArraySize * - ********************************************************/ -ArraySize::ArraySize( const std::vector& N ) -{ - d_ndim = N.size(); - d_N[0] = 0; - d_N[1] = 1; - d_N[2] = 1; - d_N[3] = 1; - d_N[4] = 1; - for ( size_t i = 0; i < d_ndim; i++ ) - d_N[i] = N[i]; - d_length = 1; - for ( unsigned long i : d_N ) - d_length *= i; - if ( d_ndim == 0 ) - d_length = 0; -} - - /******************************************************** * Explicit instantiations of Array * ********************************************************/ -template class Array; -template class Array; -template class Array; -template class Array; -template class Array; -template class Array; -template class Array; -template class Array; -template class Array; -template class Array; -template class Array; +template class Array; +template class Array; +template class Array; +template class Array; +template class Array; +template class Array; +template class Array; +template class Array; +template class Array; +template class Array; +template class Array; +template class Array; /******************************************************** * Explicit instantiations of Array * ********************************************************/ -// clang-format off -template Array::Array(); -template Array::~Array(); -template Array::Array( size_t ); -template Array::Array( size_t, size_t ); -template Array::Array( size_t, size_t, size_t ); -template Array::Array( size_t, size_t, size_t, size_t ); -template Array::Array( size_t, size_t, size_t, size_t, size_t ); -template Array::Array( const std::vector&, const bool* ); -template Array::Array( std::string ); -template Array::Array( std::initializer_list ); -template Array::Array( const Array& ); -template Array::Array( Array&& ); -template Array& Array::operator=( const Array& ); -template Array& Array::operator=( Array&& ); -template Array& Array::operator=( const std::vector& ); -template void Array::fill(bool const&); -template void Array::clear(); -template bool Array::operator==(Array const&) const; -template void Array::resize( ArraySize const& ); -// clang-format on +instantiateArrayConstructors( bool ) +template Array& Array::operator=( const std::vector& ); +template void Array::clear(); +template bool Array::operator==(Array const&) const; +template void Array::resize( ArraySize const& ); /******************************************************** * Explicit instantiations of Array * ********************************************************/ -// clang-format off -template Array, FunctionTable>::Array(); -template Array, FunctionTable>::~Array(); -template Array, FunctionTable>::Array( size_t ); -template Array, FunctionTable>::Array( size_t, size_t ); -template Array, FunctionTable>::Array( size_t, size_t, size_t ); -template Array, FunctionTable>::Array( size_t, size_t, size_t, size_t ); -template Array, FunctionTable>::Array( size_t, size_t, size_t, size_t, size_t ); -template Array, FunctionTable>::Array( const std::vector&, const std::complex* ); -template Array, FunctionTable>::Array( std::initializer_list> ); -template Array, FunctionTable>::Array( const Range>& range ); -template Array, FunctionTable>::Array( const Array, FunctionTable>& ); -template Array, FunctionTable>::Array( Array, FunctionTable>&& ); -template Array, FunctionTable>& Array, FunctionTable>::operator=( const Array, FunctionTable>& ); -template Array, FunctionTable>& Array, FunctionTable>::operator=( Array, FunctionTable>&& ); -template Array, FunctionTable>& Array, FunctionTable>::operator=( const std::vector>& ); -template void Array, FunctionTable>::resize( ArraySize const& ); -template void Array, FunctionTable>::viewRaw( ArraySize const&, std::complex*, bool, bool ); -template void Array, FunctionTable>::fill(std::complex const&); -template void Array, FunctionTable>::clear(); -template bool Array, FunctionTable>::operator==(Array, FunctionTable> const&) const; -template Array, FunctionTable> Array, FunctionTable>::repmat(std::vector > const&) const; -// clang-format on +instantiateArrayConstructors( std::complex ) +instantiateArrayConstructors( std::complex ) +template void Array,FunctionTable>::resize( ArraySize const& ); +template void Array,FunctionTable>::resize( ArraySize const& ); +template Array,FunctionTable>& Array,FunctionTable>::operator=(std::vector> const&); +template Array,FunctionTable>& Array,FunctionTable>::operator=(std::vector> const&); +template void Array,FunctionTable>::clear(); +template void Array,FunctionTable>::clear(); +template bool Array,FunctionTable>::operator==(Array,FunctionTable> const&) const; +template bool Array,FunctionTable>::operator==(Array,FunctionTable> const&) const; +template Array,FunctionTable> Array,FunctionTable>::repmat(std::vector const&) const; +template Array,FunctionTable> Array,FunctionTable>::repmat(std::vector const&) const; +template void Array,FunctionTable>::copySubset(std::vector const&, Array,FunctionTable> const&); +template void Array,FunctionTable>::copySubset(std::vector const&, Array,FunctionTable> const&); +template Array,FunctionTable> Array,FunctionTable>::subset(std::vector const&) const; +template Array,FunctionTable> Array,FunctionTable>::subset(std::vector const&) const; +template bool Array,FunctionTable>::NaNs() const; +template bool Array,FunctionTable>::NaNs() const; /******************************************************** * Explicit instantiations of Array * ********************************************************/ -// clang-format off -template Array::Array(); -template Array::~Array(); -template Array::Array( size_t ); -template Array::Array( size_t, size_t ); -template Array::Array( size_t, size_t, size_t ); -template Array::Array( size_t, size_t, size_t, size_t ); -template Array::Array( size_t, size_t, size_t, size_t, size_t ); -template Array::Array( const std::vector&, const std::string* ); -template Array::Array( std::initializer_list ); -template Array::Array( const Array& ); -template Array::Array( Array&& ); -template Array& Array::operator=( const Array& ); -template Array& Array::operator=( Array&& ); -template Array& Array::operator=( const std::vector& ); -template void Array::resize( ArraySize const& ); -// clang-format on +instantiateArrayConstructors( std::string ) +template void Array::resize( ArraySize const& ); +template void Array::clear(); +template Array &Array:: +operator=( const std::vector & ); +template bool Array::operator==(Array const&) const; + + +#if defined( USING_ICC ) +ENABLE_WARNINGS +#endif + + diff --git a/common/Array.h b/common/Array.h index 2dccfde8..ccbcbdc8 100644 --- a/common/Array.h +++ b/common/Array.h @@ -8,6 +8,7 @@ #include #include #include +#include #include @@ -74,13 +75,7 @@ public: // Constructors / assignment operators * @param N Number of elements in each dimension * @param data Optional raw array to copy the src data */ - explicit Array( const std::vector &N, const TYPE *data = NULL ); - - /*! - * Create a 1D Array with the range - * @param range Range of the data - */ - explicit Array( const Range &range ); + explicit Array( const std::vector &N, const TYPE *data = nullptr ); /*! * Create a 1D Array using a string that mimic's MATLAB @@ -94,6 +89,12 @@ public: // Constructors / assignment operators */ Array( std::initializer_list data ); + /*! + * Create a 2D Array with the given initializer lists + * @param data Input data + */ + Array( std::initializer_list> data ); + /*! * Copy constructor @@ -144,7 +145,7 @@ public: // Views/copies/subset * @param N Number of elements in each dimension * @param data Pointer to the data */ - static std::unique_ptr view( const ArraySize &N, std::shared_ptr &data ); + static std::unique_ptr view( const ArraySize &N, std::shared_ptr data ); /*! @@ -152,8 +153,8 @@ public: // Views/copies/subset * @param N Number of elements in each dimension * @param data Pointer to the data */ - static std::unique_ptr constView( - const ArraySize &N, std::shared_ptr const &data ); + static std::unique_ptr constView( const ArraySize &N, + std::shared_ptr const &data ); /*! @@ -167,7 +168,7 @@ public: // Views/copies/subset * @param N Number of elements in each dimension * @param data Pointer to the data */ - void view2( const ArraySize &N, std::shared_ptr const &data ); + void view2( const ArraySize &N, std::shared_ptr data ); /*! * Make this object a view of the raw data (expert use only). @@ -202,14 +203,30 @@ public: // Views/copies/subset */ void viewRaw( const ArraySize &N, TYPE *data, bool isCopyable = true, bool isFixedSize = true ); + /*! + * Create an array view of the given data (expert use only). + * Use view2( N, shared_ptr(data,[](TYPE*){}) ) instead. + * Note: this interface is not recommended as it does not protect from + * the src data being deleted while still being used by the Array. + * Additionally for maximum performance it does not set the internal shared_ptr + * so functions like getPtr and resize will not work correctly. + * @param N Number of elements in each dimension + * @param data Pointer to the data + */ + static inline Array staticView( const ArraySize &N, TYPE *data ) + { + Array x; + x.viewRaw( N, data, true, true ); + return x; + } /*! * Convert an array of one type to another. This may or may not allocate new memory. * @param array Input array */ template - static inline std::unique_ptr> convert( - std::shared_ptr> array ) + static inline std::unique_ptr> + convert( std::shared_ptr> array ) { auto array2 = std::make_unique>( array->size() ); array2.copy( *array ); @@ -222,8 +239,8 @@ public: // Views/copies/subset * @param array Input array */ template - static inline std::unique_ptr> convert( - std::shared_ptr> array ) + static inline std::unique_ptr> + convert( std::shared_ptr> array ) { auto array2 = std::make_unique>( array->size() ); array2.copy( *array ); @@ -235,8 +252,8 @@ public: // Views/copies/subset * Copy and convert data from another array to this array * @param array Source array */ - template - void inline copy( const Array &array ) + template + void inline copy( const Array &array ) { resize( array.size() ); copy( array.data() ); @@ -245,51 +262,55 @@ public: // Views/copies/subset /*! * Copy and convert data from a raw vector to this array. * Note: The current array must be allocated to the proper size first. - * @param array Source array + * @param data Source data */ template - void inline copy( const TYPE2 *data ) - { - for ( size_t i = 0; i < d_size.length(); i++ ) - d_data[i] = static_cast( data[i] ); - } + inline void copy( const TYPE2 *data ); /*! * Copy and convert data from this array to a raw vector. - * @param array Source array + * @param data Source data */ template - void inline copyTo( TYPE2 *data ) const - { - for ( size_t i = 0; i < d_size.length(); i++ ) - data[i] = static_cast( d_data[i] ); - } + inline void copyTo( TYPE2 *data ) const; /*! * Copy and convert data from this array to a new array */ template - Array inline cloneTo() const + Array> inline cloneTo() const { - Array dst( this->size() ); + Array> dst( this->size() ); copyTo( dst.data() ); return dst; } + /*! swap the raw data pointers for the Arrays after checking for compatibility */ void swap( Array &other ); + /*! * Fill the array with the given value - * @param value Value to fill + * @param y Value to fill */ - void fill( const TYPE &value ); + inline void fill( const TYPE &y ) + { + for ( auto &x : *this ) + x = y; + } /*! * Scale the array by the given value - * @param scale Value to scale by + * @param y Value to scale by */ - void scale( const TYPE &scale ); + template + inline void scale( const TYPE2 &y ) + { + for ( auto &x : *this ) + x *= y; + } + /*! * Set the values of this array to pow(base, exp) @@ -298,6 +319,7 @@ public: // Views/copies/subset */ void pow( const Array &base, const TYPE &exp ); + //! Destructor ~Array(); @@ -326,6 +348,10 @@ public: // Views/copies/subset inline bool empty() const { return d_size.length() == 0; } + //! Return true if the Array is not empty + inline operator bool() const { return d_size.length() != 0; } + + /*! * Resize the Array * @param N NUmber of elements @@ -371,6 +397,12 @@ public: // Views/copies/subset void reshape( const ArraySize &N ); + /*! + * Remove singleton dimensions. + */ + void squeeze(); + + /*! * Reshape the Array so that the number of dimensions is the * max of ndim and the largest dim>1. @@ -499,8 +531,8 @@ public: // Accessors * @param i3 The third index * @param i4 The fourth index */ - ARRAY_ATTRIBUTE inline const TYPE &operator()( - size_t i1, size_t i2, size_t i3, size_t i4 ) const + ARRAY_ATTRIBUTE inline const TYPE & + operator()( size_t i1, size_t i2, size_t i3, size_t i4 ) const { return d_data[d_size.index( i1, i2, i3, i4 )]; } @@ -526,8 +558,8 @@ public: // Accessors * @param i4 The fourth index * @param i5 The fifth index */ - ARRAY_ATTRIBUTE inline const TYPE &operator()( - size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) const + ARRAY_ATTRIBUTE inline const TYPE & + operator()( size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) const { return d_data[d_size.index( i1, i2, i3, i4, i5 )]; } @@ -600,6 +632,12 @@ public: // Math operations //! Concatenates the arrays along the dimension dim. static Array cat( const std::vector &x, int dim = 0 ); + //! Concatenates the arrays along the dimension dim. + static Array cat( const std::initializer_list &x, int dim = 0 ); + + //! Concatenates the arrays along the dimension dim. + static Array cat( size_t N_array, const Array *x, int dim ); + //! Concatenates a given array with the current array void cat( const Array &x, int dim = 0 ); @@ -655,20 +693,37 @@ public: // Math operations TYPE mean( const std::vector> &index ) const; //! Find all elements that match the operator - std::vector find( - const TYPE &value, std::function compare ) const; + std::vector find( const TYPE &value, + std::function compare ) const; //! Print an array - void print( - std::ostream &os, const std::string &name = "A", const std::string &prefix = "" ) const; - - //! Multiply two arrays - static Array multiply( const Array &a, const Array &b ); + void + print( std::ostream &os, const std::string &name = "A", const std::string &prefix = "" ) const; //! Transpose an array Array reverseDim() const; + /*! + * @brief Shift dimensions + * @details Shifts the dimensions of the array by N. When N is positive, + * shiftDim shifts the dimensions to the left and wraps the + * N leading dimensions to the end. When N is negative, + * shiftDim shifts the dimensions to the right and pads with singletons. + * @param N Desired shift + */ + Array shiftDim( int N ) const; + + /*! + * @brief Permute array dimensions + * @details Rearranges the dimensions of the array so that they + * are in the order specified by the vector index. + * The array produced has the same values as A but the order of the subscripts + * needed to access any particular element are rearranged as specified. + * @param index Desired order of the subscripts + */ + Array permute( const std::vector &index ) const; + //! Replicate an array a given number of times in each direction Array repmat( const std::vector &N ) const; @@ -676,8 +731,8 @@ public: // Math operations Array coarsen( const Array &filter ) const; //! Coarsen an array using the given filter - Array coarsen( - const std::vector &ratio, std::function filter ) const; + Array coarsen( const std::vector &ratio, + std::function filter ) const; /*! * Perform a element-wise operation y = f(x) @@ -692,8 +747,9 @@ public: // Math operations * @param[in] x The first array * @param[in] y The second array */ - static Array transform( - std::function fun, const Array &x, const Array &y ); + static Array transform( std::function fun, + const Array &x, + const Array &y ); /*! * axpby operation: this = alpha*x + beta*this @@ -707,7 +763,13 @@ public: // Math operations * Linear interpolation * @param[in] x Position as a decimal index */ - TYPE interp( const std::vector &x ) const; + inline TYPE interp( const std::vector &x ) const { return interp( x.data() ); } + + /*! + * Linear interpolation + * @param[in] x Position as a decimal index + */ + TYPE interp( const double *x ) const; /** * \fn equals (Array & const rhs, TYPE tol ) @@ -730,8 +792,10 @@ private: inline void checkSubsetIndex( const std::vector> &range ) const; inline std::vector> convert( const std::vector &index ) const; static inline void getSubsetArrays( const std::vector> &range, - std::array &first, std::array &last, std::array &inc, - std::array &N ); + std::array &first, + std::array &last, + std::array &inc, + std::array &N ); }; @@ -756,8 +820,8 @@ inline Array operator+( const Array &a, const Array &b ) { Array c; - const auto &fun = []( const TYPE &a, const TYPE &b ) { return a + b; }; - FUN::transform( fun, a, b, c ); + const auto &op = []( const TYPE &a, const TYPE &b ) { return a + b; }; + FUN::transform( op, a, b, c ); return c; } template @@ -765,30 +829,78 @@ inline Array operator-( const Array &a, const Array &b ) { Array c; - const auto &fun = []( const TYPE &a, const TYPE &b ) { return a - b; }; - FUN::transform( fun, a, b, c ); + const auto &op = []( const TYPE &a, const TYPE &b ) { return a - b; }; + FUN::transform( op, a, b, c ); return c; } template inline Array operator*( const Array &a, const Array &b ) { - return Array::multiply( a, b ); + Array c; + FUN::multiply( a, b, c ); + return c; } template inline Array operator*( const Array &a, const std::vector &b ) { - Array b2; + Array b2, c; b2.viewRaw( { b.size() }, const_cast( b.data() ) ); - return Array::multiply( a, b2 ); + FUN::multiply( a, b2, c ); + return c; +} +template +inline Array operator*( const TYPE &a, + const Array &b ) +{ + auto c = b; + c.scale( a ); + return c; +} +template +inline Array operator*( const Array &a, + const TYPE &b ) +{ + auto c = a; + c.scale( b ); + return c; +} + + +/******************************************************** + * Copy array * + ********************************************************/ +template +template +inline void Array::copy( const TYPE2 *data ) +{ + if ( std::is_same::value ) { + std::copy( data, data + d_size.length(), d_data ); + } else { + for ( size_t i = 0; i < d_size.length(); i++ ) + d_data[i] = static_cast( data[i] ); + } +} +template +template +inline void Array::copyTo( TYPE2 *data ) const +{ + if ( std::is_same::value ) { + std::copy( d_data, d_data + d_size.length(), data ); + } else { + for ( size_t i = 0; i < d_size.length(); i++ ) + data[i] = static_cast( d_data[i] ); + } } /******************************************************** * Convience typedefs * + * Copy array * ********************************************************/ typedef Array DoubleArray; typedef Array IntArray; + #endif diff --git a/common/Array.hpp b/common/Array.hpp index 2ad07ea4..1faaba71 100644 --- a/common/Array.hpp +++ b/common/Array.hpp @@ -17,65 +17,46 @@ /******************************************************** * External instantiations * ********************************************************/ -extern template class Array; -extern template class Array; -extern template class Array; -extern template class Array; -extern template class Array; -extern template class Array; -extern template class Array; -extern template class Array; -extern template class Array; -extern template class Array; -extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; +extern template class Array; /******************************************************** - * Helper functions * + * Macros to help instantiate functions * ********************************************************/ -template -inline typename std::enable_if::value, size_t>::type getRangeSize( - const Range &range ) -{ - return ( static_cast( range.j ) - static_cast( range.i ) ) / - static_cast( range.k ); -} -template -inline typename std::enable_if::value, size_t>::type getRangeSize( - const Range &range ) -{ - double tmp = static_cast( ( range.j - range.i ) ) / static_cast( range.k ); - return static_cast( floor( tmp + 1e-12 ) + 1 ); -} -template -inline typename std::enable_if>::value || - std::is_same>::value, - size_t>::type -getRangeSize( const Range &range ) -{ - double tmp = std::real( ( range.j - range.i ) / ( range.k ) ); - return static_cast( floor( tmp + 1e-12 ) + 1 ); -} -template -inline typename std::enable_if::value, TYPE>::type getRangeValue( - const Range &range, size_t index ) -{ - return range.i + index * range.k; -} -template -inline typename std::enable_if::value, TYPE>::type getRangeValue( - const Range &range, size_t index ) -{ - return range.k * ( range.i / range.k + index ); -} -template -inline typename std::enable_if>::value || - std::is_same>::value, - TYPE>::type -getRangeValue( const Range &range, size_t index ) -{ - return range.k * ( range.i / range.k + static_cast( index ) ); -} +// clang-format off +#define instantiateArrayConstructors( TYPE ) \ + template Array::Array(); \ + template Array::~Array(); \ + template Array::Array( const ArraySize & ); \ + template Array::Array( size_t ); \ + template Array::Array( size_t, size_t ); \ + template Array::Array( size_t, size_t, size_t ); \ + template Array::Array( size_t, size_t, size_t, size_t ); \ + template Array::Array( size_t, size_t, size_t, size_t, size_t ); \ + template Array::Array( const std::vector &, const TYPE * ); \ + template Array::Array( std::initializer_list ); \ + template Array::Array( std::initializer_list> ); \ + template Array::Array( const Array & ); \ + template Array::Array( Array && ); \ + template void Array::reshape( ArraySize const& ); \ + template void Array::squeeze(); \ + template std::unique_ptr> \ + Array::constView(ArraySize const&, std::shared_ptr const&); \ + template void Array::viewRaw( ArraySize const&, TYPE*, bool, bool ); \ + template void Array::view2(ArraySize const&, std::shared_ptr ); \ + template Array &Array::operator=( const Array & ); \ + template Array &Array::operator=( Array && ); +// clang-format on /******************************************************** @@ -126,19 +107,8 @@ Array::Array( const std::vector &N, const TYPE *da : d_isCopyable( true ), d_isFixedSize( false ) { allocate( N ); - if ( data ) { - for ( size_t i = 0; i < d_size.length(); i++ ) - d_data[i] = data[i]; - } -} -template -Array::Array( const Range &range ) - : d_isCopyable( true ), d_isFixedSize( false ) -{ - size_t N = getRangeSize( range ); - allocate( { N } ); - for ( size_t i = 0; i < N; i++ ) - d_data[i] = getRangeValue( range, i ); + if ( data ) + copy( data ); } template Array::Array( std::string str ) : d_isCopyable( true ), d_isFixedSize( false ) @@ -216,8 +186,12 @@ Array::Array( std::string str ) : d_isCopyable( true ), d_ i2 = str.length(); } allocate( data.size() ); - for ( size_t i = 0; i < data.size(); i++ ) - d_data[i] = data[i]; + if ( std::is_same::value ) { + for ( size_t i = 0; i < data.size(); i++ ) + d_data[i] = data[i]; + } else { + copy( data.data() ); + } } template Array::Array( std::initializer_list x ) @@ -229,19 +203,38 @@ Array::Array( std::initializer_list x ) d_data[i] = *it; } template +Array::Array( std::initializer_list> x ) + : d_isCopyable( true ), d_isFixedSize( false ) +{ + size_t Nx = x.size(); + size_t Ny = 0; + for ( const auto y : x ) + Ny = std::max( Ny, y.size() ); + allocate( { Nx, Ny } ); + auto itx = x.begin(); + for ( size_t i = 0; i < x.size(); ++i, ++itx ) { + auto ity = itx->begin(); + for ( size_t j = 0; j < itx->size(); ++j, ++ity ) { + d_data[i + j * Nx] = *ity; + } + } +} +template void Array::allocate( const ArraySize &N ) { if ( d_isFixedSize ) throw std::logic_error( "Array cannot be resized" ); d_size = N; auto length = d_size.length(); - if ( length == 0 ) - d_ptr.reset(); - else - d_ptr.reset( new ( std::nothrow ) TYPE[length], []( TYPE *p ) { delete[] p; } ); - d_data = d_ptr.get(); - if ( length > 0 && d_data == nullptr ) - throw std::logic_error( "Failed to allocate array" ); + d_data = nullptr; + if ( length > 0 ) { + try { + d_data = new TYPE[length]; + } catch ( ... ) { + throw std::logic_error( "Failed to allocate array" ); + } + } + d_ptr.reset( d_data, []( TYPE *p ) { delete[] p; } ); } template Array::Array( const Array &rhs ) @@ -250,18 +243,19 @@ Array::Array( const Array &rhs ) if ( !rhs.d_isCopyable ) throw std::logic_error( "Array cannot be copied" ); allocate( rhs.size() ); - for ( size_t i = 0; i < d_size.length(); i++ ) - d_data[i] = rhs.d_data[i]; + copy( rhs.d_data ); } template Array::Array( Array &&rhs ) : d_isCopyable( rhs.d_isCopyable ), d_isFixedSize( rhs.d_isFixedSize ), d_size( rhs.d_size ), - d_data( rhs.d_data ) + d_data( rhs.d_data ), + d_ptr( std::move( rhs.d_ptr ) ) { + rhs.d_size = ArraySize(); rhs.d_data = nullptr; - d_ptr = std::move( rhs.d_ptr ); + rhs.d_ptr = nullptr; } template Array &Array::operator=( const Array &rhs ) @@ -270,9 +264,8 @@ Array &Array::operator=( const Array return *this; if ( !rhs.d_isCopyable ) throw std::logic_error( "Array cannot be copied" ); - this->allocate( rhs.size() ); - for ( size_t i = 0; i < d_size.length(); i++ ) - this->d_data[i] = rhs.d_data[i]; + allocate( rhs.size() ); + copy( rhs.d_data ); return *this; } template @@ -285,15 +278,17 @@ Array &Array::operator=( Array &&rhs d_size = rhs.d_size; d_data = rhs.d_data; d_ptr = std::move( rhs.d_ptr ); + rhs.d_size = ArraySize(); rhs.d_data = nullptr; + rhs.d_ptr = nullptr; return *this; } template Array &Array::operator=( const std::vector &rhs ) { - this->allocate( ArraySize( rhs.size() ) ); + allocate( ArraySize( rhs.size() ) ); for ( size_t i = 0; i < rhs.size(); i++ ) - this->d_data[i] = rhs[i]; + d_data[i] = rhs[i]; return *this; } template @@ -331,9 +326,9 @@ static inline void moveValues( const ArraySize &N1, const ArraySize &N2, TYPE *d } } } -template -static inline typename std::enable_if::type copyValues( - const ArraySize &N1, const ArraySize &N2, const TYPE *data1, TYPE *data2 ) +template +static inline void +copyValues( const ArraySize &N1, const ArraySize &N2, const TYPE *data1, TYPE *data2 ) { for ( size_t i5 = 0; i5 < std::min( N1[4], N2[4] ); i5++ ) { for ( size_t i4 = 0; i4 < std::min( N1[3], N2[3] ); i4++ ) { @@ -349,12 +344,6 @@ static inline typename std::enable_if::type copyValues( } } } -template -static inline typename std::enable_if::type copyValues( - const ArraySize &, const ArraySize &, const TYPE *, TYPE * ) -{ - throw std::logic_error( "No copy constructor" ); -} /******************************************************** @@ -381,9 +370,11 @@ void Array::resize( const ArraySize &N ) if ( data0.use_count() <= 1 ) { // We own the data, use std:move moveValues( N0, N, data0.get(), d_data ); - } else { + } else if ( std::is_copy_constructible::value ) { // We do not own the data, copy - copyValues::value, TYPE>( N0, N, data0.get(), d_data ); + copyValues( N0, N, data0.get(), d_data ); + } else { + throw std::logic_error( "No copy constructor" ); } } } @@ -412,7 +403,7 @@ void Array::resizeDim( int dim, size_t N, const TYPE &valu /******************************************************** - * Reshape the array * + * Reshape/squeeze the array * ********************************************************/ template void Array::reshape( const ArraySize &N ) @@ -421,6 +412,85 @@ void Array::reshape( const ArraySize &N ) throw std::logic_error( "reshape is not allowed to change the array size" ); d_size = N; } +template +void Array::squeeze() +{ + d_size.squeeze(); +} + + +/******************************************************** + * Shift/permute the array * + ********************************************************/ +template +Array Array::shiftDim( int N ) const +{ + if ( N > 0 ) + N = N % d_size.ndim(); + if ( N == 0 ) { + // No shift required + return *this; + } else if ( N > 0 ) { + // Shift to the left and wrap + std::vector index( d_size.ndim() ); + size_t i = 0; + for ( size_t j=N; j +Array Array::permute( const std::vector &index ) const +{ + // Check the permutation + ASSERT( (int) index.size() == ndim() ); + for ( int i=0; i < ndim(); i++) { + ASSERT( index[i] < ndim() ); + for ( int j=0; j < i; j++) + ASSERT( index[i] != index[j] ); + } + // Create the new Array + size_t dims[5] = { 1u, 1u, 1u, 1u, 1u }; + for ( size_t i=0; ilength() ); + // Fill the data + size_t N[5] = { 1u, 1u, 1u, 1u, 1u }; + for ( int i=0; i < ndim(); i++) { + std::array ijk = { 0, 0, 0, 0, 0 }; + ijk[index[i]] = 1; + N[i] = d_size.index( ijk ); + } + size_t tmp = ( dims[0] - 1 ) * N[0] + ( dims[1] - 1 ) * N[1] + ( dims[2] - 1 ) * N[2] + ( dims[3] - 1 ) * N[3] + ( dims[4] - 1 ) * N[4] + 1; + ASSERT( tmp == length() ); + for ( size_t i4 = 0; i4 < dims[4]; i4++ ) { + for ( size_t i3 = 0; i3 < dims[3]; i3++ ) { + for ( size_t i2 = 0; i2 < dims[2]; i2++ ) { + for ( size_t i1 = 0; i1 < dims[1]; i1++ ) { + for ( size_t i0 = 0; i0 < dims[0]; i0++ ) { + size_t index2 = i0 * N[0] + i1 * N[1] + i2 * N[2] + i3 * N[3] + i4 * N[4]; + y( i0, i1, i2, i3, i4 ) = d_data[index2]; + } + } + } + } + } + return y; +} /******************************************************** @@ -428,8 +498,8 @@ void Array::reshape( const ArraySize &N ) ********************************************************/ // Helper function to check subset indices template -inline void Array::checkSubsetIndex( - const std::vector> &range ) const +inline void +Array::checkSubsetIndex( const std::vector> &range ) const { bool test = (int) range.size() == d_size.ndim(); for ( size_t d = 0; d < range.size(); d++ ) @@ -438,8 +508,8 @@ inline void Array::checkSubsetIndex( throw std::logic_error( "indices for subset are invalid" ); } template -std::vector> Array::convert( - const std::vector &index ) const +std::vector> +Array::convert( const std::vector &index ) const { std::vector> range( d_size.ndim() ); if ( index.size() % 2 != 0 || static_cast( index.size() / 2 ) < d_size.ndim() ) @@ -451,8 +521,10 @@ std::vector> Array::convert( // Helper function to return dimensions for the subset array template void Array::getSubsetArrays( const std::vector> &index, - std::array &first, std::array &last, std::array &inc, - std::array &N ) + std::array &first, + std::array &last, + std::array &inc, + std::array &N ) { first.fill( 0 ); last.fill( 0 ); @@ -467,8 +539,8 @@ void Array::getSubsetArrays( const std::vector -Array Array::subset( - const std::vector> &index ) const +Array +Array::subset( const std::vector> &index ) const { // Get the subset indicies checkSubsetIndex( index ); @@ -476,9 +548,8 @@ Array Array::subset( getSubsetArrays( index, first, last, inc, N1 ); ArraySize S1( d_size.ndim(), N1.data() ); // Create the new array - Array subset_array( S1 ); + Array subset_array( S1 ); // Fill the new array - static_assert( ArraySize::maxDim() == 5, "Not programmed for more than 5 dimensions" ); TYPE *subset_data = subset_array.data(); for ( size_t i4 = first[4], k1 = 0; i4 <= last[4]; i4 += inc[4] ) { for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) { @@ -495,22 +566,21 @@ Array Array::subset( return subset_array; } template -Array Array::subset( - const std::vector &index ) const +Array +Array::subset( const std::vector &index ) const { auto range = convert( index ); return subset( range ); } template -void Array::copySubset( - const std::vector> &index, const Array &subset ) +void Array::copySubset( const std::vector> &index, + const Array &subset ) { // Get the subset indices checkSubsetIndex( index ); std::array first, last, inc, N1; getSubsetArrays( index, first, last, inc, N1 ); // Copy the sub-array - static_assert( ArraySize::maxDim() == 5, "Not programmed for more than 5 dimensions" ); const TYPE *src_data = subset.data(); for ( size_t i4 = first[4], k1 = 0; i4 <= last[4]; i4 += inc[4] ) { for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) { @@ -526,15 +596,14 @@ void Array::copySubset( } } template -void Array::addSubset( - const std::vector> &index, const Array &subset ) +void Array::addSubset( const std::vector> &index, + const Array &subset ) { // Get the subset indices checkSubsetIndex( index ); std::array first, last, inc, N1; getSubsetArrays( index, first, last, inc, N1 ); // add the sub-array - static_assert( ArraySize::maxDim() == 5, "Not programmed for more than 5 dimensions" ); for ( size_t i4 = first[4], k1 = 0; i4 <= last[4]; i4 += inc[4] ) { for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) { for ( size_t i2 = first[2]; i2 <= last[2]; i2 += inc[2] ) { @@ -549,16 +618,16 @@ void Array::addSubset( } } template -void Array::copySubset( - const std::vector &index, const Array &subset ) +void Array::copySubset( const std::vector &index, + const Array &subset ) { auto range = convert( index ); copySubset( range, subset ); } template -void Array::addSubset( - const std::vector &index, const Array &subset ) +void Array::addSubset( const std::vector &index, + const Array &subset ) { auto range = convert( index ); addSubset( range, subset ); @@ -586,8 +655,8 @@ bool Array::operator==( const Array &rhs ) const * Get a view of an C array * ********************************************************/ template -std::unique_ptr> Array::view( - const ArraySize &N, std::shared_ptr &data ) +std::unique_ptr> +Array::view( const ArraySize &N, std::shared_ptr data ) { auto array = std::make_unique>(); array->d_size = N; @@ -596,8 +665,9 @@ std::unique_ptr> Array::view( return array; } template -std::unique_ptr> Array::constView( - const ArraySize &N, std::shared_ptr const &data ) +std::unique_ptr> +Array::constView( const ArraySize &N, + std::shared_ptr const &data ) { auto array = std::make_unique>(); array->d_size = N; @@ -612,15 +682,17 @@ void Array::view2( Array &src ) d_data = src.d_data; } template -void Array::view2( const ArraySize &N, std::shared_ptr const &data ) +void Array::view2( const ArraySize &N, std::shared_ptr data ) { d_size = N; d_ptr = data; d_data = d_ptr.get(); } template -void Array::viewRaw( - const ArraySize &N, TYPE *data, bool isCopyable, bool isFixedSize ) +void Array::viewRaw( const ArraySize &N, + TYPE *data, + bool isCopyable, + bool isFixedSize ) { d_isCopyable = isCopyable; d_isFixedSize = isFixedSize; @@ -644,20 +716,8 @@ void Array::swap( Array &other ) std::swap( d_ptr, other.d_ptr ); } template -void Array::fill( const TYPE &value ) -{ - for ( size_t i = 0; i < d_size.length(); i++ ) - d_data[i] = value; -} -template -void Array::scale( const TYPE &value ) -{ - for ( size_t i = 0; i < d_size.length(); i++ ) - d_data[i] *= value; -} -template -void Array::pow( - const Array &baseArray, const TYPE &exp ) +void Array::pow( const Array &baseArray, + const TYPE &exp ) { // not insisting on the shapes being the same // but insisting on the total size being the same @@ -674,8 +734,8 @@ void Array::pow( * Replicate the array * ********************************************************/ template -Array Array::repmat( - const std::vector &N_rep ) const +Array +Array::repmat( const std::vector &N_rep ) const { std::vector N2( d_size.begin(), d_size.end() ); if ( N2.size() < N_rep.size() ) @@ -689,7 +749,6 @@ Array Array::repmat( N2[d] *= N_rep[d]; } Array y( N2 ); - static_assert( ArraySize::maxDim() <= 5, "Not programmed for dimensions > 5" ); TYPE *y2 = y.data(); for ( size_t i4 = 0, index = 0; i4 < N1[4]; i4++ ) { for ( size_t j4 = 0; j4 < Nr[4]; j4++ ) { @@ -731,7 +790,7 @@ bool Array::NaNs() const template TYPE Array::mean( void ) const { - TYPE x = this->sum() / d_size.length(); + TYPE x = sum() / d_size.length(); return x; } template @@ -813,7 +872,6 @@ TYPE Array::min( const std::vector> &range ) checkSubsetIndex( range ); std::array first, last, inc, N1; getSubsetArrays( range, first, last, inc, N1 ); - static_assert( ArraySize::maxDim() <= 5, "Function programmed for more than 5 dimensions" ); TYPE x = std::numeric_limits::max(); for ( size_t i4 = first[4]; i4 <= last[4]; i4 += inc[4] ) { for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) { @@ -836,7 +894,6 @@ TYPE Array::max( const std::vector> &range ) checkSubsetIndex( range ); std::array first, last, inc, N1; getSubsetArrays( range, first, last, inc, N1 ); - static_assert( ArraySize::maxDim() <= 5, "Function programmed for more than 5 dimensions" ); TYPE x = std::numeric_limits::min(); for ( size_t i4 = first[4]; i4 <= last[4]; i4 += inc[4] ) { for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) { @@ -859,7 +916,6 @@ TYPE Array::sum( const std::vector> &range ) checkSubsetIndex( range ); std::array first, last, inc, N1; getSubsetArrays( range, first, last, inc, N1 ); - static_assert( ArraySize::maxDim() <= 5, "Function programmed for more than 5 dimensions" ); TYPE x = 0; for ( size_t i4 = first[4]; i4 <= last[4]; i4 += inc[4] ) { for ( size_t i3 = first[3]; i3 <= last[3]; i3 += inc[3] ) { @@ -882,7 +938,6 @@ TYPE Array::mean( const std::vector> &range checkSubsetIndex( range ); std::array first, last, inc, N1; getSubsetArrays( range, first, last, inc, N1 ); - static_assert( ArraySize::maxDim() <= 5, "Function programmed for more than 5 dimensions" ); size_t n = 1; for ( auto &d : N1 ) n *= d; @@ -919,8 +974,9 @@ TYPE Array::mean( const std::vector &index ) const * Find all elements that match the given operation * ********************************************************/ template -std::vector Array::find( - const TYPE &value, std::function compare ) const +std::vector +Array::find( const TYPE &value, + std::function compare ) const { std::vector result; result.reserve( d_size.length() ); @@ -936,8 +992,9 @@ std::vector Array::find( * Print an array to an output stream * ********************************************************/ template -void Array::print( - std::ostream &os, const std::string &name, const std::string &prefix ) const +void Array::print( std::ostream &os, + const std::string &name, + const std::string &prefix ) const { if ( d_size.ndim() == 1 ) { for ( size_t i = 0; i < d_size[0]; i++ ) @@ -961,12 +1018,11 @@ void Array::print( template Array Array::reverseDim() const { - size_t N2[ArraySize::maxDim()]; + size_t N2[5]; for ( int d = 0; d < ArraySize::maxDim(); d++ ) N2[d] = d_size[ArraySize::maxDim() - d - 1]; ArraySize S2( ArraySize::maxDim(), N2 ); Array y( S2 ); - static_assert( ArraySize::maxDim() == 5, "Not programmed for dimensions other than 5" ); TYPE *y2 = y.data(); for ( size_t i0 = 0; i0 < d_size[0]; i0++ ) { for ( size_t i1 = 0; i1 < d_size[1]; i1++ ) { @@ -991,8 +1047,8 @@ Array Array::reverseDim() const * Coarsen the array * ********************************************************/ template -Array Array::coarsen( - const Array &filter ) const +Array +Array::coarsen( const Array &filter ) const { auto S2 = size(); for ( size_t i = 0; i < S2.size(); i++ ) { @@ -1012,8 +1068,9 @@ Array Array::coarsen( for ( size_t k2 = 0; k2 < Nh[2]; k2++ ) { for ( size_t j2 = 0; j2 < Nh[1]; j2++ ) { for ( size_t i2 = 0; i2 < Nh[0]; i2++ ) { - tmp += filter( i2, j2, k2 ) * this->operator()( i1 *Nh[0] + i2, - j1 * Nh[1] + j2, k1 * Nh[2] + k2 ); + tmp += filter( i2, j2, k2 ) * operator()( i1 *Nh[0] + i2, + j1 * Nh[1] + j2, + k1 * Nh[2] + k2 ); } } } @@ -1024,7 +1081,8 @@ Array Array::coarsen( return y; } template -Array Array::coarsen( const std::vector &ratio, +Array Array::coarsen( + const std::vector &ratio, std::function & )> filter ) const { if ( ratio.size() != d_size.ndim() ) @@ -1045,7 +1103,7 @@ Array Array::coarsen( const std::vec for ( size_t k2 = 0; k2 < ratio[2]; k2++ ) { for ( size_t j2 = 0; j2 < ratio[1]; j2++ ) { for ( size_t i2 = 0; i2 < ratio[0]; i2++ ) { - tmp( i2, j2, k2 ) = this->operator()( + tmp( i2, j2, k2 ) = operator()( i1 *ratio[0] + i2, j1 * ratio[1] + j2, k1 * ratio[2] + k2 ); } } @@ -1070,13 +1128,25 @@ void Array::cat( const Array &x, int *this = cat( tmp, dim ); } template +Array Array::cat( const std::initializer_list &x, + int dim ) +{ + return cat( x.size(), x.begin(), dim ); +} +template Array Array::cat( const std::vector &x, int dim ) { - if ( x.empty() ) + return cat( x.size(), x.data(), dim ); +} +template +Array +Array::cat( size_t N_array, const Array *x, int dim ) +{ + if ( N_array == 0 ) return Array(); // Check that the dimensions match bool check = true; - for ( size_t i = 1; i < x.size(); i++ ) { + for ( size_t i = 1; i < N_array; i++ ) { check = check && x[i].ndim() == x[0].ndim(); for ( int d = 0; d < x[0].ndim(); d++ ) if ( d != dim ) @@ -1086,7 +1156,7 @@ Array Array::cat( const std::vector< throw std::logic_error( "Array dimensions do not match for concatenation" ); // Create the output array auto size = x[0].d_size; - for ( size_t i = 1; i < x.size(); i++ ) + for ( size_t i = 1; i < N_array; i++ ) size.resize( dim, size[dim] + x[i].size( dim ) ); Array out( size ); size_t N1 = 1; @@ -1097,7 +1167,7 @@ Array Array::cat( const std::vector< for ( size_t d = dim + 1; d < size.ndim(); d++ ) N3 *= size[d]; TYPE *data = out.data(); - for ( size_t i = 0, i0 = 0; i < x.size(); i++ ) { + for ( size_t i = 0, i0 = 0; i < N_array; i++ ) { const TYPE *src = x[i].data(); size_t N22 = x[i].size( dim ); for ( size_t j2 = 0; j2 < N3; j2++ ) { @@ -1117,87 +1187,82 @@ Array Array::cat( const std::vector< * Interpolate * ********************************************************/ template -struct is_compatible_double - : std::integral_constant::value || std::is_integral::value> { -}; -template -inline typename std::enable_if::value, TYPE>::type Array_interp_1D( - double x, int N, const TYPE *data ) +constexpr bool is_compatible_double() { - int i = floor( x ); - i = std::max( i, 0 ); - i = std::min( i, N - 2 ); - return ( i + 1 - x ) * data[i] + ( x - i ) * data[i + 1]; + return std::is_floating_point::value || std::is_integral::value; } template -inline typename std::enable_if::value, TYPE>::type Array_interp_2D( - double x, double y, int Nx, int Ny, const TYPE *data ) +inline TYPE Array_interp_1D( double x, int N, const TYPE *data ) { - int i = floor( x ); - i = std::max( i, 0 ); - i = std::min( i, Nx - 2 ); - double dx = x - i; - double dx2 = 1.0 - dx; - int j = floor( y ); - j = std::max( j, 0 ); - j = std::min( j, Ny - 2 ); - double dy = y - j; - double dy2 = 1.0 - dy; - double f[4] = { (double) data[i + j * Nx], (double) data[i + 1 + j * Nx], - (double) data[i + ( j + 1 ) * Nx], (double) data[i + 1 + ( j + 1 ) * Nx] }; - return ( dx * f[1] + dx2 * f[0] ) * dy2 + ( dx * f[3] + dx2 * f[2] ) * dy; + if ( is_compatible_double() ) { + int i = floor( x ); + i = std::max( i, 0 ); + i = std::min( i, N - 2 ); + return ( i + 1 - x ) * data[i] + ( x - i ) * data[i + 1]; + } else { + throw std::logic_error( "Invalid conversion" ); + } } template -inline typename std::enable_if::value, TYPE>::type Array_interp_3D( - double x, double y, double z, int Nx, int Ny, int Nz, const TYPE *data ) +inline TYPE Array_interp_2D( double x, double y, int Nx, int Ny, const TYPE *data ) { - int i = floor( x ); - i = std::max( i, 0 ); - i = std::min( i, Nx - 2 ); - double dx = x - i; - double dx2 = 1.0 - dx; - int j = floor( y ); - j = std::max( j, 0 ); - j = std::min( j, Ny - 2 ); - double dy = y - j; - double dy2 = 1.0 - dy; - int k = floor( z ); - k = std::max( k, 0 ); - k = std::min( k, Nz - 2 ); - double dz = z - k; - double dz2 = 1.0 - dz; - double f[8] = { (double) data[i + j * Nx + k * Nx * Ny], - (double) data[i + 1 + j * Nx + k * Nx * Ny], - (double) data[i + ( j + 1 ) * Nx + k * Nx * Ny], - (double) data[i + 1 + ( j + 1 ) * Nx + k * Nx * Ny], - (double) data[i + j * Nx + ( k + 1 ) * Nx * Ny], - (double) data[i + 1 + j * Nx + ( k + 1 ) * Nx * Ny], - (double) data[i + ( j + 1 ) * Nx + ( k + 1 ) * Nx * Ny], - (double) data[i + 1 + ( j + 1 ) * Nx + ( k + 1 ) * Nx * Ny] }; - double h0 = ( dx * f[1] + dx2 * f[0] ) * dy2 + ( dx * f[3] + dx2 * f[2] ) * dy; - double h1 = ( dx * f[5] + dx2 * f[4] ) * dy2 + ( dx * f[7] + dx2 * f[6] ) * dy; - return h0 * dz2 + h1 * dz; + if ( is_compatible_double() ) { + int i = floor( x ); + i = std::max( i, 0 ); + i = std::min( i, Nx - 2 ); + double dx = x - i; + double dx2 = 1.0 - dx; + int j = floor( y ); + j = std::max( j, 0 ); + j = std::min( j, Ny - 2 ); + double dy = y - j; + double dy2 = 1.0 - dy; + double f[4] = { (double) data[i + j * Nx], + (double) data[i + 1 + j * Nx], + (double) data[i + ( j + 1 ) * Nx], + (double) data[i + 1 + ( j + 1 ) * Nx] }; + return ( dx * f[1] + dx2 * f[0] ) * dy2 + ( dx * f[3] + dx2 * f[2] ) * dy; + } else { + throw std::logic_error( "Invalid conversion" ); + } } template -inline typename std::enable_if::value, TYPE>::type Array_interp_1D( - double, int, const TYPE * ) +inline TYPE +Array_interp_3D( double x, double y, double z, int Nx, int Ny, int Nz, const TYPE *data ) { - throw std::logic_error( "Invalid conversion" ); -} -template -inline typename std::enable_if::value, TYPE>::type Array_interp_2D( - double, double, int, int, const TYPE * ) -{ - throw std::logic_error( "Invalid conversion" ); -} -template -inline typename std::enable_if::value, TYPE>::type Array_interp_3D( - double, double, double, int, int, int, const TYPE * ) -{ - throw std::logic_error( "Invalid conversion" ); + if ( is_compatible_double() ) { + int i = floor( x ); + i = std::max( i, 0 ); + i = std::min( i, Nx - 2 ); + double dx = x - i; + double dx2 = 1.0 - dx; + int j = floor( y ); + j = std::max( j, 0 ); + j = std::min( j, Ny - 2 ); + double dy = y - j; + double dy2 = 1.0 - dy; + int k = floor( z ); + k = std::max( k, 0 ); + k = std::min( k, Nz - 2 ); + double dz = z - k; + double dz2 = 1.0 - dz; + double f[8] = { (double) data[i + j * Nx + k * Nx * Ny], + (double) data[i + 1 + j * Nx + k * Nx * Ny], + (double) data[i + ( j + 1 ) * Nx + k * Nx * Ny], + (double) data[i + 1 + ( j + 1 ) * Nx + k * Nx * Ny], + (double) data[i + j * Nx + ( k + 1 ) * Nx * Ny], + (double) data[i + 1 + j * Nx + ( k + 1 ) * Nx * Ny], + (double) data[i + ( j + 1 ) * Nx + ( k + 1 ) * Nx * Ny], + (double) data[i + 1 + ( j + 1 ) * Nx + ( k + 1 ) * Nx * Ny] }; + double h0 = ( dx * f[1] + dx2 * f[0] ) * dy2 + ( dx * f[3] + dx2 * f[2] ) * dy; + double h1 = ( dx * f[5] + dx2 * f[4] ) * dy2 + ( dx * f[7] + dx2 * f[6] ) * dy; + return h0 * dz2 + h1 * dz; + } else { + throw std::logic_error( "Invalid conversion" ); + } } template -TYPE Array::interp( const std::vector &x ) const +TYPE Array::interp( const double *x ) const { int ndim = 0, dim[5]; double x2[5]; @@ -1233,81 +1298,75 @@ void Array::rand() FUN::rand( *this ); } template -Array &Array::operator+=( - const Array &rhs ) +Array & +Array::operator+=( const Array &rhs ) { - const auto &fun = []( const TYPE &a, const TYPE &b ) { return a + b; }; - FUN::transform( fun, *this, rhs, *this ); + auto op = []( const TYPE &a, const TYPE &b ) { return a + b; }; + FUN::transform( op, *this, rhs, *this ); return *this; } template -Array &Array::operator-=( - const Array &rhs ) +Array & +Array::operator-=( const Array &rhs ) { - const auto &fun = []( const TYPE &a, const TYPE &b ) { return a - b; }; - FUN::transform( fun, *this, rhs, *this ); + auto op = []( const TYPE &a, const TYPE &b ) { return a - b; }; + FUN::transform( op, *this, rhs, *this ); return *this; } template Array &Array::operator+=( const TYPE &rhs ) { - const auto &fun = [rhs]( const TYPE &x ) { return x + rhs; }; - FUN::transform( fun, *this, *this ); + auto op = [rhs]( const TYPE &x ) { return x + rhs; }; + FUN::transform( op, *this, *this ); return *this; } template Array &Array::operator-=( const TYPE &rhs ) { - const auto &fun = [rhs]( const TYPE &x ) { return x - rhs; }; - FUN::transform( fun, *this, *this ); + auto op = [rhs]( const TYPE &x ) { return x - rhs; }; + FUN::transform( op, *this, *this ); return *this; } template TYPE Array::min() const { - const auto &fun = []( const TYPE &a, const TYPE &b ) { return a < b ? a : b; }; - return FUN::reduce( fun, *this ); + const auto &op = []( const TYPE &a, const TYPE &b ) { return a < b ? a : b; }; + return FUN::reduce( op, *this, d_data[0] ); } template TYPE Array::max() const { - const auto &fun = []( const TYPE &a, const TYPE &b ) { return a > b ? a : b; }; - return FUN::reduce( fun, *this ); + const auto &op = []( const TYPE &a, const TYPE &b ) { return a > b ? a : b; }; + return FUN::reduce( op, *this, d_data[0] ); } template TYPE Array::sum() const { - const auto &fun = []( const TYPE &a, const TYPE &b ) { return a + b; }; - return FUN::reduce( fun, *this ); + const auto &op = []( const TYPE &a, const TYPE &b ) { return a + b; }; + return FUN::reduce( op, *this, static_cast( 0 ) ); } template -Array Array::multiply( - const Array &a, const Array &b ) +void Array::axpby( const TYPE &alpha, + const Array &x, + const TYPE &beta ) { - Array c; - FUN::multiply( a, b, c ); - return c; + const auto &op = [alpha, beta]( const TYPE &x, const TYPE &y ) { return alpha * x + beta * y; }; + return FUN::transform( op, x, *this, *this ); } template -void Array::axpby( - const TYPE &alpha, const Array &x, const TYPE &beta ) -{ - const auto &fun = [alpha, beta]( - const TYPE &x, const TYPE &y ) { return alpha * x + beta * y; }; - return FUN::transform( fun, x, *this, *this ); -} -template -Array Array::transform( - std::function fun, const Array &x ) +Array +Array::transform( std::function fun, + const Array &x ) { Array y; FUN::transform( fun, x, y ); return y; } template -Array Array::transform( - std::function fun, const Array &x, - const Array &y ) +Array +Array::transform( std::function fun, + const Array &x, + const Array &y ) { Array z; FUN::transform( fun, x, y, z ); @@ -1319,4 +1378,5 @@ bool Array::equals( const Array &rhs, TYPE tol ) const return FUN::equals( *this, rhs, tol ); } + #endif diff --git a/common/ArraySize.h b/common/ArraySize.h index a3cb78c3..808b7863 100644 --- a/common/ArraySize.h +++ b/common/ArraySize.h @@ -1,8 +1,12 @@ #ifndef included_ArraySizeClass #define included_ArraySizeClass +#include "common/Utilities.h" #include +#include +#include +#include #include #include #include @@ -22,21 +26,22 @@ #if ( defined( DEBUG ) || defined( _DEBUG ) ) && !defined( NDEBUG ) -#define CHECK_ARRAY_LENGTH( i ) \ +#define CHECK_ARRAY_LENGTH( i, length ) \ do { \ - if ( i >= d_length ) \ + if ( i >= length ) \ throw std::out_of_range( "Index exceeds array bounds" ); \ } while ( 0 ) #else -#define CHECK_ARRAY_LENGTH( i ) \ - do { \ +#define CHECK_ARRAY_LENGTH( i, length ) \ + do { \ } while ( 0 ) #endif + // Forward declerations class FunctionTable; -template +template> class Array; @@ -46,7 +51,7 @@ class Range final { public: //! Empty constructor - constexpr Range() : i( 0 ), j( -1 ), k( 1 ) {} + Range() : i( 0 ), j( -1 ), k( 1 ) {} /*! * Create a range i:k:j (or i:j) @@ -54,8 +59,30 @@ public: * @param j_ Ending value * @param k_ Increment value */ - constexpr Range( TYPE i_, TYPE j_, TYPE k_ = 1 ) : i( i_ ), j( j_ ), k( k_ ) {} + Range( const TYPE &i_, const TYPE &j_, const TYPE &k_ = 1 ) + : i( i_ ), j( j_ ), k( k_ ) + { + } + //! Get the number of values in the range + size_t size() const + { + if ( std::is_integral::value ) { + return ( static_cast( j ) - static_cast( i ) ) / + static_cast( k ); + } else if ( std::is_floating_point::value ) { + double tmp = static_cast( ( j - i ) ) / static_cast( k ); + return static_cast( floor( tmp + 1e-12 ) + 1 ); + } else if ( std::is_same>::value || + std::is_same>::value ) { + double tmp = std::real( ( j - i ) / ( k ) ); + return static_cast( floor( tmp + 1e-12 ) + 1 ); + } else { + ERROR( "Unsupported type for range" ); + } + } + +public: TYPE i, j, k; }; @@ -65,20 +92,20 @@ class ArraySize final { public: //! Empty constructor - constexpr ArraySize() : d_ndim( 1 ), d_length( 0 ), d_N{ 0, 1, 1, 1, 1 } {} + ArraySize() : d_ndim( 1 ), d_length( 0 ), d_N{ 0, 1, 1, 1, 1 } {} /*! * Create the vector size * @param N1 Number of elements in the first dimension */ - constexpr ArraySize( size_t N1 ) : d_ndim( 1 ), d_length( N1 ), d_N{ N1, 1, 1, 1, 1 } {} + ArraySize( size_t N1 ) : d_ndim( 1 ), d_length( N1 ), d_N{ N1, 1, 1, 1, 1 } {} /*! * Create the vector size * @param N1 Number of elements in the first dimension * @param N2 Number of elements in the second dimension */ - constexpr ArraySize( size_t N1, size_t N2 ) + ArraySize( size_t N1, size_t N2 ) : d_ndim( 2 ), d_length( N1 * N2 ), d_N{ N1, N2, 1, 1, 1 } { } @@ -89,7 +116,7 @@ public: * @param N2 Number of elements in the second dimension * @param N3 Number of elements in the third dimension */ - constexpr ArraySize( size_t N1, size_t N2, size_t N3 ) + ArraySize( size_t N1, size_t N2, size_t N3 ) : d_ndim( 3 ), d_length( N1 * N2 * N3 ), d_N{ N1, N2, N3, 1, 1 } { } @@ -101,7 +128,7 @@ public: * @param N3 Number of elements in the third dimension * @param N4 Number of elements in the fourth dimension */ - constexpr ArraySize( size_t N1, size_t N2, size_t N3, size_t N4 ) + ArraySize( size_t N1, size_t N2, size_t N3, size_t N4 ) : d_ndim( 4 ), d_length( N1 * N2 * N3 * N4 ), d_N{ N1, N2, N3, N4, 1 } { } @@ -114,7 +141,7 @@ public: * @param N4 Number of elements in the fourth dimension * @param N5 Number of elements in the fifth dimension */ - constexpr ArraySize( size_t N1, size_t N2, size_t N3, size_t N4, size_t N5 ) + ArraySize( size_t N1, size_t N2, size_t N3, size_t N4, size_t N5 ) : d_ndim( 5 ), d_length( N1 * N2 * N3 * N4 * N5 ), d_N{ N1, N2, N3, N4, N5 } { } @@ -122,11 +149,14 @@ public: /*! * Create from initializer list * @param N Size of the array + * @param ndim Number of dimensions */ - constexpr ArraySize( std::initializer_list N ) + ArraySize( std::initializer_list N, int ndim = -1 ) : d_ndim( N.size() ), d_length( 0 ), d_N{ 0, 1, 1, 1, 1 } { - if ( d_ndim > maxDim() ) + if ( ndim >= 0 ) + d_ndim = ndim; + if ( d_ndim > 5 ) throw std::out_of_range( "Maximum number of dimensions exceeded" ); auto it = N.begin(); for ( size_t i = 0; i < d_ndim; i++, ++it ) @@ -144,10 +174,10 @@ public: * @param ndim Number of dimensions * @param dims Dimensions */ - constexpr ArraySize( size_t ndim, const size_t *dims ) + ArraySize( size_t ndim, const size_t *dims ) : d_ndim( ndim ), d_length( 0 ), d_N{ 0, 1, 1, 1, 1 } { - if ( d_ndim > maxDim() ) + if ( d_ndim > 5 ) throw std::out_of_range( "Maximum number of dimensions exceeded" ); for ( size_t i = 0; i < ndim; i++ ) d_N[i] = dims[i]; @@ -158,35 +188,44 @@ public: d_length = 0; } + /*! + * Create from std::array + * @param N Size of the array + */ + template + ArraySize( const std::array &N ) : ArraySize( NDIM, N.data() ) + { + } + /*! * Create from std::vector * @param N Size of the array */ - ArraySize( const std::vector &N ); + inline ArraySize( const std::vector &N ) : ArraySize( N.size(), N.data() ) {} // Copy/assignment constructors - constexpr ArraySize( ArraySize &&rhs ) = default; - constexpr ArraySize( const ArraySize &rhs ) = default; - constexpr ArraySize &operator=( ArraySize &&rhs ) = default; - constexpr ArraySize &operator=( const ArraySize &rhs ) = default; + ArraySize( ArraySize &&rhs ) = default; + ArraySize( const ArraySize &rhs ) = default; + ArraySize &operator=( ArraySize &&rhs ) = default; + ArraySize &operator=( const ArraySize &rhs ) = default; /*! * Access the ith dimension * @param i Index to access */ - constexpr ARRAY_ATTRIBUTE size_t operator[]( size_t i ) const { return d_N[i]; } + ARRAY_ATTRIBUTE size_t operator[]( size_t i ) const { return d_N[i]; } //! Return the number of dimensions - constexpr ARRAY_ATTRIBUTE uint8_t ndim() const { return d_ndim; } + ARRAY_ATTRIBUTE uint8_t ndim() const { return d_ndim; } //! Return the number of dimensions - constexpr ARRAY_ATTRIBUTE size_t size() const { return d_ndim; } + ARRAY_ATTRIBUTE size_t size() const { return d_ndim; } //! Return the total number of elements in the array - constexpr ARRAY_ATTRIBUTE size_t length() const { return d_length; } + ARRAY_ATTRIBUTE size_t length() const { return d_length; } //! Resize the dimension - constexpr void resize( uint8_t dim, size_t N ) + void resize( uint8_t dim, size_t N ) { if ( dim >= d_ndim ) throw std::out_of_range( "Invalid dimension" ); @@ -201,75 +240,141 @@ public: * max of ndim and the largest dim>1. * @param ndim Desired number of dimensions */ - constexpr void setNdim( uint8_t ndim ) { d_ndim = std::max( ndim, d_ndim ); } + void setNdim( uint8_t ndim ) { d_ndim = std::max( ndim, d_ndim ); } + + /*! + * Remove singleton dimensions + */ + void squeeze() + { + d_ndim = 0; + for ( uint8_t i = 0; i < maxDim(); i++ ) { + if ( d_N[i] != 1 ) + d_N[d_ndim++] = d_N[i]; + } + } //! Returns an iterator to the beginning - constexpr const size_t *begin() const { return d_N; } + const size_t *begin() const { return d_N; } //! Returns an iterator to the end - constexpr const size_t *end() const { return d_N + d_ndim; } + const size_t *end() const { return d_N + d_ndim; } // Check if two array sizes are equal - constexpr ARRAY_ATTRIBUTE bool operator==( const ArraySize &rhs ) const + ARRAY_ATTRIBUTE bool operator==( const ArraySize &rhs ) const { return d_ndim == rhs.d_ndim && memcmp( d_N, rhs.d_N, sizeof( d_N ) ) == 0; } // Check if two array sizes are equal (ignoring the dimension) - constexpr ARRAY_ATTRIBUTE bool approxEqual( const ArraySize &rhs ) const + ARRAY_ATTRIBUTE bool approxEqual( const ArraySize &rhs ) const { return ( length() == 0 && rhs.length() == 0 ) || memcmp( d_N, rhs.d_N, sizeof( d_N ) ) == 0; } //! Check if two matrices are not equal - constexpr ARRAY_ATTRIBUTE bool operator!=( const ArraySize &rhs ) const + ARRAY_ATTRIBUTE bool operator!=( const ArraySize &rhs ) const { return d_ndim != rhs.d_ndim || memcmp( d_N, rhs.d_N, sizeof( d_N ) ) != 0; } //! Maximum supported dimension - constexpr ARRAY_ATTRIBUTE static uint8_t maxDim() { return 5u; } + ARRAY_ATTRIBUTE static uint8_t maxDim() { return 5; } //! Get the index - constexpr ARRAY_ATTRIBUTE size_t index( size_t i ) const + ARRAY_ATTRIBUTE size_t index( size_t i ) const { - CHECK_ARRAY_LENGTH( i ); + CHECK_ARRAY_LENGTH( i, d_length ); return i; } //! Get the index - constexpr ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2 ) const + ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2 ) const { size_t index = i1 + i2 * d_N[0]; - CHECK_ARRAY_LENGTH( index ); + CHECK_ARRAY_LENGTH( index, d_length ); return index; } //! Get the index - constexpr ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2, size_t i3 ) const + ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2, size_t i3 ) const { size_t index = i1 + d_N[0] * ( i2 + d_N[1] * i3 ); - CHECK_ARRAY_LENGTH( index ); + CHECK_ARRAY_LENGTH( index, d_length ); return index; } //! Get the index - constexpr ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2, size_t i3, size_t i4 ) const + ARRAY_ATTRIBUTE size_t index( size_t i1, size_t i2, size_t i3, size_t i4 ) const { size_t index = i1 + d_N[0] * ( i2 + d_N[1] * ( i3 + d_N[2] * i4 ) ); - CHECK_ARRAY_LENGTH( index ); + CHECK_ARRAY_LENGTH( index, d_length ); return index; } //! Get the index - constexpr ARRAY_ATTRIBUTE size_t index( - size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) const + ARRAY_ATTRIBUTE size_t + index( size_t i1, size_t i2, size_t i3, size_t i4, size_t i5 ) const { size_t index = i1 + d_N[0] * ( i2 + d_N[1] * ( i3 + d_N[2] * ( i4 + d_N[3] * i5 ) ) ); - CHECK_ARRAY_LENGTH( index ); + CHECK_ARRAY_LENGTH( index, d_length ); return index; } + //! Get the index + size_t index( const std::array &i ) const + { + size_t j = 0; + for ( size_t m = 0, N = 1; m < 5; m++ ) { + j += i[m] * N; + N *= d_N[m]; + } + return j; + } + + //! Get the index + size_t index( std::initializer_list i ) const + { + size_t N = 1; + size_t j = 0; + size_t m = 0; + for ( size_t k : i ) { + j += k * N; + N *= d_N[m++]; + } + return j; + } + + //! Convert the index to ijk values + std::array ijk( size_t index ) const + { + CHECK_ARRAY_LENGTH( index, d_length ); + size_t i0 = index % d_N[0]; + index = index / d_N[0]; + size_t i1 = index % d_N[1]; + index = index / d_N[1]; + size_t i2 = index % d_N[2]; + index = index / d_N[2]; + size_t i3 = index % d_N[3]; + index = index / d_N[3]; + return { i0, i1, i2, i3, index }; + } + + //! Convert the index to ijk values + void ijk( size_t index, size_t *x ) const + { + CHECK_ARRAY_LENGTH( index, d_length ); + x[0] = index % d_N[0]; + index = index / d_N[0]; + x[1] = index % d_N[1]; + index = index / d_N[1]; + x[2] = index % d_N[2]; + index = index / d_N[2]; + x[3] = index % d_N[3]; + index = index / d_N[3]; + x[4] = index; + } + private: uint8_t d_ndim; size_t d_length; @@ -278,11 +383,11 @@ private: // Function to concatenate dimensions of two array sizes -constexpr ArraySize cat( const ArraySize &x, const ArraySize &y ) +inline ArraySize cat( const ArraySize &x, const ArraySize &y ) { - if ( x.ndim() + y.ndim() > ArraySize::maxDim() ) + if ( x.ndim() + y.ndim() > 5 ) throw std::out_of_range( "Maximum number of dimensions exceeded" ); - size_t N[ArraySize::maxDim()] = { 0 }; + size_t N[5] = { 0 }; for ( int i = 0; i < x.ndim(); i++ ) N[i] = x[i]; for ( int i = 0; i < y.ndim(); i++ ) @@ -291,4 +396,36 @@ constexpr ArraySize cat( const ArraySize &x, const ArraySize &y ) } +// Operator overloads +inline ArraySize operator*( size_t v, const ArraySize &x ) +{ + size_t N[5] = { v * x[0], v * x[1], v * x[2], v * x[3], v * x[4] }; + return ArraySize( x.ndim(), N ); +} +inline ArraySize operator*( const ArraySize &x, size_t v ) +{ + size_t N[5] = { v * x[0], v * x[1], v * x[2], v * x[3], v * x[4] }; + return ArraySize( x.ndim(), N ); +} +inline ArraySize operator-( const ArraySize &x, size_t v ) +{ + size_t N[5] = { x[0] - v, x[1] - v, x[2] - v, x[3] - v, x[4] - v }; + return ArraySize( x.ndim(), N ); +} +inline ArraySize operator+( const ArraySize &x, size_t v ) +{ + size_t N[5] = { x[0] + v, x[1] + v, x[2] + v, x[3] + v, x[4] + v }; + return ArraySize( x.ndim(), N ); +} +inline ArraySize operator+( size_t v, const ArraySize &x ) +{ + size_t N[5] = { x[0] + v, x[1] + v, x[2] + v, x[3] + v, x[4] + v }; + return ArraySize( x.ndim(), N ); +} + +#if defined( USING_ICC ) +ENABLE_WARNINGS +#endif + + #endif diff --git a/common/FunctionTable.cpp b/common/FunctionTable.cpp new file mode 100644 index 00000000..10577c0f --- /dev/null +++ b/common/FunctionTable.cpp @@ -0,0 +1,147 @@ +#include "FunctionTable.hpp" + + +/******************************************************** + * Random number generation * + ********************************************************/ +template<> char genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> int8_t genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> uint8_t genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> int16_t genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> uint16_t genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> int32_t genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> uint32_t genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> int64_t genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> uint64_t genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_int_distribution dis; + return dis( gen ); +} +template<> float genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_real_distribution dis; + return dis( gen ); +} +template<> double genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_real_distribution dis; + return dis( gen ); +} +template<> long double genRand() +{ + static std::random_device rd; + static std::mt19937 gen( rd() ); + static std::uniform_real_distribution dis; + return dis( gen ); +} + + +/******************************************************** + * axpy * + ********************************************************/ +template<> +void call_axpy( size_t N, const float alpha, const float *x, float *y ) +{ + ERROR("Not finished"); +} +template<> +void call_axpy( size_t N, const double alpha, const double *x, double *y ) +{ + ERROR("Not finished"); +} + + +/******************************************************** + * Multiply two arrays * + ********************************************************/ +template<> +void call_gemv( + size_t M, size_t N, double alpha, double beta, const double *A, const double *x, double *y ) +{ + ERROR("Not finished"); +} +template<> +void call_gemv( + size_t M, size_t N, float alpha, float beta, const float *A, const float *x, float *y ) +{ + ERROR("Not finished"); +} +template<> +void call_gemm( size_t M, + size_t N, + size_t K, + double alpha, + double beta, + const double *A, + const double *B, + double *C ) +{ + ERROR("Not finished"); +} +template<> +void call_gemm( size_t M, + size_t N, + size_t K, + float alpha, + float beta, + const float *A, + const float *B, + float *C ) +{ + ERROR("Not finished"); +} + diff --git a/common/FunctionTable.h b/common/FunctionTable.h index d0e8565d..9b162989 100644 --- a/common/FunctionTable.h +++ b/common/FunctionTable.h @@ -7,6 +7,7 @@ #include + /*! * Class FunctionTable is a serial function table class that defines * a series of operations that can be performed on the Array class. @@ -25,38 +26,55 @@ public: /*! * Perform a reduce operator y = f(x) - * @param[in] op The function operation - * Note: the operator is a template parameter - * (compared to a std::function to improve performance) - * @param[in] A The array to operate on - * @return The reduction + * @param[in] op The function operation + * Note: the operator is a template parameter to improve performance + * @param[in] A The array to operate on + * @param[in] initialValue The initial value for the reduction (0 for sum, +/- inf for min/max, + * ...) + * @return The reduction */ template - static inline TYPE reduce( LAMBDA &op, const Array &A ); + static inline TYPE reduce( LAMBDA &op, const Array &A, const TYPE &initialValue ); + + /*! + * Perform a reduce operator z = f(x,y) + * @param[in] op The function operation + * Note: the operator is a template parameter to improve performance + * @param[in] A The first array to operate on + * @param[in] B The second array to operate on + * @param[in] initialValue The initial value for the reduction (0 for sum, +/- inf for min/max, + * ...) + * @return The reduction + */ + template + static inline TYPE reduce( LAMBDA &op, + const Array &A, + const Array &B, + const TYPE &initialValue ); /*! * Perform a element-wise operation y = f(x) - * @param[in] fun The function operation - * Note: the operator is a template parameter - * (compared to a std::function to improve performance) - * @param[in] x The input array to operate on - * @param[out] y The output array + * @param[in] fun The function operation + * Note: the function is a template parameter to improve performance + * @param[in,out] x The array to operate on + * @param[out] y The output array */ template static inline void transform( LAMBDA &fun, const Array &x, Array &y ); /*! * Perform a element-wise operation z = f(x,y) - * @param[in] fun The function operation - * Note: the operator is a template parameter - * (compared to a std::function to improve performance) - * @param[in] x The first array - * @param[in] y The second array - * @param[out] z The result + * @param[in] fun The function operation + * Note: the function is a template parameter to improve performance + * @param[in] x The first array + * @param[in] y The second array + * @param[out] z The output array */ template - static inline void transform( - LAMBDA &fun, const Array &x, const Array &y, Array &z ); + static inline void transform( LAMBDA &fun, + const Array &x, + const Array &y, + Array &z ); /*! * Multiply two arrays @@ -65,8 +83,8 @@ public: * @param[out] c The output array */ template - static inline void multiply( - const Array &a, const Array &b, Array &c ); + static void + multiply( const Array &a, const Array &b, Array &c ); /*! * Perform dgemv/dgemm equavalent operation ( C = alpha*A*B + beta*C ) @@ -74,11 +92,14 @@ public: * @param[in] A The first array * @param[in] B The second array * @param[in] beta The scalar value alpha - * @param[in,out] c The output array C + * @param[in,out] C The output array C */ template - static void gemm( const TYPE alpha, const Array &A, const Array &B, - const TYPE beta, Array &C ); + static void gemm( const TYPE alpha, + const Array &A, + const Array &B, + const TYPE beta, + Array &C ); /*! * Perform axpy equavalent operation ( y = alpha*x + y ) @@ -98,9 +119,84 @@ public: template static bool equals( const Array &A, const Array &B, TYPE tol ); + template + static inline void gemmWrapper( char TRANSA, + char TRANSB, + int M, + int N, + int K, + TYPE alpha, + const TYPE *A, + int LDA, + const TYPE *B, + int LDB, + TYPE beta, + TYPE *C, + int LDC ); + + + /* Specialized Functions */ + + /*! + * Perform a element-wise operation y = max(x , 0) + * @param[in] A The input array + * @param[out] B The output array + */ + template + static void transformReLU( const Array &A, Array &B ); + + /*! + * Perform a element-wise operation B = |A| + * @param[in] A The array to operate on + * @param[out] B The output array + */ + template + static void transformAbs( const Array &A, Array &B ); + + /*! + * Perform a element-wise operation B = tanh(A) + * @param[in] A The array to operate on + * @param[out] B The output array + */ + template + static void transformTanh( const Array &A, Array &B ); + + /*! + * Perform a element-wise operation B = max(-1 , min(1 , A) ) + * @param[in] A The array to operate on + * @param[out] B The output array + */ + template + static void transformHardTanh( const Array &A, Array &B ); + + /*! + * Perform a element-wise operation B = 1 / (1 + exp(-A)) + * @param[in] A The array to operate on + * @param[out] B The output array + */ + template + static void transformSigmoid( const Array &A, Array &B ); + + /*! + * Perform a element-wise operation B = log(exp(A) + 1) + * @param[in] A The array to operate on + * @param[out] B The output array + */ + template + static void transformSoftPlus( const Array &A, Array &B ); + + /*! + * Sum the elements of the Array + * @param[in] A The array to sum + */ + template + static TYPE sum( const Array &A ); private: FunctionTable(); + + template + static inline void rand( size_t N, T *x ); }; diff --git a/common/FunctionTable.hpp b/common/FunctionTable.hpp index e5fc7786..6a611f00 100644 --- a/common/FunctionTable.hpp +++ b/common/FunctionTable.hpp @@ -2,7 +2,7 @@ #define included_FunctionTable_hpp #include "common/FunctionTable.h" -#include "common/Utilities.h" +#include "common/UtilityMacros.h" #include #include @@ -10,33 +10,16 @@ #include + /******************************************************** * Random number initialization * ********************************************************/ -template -static inline typename std::enable_if::value>::type genRand( - size_t N, TYPE* x ) -{ - std::random_device rd; - std::mt19937 gen( rd() ); - std::uniform_int_distribution dis; - for ( size_t i = 0; i < N; i++ ) - x[i] = dis( gen ); -} -template -static inline typename std::enable_if::value>::type genRand( - size_t N, TYPE* x ) -{ - std::random_device rd; - std::mt19937 gen( rd() ); - std::uniform_real_distribution dis( 0, 1 ); - for ( size_t i = 0; i < N; i++ ) - x[i] = dis( gen ); -} +template TYPE genRand(); template -inline void FunctionTable::rand( Array& x ) +inline void FunctionTable::rand( Array &x ) { - genRand( x.length(), x.data() ); + for ( size_t i = 0; i < x.length(); i++ ) + x( i ) = genRand(); } @@ -44,24 +27,39 @@ inline void FunctionTable::rand( Array& x ) * Reduction * ********************************************************/ template -inline TYPE FunctionTable::reduce( LAMBDA& op, const Array& A ) +inline TYPE FunctionTable::reduce( LAMBDA &op, const Array &A, const TYPE &initialValue ) { if ( A.length() == 0 ) return TYPE(); - const TYPE* x = A.data(); - TYPE y = x[0]; - const size_t N = A.length(); - for ( size_t i = 1; i < N; i++ ) + const TYPE *x = A.data(); + TYPE y = initialValue; + for ( size_t i = 0; i < A.length(); i++ ) y = op( x[i], y ); return y; } +template +inline TYPE FunctionTable::reduce( LAMBDA &op, + const Array &A, + const Array &B, + const TYPE &initialValue ) +{ + ARRAY_ASSERT( A.length() == B.length() ); + if ( A.length() == 0 ) + return TYPE(); + const TYPE *x = A.data(); + const TYPE *y = B.data(); + TYPE z = initialValue; + for ( size_t i = 0; i < A.length(); i++ ) + z = op( x[i], y[i], z ); + return z; +} /******************************************************** * Unary transformation * ********************************************************/ template -inline void FunctionTable::transform( LAMBDA& fun, const Array& x, Array& y ) +inline void FunctionTable::transform( LAMBDA &fun, const Array &x, Array &y ) { y.resize( x.size() ); const size_t N = x.length(); @@ -69,8 +67,10 @@ inline void FunctionTable::transform( LAMBDA& fun, const Array& x, Ar y( i ) = fun( x( i ) ); } template -inline void FunctionTable::transform( - LAMBDA& fun, const Array& x, const Array& y, Array& z ) +inline void FunctionTable::transform( LAMBDA &fun, + const Array &x, + const Array &y, + Array &z ) { if ( x.size() != y.size() ) throw std::logic_error( "Sizes of x and y do not match" ); @@ -85,25 +85,19 @@ inline void FunctionTable::transform( * axpy * ********************************************************/ template -inline void call_axpy( size_t N, const TYPE alpha, const TYPE* x, TYPE* y ); +void call_axpy( size_t N, const TYPE alpha, const TYPE *x, TYPE *y ); template<> -inline void call_axpy( size_t, const float, const float*, float* ) -{ - throw std::logic_error( "LapackWrappers not configured" ); -} +void call_axpy( size_t N, const float alpha, const float *x, float *y ); template<> -inline void call_axpy( size_t, const double, const double*, double* ) -{ - throw std::logic_error( "LapackWrappers not configured" ); -} +void call_axpy( size_t N, const double alpha, const double *x, double *y ); template -inline void call_axpy( size_t N, const TYPE alpha, const TYPE* x, TYPE* y ) +void call_axpy( size_t N, const TYPE alpha, const TYPE *x, TYPE *y ) { for ( size_t i = 0; i < N; i++ ) y[i] += alpha * x[i]; } template -void FunctionTable::axpy( const TYPE alpha, const Array& x, Array& y ) +void FunctionTable::axpy( const TYPE alpha, const Array &x, Array &y ) { if ( x.size() != y.size() ) throw std::logic_error( "Array sizes do not match" ); @@ -115,21 +109,15 @@ void FunctionTable::axpy( const TYPE alpha, const Array& x, Array -inline void call_gemv( size_t M, size_t N, TYPE alpha, TYPE beta, const TYPE* A, const TYPE* x, TYPE* y ); +void call_gemv( size_t M, size_t N, TYPE alpha, TYPE beta, const TYPE *A, const TYPE *x, TYPE *y ); template<> -inline void call_gemv( - size_t, size_t, double, double, const double*, const double*, double* ) -{ - throw std::logic_error( "LapackWrappers not configured" ); -} +void call_gemv( + size_t M, size_t N, double alpha, double beta, const double *A, const double *x, double *y ); template<> -inline void call_gemv( size_t, size_t, float, float, const float*, const float*, float* ) -{ - throw std::logic_error( "LapackWrappers not configured" ); -} +void call_gemv( + size_t M, size_t N, float alpha, float beta, const float *A, const float *x, float *y ); template -inline void call_gemv( - size_t M, size_t N, TYPE alpha, TYPE beta, const TYPE* A, const TYPE* x, TYPE* y ) +void call_gemv( size_t M, size_t N, TYPE alpha, TYPE beta, const TYPE *A, const TYPE *x, TYPE *y ) { for ( size_t i = 0; i < M; i++ ) y[i] = beta * y[i]; @@ -139,21 +127,29 @@ inline void call_gemv( } } template -inline void call_gemm( - size_t M, size_t N, size_t K, TYPE alpha, TYPE beta, const TYPE* A, const TYPE* B, TYPE* C ); +void call_gemm( + size_t M, size_t N, size_t K, TYPE alpha, TYPE beta, const TYPE *A, const TYPE *B, TYPE *C ); template<> -inline void call_gemm( size_t, size_t, size_t, double, double, const double*, const double*, double* ) -{ - throw std::logic_error( "LapackWrappers not configured" ); -} +void call_gemm( size_t M, + size_t N, + size_t K, + double alpha, + double beta, + const double *A, + const double *B, + double *C ); template<> -inline void call_gemm( size_t, size_t, size_t, float, float, const float*, const float*, float* ) -{ - throw std::logic_error( "LapackWrappers not configured" ); -} +void call_gemm( size_t M, + size_t N, + size_t K, + float alpha, + float beta, + const float *A, + const float *B, + float *C ); template -inline void call_gemm( - size_t M, size_t N, size_t K, TYPE alpha, TYPE beta, const TYPE* A, const TYPE* B, TYPE* C ) +void call_gemm( + size_t M, size_t N, size_t K, TYPE alpha, TYPE beta, const TYPE *A, const TYPE *B, TYPE *C ) { for ( size_t i = 0; i < K * M; i++ ) C[i] = beta * C[i]; @@ -165,16 +161,17 @@ inline void call_gemm( } } template -void FunctionTable::gemm( const TYPE alpha, const Array& a, const Array& b, - const TYPE beta, Array& c ) +void FunctionTable::gemm( const TYPE alpha, + const Array &a, + const Array &b, + const TYPE beta, + Array &c ) { + if ( a.size( 1 ) != b.size( 0 ) ) + throw std::logic_error( "Inner dimensions must match" ); if ( a.ndim() == 2 && b.ndim() == 1 ) { - if ( a.size( 1 ) != b.size( 0 ) ) - throw std::logic_error( "Inner dimensions must match" ); call_gemv( a.size( 0 ), a.size( 1 ), alpha, beta, a.data(), b.data(), c.data() ); } else if ( a.ndim() <= 2 && b.ndim() <= 2 ) { - if ( a.size( 1 ) != b.size( 0 ) ) - throw std::logic_error( "Inner dimensions must match" ); call_gemm( a.size( 0 ), a.size( 1 ), b.size( 1 ), alpha, beta, a.data(), b.data(), c.data() ); } else { @@ -182,17 +179,16 @@ void FunctionTable::gemm( const TYPE alpha, const Array& a, const Arr } } template -void FunctionTable::multiply( - const Array& a, const Array& b, Array& c ) +void FunctionTable::multiply( const Array &a, + const Array &b, + Array &c ) { + if ( a.size( 1 ) != b.size( 0 ) ) + throw std::logic_error( "Inner dimensions must match" ); if ( a.ndim() == 2 && b.ndim() == 1 ) { - if ( a.size( 1 ) != b.size( 0 ) ) - throw std::logic_error( "Inner dimensions must match" ); c.resize( a.size( 0 ) ); call_gemv( a.size( 0 ), a.size( 1 ), 1, 0, a.data(), b.data(), c.data() ); } else if ( a.ndim() <= 2 && b.ndim() <= 2 ) { - if ( a.size( 1 ) != b.size( 0 ) ) - throw std::logic_error( "Inner dimensions must match" ); c.resize( a.size( 0 ), b.size( 1 ) ); call_gemm( a.size( 0 ), a.size( 1 ), b.size( 1 ), 1, 0, a.data(), b.data(), c.data() ); @@ -206,8 +202,8 @@ void FunctionTable::multiply( * Check if two arrays are equal * ********************************************************/ template -inline typename std::enable_if::value, bool>::type -FunctionTableCompare( const Array& a, const Array& b, TYPE ) +inline typename std::enable_if::value, bool>::type +FunctionTableCompare( const Array &a, const Array &b, TYPE ) { bool pass = true; if ( a.size() != b.size() ) @@ -218,7 +214,7 @@ FunctionTableCompare( const Array& a, const Array& b, TYPE } template inline typename std::enable_if::value, bool>::type -FunctionTableCompare( const Array& a, const Array& b, TYPE tol ) +FunctionTableCompare( const Array &a, const Array &b, TYPE tol ) { bool pass = true; if ( a.size() != b.size() ) @@ -228,10 +224,89 @@ FunctionTableCompare( const Array& a, const Array& b, TYPE return pass; } template -bool FunctionTable::equals( const Array& a, const Array& b, TYPE tol ) +bool FunctionTable::equals( const Array &a, const Array &b, TYPE tol ) { return FunctionTableCompare( a, b, tol ); } +/******************************************************** + * Specialized Functions * + ********************************************************/ +template +void FunctionTable::transformReLU( const Array &A, Array &B ) +{ + const auto &fun = []( const TYPE &a ) { return std::max( a, static_cast( 0 ) ); }; + transform( fun, A, B ); +} + +template +void FunctionTable::transformAbs( const Array &A, Array &B ) +{ + B.resize( A.size() ); + const auto &fun = []( const TYPE &a ) { return std::abs( a ); }; + transform( fun, A, B ); +} +template +void FunctionTable::transformTanh( const Array &A, Array &B ) +{ + B.resize( A.size() ); + const auto &fun = []( const TYPE &a ) { return tanh( a ); }; + transform( fun, A, B ); +} + +template +void FunctionTable::transformHardTanh( const Array &A, + Array &B ) +{ + B.resize( A.size() ); + const auto &fun = []( const TYPE &a ) { + return std::max( -static_cast( 1.0 ), std::min( static_cast( 1.0 ), a ) ); + }; + transform( fun, A, B ); +} + +template +void FunctionTable::transformSigmoid( const Array &A, Array &B ) +{ + B.resize( A.size() ); + const auto &fun = []( const TYPE &a ) { return 1.0 / ( 1.0 + exp( -a ) ); }; + transform( fun, A, B ); +} + +template +void FunctionTable::transformSoftPlus( const Array &A, + Array &B ) +{ + B.resize( A.size() ); + const auto &fun = []( const TYPE &a ) { return log1p( exp( a ) ); }; + transform( fun, A, B ); +} + +template +TYPE FunctionTable::sum( const Array &A ) +{ + const auto &fun = []( const TYPE &a, const TYPE &b ) { return a + b; }; + return reduce( fun, A, (TYPE) 0 ); +} + +template +inline void FunctionTable::gemmWrapper( char TRANSA, + char TRANSB, + int M, + int N, + int K, + TYPE alpha, + const TYPE *A, + int LDA, + const TYPE *B, + int LDB, + TYPE beta, + TYPE *C, + int LDC ) +{ + ERROR("Not finished"); +} + + #endif diff --git a/tests/TestWriter.cpp b/tests/TestWriter.cpp index 652c3d4c..76518ac0 100644 --- a/tests/TestWriter.cpp +++ b/tests/TestWriter.cpp @@ -94,13 +94,13 @@ bool checkVar( const std::string &format, std::shared_ptr mesh, { if ( format == "new" ) IO::reformatVariable( *mesh, *variable2 ); - bool pass = true; - const IO::Variable &var1 = *variable1; - const IO::Variable &var2 = *variable2; - pass = var1.name == var2.name; - pass = pass && var1.dim == var2.dim; - pass = pass && var1.type == var2.type; - pass = pass && var1.data.length() == var2.data.length(); + bool pass = true; + const auto &var1 = *variable1; + const auto &var2 = *variable2; + pass = var1.name == var2.name; + pass = pass && var1.dim == var2.dim; + pass = pass && var1.type == var2.type; + pass = pass && var1.data.length() == var2.data.length(); if ( pass ) { for ( size_t m = 0; m < var1.data.length(); m++ ) pass = pass && approx_equal( var1.data( m ), var2.data( m ) ); @@ -133,6 +133,12 @@ void testWriter( } else if ( format == "silo-float" ) { format2 = "silo"; precision = IO::DataType::Float; + } else if ( format == "hdf5-double" ) { + format2 = "hdf5"; + precision = IO::DataType::Double; + } else if ( format == "hdf5-float" ) { + format2 = "hdf5"; + precision = IO::DataType::Float; } @@ -315,7 +321,7 @@ int main( int argc, char **argv ) set_node_vec->data = point_node_vec->data; list_node_mag->data.resize( 3 * N_tri ); list_node_vec->data.resize( 3 * N_tri, 3 ); - for ( int i = 0; i < N_points; i++ ) { + for ( int i = 0; i < N_tri; i++ ) { list_node_mag->data( 3 * i + 0 ) = distance( trilist->A[i] ); list_node_mag->data( 3 * i + 1 ) = distance( trilist->B[i] ); list_node_mag->data( 3 * i + 2 ) = distance( trilist->C[i] ); @@ -398,6 +404,8 @@ int main( int argc, char **argv ) testWriter( "new", meshData, ut ); testWriter( "silo-double", meshData, ut ); testWriter( "silo-float", meshData, ut ); + testWriter( "hdf5-double", meshData, ut ); + testWriter( "hdf5-float", meshData, ut ); // Finished ut.report();