/******************************************************************************* * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * * Copyright (C) 2012-2013 by M. Ferrero, 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 . * ******************************************************************************/ #ifndef TRIQS_GF_GFBASECLASS_H #define TRIQS_GF_GFBASECLASS_H #include #include #include #include #include #include #include "./tools.hpp" #include "./data_proxies.hpp" namespace triqs { namespace gfs { using utility::factory; using arrays::make_shape; template long get_shape(std::vector const &x) { return x.size(); } // 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... // evaluator regroup functions to evaluate the function. template struct evaluator{ static constexpr int arity = 0;}; // closest_point mechanism template struct get_closest_point; // singularity template struct singularity { typedef nothing type;}; // symmetry template struct symmetry { typedef nothing type;}; // data_proxy contains function to manipulate the data array, but no data itself. // this is used to specialize this part of the code to array of dim 3 (matrix gf), dim 1 (scalar gf) and vector (e.g. block gf, ...) template struct data_proxy; // This trait contains functions to read/write in hdf5 files. Can be specialized for some case (Cf block) 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_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); } 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. template struct factories; template struct factories,Opt> { typedef gf,Opt> gf_t; typedef tqa::mini_vector target_shape_t; typedef typename gf_t::mesh_t mesh_t; static typename gf_t::data_t make_data(mesh_t const & m, target_shape_t shape) { typename gf_t::data_t A(shape.front_append(m.size())); A() =0; return A; } 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; typedef tqa::mini_vector target_shape_t; typedef typename gf_t::mesh_t mesh_t; static typename gf_t::data_t make_data(mesh_t const & m, target_shape_t shape) { typename gf_t::data_t A(shape.front_append(m.size())); A() =0; return A; } 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 {}; typedef typename gf_t::mesh_t mesh_t; static typename gf_t::data_t make_data(mesh_t const & m, target_shape_t shape) { typename gf_t::data_t A(m.size()); A() =0; return A; } static typename gf_t::singularity_t make_singularity(mesh_t const & m, target_shape_t shape) { return typename gf_t::singularity_t {1,1} ; } }; } // gfs_implementation // make_gf and make_gf_view forward any args to them template gf make_gf(gf_mesh m, U && ... x) { return gfs_implementation::factories::make_gf(std::move(m),std::forward(x)...);} template gf make_gf(U && ... x) { return gfs_implementation::factories::make_gf(std::forward(x)...);} template gf_view make_gf_view(U && ... x) { return gfs_implementation::factories::make_gf_view(std::forward(x)...);} // The trait that "marks" the Green function TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT(ImmutableGreenFunction); template auto get_gf_data_shape(G const & g) DECL_AND_RETURN(g.get_data_shape()); // ---------------------- implementation -------------------------------- /// A common implementation class for gf and gf_view. They will only redefine contructor and = ... template class gf_impl : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction){ static_assert(!( !IsView && IsConst), "Internal error"); public : 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; typedef Target target_t; typedef Opt option_t; typedef gf_mesh mesh_t; typedef typename mesh_t::domain_t domain_t; typedef typename mesh_t::mesh_point_t mesh_point_t; typedef typename mesh_t::index_t mesh_index_t; typedef typename gfs_implementation::symmetry::type symmetry_t; typedef gfs_implementation::evaluator evaluator_t; 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 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; typedef typename std::conditional::type singularity_t; mesh_t const & mesh() const { return _mesh;} domain_t const & domain() const { return _mesh.domain();} data_t & data() { return _data;} data_t const & data() const { return _data;} singularity_t & singularity() { return _singularity;} singularity_t const & singularity() const { return _singularity;} symmetry_t const & symmetry() const { return _symmetry;} evaluator_t const & get_evaluator() const { return _evaluator;} auto get_data_shape() const DECL_AND_RETURN(get_shape(this->data())); protected: mesh_t _mesh; data_t _data; singularity_t _singularity; symmetry_t _symmetry; evaluator_t _evaluator; data_proxy_t _data_proxy; // --------------------------------Constructors ----------------------------------------------- // all protected but one, this is an implementation class, see gf/gf_view later for public one gf_impl() {} // all arrays of zero size (empty) public : //everyone can make a copy (for clef lib in particular, this one needs to be public) gf_impl(gf_impl const & x) : _mesh(x.mesh()), _data(factory(x.data())), _singularity(factory(x.singularity())), _symmetry(x.symmetry()), _evaluator(x._evaluator){} gf_impl(gf_impl &&) = default; protected: 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) : _mesh(std::forward(m)), _data(std::forward(dat)), _singularity(std::forward(sing)),_symmetry(std::forward(sy)), _evaluator(std::forward(ev)){} void operator = (gf_impl const & rhs) = delete; // done in derived class. void swap_impl (gf_impl & b) noexcept { using std::swap; swap(this->_mesh, b._mesh); swap(this->_data, b._data); swap (this->_singularity,b._singularity); swap(this->_symmetry,b._symmetry); swap(this->_evaluator,b._evaluator); } public: // ------------- All the call operators ----------------------------- // First, a simple () returns a view, like for an array... 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 ... template // match any argument list, picking out the first type : () is not permitted typename std::add_const< typename boost::lazy_disable_if< // disable the template if one the following conditions it true boost::mpl::or_< // starting condition [OR] clef::is_any_lazy // One of Args is a lazy expression , boost::mpl::bool_<(sizeof...(Args)!= evaluator_t::arity -1 ) && (evaluator_t::arity !=-1)> // if -1 : no check >, // end of OR std::result_of // what is the result type of call >::type // end of lazy_disable_if >::type // end of add_Const operator() (Arg0&& arg0, Args&&... args) const { return _evaluator(this,std::forward( arg0), std::forward(args)...); } template typename clef::_result_of::make_expr_call::type operator()(Arg0 &&arg0, Args&&... args) & { return clef::make_expr_call(*this,std::forward(arg0), std::forward(args)...); } template typename clef::_result_of::make_expr_call::type operator()(Arg0 &&arg0, Args&&... args) const & { return clef::make_expr_call(*this,std::forward(arg0), std::forward(args)...); } template typename clef::_result_of::make_expr_call::type operator()(Arg0 &&arg0, Args&&... args) && { return clef::make_expr_call(std::move(*this),std::forward(arg0), std::forward(args)...); } //// [] and access to the grid point typedef typename std::result_of::type r_type; typedef typename std::result_of::type cr_type; r_type operator[] (mesh_index_t const & arg) { return _data_proxy(_data,_mesh.index_to_linear(arg));} cr_type operator[] (mesh_index_t const & arg) const { return _data_proxy(_data,_mesh.index_to_linear(arg));} r_type operator[] (mesh_point_t const & x) { return _data_proxy(_data, x.linear_index());} cr_type operator[] (mesh_point_t const & x) const { return _data_proxy(_data, x.linear_index());} template r_type operator[] (closest_pt_wrap const & p) { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point::invoke(this,p)));} template cr_type operator[] (closest_pt_wrap const & p) const { return _data_proxy(_data, _mesh.index_to_linear( gfs_implementation::get_closest_point::invoke(this,p)));} template typename clef::_result_of::make_expr_subscript::type operator[](Arg && arg) const & { return clef::make_expr_subscript(*this,std::forward(arg));} template typename clef::_result_of::make_expr_subscript::type operator[](Arg && arg) & { return clef::make_expr_subscript(*this,std::forward(arg));} template typename clef::_result_of::make_expr_subscript::type operator[](Arg && arg) && { return clef::make_expr_subscript(std::move(*this),std::forward(arg));} /// A direct access to the grid point template r_type on_mesh (Args&&... args) { return _data_proxy(_data,_mesh.index_to_linear(mesh_index_t(std::forward(args)...)));} template cr_type on_mesh (Args&&... args) const { return _data_proxy(_data,_mesh.index_to_linear(mesh_index_t(std::forward(args)...)));} private: struct _on_mesh_wrapper_const { gf_impl const & f; _on_mesh_wrapper_const (gf_impl const & _f) : f(_f) {} template cr_type operator ()(Args && ... args) const { return f.on_mesh(std::forward(args)...);} }; struct _on_mesh_wrapper { gf_impl & f; _on_mesh_wrapper (gf_impl & _f) : f(_f) {} template r_type operator ()(Args && ... args) const { return f.on_mesh(std::forward(args)...);} }; _on_mesh_wrapper_const friend on_mesh(gf_impl const & f) { return f;} _on_mesh_wrapper friend on_mesh(gf_impl & f) { return f;} public: //----------------------------- HDF5 ----------------------------- 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; /// 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); } /// Read from HDF5 friend void h5_read (h5::group fg, std::string subgroup_name, gf_impl & g) { auto gr = fg.open_group(subgroup_name); // Check the attribute or throw auto tag_file = gr.read_triqs_hdf5_data_scheme(); 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); } //----------------------------- BOOST Serialization ----------------------------- friend class boost::serialization::access; template void serialize(Archive & ar, const unsigned int version) { ar & boost::serialization::make_nvp("data",_data); ar & boost::serialization::make_nvp("singularity",_singularity); ar & boost::serialization::make_nvp("mesh",_mesh); ar & boost::serialization::make_nvp("symmetry",_symmetry); } /// print friend std::ostream & operator << (std::ostream & out, gf_impl const & x) { return out<<(IsView ? "gf_view": "gf");} 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) { // 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); // if f is an expression, replace the placeholder with a simple tail. If f is a function callable on freq_infty, // it uses the fact that tail_non_view_t can be casted into freq_infty } // enable the writing g[om_] << .... also template friend void triqs_clef_auto_assign_subscript (gf_impl & g, RHS rhs) { triqs_clef_auto_assign(g,rhs);} private: template void triqs_clef_auto_assign_impl (RHS const & rhs, std::integral_constant) { 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) { for (auto const & w: this->mesh()) { (*this)[w] = triqs::tuple::apply(rhs,w.components_tuple()); } //for (auto w: this->mesh()) triqs::tuple::apply(*this,w.components_tuple()) = triqs::tuple::apply(rhs,w.components_tuple()); } }; // --------------------------------------------------------------------------------- ///The regular class of GF 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, bool() ){} gf(gf_const_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, typename B::data_t dat, typename B::singularity_view_t const & si, typename B::symmetry_t const & s , typename B::evaluator_t const & eval = typename B::evaluator_t () ) : B(std::move(m),std::move(dat), si,s,eval) {} typedef typename factory::target_shape_t target_shape_t; gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{} ): B(std::move(m), factory::make_data(m,shape), factory::make_singularity(m,shape), typename B::symmetry_t {}, typename B::evaluator_t{}) {} friend void swap (gf & a, gf & b) noexcept { a.swap_impl (b);} gf & operator = (gf const & rhs) { *this = gf(rhs); return *this;} // use move = gf & operator = (gf & rhs) { *this = gf(rhs); return *this;} // use move = gf & operator = (gf && rhs) noexcept { swap(*this,rhs); return *this;} template void operator = (RHS && rhs) { this->_mesh = rhs.mesh(); this->_data.resize(get_gf_data_shape(rhs)); for (auto const & w: this->mesh()) { (*this)[w] = rhs[w]; } this->_singularity = rhs.singularity(); // to be implemented : there is none in the gf_expr in particular.... //this->_symmetry = rhs.symmetry(); } }; // --------------------------------------------------------------------------------- ///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)){} 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, 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) {} friend void swap (gf_view & a, gf_view & b) noexcept { a.swap_impl (b);} 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 & rhs) { triqs_gf_view_assign_delegation(*this,rhs); return *this;} template gf_view & operator = (RHS const & rhs) { triqs_gf_view_assign_delegation(*this,rhs); return *this;} }; // class gf_view // delegate = so that I can overload it for specific RHS... template DISABLE_IF(arrays::is_scalar) triqs_gf_view_assign_delegation( gf_view g, RHS const & rhs) { if (!(g.mesh() == rhs.mesh())) TRIQS_RUNTIME_ERROR<<"Gf Assignment in View : incompatible mesh"; for (auto const & w: g.mesh()) g[w] = rhs[w]; g.singularity() = rhs.singularity(); } template ENABLE_IF(arrays::is_scalar) triqs_gf_view_assign_delegation( gf_view g, T const & x) { gf_view::data_proxy_t::assign_to_scalar(g.data(), x); g.singularity() = x; } // tool for lazy transformation 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; return {g.mesh(), g.data()(range(), std::forward(args)...), slice_target(g.singularity(), std::forward(args)...), 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; auto a = a_t {typename a_t::indexmap_type (join(g.data().shape(),make_shape(1,1))), g.data().storage()}; 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()); } /* template gf_view slice_mesh (gf_impl const & g, Args... args) { return gf_view(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), g.singularity(), g.symmetry()); }*/ }} #include "./gf_expr.hpp" #endif