From 9edda8724dcead3460ecdc0b5a1cfcad18c914ad Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Fri, 18 Oct 2013 13:39:00 +0200 Subject: [PATCH] implement gf_const_view --- test/triqs/gfs/gf_scalar.cpp | 2 +- triqs/gfs/block.hpp | 130 ++++++++------- triqs/gfs/curry.hpp | 25 +-- triqs/gfs/data_proxies.hpp | 57 +++++-- triqs/gfs/gf.hpp | 215 ++++++++++++++++--------- triqs/gfs/local/fourier_matsubara.cpp | 12 +- triqs/gfs/local/fourier_matsubara.hpp | 24 +-- triqs/gfs/local/fourier_real.cpp | 12 +- triqs/gfs/local/fourier_real.hpp | 26 +-- triqs/gfs/local/functions.cpp | 2 +- triqs/gfs/local/functions.hpp | 2 +- triqs/gfs/local/legendre_matsubara.cpp | 16 +- triqs/gfs/local/legendre_matsubara.hpp | 16 +- triqs/gfs/meshes/product.hpp | 6 +- triqs/gfs/product.hpp | 17 +- 15 files changed, 345 insertions(+), 217 deletions(-) diff --git a/test/triqs/gfs/gf_scalar.cpp b/test/triqs/gfs/gf_scalar.cpp index 25747eb1..bcc47dd7 100644 --- a/test/triqs/gfs/gf_scalar.cpp +++ b/test/triqs/gfs/gf_scalar.cpp @@ -15,7 +15,7 @@ int main() { // test hdf5 H5::H5File file("gf_scalar.h5", H5F_ACC_TRUNC); h5_write(file, "g", G); - h5_write(file, "gm", reinterpret_scalar_valued_gf_as_matrix_valued(G())); + h5_write(file, "gm", reinterpret_scalar_valued_gf_as_matrix_valued(G)); } TRIQS_CATCH_AND_ABORT; diff --git a/triqs/gfs/block.hpp b/triqs/gfs/block.hpp index 112f2c2f..d9c44f82 100644 --- a/triqs/gfs/block.hpp +++ b/triqs/gfs/block.hpp @@ -43,21 +43,24 @@ namespace triqs { namespace gfs { template struct h5_name { static std::string invoke(){ return "BlockGf";}}; - template struct h5_rw { + template struct h5_rw { - static void write (h5::group gr, gf_impl const & g) { + static void write (h5::group gr, gf_const_view g) { for (size_t i =0; i & g) { + template + static void read (h5::group gr, gf_impl &g) { // does not work : need to read the block name and remake the mesh... g._mesh = gf_mesh (gr.get_all_subgroup_names()); g._data.resize(g._mesh.size()); + //if (g._data.size() != g._mesh.size()) TRIQS_RUNTIME_ERROR << "h5 read block gf : number of block mismatch"; for (size_t i =0; i using block_gf = gf>; + template using block_gf_view = gf_view>; + template using block_gf_const_view = gf_const_view>; + // ------------------------------- Free Factories for regular type -------------------------------------------------- + // from a number and a gf to be copied + template + block_gf make_block_gf(int n, gf const & g) { + auto V = std::vector>{}; + for (int i = 0; i < n; ++i) V.push_back(g); + return {{n}, std::move(V), nothing{}, nothing{}}; + } + + // from a vector of gf (moving directly) + template + block_gf make_block_gf(std::vector> &&V) { + return {{int(V.size())}, std::move(V), nothing{}, nothing{}}; + } + + // from a vector of gf : generalized to have a different type of gf in the vector (e.g. views...) + template + block_gf make_block_gf(std::vector const &V) { + auto V2 = std::vector>{}; + for (auto const &g : V) V2.push_back(g); + return {{int(V.size())}, std::move(V2), nothing{}, nothing{}}; + } + + // from a init list of GF with the correct type + template + block_gf make_block_gf(std::initializer_list> const &V) { + return {{int(V.size())}, V, nothing{}, nothing{}}; + } + + // from vector and a gf to be copied + template + block_gf make_block_gf(std::vector const &block_names, gf const &g) { + auto V = std::vector>{}; + for (int i = 0; i < block_names.size(); ++i) V.push_back(g); + return {{block_names}, std::move(V), nothing{}, nothing{}}; + } + + // from vector, vector + template + block_gf make_block_gf(std::vector const &block_names, + std::vector> V) { + if (block_names.size() != V.size()) + TRIQS_RUNTIME_ERROR << "make_block_gf(vector, vector) : the two vectors do not have the same size !"; + return {{block_names}, std::move(V), nothing{}, nothing{}}; + } + + // from vector, init_list + template + block_gf make_block_gf(std::vector const &block_names, + std::initializer_list> const &V) { + if (block_names.size() != V.size()) TRIQS_RUNTIME_ERROR << "make_block_gf(vector, init_list) : size mismatch !"; + return {{block_names}, V, nothing{}, nothing{}}; + } + + // ------------------------------- Free Factories for view type -------------------------------------------------- + template gf_view::type::view_type> make_block_gf_view(G0 && g0, G && ... g) { auto V = std::vector::type::view_type>{std::forward(g0), std::forward(g)...}; @@ -87,54 +151,6 @@ namespace triqs { namespace gfs { //return { gf_mesh {int(V.size())}, std::move(V), nothing{}, nothing{} } ; } - // from a number and a gf to be copied - template - gf> - make_block_gf(int n, gf g) { - auto V = std::vector>{}; - for (int i =0; i - gf> - make_block_gf(std::vector> && V) { - return { {int(V.size())}, std::move(V), nothing{}, nothing{}}; - } - - // from a vector of gf : generalized to have a different type of gf in the vector (e.g. views...) - template - gf> - make_block_gf(std::vector const & V) { - auto V2 = std::vector>{}; - for (auto const & g : V) V2.push_back(g); - return { {int(V.size())}, std::move(V2), nothing{}, nothing{}}; - } - - // from a init list of GF with the correct type - template - gf> - make_block_gf(std::initializer_list> const & V) { return { {int(V.size())}, V, nothing{}, nothing{}} ; } - - // from vector, vector - template - gf> - make_block_gf(std::vector const & block_names, std::vector> V) { - if (block_names.size()!=V.size()) TRIQS_RUNTIME_ERROR << "make_block_gf(vector, vector) : the two vectors do not have the same size !"; - return { {block_names}, std::move(V), nothing{}, nothing{}}; - } - - // from vector, init_list - template - gf> - make_block_gf(std::vector const & block_names, std::initializer_list> const & V) { - if (block_names.size()!=V.size()) TRIQS_RUNTIME_ERROR << "make_block_gf(vector, init_list) : size mismatch !"; - return { {block_names}, V, nothing{}, nothing{}}; - } - - // ------------------------------- Free Factories for view type -------------------------------------------------- - template gf_view make_block_gf_view_from_vector (std::vector V) { return { {int(V.size())}, std::move(V), nothing{}, nothing{}} ; } @@ -145,8 +161,6 @@ namespace triqs { namespace gfs { template size_t n_blocks (gf const & g) { return g.mesh().size();} template size_t n_blocks (gf_view const & g) { return g.mesh().size();} - template using block_gf = gf>; - // an iterator over the block template class block_gf_iterator : @@ -165,11 +179,11 @@ namespace triqs { namespace gfs { bool at_end() const { return mesh_it.at_end();} }; - template - block_gf_iterator begin(gf_impl const & bgf) { return {bgf,false};} + template + block_gf_iterator begin(gf_impl const & bgf) { return {bgf,false};} - template - block_gf_iterator end(gf_impl const & bgf) { return {bgf,true};} + template + block_gf_iterator end(gf_impl const & bgf) { return {bgf,true};} }} #endif diff --git a/triqs/gfs/curry.hpp b/triqs/gfs/curry.hpp index 15f73151..1777c8fd 100644 --- a/triqs/gfs/curry.hpp +++ b/triqs/gfs/curry.hpp @@ -49,9 +49,9 @@ namespace triqs { namespace gfs { template auto rm_tuple_of_size_one(std::tuple const & t) DECL_AND_RETURN(t); template auto rm_tuple_of_size_one(std::tuple const & t) DECL_AND_RETURN(std::get<0>(t)); - template - gf_view< typename cart_prod_impl< triqs::tuple::filter_out_t, pos...>>::type ,Target, Opt> - partial_eval(gf_impl< cartesian_product, Target,Opt,B> const & g, IT index) { + template + gf_view< typename cart_prod_impl< triqs::tuple::filter_out_t, pos...>>::type ,Target, Opt,IsConst> + partial_eval(gf_impl< cartesian_product, Target,Opt,B,IsConst> const & g, IT index) { // meshes of the returned gf_view : just drop the mesh of the evaluated variables auto meshes_tuple_partial = triqs::tuple::filter_out(g.mesh().components()); // a view of the array of g, with the dimension sizeof...(Ms) @@ -61,7 +61,7 @@ namespace triqs { namespace gfs { // from it, we make a slice of the array of g, corresponding to the data of the returned gf_view auto arr2 = triqs::tuple::apply(arr, arr_args); // finally, we build the view on this data. - using r_t = gf_view< cart_prod< triqs::tuple::filter_out_t, pos...>> ,Target, Opt>; + using r_t = gf_view< cart_prod< triqs::tuple::filter_out_t, pos...>> ,Target, Opt,IsConst>; return r_t{ rm_tuple_of_size_one(meshes_tuple_partial), arr2, typename r_t::singularity_non_view_t{}, typename r_t::symmetry_t{} }; } @@ -78,18 +78,23 @@ namespace triqs { namespace gfs { }; // curry function ... - template - gf_view,pos...> >, - lambda_valued, Target,Opt>,pos...>>, - Opt> - curry (gf_impl, Target,Opt,B> const & g) { + template + gf_view, pos...>>, + lambda_valued, Target, Opt,IsConst>, pos...>>, Opt, IsConst> + curry(gf_view, Target, Opt, IsConst> g) { // pick up the meshed corresponding to the curryed variables auto meshes_tuple = triqs::tuple::filter(g.mesh().components()); // building the view - return {rm_tuple_of_size_one(meshes_tuple),curry_polymorphic_lambda, Target,Opt>, pos ...>{g}, nothing(), nothing()}; + return {rm_tuple_of_size_one(meshes_tuple),curry_polymorphic_lambda, Target,Opt,IsConst>, pos ...>{g}, nothing(), nothing()}; //using m_t = gf_mesh< cart_prod< triqs::tuple::filter_t,pos...>>>; //return {triqs::tuple::apply_construct(meshes_tuple),curry_polymorphic_lambda, Target,Opt>, pos ...>{g}, nothing(), nothing()}; }; + + template + auto curry(gf, Target, Opt> & g) DECL_AND_RETURN(curry(g())); + + template + auto curry(gf, Target, Opt> const & g) DECL_AND_RETURN(curry(g())); } // gf_implementation using gfs_implementation::partial_eval; diff --git a/triqs/gfs/data_proxies.hpp b/triqs/gfs/data_proxies.hpp index ca37a440..c61ed7fe 100644 --- a/triqs/gfs/data_proxies.hpp +++ b/triqs/gfs/data_proxies.hpp @@ -32,24 +32,29 @@ namespace triqs { namespace gfs { template struct data_proxy_array { /// The storage - typedef arrays::array storage_t; - typedef typename storage_t::view_type storage_view_t; + typedef arrays::array storage_t; + typedef typename storage_t::view_type storage_view_t; + typedef typename storage_t::const_view_type storage_const_view_t; /// The data access auto operator()(storage_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i)); auto operator()(storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); auto operator()(storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i)); auto operator()(storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); + auto operator()(storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); + auto operator()(storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i)); #ifdef TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS auto operator()(storage_t & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); auto operator()(storage_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); auto operator()(storage_view_t & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); auto operator()(storage_view_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); + auto operator()(storage_const_view_t & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); + auto operator()(storage_const_view_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); #endif template static void assign_to_scalar (S & data, RHS && rhs) { data() = std::forward(rhs);} - template static void rebind (storage_view_t & data, RHS && rhs) { data.rebind(rhs.data()); } + template static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); } }; //---------------------------- 3d array : returns matrices in this case ! ---------------------------------- @@ -58,15 +63,18 @@ namespace triqs { namespace gfs { /// The storage typedef arrays::array storage_t; typedef typename storage_t::view_type storage_view_t; + typedef typename storage_t::const_view_type storage_const_view_t; /// The data access auto operator()(storage_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i)); auto operator()(storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); auto operator()(storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i)); auto operator()(storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); + auto operator()(storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); + auto operator()(storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i)); #ifdef TRIQS_DATA_PROXIES_OLD_MATRIX_VIEW_PROXY - arrays::matrix_view_proxy operator()(storage_t & data, size_t i) const { return arrays::matrix_view_proxy(data,i); } + arrays::matrix_view_proxy operator()(storage_t & data, size_t i) const { return arrays::matrix_view_proxy(data,i); } arrays::const_matrix_view_proxy operator()(storage_t const & data, size_t i) const { return arrays::const_matrix_view_proxy(data,i); } arrays::matrix_view_proxy operator()(storage_view_t & data, size_t i) const { return arrays::matrix_view_proxy(data,i); } arrays::const_matrix_view_proxy operator()(storage_view_t const & data, size_t i) const { return arrays::const_matrix_view_proxy(data,i); } @@ -78,10 +86,13 @@ namespace triqs { namespace gfs { auto operator()(storage_view_t & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); auto operator()(storage_view_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); + + auto operator()(storage_const_view_t & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); + auto operator()(storage_const_view_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis())); #endif - template static void assign_to_scalar (S & data, RHS && rhs) { data() = std::forward(rhs);} - template static void rebind (storage_view_t & data, RHS && rhs) { data.rebind(rhs.data()); } + template static void assign_to_scalar(S& data, RHS&& rhs) { data() = std::forward(rhs); } + template static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); } }; //---------------------------- 1d array ---------------------------------- @@ -90,34 +101,57 @@ namespace triqs { namespace gfs { /// The storage typedef arrays::array storage_t; typedef typename storage_t::view_type storage_view_t; + typedef typename storage_t::const_view_type storage_const_view_t; /// The data access auto operator()(storage_t & data,size_t i) const -> decltype(data(i)) { return data(i);} auto operator()(storage_t const & data,size_t i) const -> decltype(data(i)) { return data(i);} auto operator()(storage_view_t & data,size_t i) const -> decltype(data(i)) { return data(i);} auto operator()(storage_view_t const & data,size_t i) const -> decltype(data(i)) { return data(i);} + auto operator()(storage_const_view_t & data,size_t i) const -> decltype(data(i)) { return data(i);} + auto operator()(storage_const_view_t const & data,size_t i) const -> decltype(data(i)) { return data(i);} template static void assign_to_scalar (S & data, RHS && rhs) { data() = std::forward(rhs);} - template static void rebind (storage_view_t & data, RHS && rhs) { data.rebind(rhs.data()); } + template static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); } }; //---------------------------- vector ---------------------------------- + template struct view_proxy : public V { + view_proxy() : V(typename V::regular_type()) {} + view_proxy(V const &v) : V(v){}; + view_proxy(view_proxy const & p) : V(p) {}; + template explicit view_proxy(Args && ... args) : V (std::forward(args)...){} + view_proxy & operator = ( view_proxy const & cp ) { this->rebind(cp); return *this;} + view_proxy & operator = ( V const & v ) { this->rebind(v); return *this;} + using V::operator=; + //template view_proxy & operator = (X && x) { V::operator=( std::forward(x) ); return *this;} + }; + template struct data_proxy_vector { typedef typename T::view_type Tv; + typedef typename T::const_view_type Tcv; /// The storage typedef std::vector storage_t; - typedef std::vector storage_view_t; + typedef std::vector> storage_view_t; + typedef std::vector> storage_const_view_t; /// The data access - T & operator()(storage_t & data, size_t i) { return data[i];} + T & operator()(storage_t & data, size_t i) { return data[i];} T const & operator()(storage_t const & data, size_t i) const { return data[i];} Tv & operator()(storage_view_t & data, size_t i) { return data[i];} Tv const & operator()(storage_view_t const & data, size_t i) const { return data[i];} + Tcv & operator()(storage_const_view_t & data, size_t i) { return data[i];} + Tcv const & operator()(storage_const_view_t const & data, size_t i) const { return data[i];} +/*Tv operator()(storage_view_t & data, size_t i) const { return data[i];} + Tv operator()(storage_view_t const & data, size_t i) const { return data[i];} + Tcv operator()(storage_const_view_t & data, size_t i) const { return data[i];} + Tcv operator()(storage_const_view_t const & data, size_t i) const { return data[i];} +*/ template static void assign_to_scalar (S & data, RHS && rhs) {for (size_t i =0; i static void rebind (storage_view_t & data, RHS && rhs) { data.clear(); for (auto & x : rhs.data()) data.push_back(x);} + template static void rebind(ST& data, RHS&& rhs) { data.clear(); for (auto & x : rhs.data()) data.push_back(x);} }; //---------------------------- lambda ---------------------------------- @@ -127,13 +161,14 @@ namespace triqs { namespace gfs { /// The storage typedef F storage_t; typedef F storage_view_t; + typedef F storage_const_view_t; /// The data access auto operator()(storage_t & data, size_t i) DECL_AND_RETURN( data(i)); auto operator()(storage_t const & data, size_t i) const DECL_AND_RETURN( data(i)); template static void assign_to_scalar (S & data, RHS && rhs) = delete; - template static void rebind (storage_view_t & data, RHS && rhs) = delete;// { data = std::forward(rhs);} + template static void rebind(ST& data, RHS&& rhs) = delete; }; diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp index cf04d737..3b45c101 100644 --- a/triqs/gfs/gf.hpp +++ b/triqs/gfs/gf.hpp @@ -33,12 +33,23 @@ namespace triqs { namespace gfs { using utility::factory; using arrays::make_shape; - template long get_shape(std::vector const & x) {return x.size();} + template long get_shape(std::vector const &x) { return x.size(); } - template struct gf_mesh; - template class gf; // the regular type - template class gf_view; // the view type - template class gf_impl; + // the gf mesh + template struct gf_mesh; + + // The regular type + template class gf; + + // The view type + template class gf_view; + + // The const view type + template + using gf_const_view = gf_view; + + // the implementation class + template class gf_impl; // various implementation traits namespace gfs_implementation { // never use using of this... @@ -63,23 +74,24 @@ namespace triqs { namespace gfs { template struct h5_name; // value is a const char * template struct h5_name { static std::string invoke(){ return h5_name::invoke() + "_s";}}; - - // the h5 write and read of gf members, so that we can specialize it e.g. for block gf - template struct h5_rw { - static void write (h5::group gr, gf_impl const & g) { + // the h5 write and read of gf members, so that we can specialize it e.g. for block gf + template struct h5_rw { + + static void write (h5::group gr, gf_const_view g) { h5_write(gr,"data",g._data); h5_write(gr,"singularity",g._singularity); h5_write(gr,"mesh",g._mesh); h5_write(gr,"symmetry",g._symmetry); } - static void read (h5::group gr, gf_impl & g) { - h5_read(gr,"data",g._data); - h5_read(gr,"singularity",g._singularity); - h5_read(gr,"mesh",g._mesh); - h5_read(gr,"symmetry",g._symmetry); - } + template + static void read (h5::group gr, gf_impl & g) { + h5_read(gr,"data",g._data); + h5_read(gr,"singularity",g._singularity); + h5_read(gr,"mesh",g._mesh); + h5_read(gr,"symmetry",g._symmetry); + } }; // factories regroup all factories (constructors..) for all types of gf. @@ -112,7 +124,7 @@ namespace triqs { namespace gfs { static typename gf_t::singularity_t make_singularity(mesh_t const & m, target_shape_t shape) { return typename gf_t::singularity_t(shape); } }; - + template struct factories { typedef gf gf_t; struct target_shape_t {}; @@ -149,10 +161,13 @@ namespace triqs { namespace gfs { // ---------------------- implementation -------------------------------- /// A common implementation class for gf and gf_view. They will only redefine contructor and = ... - template class gf_impl : + template class gf_impl : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction){ + static_assert(!( !IsView && IsConst), "Internal error"); public : - typedef gf_view view_type; + typedef gf_view mutable_view_type; + typedef gf_const_view const_view_type; + typedef typename std::conditional ::type view_type; typedef gf regular_type; typedef Variable variable_t; @@ -169,7 +184,9 @@ namespace triqs { namespace gfs { typedef gfs_implementation::data_proxy data_proxy_t; typedef typename data_proxy_t::storage_t data_regular_t; typedef typename data_proxy_t::storage_view_t data_view_t; - typedef typename std::conditional::type data_t; + typedef typename data_proxy_t::storage_const_view_t data_const_view_t; + typedef typename std::conditional::type, + data_regular_t>::type data_t; typedef typename gfs_implementation::singularity::type singularity_non_view_t; typedef typename view_type_if_exists_else_type::type singularity_view_t; @@ -207,8 +224,9 @@ namespace triqs { namespace gfs { protected: - gf_impl(gf_impl const & x): _mesh(x.mesh()), _data(factory(x.data())), - _singularity(factory(x.singularity())), _symmetry(x.symmetry()), _evaluator(x.get_evaluator()){} + template + gf_impl(G && x, bool): _mesh(x.mesh()), _data(factory(x.data())), + _singularity(factory(x.singularity())), _symmetry(x.symmetry()), _evaluator(x.get_evaluator()){} template gf_impl(M && m, D && dat, S && sing, SY && sy, EV && ev) : @@ -225,7 +243,8 @@ namespace triqs { namespace gfs { // ------------- All the call operators ----------------------------- // First, a simple () returns a view, like for an array... - view_type operator()() const { return *this;} + const_view_type operator()() const { return *this;} + view_type operator()() { return *this; } /// Calls are (perfectly) forwarded to the evaluator::operator(), except mesh_point_t and when /// there is at least one lazy argument ... @@ -241,7 +260,6 @@ namespace triqs { namespace gfs { >::type // end of add_Const operator() (Arg0&& arg0, Args&&... args) const { return _evaluator(this,std::forward( arg0), std::forward(args)...); } - // WHY SEPARATE ARg0 ? template typename clef::_result_of::make_expr_call::type operator()(Arg0 &&arg0, Args&&... args) & { @@ -260,20 +278,6 @@ namespace triqs { namespace gfs { return clef::make_expr_call(std::move(*this),std::forward(arg0), std::forward(args)...); } - /* - // on mesh component for composite meshes - // enable iif the first arg is a mesh_point_t for the first component of the mesh_t - template::value > - typename std::enable_if< MeshIsComposite && std::is_base_of< tag::mesh_point, Arg0>::value, r_type>::type - operator() (Arg0 const & arg0, Args const & ... args) - { return _data_proxy(_data, _mesh.mesh_pt_components_to_linear(arg0, args...));} - - template::value > - typename std::enable_if< MeshIsComposite && std::is_base_of< tag::mesh_point, Arg0>::value, cr_type>::type - operator() (Arg0 const & arg0, Args const & ... args) const - { return _data_proxy(_data, _mesh.mesh_pt_components_to_linear(arg0, args...));} - */ - //// [] and access to the grid point typedef typename std::result_of::type r_type; typedef typename std::result_of::type cr_type; @@ -325,13 +329,13 @@ namespace triqs { namespace gfs { friend std::string get_triqs_hdf5_data_scheme(gf_impl const & g) { return "Gf" + gfs_implementation::h5_name::invoke();} - friend class gfs_implementation::h5_rw; + friend class gfs_implementation::h5_rw; /// Write into HDF5 friend void h5_write (h5::group fg, std::string subgroup_name, gf_impl const & g) { auto gr = fg.create_group(subgroup_name); gr.write_triqs_hdf5_data_scheme(g); - gfs_implementation::h5_rw::write(gr, g); + gfs_implementation::h5_rw::write(gr, g); } /// Read from HDF5 @@ -342,7 +346,7 @@ namespace triqs { namespace gfs { auto tag_expected= get_triqs_hdf5_data_scheme(g); if (tag_file != tag_expected) TRIQS_RUNTIME_ERROR<< "h5_read : mismatch of the tag TRIQS_HDF5_data_scheme tag in the h5 group : found "<::read(gr, g); + gfs_implementation::h5_rw::read(gr, g); } //----------------------------- BOOST Serialization ----------------------------- @@ -360,7 +364,7 @@ namespace triqs { namespace gfs { friend std::ostream & triqs_nvl_formal_print(std::ostream & out, gf_impl const & x) { return out<<(IsView ? "gf_view": "gf");} // Interaction with the CLEF library : auto assignment of the gf (gf(om_) << expression fills the functions by evaluation of expression) - template friend void triqs_clef_auto_assign (gf_impl & g, RHS const & rhs) { + template friend void triqs_clef_auto_assign(gf_impl &g, RHS const &rhs) { // access to the data . Beware, we view it as a *matrix* NOT an array... (crucial for assignment to scalars !) g.triqs_clef_auto_assign_impl(rhs, typename std::is_base_of::type()); assign_from_expression(g.singularity(),rhs); @@ -376,12 +380,9 @@ namespace triqs { namespace gfs { for (auto const & w: this->mesh()) (*this)[w] = rhs(w); //for (auto const & w: this->mesh()) (*this)[w] = rhs(typename B::mesh_t::mesh_point_t::cast_t(w)); } - template void triqs_clef_auto_assign_impl (RHS const & rhs, std::integral_constant) { + template void triqs_clef_auto_assign_impl(RHS const &rhs, std::integral_constant) { for (auto const & w: this->mesh()) { - //std::cout << "abot "<mesh()) triqs::tuple::apply(*this,w.components_tuple()) = triqs::tuple::apply(rhs,w.components_tuple()); } @@ -389,15 +390,15 @@ namespace triqs { namespace gfs { // --------------------------------------------------------------------------------- ///The regular class of GF - template class gf : public gf_impl { - typedef gf_impl B; + template class gf : public gf_impl { + typedef gf_impl B; typedef gfs_implementation::factories factory; public : gf():B() {} gf(gf const & g): B(g){} gf(gf && g) noexcept : B(std::move(g)){} - gf(gf_view const & g): B(g){} + gf(gf_view const & g): B(g, bool() ){} template gf(GfType const & x,typename std::enable_if::value>::type *dummy =0 ): B() { *this = x;} gf(typename B::mesh_t m, @@ -431,14 +432,50 @@ namespace triqs { namespace gfs { }; // --------------------------------------------------------------------------------- - ///The View class of GF - template class gf_view : public gf_impl { - typedef gf_impl B; + ///The const View class of GF + template class gf_view : public gf_impl { + typedef gf_impl B; public : + gf_view() = delete; gf_view(gf_view const & g): B(g){} gf_view(gf_view && g) noexcept : B(std::move(g)){} - template gf_view(gf_impl const & g): B(g){} + gf_view(gf_impl const & g) : B(g, bool()){} // from a const_view + gf_view(gf_impl const & g): B(g, bool()){} // from a view + gf_view(gf_impl const & g): B(g, bool()){} // from a const gf + gf_view(gf_impl & g): B(g, bool()){} // from a gf & + gf_view(gf_impl && g) noexcept: B(std::move(g), bool()){} // from a gf && + + template + gf_view (typename B::mesh_t const & m, + D const & dat,typename B::singularity_view_t const & t,typename B::symmetry_t const & s, + typename B::evaluator_t const &e = typename B::evaluator_t () ) : + B(m,factory(dat),t,s,e) {} + + void rebind (gf_view const &X) noexcept { + this->_mesh = X._mesh; this->_symmetry = X._symmetry; + this->_data_proxy.rebind(this->_data,X); + this->_singularity.rebind(X._singularity); + } + gf_view & operator = (gf_view const & ) = delete; + + }; // class gf_const_view + + // --------------------------------------------------------------------------------- + ///The View class of GF + template class gf_view : public gf_impl { + typedef gf_impl B; + public : + gf_view() = delete; + gf_view(gf_view const & g): B(g){} + gf_view(gf_view && g) noexcept : B(std::move(g)){} + + gf_view(gf_impl const & g) = delete; // from a const view : impossible + gf_view(gf_impl const & g): B(g, bool()){} // from a view + gf_view(gf_impl const & g) = delete; // from a const gf : impossible + //gf_view(gf_impl const & g): B(g, bool()){} // from a gf & + gf_view(gf_impl & g): B(g, bool()){} // from a gf & + gf_view(gf_impl && g) noexcept: B(std::move(g), bool()){} // from a gf && template gf_view (typename B::mesh_t const & m, @@ -474,33 +511,67 @@ namespace triqs { namespace gfs { } // tool for lazy transformation - template struct gf_keeper{ gf_view g; gf_keeper (gf_view const & g_) : g(g_) {} }; + template struct gf_keeper { + gf_const_view g; + }; // ---------------------------------- slicing ------------------------------------ - //slice - template - gf_view slice_target (gf_view g, Args&& ... args) { - static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); - using arrays::range; - //auto sg=slice_target (g.singularity(),range(args,args+1)...); - return gf_view(g.mesh(), g.data()(range(), std::forward(args)... ), slice_target (g.singularity(), std::forward(args)...) , g.symmetry()); - } + // slice + template + gf_view slice_target(gf_view g, Args &&... args) { + static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); + using arrays::range; + return {g.mesh(), g.data()(range(), std::forward(args)...), + slice_target(g.singularity(), std::forward(args)...), g.symmetry()}; + } - template - gf_view slice_target_to_scalar (gf_view g, Args&& ... args) { - static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); - using arrays::range; - auto sg=slice_target (g.singularity(),range(args,args+1)...); - return gf_view(g.mesh(), g.data()(range(), std::forward(args)... ), sg, g.symmetry()); - } + template + gf_view slice_target(gf & g, Args &&... args) { + return slice_target(g(), std::forward(args)...); + } + + template + gf_const_view slice_target(gf const &g, Args &&... args) { + return slice_target(g(), std::forward(args)...); + } + + // slice to scalar + template + gf_view slice_target_to_scalar(gf_view g, + Args &&... args) { + static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); + using arrays::range; + return {g.mesh(), g.data()(range(), std::forward(args)...), + slice_target(g.singularity(), range(args, args + 1)...), g.symmetry()}; + } + + template + gf_view slice_target_to_scalar(gf & g, Args &&... args) { + return slice_target_to_scalar(g(), std::forward(args)...); + } + + template + gf_const_view slice_target_to_scalar(gf const &g, Args &&... args) { + return slice_target_to_scalar(g(), std::forward(args)...); + } // a scalar_valued gf can be viewed as a 1x1 matrix - template - gf_view reinterpret_scalar_valued_gf_as_matrix_valued (gf_view g) { - typedef typename gf_view::data_view_t a_t; + template + gf_view reinterpret_scalar_valued_gf_as_matrix_valued (gf_view g) { + typedef typename gf_view::data_view_t a_t; auto a = a_t {typename a_t::indexmap_type (join(g.data().shape(),make_shape(1,1))), g.data().storage()}; - return gf_view(g.mesh(), a, g.singularity(), g.symmetry()); + return {g.mesh(), a, g.singularity(), g.symmetry()}; + } + + template + gf_view reinterpret_scalar_valued_gf_as_matrix_valued (gf & g) { + return reinterpret_scalar_valued_gf_as_matrix_valued(g()); + } + + template + gf_const_view reinterpret_scalar_valued_gf_as_matrix_valued (gf const & g) { + return reinterpret_scalar_valued_gf_as_matrix_valued(g()); } /* diff --git a/triqs/gfs/local/fourier_matsubara.cpp b/triqs/gfs/local/fourier_matsubara.cpp index 92c0a043..78bc574c 100644 --- a/triqs/gfs/local/fourier_matsubara.cpp +++ b/triqs/gfs/local/fourier_matsubara.cpp @@ -44,7 +44,7 @@ namespace triqs { namespace gfs { tqa::vector g_in, g_out; - void direct (gf_view gw, gf_view const gt) { + void direct (gf_view gw, gf_const_view gt) { using namespace impl_local_matsubara; auto ta = gt(freq_infty()); //TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO @@ -81,7 +81,7 @@ namespace triqs { namespace gfs { gw.singularity() = gt.singularity();// set tail } - void inverse(gf_view gt, gf_view const gw){ + void inverse(gf_view gt, gf_const_view gw){ using namespace impl_local_matsubara; static bool Green_Function_Are_Complex_in_time = false; // If the Green function are NOT complex, then one use the symmetry property @@ -139,12 +139,12 @@ namespace triqs { namespace gfs { }; - void fourier_impl (gf_view gw , gf_view const gt, scalar_valued){ + void fourier_impl (gf_view gw , gf_const_view gt, scalar_valued){ impl_worker w; w.direct(gw,gt); } - void fourier_impl (gf_view gw , gf_view const gt, matrix_valued){ + void fourier_impl (gf_view gw , gf_const_view gt, matrix_valued){ impl_worker w; for (size_t n1=0; n1 gt , gf_view const gw, scalar_valued){ + void inverse_fourier_impl (gf_view gt , gf_const_view gw, scalar_valued){ impl_worker w; w.inverse(gt,gw); } - void inverse_fourier_impl (gf_view gt , gf_view const gw, matrix_valued){ + void inverse_fourier_impl (gf_view gt , gf_const_view gw, matrix_valued){ impl_worker w; for (size_t n1=0; n1 gw , gf_view const gt, scalar_valued); - void fourier_impl (gf_view gw , gf_view const gt, matrix_valued); - void inverse_fourier_impl (gf_view gt, gf_view const gw, scalar_valued); - void inverse_fourier_impl (gf_view gt, gf_view const gw, matrix_valued); + void fourier_impl (gf_view gw , gf_const_view gt, scalar_valued); + void fourier_impl (gf_view gw , gf_const_view gt, matrix_valued); + void inverse_fourier_impl (gf_view gt, gf_const_view gw, scalar_valued); + void inverse_fourier_impl (gf_view gt, gf_const_view gw, matrix_valued); - inline gf_view fourier (gf_view const gt) { + inline gf_view fourier (gf_const_view gt) { int L = (gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size() ); auto gw = gf{ {gt.domain(),L}, gt.data().shape().front_pop() }; auto V = gw(); fourier_impl(V, gt, matrix_valued()); return gw; } - inline gf_view fourier (gf_view const gt) { + inline gf_view fourier (gf_const_view gt) { int L = (gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size() ); auto gw = gf{ {gt.domain(),L} }; auto V = gw(); @@ -48,7 +48,7 @@ namespace triqs { namespace gfs { return gw; } - inline gf_view inverse_fourier (gf_view const gw, mesh_kind mk = half_bins) { + inline gf_view inverse_fourier (gf_const_view gw, mesh_kind mk = half_bins) { double pi = std::acos(-1); int L = (mk == full_bins ? gw.mesh().size()+1 : gw.mesh().size() ); auto gt = gf{ {gw.domain(),L}, gw.data().shape().front_pop()}; @@ -56,7 +56,7 @@ namespace triqs { namespace gfs { inverse_fourier_impl(V, gw, matrix_valued()); return gt; } - inline gf_view inverse_fourier (gf_view const gw, mesh_kind mk = half_bins) { + inline gf_view inverse_fourier (gf_const_view gw, mesh_kind mk = half_bins) { double pi = std::acos(-1); int L = (mk == full_bins ? gw.mesh().size()+1 : gw.mesh().size() ); auto gt = gf{ {gw.domain(),L} }; @@ -65,10 +65,10 @@ namespace triqs { namespace gfs { return gt; } - inline gf_keeper lazy_fourier (gf_view const & g) { return g;} - inline gf_keeper lazy_inverse_fourier (gf_view const & g) { return g;} - inline gf_keeper lazy_fourier (gf_view const & g) { return g;} - inline gf_keeper lazy_inverse_fourier (gf_view const & g) { return g;} + inline gf_keeper lazy_fourier (gf_const_view g) { return {g};} + inline gf_keeper lazy_inverse_fourier (gf_const_view g) { return {g};} + inline gf_keeper lazy_fourier (gf_const_view g) { return {g};} + inline gf_keeper lazy_inverse_fourier (gf_const_view g) { return {g};} void triqs_gf_view_assign_delegation( gf_view g, gf_keeper const & L); void triqs_gf_view_assign_delegation( gf_view g, gf_keeper const & L); diff --git a/triqs/gfs/local/fourier_real.cpp b/triqs/gfs/local/fourier_real.cpp index 504c43dc..6ebe06c7 100644 --- a/triqs/gfs/local/fourier_real.cpp +++ b/triqs/gfs/local/fourier_real.cpp @@ -37,7 +37,7 @@ namespace triqs { namespace gfs { tqa::vector g_in, g_out; - void direct (gf_view gw, gf_view const gt){ + void direct (gf_view gw, gf_const_view gt){ size_t L = gt.mesh().size(); if (gw.mesh().size() != L) TRIQS_RUNTIME_ERROR << "Meshes are different"; @@ -70,7 +70,7 @@ namespace triqs { namespace gfs { } - void inverse(gf_view gt, gf_view const gw){ + void inverse(gf_view gt, gf_const_view gw){ size_t L = gw.mesh().size(); if ( L != gt.mesh().size()) TRIQS_RUNTIME_ERROR << "Meshes are different"; @@ -108,12 +108,12 @@ namespace triqs { namespace gfs { //-------------------------------------------------------------------------------------- - void fourier_impl(gf_view gw, gf_view const gt, scalar_valued){ + void fourier_impl(gf_view gw, gf_const_view gt, scalar_valued){ impl_worker w; w.direct(gw, gt); } - void fourier_impl(gf_view gw, gf_view const gt, matrix_valued){ + void fourier_impl(gf_view gw, gf_const_view gt, matrix_valued){ impl_worker w; for (size_t n1=0; n1 gt, gf_view const gw, scalar_valued){ + void inverse_fourier_impl (gf_view gt, gf_const_view gw, scalar_valued){ impl_worker w; w.inverse(gt,gw); } - void inverse_fourier_impl (gf_view gt, gf_view const gw, matrix_valued){ + void inverse_fourier_impl (gf_view gt, gf_const_view gw, matrix_valued){ impl_worker w; for (size_t n1=0; n1 gw , gf_view const gt, scalar_valued); - void fourier_impl (gf_view gw , gf_view const gt, matrix_valued); - void inverse_fourier_impl (gf_view gt, gf_view const gw, scalar_valued); - void inverse_fourier_impl (gf_view gt, gf_view const gw, matrix_valued); + void fourier_impl (gf_view gw , gf_const_view gt, scalar_valued); + void fourier_impl (gf_view gw , gf_const_view gt, matrix_valued); + void inverse_fourier_impl (gf_view gt, gf_const_view gw, scalar_valued); + void inverse_fourier_impl (gf_view gt, gf_const_view gw, matrix_valued); - inline gf_view fourier (gf_view const gt) { + inline gf_view fourier (gf_const_view gt) { double pi = std::acos(-1); int L = gt.mesh().size(); double wmin = -pi * (L-1) / (L*gt.mesh().delta()); @@ -43,7 +43,7 @@ namespace triqs { namespace gfs { fourier_impl(V, gt, matrix_valued()); return gw; } - inline gf_view fourier (gf_view const gt) { + inline gf_view fourier (gf_const_view gt) { double pi = std::acos(-1); int L = gt.mesh().size(); double wmin = -pi * (L-1) / (L*gt.mesh().delta()); @@ -54,7 +54,7 @@ namespace triqs { namespace gfs { return gw; } - inline gf_view inverse_fourier (gf_view const gw) { + inline gf_view inverse_fourier (gf_const_view gw) { double pi = std::acos(-1); int L = gw.mesh().size(); double tmin = -pi * (L-1) / (L*gw.mesh().delta()); @@ -64,7 +64,7 @@ namespace triqs { namespace gfs { inverse_fourier_impl(V, gw, matrix_valued()); return gt; } - inline gf_view inverse_fourier (gf_view const gw) { + inline gf_view inverse_fourier (gf_const_view gw) { double pi = std::acos(-1); int L = gw.mesh().size(); double tmin = -pi * (L-1) / (L*gw.mesh().delta()); @@ -74,11 +74,11 @@ namespace triqs { namespace gfs { inverse_fourier_impl(V, gw, scalar_valued()); return gt; } - - inline gf_keeper lazy_fourier (gf_view const & g) { return g;} - inline gf_keeper lazy_inverse_fourier (gf_view const & g) { return g;} - inline gf_keeper lazy_fourier (gf_view const & g) { return g;} - inline gf_keeper lazy_inverse_fourier (gf_view const & g) { return g;} + + inline gf_keeper lazy_fourier(gf_const_view g) { return {g}; } + inline gf_keeper lazy_inverse_fourier(gf_const_view g) { return {g}; } + inline gf_keeper lazy_fourier(gf_const_view g) { return {g}; } + inline gf_keeper lazy_inverse_fourier(gf_const_view g) { return {g}; } void triqs_gf_view_assign_delegation( gf_view g, gf_keeper const & L); void triqs_gf_view_assign_delegation( gf_view g, gf_keeper const & L); diff --git a/triqs/gfs/local/functions.cpp b/triqs/gfs/local/functions.cpp index f72982b6..c1279af7 100644 --- a/triqs/gfs/local/functions.cpp +++ b/triqs/gfs/local/functions.cpp @@ -87,7 +87,7 @@ namespace triqs { namespace gfs { // compute a tail from the Legendre GF // this is Eq. 8 of our paper - local::tail_view get_tail(gf_view const & gl, int size = 10, int omin = -1) { + local::tail_view get_tail(gf_const_view gl, int size = 10, int omin = -1) { auto sh = gl.data().shape().front_pop(); local::tail t(sh, size, omin); diff --git a/triqs/gfs/local/functions.hpp b/triqs/gfs/local/functions.hpp index b768b5d0..284660c2 100644 --- a/triqs/gfs/local/functions.hpp +++ b/triqs/gfs/local/functions.hpp @@ -36,7 +36,7 @@ namespace triqs { tqa::matrix density(gf_view const & g); - local::tail_view get_tail(gf_view const & gl, int size, int omin); + local::tail_view get_tail(gf_const_view gl, int size, int omin); void enforce_discontinuity(gf_view & gl, triqs::arrays::array_view disc); diff --git a/triqs/gfs/local/legendre_matsubara.cpp b/triqs/gfs/local/legendre_matsubara.cpp index 640d270c..5fd95c1c 100644 --- a/triqs/gfs/local/legendre_matsubara.cpp +++ b/triqs/gfs/local/legendre_matsubara.cpp @@ -29,7 +29,7 @@ using namespace triqs::utility; namespace triqs { namespace gfs { -void legendre_matsubara_direct(gf_view & gw, gf_view const & gl) { +void legendre_matsubara_direct(gf_view & gw, gf_const_view gl) { gw() = 0.0; triqs::arrays::range R; @@ -45,7 +45,7 @@ void legendre_matsubara_direct(gf_view & gw, gf_view const & g } -void legendre_matsubara_inverse (gf_view & gl, gf_view const & gw) { +void legendre_matsubara_inverse (gf_view & gl, gf_const_view gw) { gl() = 0.0; @@ -62,7 +62,7 @@ void legendre_matsubara_inverse (gf_view & gl, gf_view const & } -void legendre_matsubara_direct (gf_view & gt, gf_view const & gl) { +void legendre_matsubara_direct (gf_view & gt, gf_const_view gl) { gt() = 0.0; legendre_generator L; @@ -78,7 +78,7 @@ void legendre_matsubara_direct (gf_view & gt, gf_view const & } -void legendre_matsubara_inverse (gf_view & gl, gf_view const & gt) { +void legendre_matsubara_inverse (gf_view & gl, gf_const_view gt) { gl() = 0.0; legendre_generator L; @@ -95,10 +95,10 @@ void legendre_matsubara_inverse (gf_view & gl, gf_view const & } -gf_keeper lazy_legendre_imfreq (gf_view const & gl) { return gl; } -gf_keeper lazy_legendre_imtime (gf_view const & gl) { return gl; } -gf_keeper lazy_imfreq_legendre (gf_view const & gw) { return gw; } -gf_keeper lazy_imtime_legendre (gf_view const & gt) { return gt; } +gf_keeper lazy_legendre_imfreq (gf_const_view gl) { return {gl}; } +gf_keeper lazy_legendre_imtime (gf_const_view gl) { return {gl}; } +gf_keeper lazy_imfreq_legendre (gf_const_view gw) { return {gw}; } +gf_keeper lazy_imtime_legendre (gf_const_view gt) { return {gt}; } void triqs_gf_view_assign_delegation( gf_view &gw, gf_keeper const & L) { legendre_matsubara_direct(gw, L.g); diff --git a/triqs/gfs/local/legendre_matsubara.hpp b/triqs/gfs/local/legendre_matsubara.hpp index 18b6fe3a..83a6de5f 100644 --- a/triqs/gfs/local/legendre_matsubara.hpp +++ b/triqs/gfs/local/legendre_matsubara.hpp @@ -28,18 +28,18 @@ namespace triqs { namespace gfs { - void legendre_matsubara_direct (gf_view &gw, gf_view const &gl); - void legendre_matsubara_inverse (gf_view &gl, gf_view const &gw); + void legendre_matsubara_direct (gf_view &gw, gf_const_view gl); + void legendre_matsubara_inverse (gf_view &gl, gf_const_view gw); - void legendre_matsubara_direct (gf_view >, gf_view const &gl); - void legendre_matsubara_inverse (gf_view &gl, gf_view const >); + void legendre_matsubara_direct (gf_view >, gf_const_view gl); + void legendre_matsubara_inverse (gf_view &gl, gf_const_view gt); namespace tags { struct legendre{}; } - gf_keeper lazy_legendre_imfreq (gf_view const & gl); - gf_keeper lazy_legendre_imtime (gf_view const & gl); - gf_keeper lazy_imfreq_legendre (gf_view const & gw); - gf_keeper lazy_imtime_legendre (gf_view const & gt); + gf_keeper lazy_legendre_imfreq (gf_const_view gl); + gf_keeper lazy_legendre_imtime (gf_const_view gl); + gf_keeper lazy_imfreq_legendre (gf_const_view gw); + gf_keeper lazy_imtime_legendre (gf_const_view gt); void triqs_gf_view_assign_delegation( gf_view &gw, gf_keeper const & L); void triqs_gf_view_assign_delegation( gf_view >, gf_keeper const & L); diff --git a/triqs/gfs/meshes/product.hpp b/triqs/gfs/meshes/product.hpp index f47ce443..c1703abe 100644 --- a/triqs/gfs/meshes/product.hpp +++ b/triqs/gfs/meshes/product.hpp @@ -176,9 +176,9 @@ namespace triqs { namespace gfs { // reinterpret_linear_array(m,A) returns a d-dimensionnal view of the array // with indices egal to the indices of the components of the mesh. // Very useful for slicing, currying functions. - template - arrays::array_view - reinterpret_linear_array(mesh_product const & m, arrays::array_view const & A) { + template + arrays::array_view + reinterpret_linear_array(mesh_product const & m, arrays::array_view A) { return { {join (m.all_size_as_mini_vector(), get_shape(A).front_pop())}, A.storage()}; } diff --git a/triqs/gfs/product.hpp b/triqs/gfs/product.hpp index a4a8175e..21f72f4d 100644 --- a/triqs/gfs/product.hpp +++ b/triqs/gfs/product.hpp @@ -59,29 +59,32 @@ namespace triqs { namespace gfs { template struct h5_name,tensor_valued,Opt> : h5_name,matrix_valued,Opt> {}; // a slight difference with the generic case : reinterpret the data array to avoid flattening the variables - template - struct h5_rw,tensor_valued,Opt,IsView> { - typedef gf_impl,tensor_valued,Opt,IsView> g_t; + template + struct h5_rw,tensor_valued,Opt> { + typedef gf,tensor_valued,Opt> g_t; - static void write (h5::group gr, g_t const & g) { + static void write (h5::group gr, typename g_t::const_view_type g) { h5_write(gr,"data",reinterpret_linear_array(g.mesh(), g().data())); h5_write(gr,"singularity",g._singularity); h5_write(gr,"mesh",g._mesh); h5_write(gr,"symmetry",g._symmetry); } - static void read (h5::group gr, g_t & g) { + template + static void read (h5::group gr, gf_impl,tensor_valued,Opt,IsView,false> & g) { + using G_t= gf_impl,tensor_valued,Opt,IsView,false> ; h5_read(gr,"mesh",g._mesh); - auto arr = arrays::array{}; + auto arr = arrays::array{}; h5_read(gr,"data",arr); auto sh = arr.shape(); arrays::mini_vector sh2; sh2[0] = g._mesh.size(); for (int u=1; u{sh2, std::move(arr.storage())}; + g._data = arrays::array{sh2, std::move(arr.storage())}; h5_read(gr,"singularity",g._singularity); h5_read(gr,"symmetry",g._symmetry); } + }; /// --------------------------- data access ---------------------------------