/******************************************************************************* * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * * Copyright (C) 2011 by O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software * Foundation, either version 3 of the License, or (at your option) any later * version. * * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more * details. * * You should have received a copy of the GNU General Public License along with * TRIQS. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ #ifndef TRIQS_ARRAYS_HDF5_ARRAY_PROXY_IMPL_H #define TRIQS_ARRAYS_HDF5_ARRAY_PROXY_IMPL_H //#define TRIQS_ARRAYS_DEBUG_H5_SLICE #include "../indexmaps/cuboid/domain.hpp" #include "./group_or_file.hpp" #include "../impl/sliceable_object.hpp" #include "../impl/tuple_tools.hpp" #include <boost/shared_ptr.hpp> #include <boost/make_shared.hpp> namespace triqs { namespace arrays { namespace h5 { template<int Rf> struct array_proxy_option {}; } namespace Option { template<class IM2, int Rf> struct Im_Opt_2_Opt<IM2, h5::array_proxy_option<Rf> > { typedef h5::array_proxy_option<Rf> type;}; } namespace h5 { template<int Rank, int Rank_full> class index_system { public : static const unsigned int rank= Rank, rank_full = Rank_full; typedef mini_vector<size_t, Rank> v_type; typedef mini_vector<size_t, Rank_full> v_type_full; typedef indexmaps::cuboid::domain<Rank> domain_type; domain_type const & domain() const { return mydomain;} index_system (v_type_full const & total_lengths_) { total_lens_ = total_lengths_; lens_= total_lengths_; for (size_t i =0; i<rank_full; ++i) stri_[i] = 1; for (size_t i =0; i<rank; ++i) dims[i] = lens_[i]; mydomain = domain_type (dims); } index_system (v_type const & dims_, v_type_full const & total_lengths_, v_type_full const & lengths_, v_type_full const & strides_, v_type_full const & offset_ ) { dims = dims_;total_lens_ = total_lengths_; lens_= lengths_; stri_ = strides_; off_=offset_; mydomain = domain_type (dims); } index_system (DataSpace const & ds, bool is_complex) { int rf = ds.getSimpleExtentNdims(); if ( rf != rank_full + (is_complex ? 1 : 0) ) TRIQS_RUNTIME_ERROR << "H5 : dimension error"; //int ndims = ds.getSimpleExtentDims( &lens_[0], NULL); ds.getSimpleExtentDims( &lens_[0], NULL); for (size_t i =0; i<rank; ++i) { dims[i] = lens_[i]; stri_[i] = 1; off_[i]= 0; } total_lens_=dims; mydomain = domain_type (dims); } template<bool C> DataSpace dataspace() const { return h5::dataspace_from_LS<rank_full,C> (total_lens_, lens_,stri_,off_); } size_t size() const { size_t _size = 1; for (size_t i =0; i<rank; ++i) _size *= lens_[i]; return _size;} mini_vector<size_t, rank_full> total_lengths() const { return this->total_lens_;} mini_vector<size_t, rank_full> lengths() const { return this->lens_;} mini_vector<size_t, rank_full> strides() const { return this->stri_;} private: mini_vector<hsize_t,rank_full> lens_, off_, stri_; v_type dims; v_type_full total_lens_; domain_type mydomain; }; namespace slicer_impl { template <bool BC> inline void _check_BC ( int N, int ind, size_t B) { } template <> inline void _check_BC<true> (int N, int ind, size_t B) { bool cond = (ind >= 0) && (ind < int(B)); // fix this int ... if (!cond) TRIQS_ARRAYS_KEY_ERROR << " index "<<N<<" is out of domain: \n " << ind <<" is not within [0,"<< B <<"[\n"; } template<int Rank_in, int Rank_out, int N,int P, int c, bool BoundCheck> struct slice_calc { typedef mini_vector<size_t,Rank_in> const & i_type; typedef mini_vector<size_t,Rank_in > & o_type; typedef mini_vector<size_t,Rank_out > & os_type; template< typename ArgsTuple> static void invoke(i_type li, i_type si, os_type lo_c, o_type lo, o_type so, o_type offset, ArgsTuple const & args ) { const int dP = boost::is_base_of<typename boost::tuples::element<N,ArgsTuple>::type, range >::type::value ; one_step(li[N], si[N],lo[N],so[N], offset[N] ,boost::tuples::get<N>(args)); lo_c[P] = lo[N]; slice_calc<Rank_in,Rank_out,N+1,P+dP,c-1, BoundCheck>::invoke(li,si,lo_c,lo,so, offset, args); } static void one_step(size_t li, size_t si, size_t & lo, size_t & so, size_t & offset, size_t R){ _check_BC <BoundCheck> (N, R, li); offset = R; lo =1; so= 1; } static void one_step(size_t li, size_t si, size_t & lo, size_t & so, size_t & offset, range R){ _check_BC <BoundCheck> (N, R.first(),li); _check_BC <BoundCheck> (N, (R.last()==-1 ? li : R.last()) -1 ,li); lo = ((R.last()==-1 ? li : R.last()) - R.first() + R.step()-1 )/R.step() ; // python behaviour so = R.step(); offset = R.first() ; } }; // stop the recursion template<int Ri, int Ro, int N, int P, bool BC> struct slice_calc <Ri,Ro,N,P,0,BC> : slice_calc<Ri,Ro,N,P,1,BC> { template<class T1,class T2,class T3,class T4,class T5,class T6, class T7> static void invoke(T1,T2,T3,T4,T5,T6,T7 ) {} }; }//namespace slicer_impl } // h5 namespace indexmaps { template<int R, int Rf, typename ArgsTuple> struct slicer < h5::index_system<R,Rf>, ArgsTuple> { static const unsigned int len = boost::tuples::length<ArgsTuple>::value; static_assert(len>=R, "Too few arguments in slice"); static_assert(len<=R, "Too many arguments in slice"); static_assert( (R==Rf), "Can not slice array_proxy twice (not implemented)"); static const unsigned int R2 = R - TupleTools::CountHowManyInt<ArgsTuple>::value; typedef h5::index_system< R2, Rf> return_type; static return_type invoke ( h5::index_system<R,R> const & X, ArgsTuple args) { mini_vector<size_t, R2> newdims; mini_vector<size_t,R> newoffset, newlengths, newstrides; h5::slicer_impl::slice_calc<R,R2,0,0,R,true>::invoke(X.lengths(),X.strides(),newdims, newlengths,newstrides, newoffset, args); #ifdef TRIQS_ARRAYS_DEBUG_H5_SLICE std::cerr<<"-----------------------------------------------"<<std::endl; std::cerr<<"H5 Slicing "<< X.lengths().to_string()<<X.strides().to_string()<<newlengths.to_string()<<newstrides.to_string()<< newoffset.to_string() << args<<std::endl; #endif return return_type(newdims,X.total_lengths(), newlengths,newstrides,newoffset); }; }; template<int R, int Rf> struct slicer < h5::index_system<R,Rf>, boost::tuple<> > { typedef h5::index_system< R,Rf> return_type;}; } }} #endif