LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-25 05:43:40 +01:00
Code Issues Releases Wiki Activity

[gfs] Singularity as an explicit template parameter

Browse Source 
- change the underlying data : do not flatten the linear indices of the
  mesh into a single index, keep a higher dim array.
- easier for various places, and necessary for g(nu,nu').
- work on several 2 part. containers.
- add default target (may not be always matrix_valued)
This commit is contained in:
Olivier Parcollet 2014-09-05 13:26:21 +02:00
parent 9c129cb224
commit fa07abbea9
38 changed files with 1318 additions and 813 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
test/speed/gf_s.cpp
View File

@ -1,4 +1,7 @@
#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
//#define TRIQS_ARRAYS_CHECK_WEAK_REFS
#include <triqs/gfs.hpp>
using namespace triqs::gfs;
using namespace triqs::arrays;
@ -10,26 +13,20 @@ const int N = 1000;
typedef double VALUE_TYPE;
// typedef int VALUE_TYPE;
inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8.0);}
//inline VALUE_TYPE fnt(size_t i) { return i;} //*(i+2.0)*(i-8);}
// inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8.0);}
inline VALUE_TYPE fnt(size_t i) { return i; } //*(i+2.0)*(i-8);}
// inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8);}
struct with_sliding_view {
struct with_g {
void operator()() {
double beta = 1;
auto G = gf<imfreq>{{beta, Fermion, N}, {2, 2}};
G() = 0;
//auto slv = G.data_getter.slv;
for (int u = 0; u < nl_interne; ++u)
//for (int i =0; i<N-1; ++i) { G.data_getter.slv[i](0,0) = fnt(i);}
//for (int i =0; i<N-1; ++i) { slv.set(i); slv(0,0) = fnt(i);}
//for (int i =0; i<N-1; ++i) G.data_getter(i)(0,0) = fnt(i);
// for (int i =0; i<N-1; ++i) G.on_mesh(i)(0,0) = fnt(i);
for (int i = 0; i < N - 1; ++i) G[i](0, 0) = fnt(i);
}
};
@ -43,7 +40,6 @@ struct array_code {
for (int u = 0; u < nl_interne; ++u)
for (int i = 0; i < N - 1; ++i) V(i, 0, 0) = fnt(i);
}
};
@ -52,9 +48,8 @@ struct array_code {
#include "./speed_tester.hpp"
int main() {
try {
speed_tester<array_code> (500);
speed_tester<with_sliding_view> (500);
//speed_tester<with_slices> (5000);
speed_tester<array_code>(100);
speed_tester<with_g>(100);
}
TRIQS_CATCH_AND_ABORT;
}

8
test/triqs/gfs/curry_and_fourier.cpp
View File

@ -33,11 +33,9 @@ try {
auto G_w_wn_curry0 = curry<0>(G_w_wn);
auto G_w_tau_curry0 = curry<0>(G_w_tau);
for (auto const & w : G_w_wn_curry0.mesh()) G_w_wn_curry0[w] = fourier(G_w_tau_curry0[w]);
G_w_wn_curry0[w_] << fourier(G_w_tau_curry0[w_]);
curry<0>(G_w_wn) [w_] << fourier(curry<0>(G_w_tau)[w_]);
//for (auto const & w : G_w_wn_curry0.mesh()) G_w_wn_curry0[w] = fourier(G_w_tau_curry0[w]);
//G_w_wn_curry0[w_] << fourier(G_w_tau_curry0[w_]);
//curry<0>(G_w_wn) [w_] << fourier(curry<0>(G_w_tau)[w_]);
}
// temp fix : THE TEST DOES NOT RUN !!
//TRIQS_CATCH_AND_ABORT;

2
test/triqs/gfs/g_k_om.cpp
View File

@ -16,7 +16,7 @@ int main() {
auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix<double>(2)}};
auto g_eps = gf<bz>{{bz_, 1000}, {1, 1}};
auto g_eps = gf<bz>{{bz_, 100}, {1, 1}};
auto G = gf<cartesian_product<bz, imfreq>>{{{bz_, 100}, {beta, Fermion, 100}}, {1, 1}};

40
test/triqs/gfs/g_k_tail.cpp Normal file
View File

@ -0,0 +1,40 @@
#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
#include <triqs/gfs.hpp>
#include <triqs/gfs/bz.hpp>
#include <triqs/gfs/m_tail.hpp>
namespace h5 = triqs::h5;
using namespace triqs::gfs;
using namespace triqs::clef;
using namespace triqs::arrays;
using namespace triqs::lattice;
using local::tail;
#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
// THE NAME bz is TOO SHORT
int main() {
try {
double beta = 1;
auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix<double>(2)}};
auto t = tail{1,1};
auto G = gf<bz, tail>{{bz_, 100}};
// try to assign to expression
placeholder<0> k_;
placeholder<1> w_;
auto eps = make_expr( [](k_t const& k) { return -2 * (cos(k(0)) + cos(k(1))); }) ;
G(k_) << eps(k_) * t;
// hdf5
h5::file file("ess_g_k_tail.h5", H5F_ACC_TRUNC );
h5_write(file, "g", G);
}
TRIQS_CATCH_AND_ABORT;
}

40
test/triqs/gfs/g_nu_nup.cpp Normal file
View File

@ -0,0 +1,40 @@
#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
#include <triqs/gfs.hpp>
namespace h5 = triqs::h5;
using namespace triqs::gfs;
using namespace triqs::clef;
using namespace triqs::arrays;
#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
int main() {
try {
double beta = 1;
// auto G = gf<cartesian_product<imfreq, imfreq>, tensor_valued<4>, matrix_operations>{
// {{beta, Fermion, 100}, {beta, Fermion, 100}}, {2, 2, 2, 2}};
/*
auto G = gf<imfreq_x2>{{beta, Fermion, 100}, {2, 2}};
// try to assign to expression
placeholder<0> nu_;
placeholder<1> nup_;
G(nu_, nup_) << 1 / (nu_ + nup_ + 1);
for (auto w : G.mesh())
std::cerr << G[w](0,0) <<std::endl;
// inverse
auto Ginv = inverse(G);
// hdf5
h5::file file("ess_g_nu_nup.h5", H5F_ACC_TRUNC );
h5_write(file, "g", G);
*/
}
TRIQS_CATCH_AND_ABORT;
}

44
test/triqs/gfs/g_om_nu_nup.cpp Normal file
View File

@ -0,0 +1,44 @@
#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
#include <triqs/gfs.hpp>
namespace h5 = triqs::h5;
using namespace triqs::gfs;
using namespace triqs::clef;
using namespace triqs::arrays;
#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
int main() {
try {
double beta = 1;
/*
auto G = gf<imfreq_bff, tensor_valued<4>>{{{beta, Boson, 3}, {beta, Fermion, 10}, {beta, Fermion, 10}}, {2, 2, 2, 2}};
// try to assign to expression
placeholder<2> om_;
placeholder<0> nu_;
placeholder<1> nup_;
G(om_, nu_, nup_) << 1 / (nu_ + nup_ + om_ + 1);
for (auto w : G.mesh())
std::cerr << G[w] <<std::endl;
// hdf5
h5::file file("ess_g_om_nu_nup.h5", H5F_ACC_TRUNC );
h5_write(file, "g", G);
// Curry
//auto g_cur = curry<0>(G);
auto g0 = partial_eval<0>(G,0);
//auto gg = g_cur[0];
//h5_write(file, "g[0]c", gg);
h5_write(file, "g[0]", g0);
*/
}
TRIQS_CATCH_AND_ABORT;
}

2
test/triqs/gfs/gf_scalar.cpp
View File

@ -1,4 +1,4 @@
//#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
#include <triqs/gfs.hpp>
#include <triqs/gfs/local/functions.hpp>
using namespace triqs::gfs;

2
test/triqs/gfs/gfv2.cpp
View File

@ -76,6 +76,8 @@ int main() {
TEST( G( 0) ) ;
TEST( Gc( 0) ) ;
// error
//TEST (G(1.0_j));
// operations on gf
G3 = G +2* Gc;

8
triqs/arrays/linalg/det_and_inverse.hpp
View File

@ -41,10 +41,16 @@ namespace arrays {
*/
template <typename A> struct inverse_lazy;
template <class A> inverse_lazy<typename utility::remove_rvalue_ref<A>::type> inverse(A &&a) {
template <class A>
std14::enable_if_t<ImmutableMatrix<std14::remove_reference_t<A>>::value,
inverse_lazy<typename utility::remove_rvalue_ref<A>::type>>
inverse(A &&a) {
return {std::forward<A>(a)};
}
template <class A>
std14::enable_if_t<ImmutableArray<std14::remove_reference_t<A>>::value> inverse(A &&a) = delete; // arrays can not be inverted.
// ----------------- implementation -----------------------------------------
// worker takes a contiguous view and compute the det and inverse in two steps.

2
triqs/arrays/matrix.hpp
View File

@ -217,7 +217,7 @@ namespace triqs { namespace arrays {
}
template<typename ArrayType>
matrix_view<typename ArrayType::value_type, ArrayType::opt_flags, ArrayType::traversal_order>
matrix_view<typename ArrayType::value_type, ArrayType::opt_flags, ArrayType::traversal_order, true>
make_matrix_view(ArrayType const & a) {
static_assert(ArrayType::rank ==2, "make_matrix_view only works for array of rank 2");
return a;

135
triqs/gfs/block.hpp
View File

@ -18,8 +18,7 @@
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef TRIQS_GF_BLOCK_H
#define TRIQS_GF_BLOCK_H
#pragma once
#include "./tools.hpp"
#include "./gf.hpp"
#include "./local/tail.hpp"
@ -31,7 +30,7 @@ namespace gfs {
struct block_index {};
template <typename Opt> struct gf_mesh<block_index, Opt> : discrete_mesh<discrete_domain> {
typedef discrete_mesh<discrete_domain> B;
using B = discrete_mesh<discrete_domain>;
gf_mesh() = default;
gf_mesh(int s) : B(s) {}
gf_mesh(discrete_domain const &d) : B(d) {}
@ -42,20 +41,20 @@ namespace gfs {
/// --------------------------- hdf5 ---------------------------------
template <typename Target, typename Opt> struct h5_name<block_index, Target, Opt> {
template <typename Target> struct h5_name<block_index, Target, nothing, void> {
static std::string invoke() { return "BlockGf"; }
};
template <typename Target, typename Opt> struct h5_rw<block_index, Target, Opt> {
template <typename Target> struct h5_rw<block_index, Target, nothing, void> {
static void write(h5::group gr, gf_const_view<block_index, Target, Opt> g) {
static void write(h5::group gr, gf_const_view<block_index, Target, nothing, void> g) {
for (size_t i = 0; i < g.mesh().size(); ++i) h5_write(gr, g.mesh().domain().names()[i], g._data[i]);
// h5_write(gr,"symmetry",g._symmetry);
}
template <bool IsView> static void read(h5::group gr, gf_impl<block_index, Target, Opt, IsView, false> &g) {
template <bool IsView> static void read(h5::group gr, gf_impl<block_index, Target, nothing, void, IsView, false> &g) {
// does not work : need to read the block name and remake the mesh...
g._mesh = gf_mesh<block_index, Opt>(gr.get_all_subgroup_names());
g._mesh = gf_mesh<block_index>(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 < g.mesh().size(); ++i) h5_read(gr, g.mesh().domain().names()[i], g._data[i]);
@ -65,15 +64,15 @@ namespace gfs {
/// --------------------------- data access ---------------------------------
template <typename Target, typename Opt>
struct data_proxy<block_index, Target, Opt> : data_proxy_vector<typename regular_type_if_exists_else_type<Target>::type> {};
template <typename Target>
struct data_proxy<block_index, Target, void> : data_proxy_vector<typename regular_type_if_exists_else_type<Target>::type> {};
// ------------------------------- Factories --------------------------------------------------
template <typename Target, typename Opt> struct factories<block_index, Target, Opt> {
typedef gf_mesh<block_index, Opt> mesh_t;
typedef gf<block_index, Target> gf_t;
typedef gf_view<block_index, Target> gf_view_t;
template <typename Target> struct factories<block_index, Target, nothing, void> {
using mesh_t = gf_mesh<block_index>;
using gf_t = gf<block_index, Target>;
using gf_view_t = gf_view<block_index, Target>;
struct target_shape_t {};
static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t) { return std::vector<Target>(m.size()); }
@ -93,54 +92,50 @@ namespace gfs {
// ------------------------------- Free Factories for regular type --------------------------------------------------
// from a number and a gf to be copied
template <typename Variable, typename Target, typename Opt>
block_gf<Variable, Target, Opt> make_block_gf(int n, gf<Variable, Target, Opt> const &g) {
auto V = std::vector<gf<Variable, Target, Opt>>{};
template <typename... A> block_gf<A...> make_block_gf(int n, gf<A...> const &g) {
auto V = std::vector<gf<A...>>{};
for (int i = 0; i < n; ++i) V.push_back(g);
return {{n}, std::move(V), nothing{}, nothing{}, nothing{}};
}
// from a vector of gf (moving directly)
template <typename Variable, typename Target, typename Opt>
block_gf<Variable, Target, Opt> make_block_gf(std::vector<gf<Variable, Target, Opt>> V) {
template <typename... A> block_gf<A...> make_block_gf(std::vector<gf<A...>> V) {
int s = V.size(); // DO NOT use V.size in next statement, the V is moved and the order of arg. evaluation is undefined.
return {{s}, std::move(V), nothing{}, nothing{}, nothing{}};
}
/*
// from a vector of gf : generalized to have a different type of gf in the vector (e.g. views...)
template <typename Variable, typename Target, typename Opt, typename GF2>
block_gf<Variable, Target, Opt> make_block_gf(std::vector<GF2> const &V) {
auto V2 = std::vector<gf<Variable, Target, Opt>>{};
template <typename Variable, typename Target, typename Singularity, typename Opt, typename GF2>
block_gf<Variable, Target, Singularity, Opt> make_block_gf(std::vector<GF2> const &V) {
auto V2 = std::vector<gf<Variable, Target, Singularity, Opt>>{};
for (auto const &g : V) V2.push_back(g);
return {{int(V.size())}, std::move(V2), nothing{}, nothing{}, nothing{}};
}
*/
// from a init list of GF with the correct type
template <typename Variable, typename Target, typename Opt>
block_gf<Variable, Target, Opt> make_block_gf(std::initializer_list<gf<Variable, Target, Opt>> const &V) {
template <typename... A> block_gf<A...> make_block_gf(std::initializer_list<gf<A...>> const &V) {
return {{int(V.size())}, V, nothing{}, nothing{}, nothing{}};
}
// from vector<string> and a gf to be copied
template <typename Variable, typename Target, typename Opt>
block_gf<Variable, Target, Opt> make_block_gf(std::vector<std::string> block_names, gf<Variable, Target, Opt> const &g) {
auto V = std::vector<gf<Variable, Target, Opt>>{};
template <typename... A> block_gf<A...> make_block_gf(std::vector<std::string> block_names, gf<A...> const &g) {
auto V = std::vector<gf<A...>>{};
for (int i = 0; i < block_names.size(); ++i) V.push_back(g);
return {{block_names}, std::move(V), nothing{}, nothing{}, nothing{}};
}
// from vector<string>, vector<gf>
template <typename Variable, typename Target, typename Opt>
block_gf<Variable, Target, Opt> make_block_gf(std::vector<std::string> block_names, std::vector<gf<Variable, Target, Opt>> V) {
template <typename... A> block_gf<A...> make_block_gf(std::vector<std::string> block_names, std::vector<gf<A...>> V) {
if (block_names.size() != V.size())
TRIQS_RUNTIME_ERROR << "make_block_gf(vector<string>, vector<gf>) : the two vectors do not have the same size !";
return {{block_names}, std::move(V), nothing{}, nothing{}, nothing{}};
}
// from vector<string>, init_list<GF>
template <typename Variable, typename Target, typename Opt>
block_gf<Variable, Target, Opt> make_block_gf(std::vector<std::string> block_names,
std::initializer_list<gf<Variable, Target, Opt>> const &V) {
template <typename... A>
block_gf<A...> make_block_gf(std::vector<std::string> block_names, std::initializer_list<gf<A...>> const &V) {
if (block_names.size() != V.size()) TRIQS_RUNTIME_ERROR << "make_block_gf(vector<string>, init_list) : size mismatch !";
return {{block_names}, V, nothing{}, nothing{}, nothing{}};
}
@ -148,10 +143,9 @@ namespace gfs {
// ------------------------------- Free Factories for view type --------------------------------------------------
template <typename G0, typename... G>
gf_view<block_index, typename std::remove_reference<G0>::type::view_type> make_block_gf_view(G0 &&g0, G &&... g) {
auto V = std::vector<typename std::remove_reference<G0>::type::view_type>{std::forward<G0>(g0), std::forward<G>(g)...};
gf_view<block_index, typename std14::decay_t<G0>::view_type> make_block_gf_view(G0 &&g0, G &&... g) {
auto V = std::vector<typename std14::decay_t<G0>::view_type>{std::forward<G0>(g0), std::forward<G>(g)...};
return {{int(V.size())}, std::move(V), nothing{}, nothing{}, nothing{}};
// return { gf_mesh<block_index, Opt> {int(V.size())}, std::move(V), nothing{}, nothing{} } ;
}
template <typename GF> gf_view<block_index, typename GF::regular_type> make_block_gf_view_from_vector(std::vector<GF> V) {
@ -167,44 +161,46 @@ namespace gfs {
// ------------------------------- Extend reinterpret_scalar_valued_gf_as_matrix_valued for block gf ------
template <typename Variable, typename Opt, typename Opt2, bool IsConst>
gf_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2, IsConst>
reinterpret_scalar_valued_gf_as_matrix_valued(gf_view<block_index, gf<Variable, scalar_valued, Opt>, Opt2, IsConst> bg) {
std::vector<gf_view<Variable, matrix_valued, Opt>> V;
// TODO simplify ?
template <typename Variable, typename Singularity, typename Opt, bool IsConst>
gf_view<block_index, gf<Variable, matrix_valued, Singularity, Opt>, nothing, void, IsConst>
reinterpret_scalar_valued_gf_as_matrix_valued(
gf_view<block_index, gf<Variable, scalar_valued, Singularity, Opt>, nothing, void, IsConst> bg) {
std::vector<gf_view<Variable, matrix_valued, Singularity, Opt>> V;
for (auto &g : bg) V.push_back(reinterpret_scalar_valued_gf_as_matrix_valued(g));
return make_block_gf_view_from_vector(std::move(V));
}
template <typename Variable, typename Opt, typename Opt2>
gf_const_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2>
reinterpret_scalar_valued_gf_as_matrix_valued(gf<block_index, gf<Variable, scalar_valued, Opt>, Opt2> const &bg) {
template <typename Variable, typename Singularity, typename Opt>
block_gf_const_view<Variable, matrix_valued, Singularity, Opt>
reinterpret_scalar_valued_gf_as_matrix_valued(block_gf<Variable, scalar_valued, Singularity, Opt> const &bg) {
return reinterpret_scalar_valued_gf_as_matrix_valued(bg());
}
template <typename Variable, typename Opt, typename Opt2>
gf_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2>
reinterpret_scalar_valued_gf_as_matrix_valued(gf<block_index, gf<Variable, scalar_valued, Opt>, Opt2> &bg) {
template <typename Variable, typename Singularity, typename Opt>
block_gf_view<Variable, matrix_valued, Singularity, Opt>
reinterpret_scalar_valued_gf_as_matrix_valued(block_gf<Variable, scalar_valued, Singularity, Opt> &bg) {
return reinterpret_scalar_valued_gf_as_matrix_valued(bg());
}
// ------------------------------- Free functions --------------------------------------------------
// a simple function to get the number of blocks
template <typename T> size_t n_blocks(gf<block_index, T> const &g) { return g.mesh().size(); }
template <typename T> size_t n_blocks(gf_view<block_index, T> const &g) { return g.mesh().size(); }
template <typename... T> size_t n_blocks(gf<block_index, T...> const &g) { return g.mesh().size(); }
template <typename... T> size_t n_blocks(gf_view<block_index, T...> const &g) { return g.mesh().size(); }
// ------------------------------- an iterator over the blocks --------------------------------------------------
template<typename T> using __get_target = typename std::remove_reference<decltype(std::declval<T>()[0])>::type;
template <typename T> using __get_target = std14::remove_reference_t<decltype(std::declval<T>()[0])>;
// iterator
template <typename G>
class block_gf_iterator
: public boost::iterator_facade<block_gf_iterator<G>, __get_target<G>, boost::forward_traversal_tag, __get_target<G> &> {
friend class boost::iterator_core_access;
typedef typename std::remove_reference<G>::type big_gf_t;
using big_gf_t = typename std::remove_reference<G>::type;
big_gf_t &big_gf;
typedef typename big_gf_t::mesh_t::const_iterator mesh_iterator_t;
using mesh_iterator_t = typename big_gf_t::mesh_t::const_iterator;
mesh_iterator_t mesh_it;
__get_target<G> &dereference() const { return big_gf[*mesh_it]; }
@ -217,25 +213,27 @@ namespace gfs {
};
//------------
template <typename Target, typename Opt, bool B, bool C>
block_gf_iterator<gf_impl<block_index, Target, Opt, B, C>> begin(gf_impl<block_index, Target, Opt, B, C> &bgf) {
template <typename Target, typename Singularity, typename Opt, bool B, bool C>
block_gf_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
begin(gf_impl<block_index, Target, Singularity, Opt, B, C> &bgf) {
return {bgf, false};
}
//------------
template <typename Target, typename Opt, bool B, bool C>
block_gf_iterator<gf_impl<block_index, Target, Opt, B, C>> end(gf_impl<block_index, Target, Opt, B, C> &bgf) {
template <typename Target, typename Singularity, typename Opt, bool B, bool C>
block_gf_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
end(gf_impl<block_index, Target, Singularity, Opt, B, C> &bgf) {
return {bgf, true};
}
//----- const iterator
template <typename G>
class block_gf_const_iterator
: public boost::iterator_facade<block_gf_const_iterator<G>, __get_target<G>, boost::forward_traversal_tag, __get_target<G> const &> {
class block_gf_const_iterator : public boost::iterator_facade<block_gf_const_iterator<G>, __get_target<G>,
boost::forward_traversal_tag, __get_target<G> const &> {
friend class boost::iterator_core_access;
typedef typename std::remove_reference<G>::type big_gf_t;
using big_gf_t = std14::remove_reference_t<G>;
big_gf_t const &big_gf;
typedef typename big_gf_t::mesh_t::const_iterator mesh_iterator_t;
using mesh_iterator_t = typename big_gf_t::mesh_t::const_iterator;
mesh_iterator_t mesh_it;
__get_target<G> const &dereference() const { return big_gf[*mesh_it]; }
@ -247,26 +245,29 @@ namespace gfs {
bool at_end() const { return mesh_it.at_end(); }
};
template <typename Target, typename Opt, bool B, bool C>
block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> begin(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
template <typename Target, typename Singularity, typename Opt, bool B, bool C>
block_gf_const_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
begin(gf_impl<block_index, Target, Singularity, Opt, B, C> const &bgf) {
return {bgf, false};
}
template <typename Target, typename Opt, bool B, bool C>
block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> end(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
template <typename Target, typename Singularity, typename Opt, bool B, bool C>
block_gf_const_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
end(gf_impl<block_index, Target, Singularity, Opt, B, C> const &bgf) {
return {bgf, true};
}
template <typename Target, typename Opt, bool B, bool C>
block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> cbegin(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
template <typename Target, typename Singularity, typename Opt, bool B, bool C>
block_gf_const_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
cbegin(gf_impl<block_index, Target, Singularity, Opt, B, C> const &bgf) {
return {bgf, false};
}
template <typename Target, typename Opt, bool B, bool C>
block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> cend(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
template <typename Target, typename Singularity, typename Opt, bool B, bool C>
block_gf_const_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
cend(gf_impl<block_index, Target, Singularity, Opt, B, C> const &bgf) {
return {bgf, true};
}
}
}
#endif

4
triqs/gfs/bz.hpp
View File

@ -37,7 +37,7 @@ namespace gfs {
namespace gfs_implementation {
// h5 name
template <typename Opt> struct h5_name<bz, matrix_valued, Opt> {
template <typename Singularity, typename Opt> struct h5_name<bz, matrix_valued, Singularity, Opt> {
static std::string invoke() { return "BZ"; }
};
@ -60,7 +60,7 @@ namespace gfs {
// ------------- evaluator -------------------
// handle the case where the matsu. freq is out of grid...
template <typename Target, typename Opt> struct evaluator<bz, Target, Opt> {
template <typename Target, typename Singularity, typename Opt> struct evaluator<bz, Target, Singularity, Opt> {
static constexpr int arity = 1;
template <typename G>

118
triqs/gfs/curry.hpp
View File

@ -18,13 +18,16 @@
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef TRIQS_GF_CURRY_H
#define TRIQS_GF_CURRY_H
#pragma once
#include "./product.hpp"
namespace triqs { namespace gfs {
template<typename F> struct lambda_valued {};
template <typename Var, typename M, typename L> gf_view<Var, lambda_valued<L>> make_gf_view_lambda_valued(M m, L l) {
return {std::move(m), l, nothing(), nothing(), {}};
}
namespace gfs_implementation {
/// --------------------------- data access ---------------------------------
@ -34,7 +37,7 @@ namespace triqs { namespace gfs {
/// --------------------------- Factories ---------------------------------
template<typename F, typename Opt, typename ... Ms>
struct factories<cartesian_product<Ms...>, lambda_valued<F>, Opt> {};
struct factories<cartesian_product<Ms...>, lambda_valued<F>, nothing, Opt> {};
/// --------------------------- partial_eval ---------------------------------
// partial_eval<0> (g, 1) returns : x -> g(1,x)
@ -42,82 +45,85 @@ namespace triqs { namespace gfs {
// a technical trait: from a tuple of mesh, return the mesh (either M if it is a tuple of size 1, or the corresponding cartesian_product<M..>).
template<typename ... Ms> struct cart_prod_impl;
template<typename ... Ms> using cart_prod = typename cart_prod_impl<Ms...>::type;
template<typename ... Ms> struct cart_prod_impl<std::tuple<Ms...>> { using type = cartesian_product<Ms...>;};
template<typename M> struct cart_prod_impl<std::tuple<M>> { using type = M;};
template<typename ... Ms> using cart_prod = typename cart_prod_impl<Ms...>::type;
template<typename M0, typename M1, typename ...M> auto rm_tuple_of_size_one(std::tuple<M0,M1,M...> const & t) DECL_AND_RETURN(t);
template<typename M> auto rm_tuple_of_size_one(std::tuple<M> const & t) DECL_AND_RETURN(std::get<0>(t));
// The implementation (can be overloaded for some types), so put in a struct to have partial specialization
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst> struct partial_eval_impl;
// as_tuple leaves a tuple intact and wrap everything else in a tuple...
template<typename T> std::tuple<T> as_tuple(T && x) { return std::tuple<T> {std::forward<T>(x)};}
template<typename ... T> std::tuple<T...> as_tuple(std::tuple<T...> && x) { return std::forward<T...>(x);}
template<typename ... T> std::tuple<T...> const & as_tuple(std::tuple<T...> const & x) { return x;}
template<typename ... T> std::tuple<T...> & as_tuple(std::tuple<T...> & x) { return x;}
template <int... pos, typename Opt, typename Target, bool IsConst, typename IT, typename... Ms>
gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, IsConst>
partial_eval(gf_view<cartesian_product<Ms...>, Target, Opt, IsConst> 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<pos...>(g.mesh().components());
// a view of the array of g, with the dimension sizeof...(Ms)
auto arr = reinterpret_linear_array(g.mesh(),g.data()); // NO the second () forces a view
// now rebuild a tuple of the size sizeof...(Ms), containing the indices and range at the position of evaluated variables.
auto arr_args = triqs::tuple::inverse_filter<sizeof...(Ms),pos...>(as_tuple(index), arrays::range());
// 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, std::tuple_cat(arr_args, std::make_tuple(arrays::ellipsis{})));
// finally, we build the view on this data.
using r_t = gf_view< cart_prod< triqs::tuple::filter_out_t<std::tuple<Ms...>, 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{} };
// The user function
template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt, bool C, typename... T>
auto partial_eval(gf_view<Variable, Target, Singularity, Opt, C> g, T&&... x) {
return partial_eval_impl<Variable, Target, Singularity, Opt, C>::template invoke<pos...>(g(), std::forward<T>(x)...);
}
template <int... pos, typename Opt, typename Target, typename IT, typename... Ms>
gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, false>
partial_eval(gf<cartesian_product<Ms...>, Target, Opt> & g, IT index) {
return partial_eval<pos...>(g(),index);
template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt, typename... T>
auto partial_eval(gf<Variable, Target, Singularity, Opt>& g, T&&... x) {
return partial_eval_impl<Variable, Target, Singularity, Opt, false>::template invoke<pos...>(g(), std::forward<T>(x)...);
}
template <int... pos, typename Opt, typename Target, typename IT, typename... Ms>
gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, true>
partial_eval(gf<cartesian_product<Ms...>, Target, Opt> const& g, IT index) {
return partial_eval<pos...>(g(),index);
template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt, typename... T>
auto partial_eval(gf<Variable, Target, Singularity, Opt> const& g, T&&... x) {
return partial_eval_impl<Variable, Target, Singularity, Opt, true>::template invoke<pos...>(g(), std::forward<T>(x)...);
}
/// --------------------------- curry ---------------------------------
// curry<0>(g) returns : x-> y... -> g(x,y...)
// curry<1>(g) returns : y-> x,z... -> g(x,y,z...)
// to adapt the partial_eval as a polymorphic lambda (replace by a lambda in c++14)
template<typename Gview, int ... pos> struct curry_polymorphic_lambda {
Gview g;
template<typename ...I> auto operator()(I ... i) const DECL_AND_RETURN(partial_eval<pos...>(g,std::make_tuple(i...)));
friend int get_shape(curry_polymorphic_lambda const&) { return 0;}// no shape here, but needed for compilation
//void resize(int){}
};
// curry function ...
template <int... pos, typename Target, typename Opt, bool IsConst, typename... Ms>
gf_view<cart_prod<triqs::tuple::filter_t<std::tuple<Ms...>, pos...>>,
lambda_valued<curry_polymorphic_lambda<gf_view<cartesian_product<Ms...>, Target, Opt,IsConst>, pos...>>, Opt, IsConst>
curry(gf_view<cartesian_product<Ms...>, Target, Opt, IsConst> g) {
// The implementation (can be overloaded for some types)
template <int... pos, typename Target, typename Singularity, typename Opt, bool IsConst, typename... Ms>
auto curry_impl(gf_view<cartesian_product<Ms...>, Target, Singularity, Opt, IsConst> g) {
// pick up the meshed corresponding to the curryed variables
auto meshes_tuple = triqs::tuple::filter<pos...>(g.mesh().components());
// building the view
return {rm_tuple_of_size_one(meshes_tuple),curry_polymorphic_lambda<gf_view<cartesian_product<Ms...>, Target,Opt,IsConst>, pos ...>{g}, nothing(), nothing()};
//using m_t = gf_mesh< cart_prod< triqs::tuple::filter_t<std::tuple<Ms...>,pos...>>>;
//return {triqs::tuple::apply_construct<m_t>(meshes_tuple),curry_polymorphic_lambda<gf_view<cartesian_product<Ms...>, Target,Opt>, pos ...>{g}, nothing(), nothing()};
using var_t = cart_prod<triqs::tuple::filter_t<std::tuple<Ms...>, pos...>>;
auto m = triqs::tuple::apply_construct<gf_mesh<var_t>>(meshes_tuple);
auto l = [g](auto&&... x) { return partial_eval<pos...>(g, x...); };
return make_gf_view_lambda_valued<var_t>(m, l);
};
template <int... pos, typename Target, typename Opt, typename... Ms>
auto curry(gf<cartesian_product<Ms...>, Target, Opt> & g) DECL_AND_RETURN(curry<pos...>(g()));
// The user function
template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst>
auto curry(gf_view<Variable, Target, Singularity, Opt, IsConst> g) {
return curry_impl<pos...>(g());
}
template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt>
auto curry(gf<Variable, Target, Singularity, Opt>& g) {
return curry_impl<pos...>(g());
}
template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt>
auto curry(gf<Variable, Target, Singularity, Opt> const& g) {
return curry_impl<pos...>(g());
}
template <int... pos, typename Target, typename Opt, typename... Ms>
auto curry(gf<cartesian_product<Ms...>, Target, Opt> const & g) DECL_AND_RETURN(curry<pos...>(g()));
//---------------------------------------------
// A generic impl. for cartesian product
template <typename Target, typename Singularity, typename Opt, bool IsConst, typename... Ms>
struct partial_eval_impl<cartesian_product<Ms...>, Target, Singularity, Opt, IsConst> {
template <int... pos, typename... T>
static auto invoke(gf_view<cartesian_product<Ms...>, Target, Singularity, Opt, IsConst> g, T const&... x) {
using var_t = cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>;
// meshes of the returned gf_view : just drop the mesh of the evaluated variables
auto meshes_tuple_partial = triqs::tuple::filter_out<pos...>(g.mesh().components());
auto m = triqs::tuple::apply_construct<gf_mesh<var_t>>(meshes_tuple_partial);
// now rebuild a tuple of the size sizeof...(Ms), containing the indices and range at the position of evaluated variables.
auto arr_args = triqs::tuple::inverse_filter<sizeof...(Ms), pos...>(std::make_tuple(x...), arrays::range());
// 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(g.data(), std::tuple_cat(arr_args, std::make_tuple(arrays::ellipsis{})));
auto singv = partial_eval<pos...>(g.singularity(), x...);
using r_sing_t = typename decltype(singv)::regular_type;
// finally, we build the view on this data.
using r_t = gf_view<var_t, Target, r_sing_t, Opt, IsConst>;
return r_t{m, arr2, singv, {}, {}};
}
};
} // gf_implementation
using gfs_implementation::partial_eval;
using gfs_implementation::curry;
}}
#endif

164
triqs/gfs/data_proxies.hpp
View File

@ -23,96 +23,98 @@
#include <triqs/utility/first_include.hpp>
#include <utility>
#include <triqs/arrays.hpp>
//#include "./matrix_view_proxy.hpp"
#include "../arrays/matrix_tensor_proxy.hpp"
#define TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS
namespace triqs { namespace gfs {
//---------------------------- generic case array of dim R----------------------------------
//---------------------------- common stuff for array proxies ----------------------------------
template<typename T, int R> struct data_proxy_array {
/// The storage
typedef arrays::array<T, R> storage_t;
typedef typename storage_t::view_type storage_view_t;
typedef typename storage_t::const_view_type storage_const_view_t;
template <typename T, int D> struct data_proxy_array_common {
using storage_t = arrays::array<T, D>;
using storage_view_t = typename storage_t::view_type;
using storage_const_view_t = typename storage_t::const_view_type;
/// 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
// from the shape of the mesh and the target, make the shape of the array. default is to glue them
template <typename S1, typename S2> static auto join_shape(S1 const& s1, S2 const& s2) RETURN(join(s1, s2));
template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs) { data() = std::forward<RHS>(rhs);}
template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); }
};
//---------------------------- generic case array of dim R----------------------------------
template <typename T, int R> struct data_proxy_array : data_proxy_array_common<T, R> {
using B = data_proxy_array_common<T, R>;
/// The data access
#ifdef TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS
template <typename S> auto operator()(S& data, long i) const DECL_AND_RETURN(data(i, arrays::ellipsis()));
#else
auto operator()(B::storage_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i));
auto operator()(B::storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
auto operator()(B::storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i));
auto operator()(B::storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
auto operator()(B::storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
auto operator()(B::storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
#endif
};
//---------------------------- 3d array : returns matrices in this case ! ----------------------------------
template<typename T> struct data_proxy_array<T,3> {
/// The storage
typedef arrays::array<T,3> 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<storage_t,0> operator()(storage_t & data, size_t i) const { return arrays::matrix_view_proxy<storage_t,0>(data,i); }
arrays::const_matrix_view_proxy<storage_t,0> operator()(storage_t const & data, size_t i) const { return arrays::const_matrix_view_proxy<storage_t,0>(data,i); }
arrays::matrix_view_proxy<storage_view_t,0> operator()(storage_view_t & data, size_t i) const { return arrays::matrix_view_proxy<storage_view_t,0>(data,i); }
arrays::const_matrix_view_proxy<storage_view_t,0> operator()(storage_view_t const & data, size_t i) const { return arrays::const_matrix_view_proxy<storage_view_t,0>(data,i); }
#endif
template <typename T> struct data_proxy_array<T, 3> : data_proxy_array_common<T, 3> {
using B = data_proxy_array_common<T, 3>;
#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()));
template <typename S> auto operator()(S & data, long i) const RETURN(make_matrix_view(data(i, arrays::ellipsis())));
#else
/// The data access
auto operator()(B::storage_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i));
auto operator()(B::storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
auto operator()(B::storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i));
auto operator()(B::storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
auto operator()(B::storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
auto operator()(B::storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
#endif
template <typename S, typename RHS> static void assign_to_scalar(S& data, RHS&& rhs) { data() = std::forward<RHS>(rhs); }
template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); }
};
//---------------------------- 1d array ----------------------------------
template<typename T> struct data_proxy_array<T,1>{
/// The storage
typedef arrays::array<T,1> storage_t;
typedef typename storage_t::view_type storage_view_t;
typedef typename storage_t::const_view_type storage_const_view_t;
template <typename T> struct data_proxy_array<T, 1> : data_proxy_array_common<T, 1> {
template <typename S> AUTO_DECL operator()(S& data, long i) const RETURN(data(i));
};
/// 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);}
//---------------------------- multi variable ----------------------------------
template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs) { data() = std::forward<RHS>(rhs);}
template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); }
template <typename T, int TotalDim> struct data_proxy_array_multivar : data_proxy_array_common<T, TotalDim> {
// using the standard technique from tuple::apply with a sequence
template <typename S, typename Tu, size_t... Is>
AUTO_DECL _impl(S& data, Tu const& tu, std14::index_sequence<Is...>) const RETURN(data(std::get<Is>(tu)..., arrays::ellipsis()));
template <typename S, typename Tu>
AUTO_DECL operator()(S& data, Tu const& tu) const RETURN(_impl(data, tu, triqs::tuple::_get_seq<Tu>()));
};
//---------------------------- multi variable ----------------------------------
template <typename T, int TotalDim> struct data_proxy_array_multivar_matrix_valued : data_proxy_array_common<T, TotalDim> {
// using the standard technique from tuple::apply with a sequence
template <typename S, typename Tu, size_t... Is>
AUTO_DECL _impl(S& data, Tu const& tu, std14::index_sequence<Is...>) const RETURN(make_matrix_view(data(std::get<Is>(tu)..., arrays::range(), arrays::range())));
template <typename S, typename Tu>
AUTO_DECL operator()(S& data, Tu const& tu) const RETURN(_impl(data, tu, triqs::tuple::_get_seq<Tu>()));
};
//---------------------------- multi variable with index mixer----------------------------------
template <typename T, int Nvar, int TargetDim, typename IndexMixer>
struct data_proxy_array_index_mixer : data_proxy_array_common<T, Nvar + TargetDim> {
template <typename S1, typename S2> static utility::mini_vector<int, Nvar + TargetDim> join_shape(S1 const& s1, S2 const& s2) {
return tuple::apply_construct_parenthesis<utility::mini_vector<int, Nvar + TargetDim>>(IndexMixer::invoke(s1, s2));
}
template <typename S, typename Tu> auto operator()(S& data, Tu const& tu) const {
return tuple::apply(data, IndexMixer::invoke(tu, tuple::make_tuple_repeat<TargetDim>(arrays::range())));
}
};
//---------------------------- vector ----------------------------------
@ -129,25 +131,24 @@ namespace triqs { namespace gfs {
};
template <typename T> struct data_proxy_vector {
typedef typename T::view_type Tv;
typedef typename T::const_view_type Tcv;
using Tv = typename T::view_type;
using Tcv = typename T::const_view_type;
/// The storage
typedef std::vector<T> storage_t;
typedef std::vector<view_proxy<Tv>> storage_view_t;
typedef std::vector<view_proxy<Tcv>> storage_const_view_t;
using storage_t = std::vector<T>;
using storage_view_t = std::vector<view_proxy<Tv>>;
using storage_const_view_t = std::vector<view_proxy<Tcv>>;
/// The data access
template <typename S> AUTO_DECL operator()(S& data, size_t i) const 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<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs) {for (size_t i =0; i<data.size(); ++i) data[i] = rhs;}
@ -159,13 +160,12 @@ namespace triqs { namespace gfs {
template <typename F> struct data_proxy_lambda {
/// The storage
typedef F storage_t;
typedef F storage_view_t;
typedef F storage_const_view_t;
using storage_t = F;
using storage_view_t = F;
using storage_const_view_t = F;
/// 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 <typename S, typename ... I> AUTO_DECL operator()(S& data, I const& ...i) const RETURN(data(i...));
template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs) = delete;
template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) = delete;

400
triqs/gfs/gf.hpp
View File

@ -2,7 +2,7 @@
*
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
*
* Copyright (C) 2012-2013 by M. Ferrero, O. Parcollet
* Copyright (C) 2012-2014 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
@ -42,40 +42,58 @@ namespace gfs {
using arrays::matrix_view;
using triqs::make_clone;
// The default target, for each Variable.
template <typename Variable> struct gf_default_target {
using type = matrix_valued;
};
template <typename Variable> using gf_default_target_t = typename gf_default_target<Variable>::type;
// The default singularity, for each Variable.
template <typename Variable, typename Target> struct gf_default_singularity {
using type = nothing;
};
template <typename Variable, typename Target>
using gf_default_singularity_t = typename gf_default_singularity<Variable, Target>::type;
// the gf mesh
template <typename Variable, typename Opt = void> struct gf_mesh;
// The regular type
template <typename Variable, typename Target = matrix_valued, typename Opt = void> class gf;
template <typename Variable, typename Target = gf_default_target_t<Variable>,
typename Singularity = gf_default_singularity_t<Variable, Target>, typename Opt = void>
class gf;
// The view type
template <typename Variable, typename Target = matrix_valued, typename Opt = void, bool IsConst = false> class gf_view;
template <typename Variable, typename Target = gf_default_target_t<Variable>,
typename Singularity = gf_default_singularity_t<Variable, Target>, typename Opt = void, bool IsConst = false>
class gf_view;
// The const view type
template <typename Variable, typename Target = matrix_valued, typename Opt = void>
using gf_const_view = gf_view<Variable, Target, Opt, true>;
template <typename Variable, typename Target = gf_default_target_t<Variable>,
typename Singularity = gf_default_singularity_t<Variable, Target>, typename Opt = void>
using gf_const_view = gf_view<Variable, Target, Singularity, Opt, true>;
// the implementation class
template <typename Variable, typename Target, typename Opt, bool IsView, bool IsConst> class gf_impl;
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst> class gf_impl;
// various implementation traits
namespace gfs_implementation { // never use using of this...
// evaluator regroup functions to evaluate the function.
template <typename Variable, typename Target, typename Opt> struct evaluator {
template <typename Variable, typename Target, typename Singularity, typename Opt> struct evaluator {
static constexpr int arity = 0;
};
// closest_point mechanism
template <typename Variable, typename Target, typename Opt> struct get_closest_point;
template <typename Variable, typename Target, typename Singularity, typename Opt> struct get_closest_point;
// singularity
template <typename Variable, typename Target, typename Opt> struct singularity {
using type = nothing;
};
//template <typename Variable, typename Target, typename Opt> struct singularity {
// using type = nothing;
//};
// symmetry
template <typename Variable, typename Target, typename Opt> struct symmetry {
template <typename Variable, typename Target, typename Singularity, typename Opt> struct symmetry {
using type = nothing;
};
@ -94,21 +112,29 @@ namespace gfs {
template <typename Variable, typename Target, typename Opt, typename Enable = void> struct data_proxy;
// Traits to read/write in hdf5 files. Can be specialized for some case (Cf block). Defined below
template <typename Variable, typename Target, typename Opt> struct h5_name; // value is a const char
template <typename Variable, typename Target, typename Opt> struct h5_rw;
template <typename Variable, typename Target, typename Singularity, typename Opt> struct h5_name; // value is a const char
template <typename Variable, typename Target, typename Singularity, typename Opt> struct h5_rw;
// factories regroup all factories (constructors..) for all types of gf. Defaults implemented below.
template <typename Variable, typename Target, typename Opt> struct factories;
template <typename Variable, typename Target, typename Singularity, typename Opt> struct factories;
} // gfs_implementation
// The trait that "marks" the Green function
TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT(ImmutableGreenFunction);
template <typename G> auto get_gf_data_shape(G const &g) DECL_AND_RETURN(g.get_data_shape());
template <typename Variable, typename Target, typename Singularity, typename Opt>
auto get_gf_data_shape(gf<Variable, Target, Singularity, Opt> const &g) RETURN(get_shape(g.data()));
template <typename Variable, typename Target, typename Opt, bool IsView, bool IsConst>
auto get_target_shape(gf_impl<Variable, Target, Opt, IsView, IsConst> const &g) DECL_AND_RETURN(g.data().shape().front_pop());
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst>
auto get_gf_data_shape(gf_view<Variable, Target, Singularity, Opt, IsConst> const &g) RETURN(get_shape(g.data()));
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst>
auto get_gf_data_shape(gf_impl<Variable, Target, Singularity, Opt, IsView, IsConst> const &g) RETURN(get_shape(g.data()));
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst>
auto get_target_shape(gf_impl<Variable, Target, Singularity, Opt, IsView, IsConst> const &g)
DECL_AND_RETURN(g.data().shape().front_pop());
// ---------------------- implementation --------------------------------
@ -116,15 +142,15 @@ namespace gfs {
template <typename T> long get_shape(std::vector<T> const &x) { return x.size(); }
/// A common implementation class for gf and gf_view. They will only redefine contructor and = ...
template <typename Variable, typename Target, typename Opt, bool IsView, bool IsConst>
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst>
class gf_impl : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction) {
static_assert(!(!IsView && IsConst), "Internal error");
public:
using mutable_view_type = gf_view<Variable, Target, Opt>;
using const_view_type = gf_const_view<Variable, Target, Opt>;
using mutable_view_type = gf_view<Variable, Target, Singularity, Opt>;
using const_view_type = gf_const_view<Variable, Target, Singularity, Opt>;
using view_type = typename std::conditional<IsConst, const_view_type, mutable_view_type>::type;
using regular_type = gf<Variable, Target, Opt>;
using regular_type = gf<Variable, Target, Singularity, Opt>;
using variable_t = Variable;
using target_t = Target;
@ -134,9 +160,9 @@ namespace gfs {
using domain_t = typename mesh_t::domain_t;
using mesh_point_t = typename mesh_t::mesh_point_t;
using mesh_index_t = typename mesh_t::index_t;
using symmetry_t = typename gfs_implementation::symmetry<Variable, Target, Opt>::type;
using symmetry_t = typename gfs_implementation::symmetry<Variable, Target, Singularity, Opt>::type;
using indices_t = typename gfs_implementation::indices<Target, Opt>::type;
using evaluator_t = gfs_implementation::evaluator<Variable, Target, Opt>;
using evaluator_t = gfs_implementation::evaluator<Variable, Target, Singularity, Opt>;
using data_proxy_t = gfs_implementation::data_proxy<Variable, Target, Opt>;
using data_regular_t = typename data_proxy_t::storage_t;
@ -144,10 +170,13 @@ namespace gfs {
using data_const_view_t = typename data_proxy_t::storage_const_view_t;
using data_t = std14::conditional_t<IsView, std14::conditional_t<IsConst, data_const_view_t, data_view_t>, data_regular_t>;
using singularity_non_view_t = typename gfs_implementation::singularity<Variable, Target, Opt>::type;
//using singularity_non_view_t = typename gfs_implementation::singularity<Variable, Target, Opt>::type;
using singularity_non_view_t = Singularity;
using singularity_view_t = typename view_type_if_exists_else_type<singularity_non_view_t>::type;
using singularity_t = std14::conditional_t<IsView, singularity_view_t, singularity_non_view_t>;
//using is_container_t = void; // a tag to recognize the container
mesh_t const &mesh() const { return _mesh; }
domain_t const &domain() const { return _mesh.domain(); }
data_t &data() { return _data; }
@ -158,8 +187,6 @@ namespace gfs {
indices_t const &indices() const { return _indices; }
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;
@ -259,8 +286,8 @@ namespace gfs {
// ------------- All the [] operators -----------------------------
// [] and access to the grid point
using r_type = typename std::result_of<data_proxy_t(data_t &, size_t)>::type;
using cr_type = typename std::result_of<data_proxy_t(data_t const &, size_t)>::type;
using r_type = std14::result_of_t<data_proxy_t(data_t &, typename mesh_t::linear_index_t)>;
using cr_type = std14::result_of_t<data_proxy_t(data_t const &, typename mesh_t::linear_index_t)>;
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)); }
@ -269,12 +296,12 @@ namespace gfs {
cr_type operator[](mesh_point_t const &x) const { return _data_proxy(_data, x.linear_index()); }
template <typename... U> r_type operator[](closest_pt_wrap<U...> const &p) {
return _data_proxy(_data,
_mesh.index_to_linear(gfs_implementation::get_closest_point<Variable, Target, Opt>::invoke(this, p)));
return _data_proxy(
_data, _mesh.index_to_linear(gfs_implementation::get_closest_point<Variable, Target, Singularity, Opt>::invoke(this, p)));
}
template <typename... U> cr_type operator[](closest_pt_wrap<U...> const &p) const {
return _data_proxy(_data,
_mesh.index_to_linear(gfs_implementation::get_closest_point<Variable, Target, Opt>::invoke(this, p)));
return _data_proxy(
_data, _mesh.index_to_linear(gfs_implementation::get_closest_point<Variable, Target, Singularity, Opt>::invoke(this, p)));
}
template <typename Arg>
@ -321,17 +348,17 @@ namespace gfs {
//----------------------------- HDF5 -----------------------------
friend std::string get_triqs_hdf5_data_scheme(gf_impl const &g) {
auto s = gfs_implementation::h5_name<Variable, Target, Opt>::invoke();
auto s = gfs_implementation::h5_name<Variable, Target, Singularity, Opt>::invoke();
return (s == "BlockGf" ? s : "Gf" + s);
}
friend struct gfs_implementation::h5_rw<Variable, Target, Opt>;
friend struct gfs_implementation::h5_rw<Variable, Target, Singularity, Opt>;
/// 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<Variable, Target, Opt>::write(gr, g);
gfs_implementation::h5_rw<Variable, Target, Singularity, Opt>::write(gr, g);
}
/// Read from HDF5
@ -343,7 +370,7 @@ namespace gfs {
if (tag_file != tag_expected)
TRIQS_RUNTIME_ERROR << "h5_read : mismatch of the tag TRIQS_HDF5_data_scheme tag in the h5 group : found " << tag_file
<< " while I expected " << tag_expected;
gfs_implementation::h5_rw<Variable, Target, Opt>::read(gr, g);
gfs_implementation::h5_rw<Variable, Target, Singularity, Opt>::read(gr, g);
}
//----------------------------- BOOST Serialization -----------------------------
@ -365,8 +392,8 @@ namespace gfs {
// -------------------------Interaction with the CLEF library : auto assignement implementation -----------------
// auto assignment of the gf (gf(om_) << expression fills the functions by evaluation of expression)
template <typename RHS, typename Variable, typename Target, typename Opt, bool IsView>
void triqs_clef_auto_assign(gf_impl<Variable, Target, Opt, IsView, false> &g, RHS const &rhs) {
template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt, bool IsView>
void triqs_clef_auto_assign(gf_impl<Variable, Target, Singularity, Opt, IsView, false> &g, RHS const &rhs) {
triqs_clef_auto_assign_impl(g, rhs, typename std::is_base_of<tag::composite, gf_mesh<Variable, Opt>>::type());
assign_from_expression(g.singularity(), rhs);
// access to the data . Beware, we view it as a *matrix* NOT an array... (crucial for assignment to scalars !)
@ -375,8 +402,8 @@ namespace gfs {
}
// enable the writing g[om_] << .... also
template <typename RHS, typename Variable, typename Target, typename Opt, bool IsView>
void triqs_clef_auto_assign_subscript(gf_impl<Variable, Target, Opt, IsView, false> &g, RHS const &rhs) {
template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt, bool IsView>
void triqs_clef_auto_assign_subscript(gf_impl<Variable, Target, Singularity, Opt, IsView, false> &g, RHS const &rhs) {
triqs_clef_auto_assign(g, rhs);
}
@ -389,8 +416,8 @@ namespace gfs {
triqs_clef_auto_assign(std::forward<G>(g), std::forward<clef::make_fun_impl<Expr, Is...>>(rhs));
}
template <typename RHS, typename Variable, typename Target, typename Opt, bool IsView>
void triqs_clef_auto_assign_impl(gf_impl<Variable, Target, Opt, IsView, false> &g, RHS const &rhs,
template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt, bool IsView>
void triqs_clef_auto_assign_impl(gf_impl<Variable, Target, Singularity, Opt, IsView, false> &g, RHS const &rhs,
std::integral_constant<bool, false>) {
for (auto const &w : g.mesh()) {
triqs_gf_clef_auto_assign_impl_aux_assign(g[w], rhs(w));
@ -398,8 +425,8 @@ namespace gfs {
}
}
template <typename RHS, typename Variable, typename Target, typename Opt, bool IsView>
void triqs_clef_auto_assign_impl(gf_impl<Variable, Target, Opt, IsView, false> &g, RHS const &rhs,
template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt, bool IsView>
void triqs_clef_auto_assign_impl(gf_impl<Variable, Target, Singularity, Opt, IsView, false> &g, RHS const &rhs,
std::integral_constant<bool, true>) {
for (auto const &w : g.mesh()) {
triqs_gf_clef_auto_assign_impl_aux_assign(g[w], triqs::tuple::apply(rhs, w.components_tuple()));
@ -409,19 +436,19 @@ namespace gfs {
// -------------------------The regular class of GF --------------------------------------------------------
template <typename Variable, typename Target, typename Opt> class gf : public gf_impl<Variable, Target, Opt, false, false> {
using B = gf_impl<Variable, Target, Opt, false, false>;
using factory = gfs_implementation::factories<Variable, Target, Opt>;
template <typename Variable, typename Target, typename Singularity, typename Opt> class gf : public gf_impl<Variable, Target, Singularity, Opt, false, false> {
using B = gf_impl<Variable, Target, Singularity, Opt, false, false>;
using factory = gfs_implementation::factories<Variable, Target, Singularity, Opt>;
public:
gf() : B() {}
gf(gf const &g) : B(g) {}
gf(gf &&g) noexcept : B(std::move(g)) {}
gf(gf_view<Variable, Target, Opt> const &g) : B(g, bool {}) {}
gf(gf_const_view<Variable, Target, Opt> const &g) : B(g, bool {}) {}
gf(gf_view<Variable, Target, Singularity, Opt> const &g) : B(g, bool {}) {}
gf(gf_const_view<Variable, Target, Singularity, Opt> const &g) : B(g, bool {}) {}
template <typename GfType>
gf(GfType const &x, typename std::enable_if<ImmutableGreenFunction<GfType>::value>::type *dummy = 0)
gf(GfType const &x, std14::enable_if_t<ImmutableGreenFunction<GfType>::value> *dummy = 0)
: B() {
*this = x;
}
@ -481,27 +508,27 @@ namespace gfs {
// But in chained clef expression, A(i_)(om_) where A is an array of gf
// we need to call it with the gf, not gf_impl (or the template resolution find the deleted funciton in clef).
// Another fix is to make gf, gf_view in the expression tree, but this requires using CRPT in gf_impl...
template <typename RHS, typename Variable, typename Target, typename Opt>
void triqs_clef_auto_assign(gf<Variable, Target, Opt> &g, RHS const &rhs) {
triqs_clef_auto_assign( static_cast<gf_impl<Variable, Target, Opt, false, false>&>(g), rhs);
template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt>
void triqs_clef_auto_assign(gf<Variable, Target, Singularity, Opt> &g, RHS const &rhs) {
triqs_clef_auto_assign( static_cast<gf_impl<Variable, Target, Singularity, Opt, false, false>&>(g), rhs);
}
// --------------------------The const View class of GF -------------------------------------------------------
template <typename Variable, typename Target, typename Opt>
class gf_view<Variable, Target, Opt, true> : public gf_impl<Variable, Target, Opt, true, true> {
using B = gf_impl<Variable, Target, Opt, true, true>;
template <typename Variable, typename Target, typename Singularity, typename Opt>
class gf_view<Variable, Target, Singularity, Opt, true> : public gf_impl<Variable, Target, Singularity, Opt, true, true> {
using B = gf_impl<Variable, Target, Singularity, Opt, true, true>;
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<Variable, Target, Opt, true, true> const &g) : B(g, bool {}) {} // from a const_view
gf_view(gf_impl<Variable, Target, Opt, true, false> const &g) : B(g, bool {}) {} // from a view
gf_view(gf_impl<Variable, Target, Opt, false, false> const &g) : B(g, bool {}) {} // from a const gf
gf_view(gf_impl<Variable, Target, Opt, false, false> &g) : B(g, bool {}) {} // from a gf &
gf_view(gf_impl<Variable, Target, Opt, false, false> &&g) noexcept : B(std::move(g), bool {}) {} // from a gf &&
gf_view(gf_impl<Variable, Target, Singularity, Opt, true, true> const &g) : B(g, bool {}) {} // from a const_view
gf_view(gf_impl<Variable, Target, Singularity, Opt, true, false> const &g) : B(g, bool {}) {} // from a view
gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> const &g) : B(g, bool {}) {} // from a const gf
gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> &g) : B(g, bool {}) {} // from a gf &
gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> &&g) noexcept : B(std::move(g), bool {}) {} // from a gf &&
template <typename D>
gf_view(typename B::mesh_t const &m, D const &dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s,
@ -517,7 +544,7 @@ namespace gfs {
this->name = X.name;
}
void rebind(gf_view<Variable, Target, Opt, false> const &X) noexcept {
void rebind(gf_view<Variable, Target, Singularity, Opt, false> const &X) noexcept {
rebind(gf_view{X});
}
@ -526,20 +553,20 @@ namespace gfs {
// ------------------------- The View class of GF -------------------------------------------------------
template <typename Variable, typename Target, typename Opt>
class gf_view<Variable, Target, Opt, false> : public gf_impl<Variable, Target, Opt, true, false> {
using B = gf_impl<Variable, Target, Opt, true, false>;
template <typename Variable, typename Target, typename Singularity, typename Opt>
class gf_view<Variable, Target, Singularity, Opt, false> : public gf_impl<Variable, Target, Singularity, Opt, true, false> {
using B = gf_impl<Variable, Target, Singularity, Opt, true, false>;
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<Variable, Target, Opt, true, true> const &g) = delete; // from a const view : impossible
gf_view(gf_impl<Variable, Target, Opt, true, false> const &g) : B(g, bool {}) {} // from a view
gf_view(gf_impl<Variable, Target, Opt, false, false> const &g) = delete; // from a const gf : impossible
gf_view(gf_impl<Variable, Target, Opt, false, false> &g) : B(g, bool {}) {} // from a gf &
gf_view(gf_impl<Variable, Target, Opt, false, false> &&g) noexcept : B(std::move(g), bool {}) {} // from a gf &&
gf_view(gf_impl<Variable, Target, Singularity, Opt, true, true> const &g) = delete; // from a const view : impossible
gf_view(gf_impl<Variable, Target, Singularity, Opt, true, false> const &g) : B(g, bool {}) {} // from a view
gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> const &g) = delete; // from a const gf : impossible
gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> &g) : B(g, bool {}) {} // from a gf &
gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> &&g) noexcept : B(std::move(g), bool {}) {} // from a gf &&
template <typename D>
gf_view(typename B::mesh_t const &m, D &&dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s,
@ -569,69 +596,81 @@ namespace gfs {
}; // class gf_view
// delegate = so that I can overload it for specific RHS...
template <typename Variable, typename Target, typename Opt, typename RHS>
DISABLE_IF(arrays::is_scalar<RHS>) triqs_gf_view_assign_delegation(gf_view<Variable, Target, Opt> g, RHS const &rhs) {
template <typename Variable, typename Target, typename Singularity, typename Opt, typename RHS>
DISABLE_IF(arrays::is_scalar<RHS>) triqs_gf_view_assign_delegation(gf_view<Variable, Target, Singularity, Opt> g, RHS const &rhs) {
if (!(g.mesh() == rhs.mesh()))
TRIQS_RUNTIME_ERROR << "Gf Assignment in View : incompatible mesh" << g.mesh() << " vs " << rhs.mesh();
for (auto const &w : g.mesh()) g[w] = rhs[w];
g.singularity() = rhs.singularity();
}
template <typename Variable, typename Target, typename Opt, typename T>
ENABLE_IF(arrays::is_scalar<T>) triqs_gf_view_assign_delegation(gf_view<Variable, Target, Opt> g, T const &x) {
gf_view<Variable, Target, Opt>::data_proxy_t::assign_to_scalar(g.data(), x);
template <typename Variable, typename Target, typename Singularity, typename Opt, typename T>
ENABLE_IF(arrays::is_scalar<T>) triqs_gf_view_assign_delegation(gf_view<Variable, Target, Singularity, Opt> g, T const &x) {
gf_view<Variable, Target, Singularity, Opt>::data_proxy_t::assign_to_scalar(g.data(), x);
g.singularity() = x;
}
// tool for lazy transformation
template <typename Tag, typename D, typename Target = matrix_valued, typename Opt = void> struct gf_keeper {
gf_const_view<D, Target, Opt> g;
// TODO : clean this : typename G
template <typename Tag, typename D, typename Target = matrix_valued,
typename Singularity = gf_default_singularity_t<D, Target>, typename Opt = void>
struct gf_keeper {
gf_const_view<D, Target, Singularity, Opt> g;
};
// Cf gf
template <typename RHS, typename Variable, typename Target, typename Opt>
void triqs_clef_auto_assign(gf_view<Variable, Target, Opt> &g, RHS const &rhs) {
triqs_clef_auto_assign( static_cast<gf_impl<Variable, Target, Opt, true, false>&>(g), rhs);
template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt>
void triqs_clef_auto_assign(gf_view<Variable, Target, Singularity, Opt> &g, RHS const &rhs) {
triqs_clef_auto_assign( static_cast<gf_impl<Variable, Target, Singularity, Opt, true, false>&>(g), rhs);
}
// ---------------------------------- slicing ------------------------------------
// slice
template <typename Variable, typename Target, typename Opt, bool IsConst, typename... Args>
gf_view<Variable, matrix_valued, Opt, IsConst> slice_target(gf_view<Variable, Target, Opt, IsConst> g, Args &&... args) {
static_assert(std::is_same<Target, matrix_valued>::value, "slice_target only for matrix_valued GF's");
using arrays::range;
return {g.mesh(), g.data()(range(), std::forward<Args>(args)...),
slice_target(g.singularity(), std::forward<Args>(args)...), g.symmetry(), slice(g.indices(),std::forward<Args>(args)...), g.name };
}
template <typename Variable, typename Target, typename Opt, typename... Args>
gf_view<Variable, matrix_valued, Opt> slice_target(gf<Variable, Target, Opt> &g, Args &&... args) {
return slice_target(g(), std::forward<Args>(args)...);
}
template <typename Variable, typename Target, typename Opt, typename... Args>
gf_const_view<Variable, matrix_valued, Opt> slice_target(gf<Variable, Target, Opt> const &g, Args &&... args) {
return slice_target(g(), std::forward<Args>(args)...);
}
// slice to scalar
template <typename Variable, typename Target, typename Opt, bool IsConst, typename... Args>
gf_view<Variable, scalar_valued, Opt, IsConst> slice_target_to_scalar(gf_view<Variable, Target, Opt, IsConst> g,
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst, typename... Args>
gf_view<Variable, matrix_valued, Singularity, Opt, IsConst> slice_target(gf_view<Variable, Target, Singularity, Opt, IsConst> g,
Args &&... args) {
static_assert(std::is_same<Target, matrix_valued>::value, "slice_target only for matrix_valued GF's");
using arrays::range;
return {g.mesh(), g.data()(range(), std::forward<Args>(args)...),
slice_target(g.singularity(), range(args, args + 1)...), g.symmetry(), {}, g.name};
slice_target(g.singularity(), std::forward<Args>(args)...), g.symmetry(),
slice(g.indices(), std::forward<Args>(args)...), g.name};
}
template <typename Variable, typename Target, typename Opt, typename... Args>
gf_view<Variable, scalar_valued, Opt> slice_target_to_scalar(gf<Variable, Target, Opt> &g, Args &&... args) {
template <typename Variable, typename Target, typename Singularity, typename Opt, typename... Args>
gf_view<Variable, matrix_valued, Singularity, Opt> slice_target(gf<Variable, Target, Singularity, Opt> &g, Args &&... args) {
return slice_target(g(), std::forward<Args>(args)...);
}
template <typename Variable, typename Target, typename Singularity, typename Opt, typename... Args>
gf_const_view<Variable, matrix_valued, Singularity, Opt> slice_target(gf<Variable, Target, Singularity, Opt> const &g,
Args &&... args) {
return slice_target(g(), std::forward<Args>(args)...);
}
// slice to scalar
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst, typename... Args>
gf_view<Variable, scalar_valued, Singularity, Opt, IsConst>
slice_target_to_scalar(gf_view<Variable, Target, Singularity, Opt, IsConst> g, Args &&... args) {
static_assert(std::is_same<Target, matrix_valued>::value, "slice_target only for matrix_valued GF's");
using arrays::range;
return {g.mesh(),
g.data()(range(), std::forward<Args>(args)...),
slice_target(g.singularity(), range(args, args + 1)...),
g.symmetry(),
{},
g.name};
}
template <typename Variable, typename Target, typename Singularity, typename Opt, typename... Args>
gf_view<Variable, scalar_valued, Singularity, Opt> slice_target_to_scalar(gf<Variable, Target, Singularity, Opt> &g,
Args &&... args) {
return slice_target_to_scalar(g(), std::forward<Args>(args)...);
}
template <typename Variable, typename Target, typename Opt, typename... Args>
gf_const_view<Variable, scalar_valued, Opt> slice_target_to_scalar(gf<Variable, Target, Opt> const &g, Args &&... args) {
template <typename Variable, typename Target, typename Singularity, typename Opt, typename... Args>
gf_const_view<Variable, scalar_valued, Singularity, Opt> slice_target_to_scalar(gf<Variable, Target, Singularity, Opt> const &g,
Args &&... args) {
return slice_target_to_scalar(g(), std::forward<Args>(args)...);
}
@ -647,32 +686,26 @@ namespace gfs {
// ---------------------------------- target reinterpretation ------------------------------------
// a scalar_valued gf can be viewed as a 1x1 matrix
template <typename Variable, typename Opt, bool IsConst>
gf_view<Variable, matrix_valued, Opt, IsConst>
reinterpret_scalar_valued_gf_as_matrix_valued(gf_view<Variable, scalar_valued, Opt, IsConst> g) {
using a_t = typename gf_view<Variable, matrix_valued, Opt, IsConst>::data_view_t;
template <typename Variable, typename Singularity, typename Opt, bool IsConst>
gf_view<Variable, matrix_valued, Singularity, Opt, IsConst>
reinterpret_scalar_valued_gf_as_matrix_valued(gf_view<Variable, scalar_valued, Singularity, Opt, IsConst> g) {
using a_t = typename gf_view<Variable, matrix_valued, Singularity, Opt, IsConst>::data_view_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(), {}, g.name};
}
template <typename Variable, typename Opt>
gf_view<Variable, matrix_valued, Opt> reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Opt> &g) {
template <typename Variable, typename Singularity, typename Opt>
gf_view<Variable, matrix_valued, Singularity, Opt>
reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Singularity, Opt> &g) {
return reinterpret_scalar_valued_gf_as_matrix_valued(g());
}
template <typename Variable, typename Opt>
gf_const_view<Variable, matrix_valued, Opt>
reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Opt> const &g) {
template <typename Variable, typename Singularity, typename Opt>
gf_const_view<Variable, matrix_valued, Singularity, Opt>
reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Singularity, Opt> const &g) {
return reinterpret_scalar_valued_gf_as_matrix_valued(g());
}
/*
template<typename Variable1,typename Variable2, typename Target, typename Opt, bool V, typename... Args>
gf_view<Variable2,Target,Opt> slice_mesh (gf_impl<Variable1,Target,Opt,V> const & g, Args... args) {
return gf_view<Variable2,Target,Opt>(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()),
g.singularity(), g.symmetry());
}*/
// ---------------------------------- some functions : invert, conjugate, transpose, ... ------------------------------------
// ---- inversion
@ -684,33 +717,36 @@ namespace gfs {
}
}
template <typename Variable, typename Opt>
void invert_in_place(gf_view<Variable, matrix_valued, Opt> g) {
template <typename Variable, typename Singularity, typename Opt>
void invert_in_place(gf_view<Variable, matrix_valued, Singularity, Opt> g) {
_gf_invert_data_in_place(g.data());
g.singularity() = inverse(g.singularity());
}
template <typename Variable, typename Opt> gf<Variable, matrix_valued, Opt> inverse(gf<Variable, matrix_valued, Opt> const & g) {
template <typename Variable, typename Singularity, typename Opt>
gf<Variable, matrix_valued, Singularity, Opt> inverse(gf<Variable, matrix_valued, Singularity, Opt> const &g) {
auto res = g;
gf_view<Variable, matrix_valued, Opt> v = res;
gf_view<Variable, matrix_valued, Singularity, Opt> v = res;
invert_in_place(v);
return res;
}
template <typename Variable, typename Opt, bool B> gf<Variable, matrix_valued, Opt> inverse(gf_view<Variable, matrix_valued, Opt, B> g) {
return inverse(gf<Variable, matrix_valued, Opt>(g));
template <typename Variable, typename Singularity, typename Opt, bool B>
gf<Variable, matrix_valued, Singularity, Opt> inverse(gf_view<Variable, matrix_valued, Singularity, Opt, B> g) {
return inverse(gf<Variable, matrix_valued, Singularity, Opt>(g));
}
// ---- transpose : a new gf
template <typename Variable, typename Opt>
gf<Variable, matrix_valued, Opt> transpose(gf_view<Variable, matrix_valued, Opt> g) {
template <typename Variable, typename Singularity, typename Opt>
gf<Variable, matrix_valued, Singularity, Opt> transpose(gf_view<Variable, matrix_valued, Singularity, Opt> g) {
return {g.mesh(), transposed_view(g.data(), 0, 2, 1), transpose(g.singularity()), g.symmetry(), transpose(g.indices()), g.name};
}
// ---- conjugate : always a new function -> changelog
template <typename Variable, typename Opt> gf<Variable, matrix_valued, Opt> conj(gf_view<Variable, matrix_valued, Opt> g) {
template <typename Variable, typename Singularity, typename Opt>
gf<Variable, matrix_valued, Singularity, Opt> conj(gf_view<Variable, matrix_valued, Singularity, Opt> g) {
return {g.mesh(), conj(g.data()), conj(g.singularity()), g.symmetry(), g.indices(), g.name};
}
@ -739,100 +775,67 @@ namespace gfs {
return res;
}
template <typename Variable, typename Opt, typename T>
gf<Variable, matrix_valued, Opt> operator*(gf<Variable, matrix_valued, Opt> g, matrix<T> r) {
template <typename Variable, typename Singularity, typename Opt, typename T>
gf<Variable, matrix_valued, Singularity, Opt> operator*(gf<Variable, matrix_valued, Singularity, Opt> g, matrix<T> r) {
_gf_data_mul_R(g.data(), r);
g.singularity() = g.singularity() * r;
return g;
}
template <typename Variable, typename Opt, typename T>
gf<Variable, matrix_valued, Opt> operator*(matrix<T> l, gf<Variable, matrix_valued, Opt> g) {
template <typename Variable, typename Singularity, typename Opt, typename T>
gf<Variable, matrix_valued, Singularity, Opt> operator*(matrix<T> l, gf<Variable, matrix_valued, Singularity, Opt> g) {
_gf_data_mul_L(l, g.data());
g.singularity() = l * g.singularity();
return g;
}
template <typename Variable, typename Opt, typename T>
gf<Variable, matrix_valued, Opt> L_G_R(matrix<T> l, gf<Variable, matrix_valued, Opt> g, matrix<T> r) {
auto res = gf<Variable, matrix_valued, Opt>{g.mesh(), {int(first_dim(l)), int(second_dim(r))}};
template <typename Variable, typename Singularity, typename Opt, typename T>
gf<Variable, matrix_valued, Singularity, Opt> L_G_R(matrix<T> l, gf<Variable, matrix_valued, Singularity, Opt> g, matrix<T> r) {
auto res = gf<Variable, matrix_valued, Singularity, Opt>{g.mesh(), {int(first_dim(l)), int(second_dim(r))}};
res.data() = _gf_data_mul_LR(l, g.data(), r);
res.singularity() = l * g.singularity() * r;
return res;
}
namespace gfs_implementation { // implement some default traits
// ------------------------- default factories ---------------------
// ----- tensor_valued
template <int R, typename Var, typename Opt> struct factories<Var, tensor_valued<R>, Opt> {
using gf_t = gf<Var, tensor_valued<R>, Opt>;
using target_shape_t = arrays::mini_vector<int, R>;
template <typename Var, typename Target, int VarDim, typename Singularity, typename Opt> struct factories_default_impl {
using gf_t = gf<Var, Target, Singularity, Opt>;
using target_shape_t = arrays::mini_vector<int, VarDim>;
using mesh_t = typename gf_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()));
typename gf_t::data_t A(gf_t::data_proxy_t::join_shape(m.size_of_components(), shape));
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);
return typename gf_t::singularity_t{shape};
}
};
// ----- matrix_valued
template <typename Var, typename Opt> struct factories<Var, matrix_valued, Opt> {
using gf_t = gf<Var, matrix_valued, Opt>;
using target_shape_t = arrays::mini_vector<int, 2>;
using mesh_t = typename gf_t::mesh_t;
template <int R, typename Var, typename Singularity, typename Opt>
struct factories<Var, tensor_valued<R>, Singularity, Opt> : factories_default_impl<Var, tensor_valued<R>, R, Singularity, Opt> {};
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;
}
template <typename Var, typename Singularity, typename Opt>
struct factories<Var, matrix_valued, Singularity, Opt> : factories_default_impl<Var, matrix_valued, 2, Singularity, Opt> {};
static typename gf_t::singularity_t make_singularity(mesh_t const &m, target_shape_t shape) {
return typename gf_t::singularity_t(shape);
}
};
// ----- scalar_valued
template <typename Var, typename Opt> struct factories<Var, scalar_valued, Opt> {
using gf_t = gf<Var, scalar_valued, Opt>;
struct target_shape_t {
target_shape_t front_pop() const { // this make the get_target_shape function works in this case...
return {};
}
target_shape_t() = default;
template <typename T> target_shape_t(utility::mini_vector<T, 0>) {}
};
using mesh_t = typename gf_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};
}
};
template <typename Var, typename Singularity, typename Opt>
struct factories<Var, scalar_valued, Singularity, Opt> : factories_default_impl<Var, scalar_valued, 0, Singularity, Opt> {};
// --------------------- hdf5 ---------------------------------------
// scalar function : just add a _s
template <typename Variable, typename Opt> struct h5_name<Variable, scalar_valued, Opt> {
static std::string invoke() { return h5_name<Variable, matrix_valued, Opt>::invoke() + "_s"; }
template <typename Variable, typename Singularity, typename Opt> struct h5_name<Variable, scalar_valued, Singularity, Opt> {
static std::string invoke() { return h5_name<Variable, matrix_valued, Singularity, Opt>::invoke() + "_s"; }
};
// the h5 write and read of gf members, so that we can specialize it e.g. for block gf
template <typename Variable, typename Target, typename Opt> struct h5_rw {
template <typename Variable, typename Target, typename Singularity, typename Opt> struct h5_rw {
static void write(h5::group gr, gf_const_view<Variable, Target, Opt> g) {
static void write(h5::group gr, gf_const_view<Variable, Target, Singularity, Opt> g) {
h5_write(gr, "data", g._data);
h5_write(gr, "singularity", g._singularity);
h5_write(gr, "mesh", g._mesh);
@ -841,7 +844,7 @@ namespace gfs {
// h5_write(gr, "name", g.name);
}
template <bool B> static void read(h5::group gr, gf_impl<Variable, Target, Opt, B, false> &g) {
template <bool B> static void read(h5::group gr, gf_impl<Variable, Target, Singularity, Opt, B, false> &g) {
h5_read(gr, "data", g._data);
h5_read(gr, "singularity", g._singularity);
h5_read(gr, "mesh", g._mesh);
@ -855,19 +858,20 @@ namespace gfs {
namespace mpi {
template <typename Variable, typename Target, typename Opt>
struct mpi_impl<gfs::gf<Variable, Target, Opt>, void> : mpi_impl_triqs_gfs<gfs::gf<Variable, Target, Opt>> {};
template <typename Variable, typename Target, typename Opt, bool IsConst>
struct mpi_impl<gfs::gf_view<Variable, Target, Opt, IsConst>, void> : mpi_impl_triqs_gfs<gfs::gf_view<Variable, Target, Opt, IsConst>> {};
template <typename Variable, typename Target, typename Singularity, typename Opt>
struct mpi_impl<gfs::gf<Variable, Target, Singularity, Opt>,
void> : mpi_impl_triqs_gfs<gfs::gf<Variable, Target, Singularity, Opt>> {};
template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst>
struct mpi_impl<gfs::gf_view<Variable, Target, Singularity, Opt, IsConst>,
void> : mpi_impl_triqs_gfs<gfs::gf_view<Variable, Target, Singularity, Opt, IsConst>> {};
}
}
// same as for arrays : views cannot be swapped by the std::swap. Delete it
namespace std {
template <typename Variable, typename Target, typename Opt, bool C1, bool C2>
void swap(triqs::gfs::gf_view<Variable, Target, Opt, C1> &a, triqs::gfs::gf_view<Variable, Target, Opt, C2> &b) = delete;
template <typename Variable, typename Target, typename Singularity, typename Opt, bool C1, bool C2>
void swap(triqs::gfs::gf_view<Variable, Target, Singularity, Opt, C1> &a,
triqs::gfs::gf_view<Variable, Target, Singularity, Opt, C2> &b) = delete;
}
#include "./gf_expr.hpp"

94
triqs/gfs/gf_expr.hpp
View File

@ -30,9 +30,9 @@ namespace triqs { namespace gfs {
// a wrapper for scalars
template <typename S> struct scalar_wrap {
typedef void variable_t;
typedef void target_t;
typedef void option_t;
using variable_t = void;
using target_t = void;
using option_t = void;
S s;
template <typename T> scalar_wrap(T &&x) : s(std::forward<T>(x)) {}
S singularity() const { return s; }
@ -43,85 +43,99 @@ namespace triqs { namespace gfs {
// Combine the two meshes of LHS and RHS : need to specialize where there is a scalar
struct combine_mesh {
template<typename L, typename R>
auto operator() (L && l, R && r) const -> decltype(std::forward<L>(l).mesh()) {
if (!(l.mesh() == r.mesh())) TRIQS_RUNTIME_ERROR << "Mesh mismatch : in Green Function Expression "<< l.mesh()<<" vs" <<r.mesh();
template <typename L, typename R> auto operator()(L &&l, R &&r) const -> decltype(std::forward<L>(l).mesh()) {
if (!(l.mesh() == r.mesh()))
TRIQS_RUNTIME_ERROR << "Mesh mismatch : in Green Function Expression " << l.mesh() << " vs" << r.mesh();
return std::forward<L>(l).mesh();
}
template<typename S, typename R> auto operator() (scalar_wrap<S> const &, R && r) const DECL_AND_RETURN(std::forward<R>(r).mesh());
template<typename S, typename L> auto operator() (L && l, scalar_wrap<S> const &) const DECL_AND_RETURN(std::forward<L>(l).mesh());
template <typename S, typename R>
auto operator()(scalar_wrap<S> const &, R &&r) const DECL_AND_RETURN(std::forward<R>(r).mesh());
template <typename S, typename L>
auto operator()(L &&l, scalar_wrap<S> const &) const DECL_AND_RETURN(std::forward<L>(l).mesh());
};
// Same thing to get the data shape
// NB : could be unified to one combine<F>, where F is a functor, but an easy usage requires polymorphic lambda ...
struct combine_shape {
template<typename L, typename R>
auto operator() (L && l, R && r) const -> decltype(get_gf_data_shape(std::forward<L>(l))) {
template <typename L, typename R> auto operator()(L &&l, R &&r) const -> decltype(get_gf_data_shape(std::forward<L>(l))) {
if (!(get_gf_data_shape(l) == get_gf_data_shape(r)))
TRIQS_RUNTIME_ERROR << "Shape mismatch in Green Function Expression: " << get_gf_data_shape(l) << " vs "<< get_gf_data_shape(r);
TRIQS_RUNTIME_ERROR << "Shape mismatch in Green Function Expression: " << get_gf_data_shape(l) << " vs "
<< get_gf_data_shape(r);
return get_gf_data_shape(std::forward<L>(l));
}
template<typename S, typename R> auto operator() (scalar_wrap<S> const &, R && r) const DECL_AND_RETURN(get_gf_data_shape(std::forward<R>(r)));
template<typename S, typename L> auto operator() (L && l, scalar_wrap<S> const &) const DECL_AND_RETURN(get_gf_data_shape(std::forward<L>(l)));
template <typename S, typename R>
auto operator()(scalar_wrap<S> const &, R &&r) const DECL_AND_RETURN(get_gf_data_shape(std::forward<R>(r)));
template <typename S, typename L>
auto operator()(L &&l, scalar_wrap<S> const &) const DECL_AND_RETURN(get_gf_data_shape(std::forward<L>(l)));
};
template<typename T> struct node_t : std::conditional<utility::is_in_ZRC<T>::value, scalar_wrap<T>, typename remove_rvalue_ref<T>::type> {};
template <typename T> using node_t = std14::conditional_t<utility::is_in_ZRC<T>::value, scalar_wrap<T>, typename remove_rvalue_ref<T>::type>;
template <typename A, typename B> struct _or_ {typedef void type;};
template <typename A> struct _or_<A,A> {typedef A type;};
template <typename A> struct _or_<void,A> {typedef A type;};
template <typename A> struct _or_<A,void> {typedef A type;};
template <> struct _or_<void,void> {typedef void type;};
template <typename A, typename B> struct _or_ {using type=void ;};
template <typename A> struct _or_<A,A> {using type=A ;};
template <typename A> struct _or_<void,A> {using type=A ;};
template <typename A> struct _or_<A,void> {using type=A ;};
template <> struct _or_<void,void> {using type=void ;};
}// gfs_expr_tools
template <typename Tag, typename L, typename R> struct gf_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction) {
typedef typename std::remove_reference<L>::type L_t;
typedef typename std::remove_reference<R>::type R_t;
typedef typename gfs_expr_tools::_or_<typename L_t::variable_t,typename R_t::variable_t>::type variable_t;
typedef typename gfs_expr_tools::_or_<typename L_t::target_t,typename R_t::target_t>::type target_t;
typedef typename gfs_expr_tools::_or_<typename L_t::option_t,typename R_t::option_t>::type option_t;
using L_t = typename std::remove_reference<L>::type;
using R_t = typename std::remove_reference<R>::type;
using variable_t = typename gfs_expr_tools::_or_<typename L_t::variable_t, typename R_t::variable_t>::type;
using target_t = typename gfs_expr_tools::_or_<typename L_t::target_t, typename R_t::target_t>::type;
using option_t = typename gfs_expr_tools::_or_<typename L_t::option_t, typename R_t::option_t>::type;
static_assert(!std::is_same<variable_t, void>::value, "Cannot combine two gf expressions with different variables");
static_assert(!std::is_same<target_t, void>::value, "Cannot combine two gf expressions with different target");
L l; R r;
L l;
R r;
template <typename LL, typename RR> gf_expr(LL &&l_, RR &&r_) : l(std::forward<LL>(l_)), r(std::forward<RR>(r_)) {}
auto mesh() const DECL_AND_RETURN(gfs_expr_tools::combine_mesh()(l, r));
auto singularity() const DECL_AND_RETURN(utility::operation<Tag>()(l.singularity(), r.singularity()));
auto get_data_shape() const DECL_AND_RETURN (gfs_expr_tools::combine_shape()(l,r));
template<typename KeyType> auto operator[](KeyType && key) const DECL_AND_RETURN(utility::operation<Tag>()(l[std::forward<KeyType>(key)] , r[std::forward<KeyType>(key)]));
template<typename ... Args> auto operator()(Args && ... args) const DECL_AND_RETURN(utility::operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...)));
friend std::ostream &operator <<(std::ostream &sout, gf_expr const &expr){return sout << "("<<expr.l << " "<<utility::operation<Tag>::name << " "<<expr.r<<")" ; }
template <typename KeyType>
auto operator[](KeyType &&key) const
DECL_AND_RETURN(utility::operation<Tag>()(l[std::forward<KeyType>(key)], r[std::forward<KeyType>(key)]));
template <typename... Args>
auto operator()(Args &&... args) const
DECL_AND_RETURN(utility::operation<Tag>()(l(std::forward<Args>(args)...), r(std::forward<Args>(args)...)));
friend std::ostream &operator<<(std::ostream &sout, gf_expr const &expr) {
return sout << "(" << expr.l << " " << utility::operation<Tag>::name << " " << expr.r << ")";
}
};
template <typename Tag, typename L, typename R>
auto get_gf_data_shape(gf_expr<Tag, L, R> const &g) RETURN(gfs_expr_tools::combine_shape()(g.l, g.r));
// -------------------------------------------------------------------
// a special case : the unary operator !
template <typename L> struct gf_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction) {
typedef typename std::remove_reference<L>::type L_t;
typedef typename L_t::variable_t variable_t;
typedef typename L_t::target_t target_t;
typedef typename L_t::option_t option_t;
using L_t = typename std::remove_reference<L>::type;
using variable_t = typename L_t::variable_t;
using target_t = typename L_t::target_t;
using option_t = typename L_t::option_t;
L l;
template<typename LL> gf_unary_m_expr(LL && l_) : l(std::forward<LL>(l_)) {}
auto mesh() const DECL_AND_RETURN(l.mesh());
auto singularity() const DECL_AND_RETURN(l.singularity());
AUTO_DECL get_data_shape() const RETURN (get_gf_data_shape(l));
template<typename KeyType> auto operator[](KeyType&& key) const DECL_AND_RETURN( -l[key]);
template<typename ... Args> auto operator()(Args && ... args) const DECL_AND_RETURN( -l(std::forward<Args>(args)...));
friend std::ostream &operator <<(std::ostream &sout, gf_unary_m_expr const &expr){return sout << '-'<<expr.l; }
};
template <typename L> AUTO_DECL get_gf_data_shape(gf_unary_m_expr<L> const &g) RETURN(get_gf_data_shape(g.l));
// -------------------------------------------------------------------
// Now we can define all the C++ operators ...
#define DEFINE_OPERATOR(TAG, OP, TRAIT1, TRAIT2) \
template<typename A1, typename A2>\
typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, \
gf_expr<utility::tags::TAG, typename gfs_expr_tools::node_t<A1>::type, typename gfs_expr_tools::node_t<A2>::type>>::type\
std14::enable_if_t<TRAIT1<A1>::value && TRAIT2 <A2>::value, \
gf_expr<utility::tags::TAG, gfs_expr_tools::node_t<A1>, gfs_expr_tools::node_t<A2>>>\
operator OP (A1 && a1, A2 && a2) { return {std::forward<A1>(a1),std::forward<A2>(a2)};}
DEFINE_OPERATOR(plus, +, ImmutableGreenFunction,ImmutableGreenFunction);
@ -136,11 +150,9 @@ namespace triqs { namespace gfs {
// the unary is special
template <typename A1>
typename std::enable_if<
ImmutableGreenFunction<A1>::value,
gf_unary_m_expr<typename gfs_expr_tools::node_t<A1>::type >
>::type
operator - (A1 && a1) { return {std::forward<A1>(a1)};}
std14::enable_if_t<ImmutableGreenFunction<A1>::value, gf_unary_m_expr<gfs_expr_tools::node_t<A1>>> operator-(A1 &&a1) {
return {std::forward<A1>(a1)};
}
// Now the inplace operator. Because of expression template, there are useless for speed
// we implement them trivially.

136
triqs/gfs/imfreq.hpp
View File

@ -35,21 +35,27 @@ namespace gfs {
// using matsubara_freq_mesh::matsubara_freq_mesh;
};
// singularity
template <> struct gf_default_singularity<imfreq, matrix_valued> {
using type = local::tail;
};
template <> struct gf_default_singularity<imfreq, scalar_valued> {
using type = local::tail;
};
namespace gfs_implementation {
// singularity
template <> struct singularity<imfreq, matrix_valued, void> {
using type = local::tail;
};
template <> struct singularity<imfreq, scalar_valued, void> {
using type = local::tail;
};
/// --------------------------- hdf5 ---------------------------------
// h5 name
template <typename Opt> struct h5_name<imfreq, matrix_valued, Opt> {
template <typename S, typename Opt> struct h5_name<imfreq, matrix_valued, S, Opt> {
static std::string invoke() { return "ImFreq"; }
};
/// --------------------------- data access ---------------------------------
template <typename Opt> struct data_proxy<imfreq, matrix_valued, Opt> : data_proxy_array<std::complex<double>, 3> {};
template <typename Opt> struct data_proxy<imfreq, scalar_valued, Opt> : data_proxy_array<std::complex<double>, 1> {};
/// --------------------------- evaluator ---------------------------------
// simple evaluation : take the point on the grid...
@ -67,104 +73,90 @@ namespace gfs {
// ------------- evaluator -------------------
// handle the case where the matsu. freq is out of grid...
template <typename Target, typename Opt> struct evaluator<imfreq, Target, Opt> {
static constexpr int arity = 1;
private:
struct _eval_imfreq_base_impl {
static constexpr int arity = 1;
template <typename G> int sh(G const * g) const { return (g->mesh().domain().statistic == Fermion ? 1 : 0);}
// dispatch for 2x2 cases : matrix/scalar and tail/no_tail ( true means no_tail)
// int -> replace by matsubara_freq
template <typename G>
std::complex<double> _call_impl(G const *g, matsubara_freq const &f, scalar_valued, std::false_type) const {
AUTO_DECL operator()(G const *g, int n) const
RETURN((*g)(matsubara_freq(n, g->mesh().domain().beta, g->mesh().domain().statistic)));
template <typename G> typename G::singularity_t const &operator()(G const *g, freq_infty const &) const {
return g->singularity();
}
};
// --- various 4 specializations
// scalar_valued, tail
template <typename Opt> struct evaluator<imfreq, scalar_valued, local::tail, Opt> : _eval_imfreq_base_impl {
using _eval_imfreq_base_impl::operator();
template <typename G> std::complex<double> operator()(G const *g, matsubara_freq const &f) const {
if (g->mesh().positive_only()) { // only positive Matsubara frequencies
if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n];
if ((f.n < 0) && ((-f.n-sh(g)) < g->mesh().size())) return conj((*g)[-f.n-sh(g)]);
}
else{
if ((f.n < 0) && ((-f.n - this->sh(g)) < g->mesh().size())) return conj((*g)[-f.n - this->sh(g)]);
} else {
if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size() + g->mesh().first_index())) return (*g)[f.n];
}
return g->singularity().evaluate(f)(0, 0);
}
};
template <typename G>
std::complex<double> _call_impl(G const *g, matsubara_freq const &f, scalar_valued, std::true_type) const {
// scalar_valued, no tail
template <typename Opt> struct evaluator<imfreq, scalar_valued, nothing, Opt> : _eval_imfreq_base_impl {
using _eval_imfreq_base_impl::operator();
template <typename G> std::complex<double> operator()(G const *g, matsubara_freq const &f) const {
if (g->mesh().positive_only()) { // only positive Matsubara frequencies
if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n];
if ((f.n < 0) && ((-f.n-sh(g)) < g->mesh().size())) return conj((*g)[-f.n-sh(g)]);
}
else{
if ((f.n < 0) && ((-f.n - this->sh(g)) < g->mesh().size())) return conj((*g)[-f.n - this->sh(g)]);
} else {
if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size() + g->mesh().first_index())) return (*g)[f.n];
}
return 0;
}
};
template <typename G>
arrays::matrix_const_view<std::complex<double>> _call_impl(G const *g, matsubara_freq const &f, matrix_valued,
std::false_type) const {
// matrix_valued, tail
template <typename Opt> struct evaluator<imfreq, matrix_valued, local::tail, Opt> : _eval_imfreq_base_impl {
using _eval_imfreq_base_impl::operator();
template <typename G> arrays::matrix_const_view<std::complex<double>> operator()(G const *g, matsubara_freq const &f) const {
if (g->mesh().positive_only()) { // only positive Matsubara frequencies
if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n]();
if ((f.n < 0) && ((-f.n-sh(g)) < g->mesh().size()))
return arrays::matrix<std::complex<double>>{conj((*g)[-f.n-sh(g)]())};
}
else{
if ((f.n < 0) && ((-f.n - this->sh(g)) < g->mesh().size()))
return arrays::matrix<std::complex<double>>{conj((*g)[-f.n - this->sh(g)]())};
} else {
if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size() + g->mesh().first_index())) return (*g)[f.n];
}
return g->singularity().evaluate(f);
}
};
template <typename G>
arrays::matrix_const_view<std::complex<double>> _call_impl(G const *g, matsubara_freq const &f, matrix_valued,
std::true_type) const {
// matrix_valued, no tail
template <typename Opt> struct evaluator<imfreq, matrix_valued, nothing, Opt> : _eval_imfreq_base_impl {
using _eval_imfreq_base_impl::operator();
template <typename G> arrays::matrix_const_view<std::complex<double>> operator()(G const *g, matsubara_freq const &f) const {
if (g->mesh().positive_only()) { // only positive Matsubara frequencies
if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n]();
if ((f.n < 0) && ((-f.n-sh(g)) < g->mesh().size()))
return arrays::matrix<std::complex<double>>{conj((*g)[-f.n-sh(g)]())};
}
else{
if ((f.n < 0) && ((-f.n - this->sh(g)) < g->mesh().size()))
return arrays::matrix<std::complex<double>>{conj((*g)[-f.n - this->sh(g)]())};
} else {
if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size() + g->mesh().first_index())) return (*g)[f.n];
}
auto r = arrays::matrix<std::complex<double>>{get_target_shape(*g)};
r() = 0;
return r;
}
// does not work on gcc 4.8.1 ???
/* template <typename G>
auto operator()(G const *g, matsubara_freq const &f) const
DECL_AND_RETURN(_call_impl(g, f, Target{}, std::integral_constant<bool, std::is_same<Opt, no_tail>::value>{}));
*/
public:
template <typename G>
typename std::conditional<std::is_same<Target, matrix_valued>::value, arrays::matrix_const_view<std::complex<double>>,
std::complex<double>>::type
operator()(G const *g, matsubara_freq const &f) const {
return _call_impl(g, f, Target{}, std::integral_constant<bool, std::is_same<Opt, no_tail>::value>{});
}
// int -> replace by matsubara_freq
template <typename G> auto operator()(G const *g, int n) const DECL_AND_RETURN((*g)(matsubara_freq(n,g->mesh().domain().beta,g->mesh().domain().statistic)));
#ifdef __clang__
// to generate a clearer error message ? . Only ok on clang ?
template <int n> struct error {
static_assert(n > 0, "Green function cannot be evaluated on a complex number !");
};
template <typename G> error<0> operator()(G const *g, std::complex<double>) const {
return {};
}
#endif
template <typename G> typename G::singularity_t const &operator()(G const *g, freq_infty const &) const {
return g->singularity();
}
};
/// --------------------------- data access ---------------------------------
template <typename Opt> struct data_proxy<imfreq, matrix_valued, Opt> : data_proxy_array<std::complex<double>, 3> {};
template <typename Opt> struct data_proxy<imfreq, scalar_valued, Opt> : data_proxy_array<std::complex<double>, 1> {};
} // gfs_implementation
// specific operations (for legacy python code).

122
triqs/gfs/imfreq_bff.hpp Normal file
View File

@ -0,0 +1,122 @@
/*******************************************************************************
*
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
*
* Copyright (C) 2014 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/>.
*
******************************************************************************/
#pragma once
#include "./product.hpp"
#include "./meshes/matsubara_freq.hpp"
namespace triqs {
namespace gfs {
// short cut. Here only the change compare to default multi var implementation
using imfreq_bff = cartesian_product<imfreq, imfreq, imfreq>;
//struct imfreq_bff {}; // = cartesian_product<imfreq, imfreq, imfreq>;
using imfreq_mesh_3 = mesh_product<matsubara_freq_mesh, matsubara_freq_mesh, matsubara_freq_mesh>;
// The default target for this mesh
template <> struct gf_default_target<imfreq_bff> {
using type = tensor_valued<4>;
};
// reimplement a simpler constructor, which enforces that the 2 meshes are equal.
template <typename Opt>
struct gf_mesh<imfreq_bff, Opt> : imfreq_mesh_3 {
gf_mesh() = default;
gf_mesh(matsubara_freq_mesh const& m1, matsubara_freq_mesh const& m2) : imfreq_mesh_3{m1, m2, m2} {}
};
namespace gfs_implementation {
/// --------------------------- data access ---------------------------------
struct imfreq_bff_indices_mixer {
template <typename MI, typename TI> static auto invoke(MI const& m, TI const& t) {
return std::make_tuple(std::get<0>(m), std::get<1>(m), std::get<0>(t), std::get<1>(t), std::get<2>(m), std::get<2>(t),
std::get<3>(t));
}
};
template <>
struct data_proxy<imfreq_bff, tensor_valued<4>, void> : data_proxy_array_index_mixer<std::complex<double>, 3, 4,
imfreq_bff_indices_mixer> {
//template <typename S, typename Tu> auto operator()(S& data, Tu const& tu) const {
// return data(std::get<0>(tu), arrays::range(), std::get<1>(tu), arrays::range(), std::get<2>(tu), arrays::range());
// }
};
// ------------------------------- evaluator --------------------------------------------------
//template <typename Target, typename Opt>
//struct evaluator<imfreq_bff, Target, nothing, Opt> : evaluator<cartesian_product<imfreq, imfreq, imfreq>, Target, nothing, Opt> {};
/// --------------------------- hdf5 ---------------------------------
template <typename Opt> struct h5_name<imfreq_bff, tensor_valued<4>, nothing, Opt> {
static std::string invoke() { return "imfreq_bff"; }
};
/// --------------------------- partial eval ---------------------------------
template <bool IsConst> struct partial_eval_impl<imfreq_x2, matrix_valued, nothing, void, IsConst> {
template <int pos, typename Omega>
gf_view<imfreq_x2, matrix_valued, nothing, void, IsConst>
partial_eval(gf_view<imfreq_bff, tensor_valued<4>, nothing, void, IsConst> g, Omega const& omega) {
static_assert(pos == 0, "EE");
auto& m = g.mesh().components();
auto av = g.data()(std::get<0>(m).index_to_linear(omega), arrays::ellipsis{});
return {{std::get<1>(m)}, av, {}, {}, {}};
}
};
} // gfs_implementation
gf_mesh<imfreq> get_bosonic_mesh(gf_view<imfreq_bff, tensor_valued<4>> g) { return std::get<0>(g.mesh().components()); }
/* template <typename Omega>
auto partial_eval(gf<imfreq_bff, tensor_valued<4>>& g, Omega const& omega) RETURN(partial_eval(g(), omega));
template <typename Omega>
auto partial_eval(gf<imfreq_bff, tensor_valued<4>> const& g, Omega const& omega) RETURN(partial_eval(g(), omega));
/// --------------------------- curry ---------------------------------
template <bool IsConst> auto curry_on_bosonic_freq(gf_view<imfreq_bff, tensor_valued<4>, nothing, void, IsConst> g) {
return make_gf_view_lambda_valued<imfreq>(std::get<0>(g.mesh().components()),
[g](auto&& x) { return partial_eval(g, x); });
}
auto curry_on_bosonic_freq(gf<imfreq_bff, tensor_valued<4>>& g) RETURN(curry_on_bosonic_freq(g()));
auto curry_on_bosonic_freq(gf<imfreq_bff, tensor_valued<4>> const& g) RETURN(curry_on_bosonic_freq(g()));
*/
}
}
/* template<typename G>
struct curry_polymorphic_lambda2 {
G g;
template <typename... I> AUTO_DECL operator()(I &&... i) const RETURN(partial_eval(g, i...));
};
*/

97
triqs/gfs/imfreq_x2.hpp Normal file
View File

@ -0,0 +1,97 @@
/*******************************************************************************
*
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
*
* Copyright (C) 2014 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/>.
*
******************************************************************************/
#pragma once
#include "./product.hpp"
#include "./meshes/matsubara_freq.hpp"
namespace triqs {
namespace gfs {
struct matrix_operations {};
// in this version, we store all frequencies, to be able to write matrix operation simply.
// memory usage is x2 too large in principle.
// short cut. Here only the change compare to default multi var implementation
using imfreq_x2 = cartesian_product<imfreq, imfreq>;
// reimplement a simpler constructor, which enforces that the 2 meshes are equal.
template <typename Opt> struct gf_mesh<imfreq_x2, Opt> : mesh_product<matsubara_freq_mesh, matsubara_freq_mesh> {
using B = mesh_product<matsubara_freq_mesh, matsubara_freq_mesh>;
gf_mesh() = default;
gf_mesh(matsubara_freq_mesh const& m) : B{m, m} {}
//gf_mesh(matsubara_freq_mesh const & m1, matsubara_freq_mesh const& m2) : B{m1,m2} {} //needed for curry. DO NOT document.
template <typename... T>
gf_mesh(T&&... x)
: B{matsubara_freq_mesh(std::forward<T>(x)...), matsubara_freq_mesh(std::forward<T>(x)...)} {}
};
// The default target for this mesh
template <> struct gf_default_target<imfreq_x2> {
using type = tensor_valued<4>;
};
namespace gfs_implementation {
/// --------------------------- data access ---------------------------------
struct imfreq_x2_indices_mixer {
template <typename MI, typename TI> static auto invoke(MI const& m, TI const& t) {
return std::make_tuple(std::get<0>(m), std::get<0>(t), std::get<1>(t), std::get<1>(m), std::get<2>(t), std::get<3>(t));
}
};
template <>
struct data_proxy<imfreq_x2, tensor_valued<4>, void> : data_proxy_array_index_mixer<std::complex<double>, 2, 4,
imfreq_x2_indices_mixer> {
//template <typename S, typename Tu> auto operator()(S& data, Tu const& tu) const {
// return data(std::get<0>(tu), arrays::range(), std::get<1>(tu), arrays::range(), std::get<2>(tu), arrays::range());
// }
};
/*
// Change the ordering of the indices in the array to allow matrix operations
template <typename Opt> struct data_proxy<imfreq_x2, tensor_valued<4>, Opt> : data_proxy_array_common<std::complex<double>, 6> {
template <typename S1, typename S2> static utility::mini_vector<int, 4> join_shape(S1 const& s1, S2 const& s2) {
return {int(s1[0]), s2[0], int(s1[1]), s2[1]};
// TODO : clean this size_t....
}
template <typename S, typename Tu>
AUTO_DECL operator()(S& data, Tu const& tu) const
RETURN(make_matrix_view(data(std::get<0>(tu), arrays::range(), std::get<1>(tu), arrays::range())));
};
*/
} // gfs_implementation
/// --------------------------- inverse ---------------------------------
// the generic inverse is fine. We only need to redo the invert_in_place.
template <typename Opt> void invert_in_place(gf_view<imfreq_x2, matrix_valued, nothing, Opt> g) {
auto arr = make_matrix_view(group_indices_view(g.data(), {0, 1, 2}, {3, 4, 5})); // a view of the array properly regrouped
arr = inverse(arr); // inverse as a big matrix (nu,n) x (nu', n')
// no singularity
}
}
}

23
triqs/gfs/imtime.hpp
View File

@ -37,18 +37,18 @@ namespace gfs {
// using matsubara_time_mesh::matsubara_time_mesh;
};
// singularity
template <> struct gf_default_singularity<imtime, matrix_valued> {
using type = local::tail;
};
template <> struct gf_default_singularity<imtime, scalar_valued> {
using type = local::tail;
};
namespace gfs_implementation {
// singularity. If no_tail is given, then it is the default (nothing)
template <> struct singularity<imtime, matrix_valued, void> {
using type = local::tail;
};
template <> struct singularity<imtime, scalar_valued, void> {
using type = local::tail;
};
// h5 name
template <typename Opt> struct h5_name<imtime, matrix_valued, Opt> {
template <typename Singularity, typename Opt> struct h5_name<imtime, matrix_valued, Singularity, Opt> {
static std::string invoke() { return "ImTime"; }
};
@ -59,7 +59,7 @@ namespace gfs {
/// --------------------------- closest mesh point on the grid ---------------------------------
template <typename Opt, typename Target> struct get_closest_point<imtime, Target, Opt> {
template <typename Singularity, typename Opt, typename Target> struct get_closest_point<imtime, Target, Singularity, Opt> {
// index_t is int
template <typename G, typename T> static int invoke(G const *g, closest_pt_wrap<T> const &p) {
double x = (g->mesh().kind() == half_bins ? double(p.value) : double(p.value) + 0.5 * g->mesh().delta());
@ -98,7 +98,8 @@ namespace gfs {
};
// now evaluator
template <typename Opt, typename Target> struct evaluator<imtime, Target, Opt> : evaluator_one_var<imtime> {};
template <typename Singularity, typename Opt, typename Target>
struct evaluator<imtime, Target, Singularity, Opt> : evaluator_one_var<imtime> {};
} // gfs_implementation.
}

18
triqs/gfs/legendre.hpp
View File

@ -26,7 +26,8 @@
#include "./domains/legendre.hpp"
#include "./meshes/discrete.hpp"
namespace triqs { namespace gfs {
namespace triqs {
namespace gfs {
struct legendre {};
@ -40,17 +41,19 @@ namespace triqs { namespace gfs {
namespace gfs_implementation {
// h5 name
template<typename Opt> struct h5_name<legendre,matrix_valued,Opt> { static std::string invoke(){ return "Legendre";}};
template <typename Opt> struct h5_name<legendre, matrix_valued, nothing, Opt> {
static std::string invoke() { return "Legendre"; }
};
/// --------------------------- evaluator ---------------------------------
// Not finished, not tested
template<typename Opt>
struct evaluator<legendre,matrix_valued,Opt> {
template <typename Opt> struct evaluator<legendre, matrix_valued, nothing, Opt> {
static constexpr int arity = 1;
// ERROR : give a double and interpolate
template<typename G>
arrays::matrix_view<double > operator() (G const * g,long n) const {return g->data()(n, arrays::range(), arrays::range()); }
template <typename G> arrays::matrix_view<double> operator()(G const* g, long n) const {
return g->data()(n, arrays::range(), arrays::range());
}
};
/// --------------------------- data access ---------------------------------
@ -59,6 +62,7 @@ namespace triqs { namespace gfs {
template <typename Opt> struct data_proxy<legendre, scalar_valued, Opt> : data_proxy_array<double, 1> {};
} // gfs_implementation
}}
}
}
#endif

4
triqs/gfs/local/fit_tail.cpp
View File

@ -25,8 +25,8 @@ namespace triqs { namespace gfs { namespace local {
int size1 = n_max - n_min + 1;
// size2 is the number of moments
arrays::matrix<double, 2> A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT);
arrays::matrix<double, 2> B(size1, 1, FORTRAN_LAYOUT);
arrays::matrix<double> A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT);
arrays::matrix<double> B(size1, 1, FORTRAN_LAYOUT);
arrays::vector<double> S(std::max(size_even, size_odd));
const double rcond = 0.0;
int rank;

32
triqs/gfs/local/fourier_matsubara.hpp
View File

@ -18,8 +18,7 @@
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef TRIQS_GF_LOCAL_FOURIER_MATSU_H
#define TRIQS_GF_LOCAL_FOURIER_MATSU_H
#pragma once
#include "fourier_base.hpp"
#include <triqs/gfs/imfreq.hpp>
@ -28,12 +27,12 @@
namespace triqs {
namespace gfs {
template <typename Target, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, imtime, Target, Opt> fourier(gf_impl<imtime, Target, Opt, V, C> const& g) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, imtime, Target, Singularity, Opt> fourier(gf_impl<imtime, Target, Singularity, Opt, V, C> const& g) {
return {g};
}
template <typename Target, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, imfreq, Target, Opt> inverse_fourier(gf_impl<imfreq, Target, Opt, V, C> const& g) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, imfreq, Target, Singularity, Opt> inverse_fourier(gf_impl<imfreq, Target, Singularity, Opt, V, C> const& g) {
return {g};
}
@ -49,30 +48,30 @@ namespace gfs {
void triqs_gf_view_assign_delegation(gf_view<imfreq, matrix_valued, no_tail> g,
gf_keeper<tags::fourier, imtime, matrix_valued, no_tail> const& L);
template <typename Target, typename Opt, bool V, bool C>
gf<imfreq, Target, Opt> make_gf_from_fourier(gf_impl<imtime, Target, Opt, V, C> const& gt, int n_iw) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf<imfreq, Target, Singularity, Opt> make_gf_from_fourier(gf_impl<imtime, Target, Singularity, Opt, V, C> const& gt, int n_iw) {
auto m = gf_mesh<imfreq, Opt>{gt.mesh().domain(), n_iw};
auto gw = gf<imfreq, Target, Opt>{m, get_target_shape(gt)};
auto gw = gf<imfreq, Target, Singularity, Opt>{m, get_target_shape(gt)};
gw() = fourier(gt);
return gw;
}
template <typename Target, typename Opt, bool V, bool C>
gf<imfreq, Target, Opt> make_gf_from_fourier(gf_impl<imtime, Target, Opt, V, C> const& gt) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf<imfreq, Target, Singularity, Opt> make_gf_from_fourier(gf_impl<imtime, Target, Singularity, Opt, V, C> const& gt) {
return make_gf_from_fourier(gt, (gt.mesh().size() - (gt.mesh().kind() == full_bins ? 1 : 0)) / 2);
}
template <typename Target, typename Opt, bool V, bool C>
gf<imtime, Target, Opt> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Opt, V, C> const& gw, int n_tau,
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf<imtime, Target, Singularity, Opt> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Singularity, Opt, V, C> const& gw, int n_tau,
mesh_kind mk = full_bins) {
auto m = gf_mesh<imtime, Opt>{gw.mesh().domain(), n_tau};
auto gt = gf<imtime, Target, Opt>{m, get_target_shape(gw)};
auto gt = gf<imtime, Target, Singularity, Opt>{m, get_target_shape(gw)};
gt() = inverse_fourier(gw);
return gt;
}
template <typename Target, typename Opt, bool V, bool C>
gf<imtime, Target, Opt> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Opt, V, C> const& gw, mesh_kind mk = full_bins) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf<imtime, Target, Singularity, Opt> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Singularity, Opt, V, C> const& gw, mesh_kind mk = full_bins) {
return make_gf_from_inverse_fourier(gw, 2 * gw.mesh().size() + (mk == full_bins ? 1 : 0), mk);
}
}
@ -82,5 +81,4 @@ namespace clef {
TRIQS_CLEF_MAKE_FNT_LAZY(inverse_fourier);
}
}
#endif

21
triqs/gfs/local/fourier_real.hpp
View File

@ -18,21 +18,19 @@
* TRIQS. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef TRIQS_GF_LOCAL_FOURIER_REAL_H
#define TRIQS_GF_LOCAL_FOURIER_REAL_H
#pragma once
#include "fourier_base.hpp"
#include <triqs/gfs/refreq.hpp>
#include <triqs/gfs/retime.hpp>
namespace triqs { namespace gfs {
template <typename Target, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, retime, Target> fourier(gf_impl<retime, Target, Opt, V, C> const& g) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, retime, Target, Singularity> fourier(gf_impl<retime, Target, Singularity, Opt, V, C> const& g) {
return {g};
}
template <typename Target, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, refreq, Target> inverse_fourier(gf_impl<refreq, Target, Opt, V, C> const& g) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf_keeper<tags::fourier, refreq, Target, Singularity> inverse_fourier(gf_impl<refreq, Target, Singularity, Opt, V, C> const& g) {
return {g};
}
@ -58,20 +56,19 @@ namespace triqs { namespace gfs {
return {tmin, tmax, L};
}
template <typename Target, typename Opt, bool V, bool C>
gf_view<refreq, Target> make_gf_from_fourier(gf_impl<retime, Target, Opt, V, C> const& gt) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf_view<refreq, Target> make_gf_from_fourier(gf_impl<retime, Target, Singularity, Opt, V, C> const& gt) {
auto gw = gf<refreq, Target>{make_mesh_fourier_compatible(gt.mesh()), get_target_shape(gt)};
gw() = fourier(gt);
return gw;
}
template <typename Target, typename Opt, bool V, bool C>
gf_view<retime, Target> make_gf_from_inverse_fourier(gf_impl<refreq, Target, Opt, V, C> const& gw) {
template <typename Target, typename Singularity, typename Opt, bool V, bool C>
gf_view<retime, Target> make_gf_from_inverse_fourier(gf_impl<refreq, Target, Singularity, Opt, V, C> const& gw) {
auto gt = gf<retime, Target>{make_mesh_fourier_compatible(gw.mesh()), get_target_shape(gw)};
gt() = inverse_fourier(gw);
return gt;
}
}}
#endif

15
triqs/gfs/local/no_tail.hpp
View File

@ -26,10 +26,11 @@
namespace triqs {
namespace gfs {
struct no_tail {};
//struct no_tail {};
using no_tail = nothing;
template <typename Variable, typename Target, bool V, bool C>
gf_view<Variable, Target, no_tail, C> make_gf_view_without_tail(gf_impl<Variable, Target, void, V, C> const &g) {
template <typename Variable, typename Target, typename S, typename Opt, bool V, bool C>
gf_view<Variable, Target, no_tail, Opt, C> make_gf_view_without_tail(gf_impl<Variable, Target, S, Opt, V, C> const &g) {
return {g.mesh(), g.data(), {}, g.symmetry(), g.indices(), g.name};
}
@ -46,15 +47,15 @@ namespace gfs {
}
}
template <typename Variable, typename Target, bool V, bool C>
gf_view<Variable, Target> make_gf_from_g_and_tail(gf_impl<Variable, Target, no_tail, V, C> const &g, local::tail t) {
template <typename Variable, typename Target, typename S, typename Opt, bool V, bool C>
gf_view<Variable, Target> make_gf_from_g_and_tail(gf_impl<Variable, Target, S, Opt, V, C> const &g, local::tail t) {
details::_equal_or_throw(t.shape(), get_target_shape(g));
auto g2 = gf<Variable, Target, no_tail>{g}; // copy the function without tail
return {std::move(g2.mesh()), std::move(g2.data()), std::move(t), g2.symmetry()};
}
template <typename Variable, typename Target, bool V, bool C>
gf_view<Variable, Target, void, C> make_gf_view_from_g_and_tail(gf_impl<Variable, Target, no_tail, V, C> const &g,
template <typename Variable, typename Target, typename S, typename Opt, bool V, bool C>
gf_view<Variable, Target, local::tail, Opt, C> make_gf_view_from_g_and_tail(gf_impl<Variable, Target, no_tail, Opt, V, C> const &g,
local::tail_view t) {
details::_equal_or_throw(t.shape(), get_target_shape(g));
return {g.mesh(), g.data(), t, g.symmetry()};

2
triqs/gfs/local/tail.hpp
View File

@ -53,6 +53,7 @@ namespace triqs { namespace gfs { namespace local {
public:
TRIQS_MPI_IMPLEMENTED_VIA_BOOST;
typedef tail_view view_type;
typedef tail_view const_view_type; // not nice
typedef tail regular_type;
typedef arrays::array <dcomplex,3> data_regular_type;
@ -232,6 +233,7 @@ namespace triqs { namespace gfs { namespace local {
typedef tqa::mini_vector<size_t,2> shape_type;
tail(size_t N1, size_t N2, size_t size_=10, long order_min=-1): B(N1,N2,size_,order_min) {}
tail(shape_type const & sh, size_t size_=10, long order_min=-1): B(sh[0],sh[1],size_,order_min) {}
tail(tqa::mini_vector<int,0>) : tail(1,1) {}
tail(B::data_type const &d, B::mask_type const &m, long order_min): B(d, m, order_min) {}
tail(tail const & g): B(g) {}
tail(tail_view const & g): B(g) {}

103
triqs/gfs/m_tail.hpp Normal file
View File

@ -0,0 +1,103 @@
/*******************************************************************************
*
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
*
* Copyright (C) 2012-2014 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/>.
*
******************************************************************************/
#pragma once
#include "./tools.hpp"
#include "./gf.hpp"
#include "./local/tail.hpp"
#include "./meshes/discrete.hpp"
namespace triqs {
namespace gfs {
namespace gfs_implementation {
/*template <typename T> constexpr bool has_a_singularity() {
return std::is_same<typename singularity<V1, Target>::type, local::tail>::value;
}
/// --------------------------- singularity ---------------------------------
template <typename V1, typename V2, typename Target, typename Opt> struct singularity<cartesian_product<V1,V2>, Target, Opt> {
using type = std14::conditional_t < (!has_a_singularity<V1>()) && has_a_singularity<V2>(), gf<V1, local::tail>, nothing > ;
};
*/
/// --------------------------- hdf5 ---------------------------------
template <typename Variable, typename Opt> struct h5_name<Variable, local::tail, nothing, Opt> {
static std::string invoke() { return "xxxxx"; }
};
template <typename Variable, typename Opt> struct h5_rw<Variable, local::tail, nothing, Opt> {
static void write(h5::group gr, gf_const_view<Variable, local::tail, nothing, Opt> g) {
for (size_t i = 0; i < g.mesh().size(); ++i) h5_write(gr, std::to_string(i), g._data[i]);
// h5_write(gr,"symmetry",g._symmetry);
}
template <bool IsView> static void read(h5::group gr, gf_impl<Variable, local::tail, nothing, Opt, IsView, false> &g) {
// does not work : need to read the block name and remake the mesh...
//g._mesh = gf_mesh<block_index, Opt>(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 < g.mesh().size(); ++i) h5_read(gr, g.mesh().domain().names()[i], g._data[i]);
// h5_read(gr,"symmetry",g._symmetry);
}
};
/// --------------------------- data access ---------------------------------
template <typename Variable, typename Opt> struct data_proxy<Variable, local::tail, Opt> : data_proxy_vector<local::tail> {};
// ------------------------------- Factories --------------------------------------------------
template <typename Variable, typename Opt> struct factories<Variable, local::tail, nothing, Opt> {
using mesh_t=gf_mesh<Variable, Opt> ;
using gf_t=gf<Variable, local::tail> ;
//using gf_view_t=gf_view<block_index, Target> ;
struct target_shape_t {};
static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t) { return std::vector<local::tail>(m.size()); }
static nothing make_singularity(mesh_t const &m, target_shape_t) {
return {};
}
};
// partial_eval
template <typename Variable, typename Opt, bool IsConst> struct partial_eval_impl<Variable, local::tail, nothing, Opt, IsConst> {
using gv_t = gf_view<Variable, local::tail, nothing, Opt, IsConst>;
template <int... pos, typename... T> static auto invoke(gv_t g, T const &... x) {
return invoke_impl(g, std14::index_sequence<pos...>(), x...);
}
template <typename T> static local::tail_view invoke_impl(gv_t g, std14::index_sequence<0>, T const &x) {
return g[g.mesh().index_to_linear(x)];
}
template <typename T> static nothing invoke_impl(gv_t g, std14::index_sequence<1>, T const &x) {
return {};
}
};
} // gfs_implementation
}}

8
triqs/gfs/meshes/discrete.hpp
View File

@ -29,12 +29,18 @@ namespace gfs {
using domain_t = Domain;
using index_t = long;
using linear_index_t = long;
discrete_mesh(domain_t dom) : _dom(std::move(dom)) {}
discrete_mesh() = default;
domain_t const &domain() const { return _dom; }
long size() const { return _dom.size(); }
size_t size() const { return _dom.size(); }
///
utility::mini_vector<size_t, 1> size_of_components() const {
return {size()};
}
/// Conversions point <-> index <-> discrete_index
long index_to_point(index_t ind) const { return ind; }

8
triqs/gfs/meshes/linear.hpp
View File

@ -37,6 +37,7 @@ namespace gfs {
using domain_t = Domain;
using index_t = long;
using linear_index_t = long;
using domain_pt_t = typename domain_t::point_t;
static_assert(!std::is_base_of<std::complex<double>, domain_pt_t>::value,
@ -64,7 +65,12 @@ namespace gfs {
}
domain_t const &domain() const { return _dom; }
long size() const { return L; }
size_t size() const { return L; }
utility::mini_vector<size_t, 1> size_of_components() const {
return {size()};
}
double delta() const { return del; }
double x_max() const { return xmax; }
double x_min() const { return xmin; }

7
triqs/gfs/meshes/matsubara_freq.hpp
View File

@ -21,6 +21,7 @@
#pragma once
#include "./mesh_tools.hpp"
#include "../domains/matsubara.hpp"
#include <triqs/mpi.hpp>
namespace triqs {
namespace gfs {
@ -31,6 +32,7 @@ namespace gfs {
using domain_t = matsubara_domain<true>;
///
using index_t = long;
using linear_index_t = long;
///
using domain_pt_t = typename domain_t::point_t;
@ -98,6 +100,11 @@ namespace gfs {
/// Size (linear) of the mesh of the window
long size() const { return _last_index_window - _first_index_window + 1; }
///
utility::mini_vector<size_t, 1> size_of_components() const {
return {size_t(size())};
}
/// From an index of a point in the mesh, returns the corresponding point in the domain
domain_pt_t index_to_point(index_t ind) const { return 1_j * M_PI * (2 * ind + (_dom.statistic == Fermion)) / _dom.beta; }

53
triqs/gfs/meshes/product.c14.hpp
View File

@ -2,7 +2,7 @@
*
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
*
* Copyright (C) 2012-2013 by O. Parcollet
* Copyright (C) 2012-2014 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
@ -34,6 +34,7 @@ namespace gfs {
using m_tuple_t = std::tuple<Meshes...>;
using m_pt_tuple_t = std::tuple<typename Meshes::mesh_point_t...>;
using domain_pt_t = typename domain_t::point_t;
using linear_index_t = std::tuple<typename Meshes::linear_index_t...>;
static constexpr int dim = sizeof...(Meshes);
@ -49,7 +50,7 @@ namespace gfs {
size_t size() const {
return triqs::tuple::fold([](auto const &m, size_t R) { return R * m.size(); }, m_tuple, 1);
}
/*
/// Scatter the first mesh over the communicator c
friend mesh_product mpi_scatter(mesh_product const &m, mpi::communicator c, int root) {
auto r = m; // same domain, but mesh with a window. Ok ?
@ -58,10 +59,18 @@ namespace gfs {
}
/// Opposite of scatter : rebuild the original mesh, without a window
friend matsubara_freq_mesh mpi_gather(matsubara_freq_mesh m, mpi::communicator c, int root) {
friend mesh_product mpi_gather(mesh_product m, mpi::communicator c, int root) {
auto r = m; // same domain, but mesh with a window. Ok ?
std::get<0>(r.m_tuple) = mpi_gather(std::get<0>(r.m_tuple), c, root);
return r;
*/
/// The sizes of all mesh components
utility::mini_vector<size_t, dim> size_of_components() const {
utility::mini_vector<size_t, dim> res;
auto l = [&res](int i, auto const &m) mutable { res[i] = m.size(); };
triqs::tuple::for_each_enumerate(m_tuple, l);
return res;
}
/// Conversions point <-> index <-> linear_index
@ -72,32 +81,23 @@ namespace gfs {
return res;
}
/// Flattening index to linear : index[0] + component[0].size * (index[1] + component[1].size* (index[2] + ....))
size_t index_to_linear(index_t const &ii) const {
auto l = [](auto const &m, auto const &i, size_t R) { return m.index_to_linear(i) + R * m.size(); };
return triqs::tuple::fold(l, reverse(m_tuple), reverse(ii), size_t(0));
/// The linear_index is the tuple of the linear_index of the components
linear_index_t index_to_linear(index_t const &ind) const {
auto l = [](auto const &m, auto const &i) { return m.index_to_linear(i); };
return triqs::tuple::map_on_zip(l, m_tuple, ind);
}
/// Flattening index to linear : index[0] + component[0].size * (index[1] + component[1].size* (index[2] + ....))
size_t mp_to_linear(m_pt_tuple_t const &mp) const {
auto l = [](auto const &m, auto const &p, size_t R) { return p.linear_index() + R * m.size(); };
return triqs::tuple::fold(l, reverse(m_tuple), reverse(mp), size_t(0));
}
utility::mini_vector<size_t, dim> shape() const {
utility::mini_vector<size_t, dim> res;
auto l = [&res](int i, auto const &m) mutable { res[i] = m.size(); };
triqs::tuple::for_each_enumerate(m_tuple, l);
return res;
///
linear_index_t mp_to_linear(m_pt_tuple_t const &mp) const {
auto l = [](auto const &p) { return p.linear_index(); };
return triqs::tuple::map(l, mp);
}
// Same but a variadic list of mesh_point_t
template <typename... MP> size_t mesh_pt_components_to_linear(MP const &... mp) const {
static_assert(std::is_same<std::tuple<MP...>, m_pt_tuple_t>::value, "Call incorrect ");
// static_assert(std::is_same< std::tuple<typename std::remove_cv<typename std::remove_reference<MP>::type>::type...>,
// m_pt_tuple_t>::value, "Call incorrect ");
return mp_to_linear(std::forward_as_tuple(mp...));
} // speed test ? or make a variadic fold...
}
/// The wrapper for the mesh point
class mesh_point_t : tag::mesh_point {
@ -116,7 +116,7 @@ namespace gfs {
, _atend(false) {}
mesh_point_t(mesh_product const &m_) : m(&m_), _c(triqs::tuple::map(F1(), m_.m_tuple)), _atend(false) {}
m_pt_tuple_t const &components_tuple() const { return _c; }
size_t linear_index() const { return m->mp_to_linear(_c); }
linear_index_t linear_index() const { return m->mp_to_linear(_c); }
const mesh_product *mesh() const { return m; }
using cast_t = domain_pt_t;
@ -193,14 +193,5 @@ namespace gfs {
template <int pos, typename P> decltype(auto) get_component(P const &p) { return std::get<pos>(p.components_tuple()); }
// Given a composite mesh m , and a linear array of storage A
// 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 <typename... Meshes, typename T, ull_t OptionsFlags, ull_t To, int R, bool B, bool C>
arrays::array_view<T, sizeof...(Meshes) + R - 1, OptionsFlags, To, true, C>
reinterpret_linear_array(mesh_product<Meshes...> const &m, arrays::array_view<T, R, OptionsFlags, To, B, C> A) {
return {{join(m.shape(), get_shape(A).front_pop())}, A.storage()};
}
}
}

73
triqs/gfs/product.hpp
View File

@ -35,6 +35,7 @@ namespace gfs {
// use alias
template <typename... Ms> struct cartesian_product<std::tuple<Ms...>> : cartesian_product<Ms...> {};
/// TODO : Put inheriting constructor, simpler...
// the mesh is simply a cartesian product
template <typename Opt, typename... Ms> struct gf_mesh<cartesian_product<Ms...>, Opt> : mesh_product<gf_mesh<Ms, Opt>...> {
// using mesh_product< gf_mesh<Ms,Opt> ... >::mesh_product< gf_mesh<Ms,Opt> ... > ;
@ -45,56 +46,36 @@ namespace gfs {
namespace gfs_implementation {
/// --------------------------- hdf5 ---------------------------------
// h5 name : name1_x_name2_.....
template <typename Opt, typename... Ms> struct h5_name<cartesian_product<Ms...>, matrix_valued, Opt> {
static std::string invoke() {
return triqs::tuple::fold([](std::string a, std::string b) { return a + std::string(b.empty() ? "" : "_x_") + b; },
std::make_tuple(h5_name<Ms, matrix_valued, Opt>::invoke()...), std::string());
}
};
template <typename Opt, int R, typename... Ms>
struct h5_name<cartesian_product<Ms...>, tensor_valued<R>, Opt> : h5_name<cartesian_product<Ms...>, matrix_valued, Opt> {};
// a slight difference with the generic case : reinterpret the data array to avoid flattening the variables
template <typename Opt, int R, typename... Ms> struct h5_rw<cartesian_product<Ms...>, tensor_valued<R>, Opt> {
using g_t = gf<cartesian_product<Ms...>, tensor_valued<R>, Opt>;
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);
}
template <bool IsView>
static void read(h5::group gr, gf_impl<cartesian_product<Ms...>, tensor_valued<R>, Opt, IsView, false> &g) {
using G_t = gf_impl<cartesian_product<Ms...>, tensor_valued<R>, Opt, IsView, false>;
h5_read(gr, "mesh", g._mesh);
auto arr = arrays::array<typename G_t::data_t::value_type, sizeof...(Ms) + R>{};
h5_read(gr, "data", arr);
auto sh = arr.shape();
arrays::mini_vector<size_t, R + 1> sh2;
sh2[0] = g._mesh.size();
for (int u = 1; u < R + 1; ++u) sh2[u] = sh[sizeof...(Ms) - 1 + u];
g._data = arrays::array<typename G_t::data_t::value_type, R + 1>{sh2, std::move(arr.storage())};
h5_read(gr, "singularity", g._singularity);
h5_read(gr, "symmetry", g._symmetry);
}
};
/// --------------------------- data access ---------------------------------
template <typename Opt, typename... Ms>
struct data_proxy<cartesian_product<Ms...>, scalar_valued, Opt> : data_proxy_array<std::complex<double>, 1> {};
template <typename Opt, typename... Ms>
struct data_proxy<cartesian_product<Ms...>, matrix_valued, Opt> : data_proxy_array<std::complex<double>, 3> {};
template <int R, typename Opt, typename... Ms>
struct data_proxy<cartesian_product<Ms...>, tensor_valued<R>, Opt> : data_proxy_array<std::complex<double>, R + 1> {};
struct data_proxy<cartesian_product<Ms...>, scalar_valued, Opt> : data_proxy_array_multivar<std::complex<double>,
sizeof...(Ms)> {};
template <typename Opt, typename... Ms>
struct data_proxy<cartesian_product<Ms...>, matrix_valued, Opt> : data_proxy_array_multivar_matrix_valued<std::complex<double>,
2 + sizeof...(Ms)> {};
template <int R, typename Opt, typename... Ms>
struct data_proxy<cartesian_product<Ms...>, tensor_valued<R>, Opt> : data_proxy_array_multivar<std::complex<double>,
R + sizeof...(Ms)> {};
// special case ? Or make a specific container....
template <typename Opt, typename M0>
struct data_proxy<cartesian_product<M0,imtime>, matrix_valued, Opt> : data_proxy_array<double, 3> {};
struct data_proxy<cartesian_product<M0, imtime>, matrix_valued, Opt> : data_proxy_array_multivar_matrix_valued<double, 2 + 2> {
};
/// --------------------------- hdf5 ---------------------------------
// h5 name : name1_x_name2_.....
template <typename Opt, typename... Ms> struct h5_name<cartesian_product<Ms...>, matrix_valued, nothing, Opt> {
static std::string invoke() {
return triqs::tuple::fold([](std::string a, std::string b) { return a + std::string(b.empty() ? "" : "_x_") + b; },
std::make_tuple(h5_name<Ms, matrix_valued, nothing, Opt>::invoke()...), std::string());
}
};
template <typename Opt, int R, typename... Ms>
struct h5_name<cartesian_product<Ms...>, tensor_valued<R>, nothing, Opt> : h5_name<cartesian_product<Ms...>, matrix_valued, nothing, Opt> {};
/// --------------------------- evaluator ---------------------------------
@ -146,7 +127,7 @@ namespace gfs {
};
// now the multi d evaluator itself.
template <typename Target, typename Opt, typename... Ms> struct evaluator<cartesian_product<Ms...>, Target, Opt> {
template <typename Target, typename Opt, typename... Ms> struct evaluator<cartesian_product<Ms...>, Target, nothing, Opt> {
static constexpr int arity = sizeof...(Ms);
mutable std::tuple<evaluator_fnt_on_mesh<Ms>...> evals;

20
triqs/gfs/refreq.hpp
View File

@ -35,18 +35,18 @@ namespace gfs {
//using segment_mesh::segment_mesh;
};
// singularity
template <> struct gf_default_singularity<refreq, matrix_valued> {
using type = local::tail;
};
template <> struct gf_default_singularity<refreq, scalar_valued> {
using type = local::tail;
};
namespace gfs_implementation {
// singularity
template <typename Opt> struct singularity<refreq, matrix_valued, Opt> {
using type = local::tail;
};
template <typename Opt> struct singularity<refreq, scalar_valued, Opt> {
using type = local::tail;
};
// h5 name
template <typename Opt> struct h5_name<refreq, matrix_valued, Opt> {
template <typename Singularity, typename Opt> struct h5_name<refreq, matrix_valued, Singularity, Opt> {
static std::string invoke() { return "ReFreq"; }
};
@ -56,7 +56,7 @@ namespace gfs {
struct evaluator_fnt_on_mesh<refreq> TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh,
evaluator_grid_linear_interpolation);
template <typename Opt, typename Target> struct evaluator<refreq, Target, Opt> : evaluator_one_var<refreq> {};
template <typename Singularity, typename Opt, typename Target> struct evaluator<refreq, Target, Singularity, Opt> : evaluator_one_var<refreq> {};
/// --------------------------- data access ---------------------------------

20
triqs/gfs/retime.hpp
View File

@ -35,18 +35,18 @@ namespace gfs {
//using segment_mesh::segment_mesh;
};
// singularity
template <> struct gf_default_singularity<retime, matrix_valued> {
using type = local::tail;
};
template <> struct gf_default_singularity<retime, scalar_valued> {
using type = local::tail;
};
namespace gfs_implementation {
// singularity
template <typename Opt> struct singularity<retime, matrix_valued, Opt> {
using type = local::tail;
};
template <typename Opt> struct singularity<retime, scalar_valued, Opt> {
using type = local::tail;
};
// h5 name
template <typename Opt> struct h5_name<retime, matrix_valued, Opt> {
template <typename Singularity, typename Opt> struct h5_name<retime, matrix_valued, Singularity, Opt> {
static std::string invoke() { return "ReTime"; }
};
@ -55,7 +55,7 @@ namespace gfs {
struct evaluator_fnt_on_mesh<retime> TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh,
evaluator_grid_linear_interpolation);
template <typename Opt, typename Target> struct evaluator<retime, Target, Opt> : evaluator_one_var<retime> {};
template <typename Singularity, typename Opt, typename Target> struct evaluator<retime, Target, Singularity, Opt> : evaluator_one_var<retime> {};
/// --------------------------- data access ---------------------------------
template <typename Opt> struct data_proxy<retime, matrix_valued, Opt> : data_proxy_array<std::complex<double>, 3> {};

1
triqs/gfs/tools.hpp
View File

@ -186,6 +186,7 @@ namespace gfs {
bool is_empty() const { return false;}
};
template<int ... pos, typename ...T> nothing partial_eval(nothing, T&&...) { return {};}
inline nothing transpose(nothing) { return {};}
inline nothing inverse(nothing) { return {};}
inline nothing conj(nothing) { return {};}

8
triqs/lattice/bz_mesh.hpp
View File

@ -30,6 +30,7 @@ namespace lattice {
using domain_t = brillouin_zone;
using index_t = long;
using linear_index_t = long;
using domain_pt_t = typename domain_t::point_t;
bz_mesh() = default;
@ -37,7 +38,12 @@ namespace lattice {
bz_mesh(brillouin_zone const &bz, std::vector<k_t> k_pt_stack) : bz(bz), k_pt_stack(std::move(k_pt_stack)) {}
domain_t const &domain() const { return bz; }
long size() const { return k_pt_stack.size(); }
size_t size() const { return k_pt_stack.size(); }
///
utility::mini_vector<size_t, 1> size_of_components() const {
return {size()};
}
/// Conversions point <-> index <-> linear_index
domain_pt_t const &index_to_point(index_t i) const { return k_pt_stack[i]; }

52
triqs/utility/mini_vector.hpp
View File

@ -38,7 +38,7 @@ namespace triqs { namespace utility {
template <typename T, int Rank>
class mini_vector {
T _data[Rank];
T _data[(Rank!=0 ? Rank : 1)];
friend class boost::serialization::access;
template<class Archive> void serialize(Archive & ar, const unsigned int version) { ar & TRIQS_MAKE_NVP("_data",_data); }
void init() { for (int i=0;i<Rank; ++i) _data[i] = 0;}
@ -126,6 +126,16 @@ namespace triqs { namespace utility {
friend std::stringstream & operator << ( std::stringstream & out, mini_vector const & v ) { out<<v.to_string(); return out;}
}; // class mini_vector
// ------------ specialization for size 0 -------------
template <typename T> class mini_vector<T, 0> {
T _data[1];
public:
T & operator[](size_t i) { return _data[i];}
const T & operator[](size_t i) const { return _data[i];}
};
// ------------- Comparison --------------------------------------
template <typename T, int R>
bool operator ==(mini_vector<T,R> const & v1, mini_vector<T,R> const & v2) {
for (int i=0;i<R; ++i) { if (v1[i]!=v2[i]) return false;}
@ -134,14 +144,16 @@ namespace triqs { namespace utility {
template <typename T, int R> bool operator !=(mini_vector<T,R> const & v1, mini_vector<T,R> const & v2) { return (!(v1==v2));}
template <typename T, int R1, int R2>
mini_vector<T, R1+R2> join (mini_vector<T,R1> const & v1, mini_vector<T,R2> const & v2) {
mini_vector<T, R1+R2> res;
// ------------- join --------------------------------------
template <typename T1, typename T2, int R1, int R2>
mini_vector<T1, R1+R2> join (mini_vector<T1,R1> const & v1, mini_vector<T2,R2> const & v2) {
mini_vector<T1, R1+R2> res;
for (int i=0;i<R1; ++i) res[i]=v1[i];
for (int i=0;i<R2; ++i) res[R1+i]=v2[i];
return res;
}
// ------------- dot --------------------------------------
template <typename T1, typename T2, int Rank>
T1 dot_product(mini_vector<T1,Rank> const & v1, mini_vector<T2,Rank> const & v2) {
T1 res=0;
@ -149,25 +161,21 @@ namespace triqs { namespace utility {
return res;
}
struct tuple_to_mini_vector_aux { template<typename M, typename V> V * operator()(M const & m, V * v) { *v = m; return ++v;}};
#ifndef TRIQS_C11
// ------------- transform a tuple into a minivector --------------------------------------
// change : the first version crash clang 3.3, but not svn version.
// must be a bug, corrected since then
/*
template<typename T, typename ... U>
mini_vector<T,sizeof...(U)> tuple_to_mini_vector(std::tuple<U...> const & t) {
mini_vector<T,sizeof...(U)> res;
triqs::tuple::fold(tuple_to_mini_vector_aux(),t,&res[0]);
return res;
template <typename T, typename TU> mini_vector<T, std::tuple_size<TU>::value> tuple_to_mini_vector(TU const &t) {
return tuple::apply_construct_parenthesis<mini_vector<T, std::tuple_size<TU>::value>>(t);
}
*/
template<typename T, typename TU>
mini_vector<T,std::tuple_size<TU>::value> tuple_to_mini_vector(TU const & t) {
mini_vector<T,std::tuple_size<TU>::value> res;
triqs::tuple::fold(tuple_to_mini_vector_aux(),t,&res[0]);
return res;
}
}}//namespace triqs::arrays
#endif
}}//namespace triqs::arrays
namespace std { // overload std::get to work with it
template <int i, typename T, int R> AUTO_DECL get(triqs::utility::mini_vector<T, R> &v) RETURN(v[i]);
template <int i, typename T, int R> AUTO_DECL get(triqs::utility::mini_vector<T, R> const &v) RETURN(v[i]);
}
#endif

38
triqs/utility/tuple_tools.hpp
View File

@ -68,6 +68,25 @@ namespace std {
namespace triqs { namespace tuple {
/// Repeat an element
template<typename T, int R> struct make_tuple_repeat_impl;
template <typename T> struct make_tuple_repeat_impl<T,1> {
static std::tuple<T> invoke(T&&x) { return {x};}
};
template <typename T> struct make_tuple_repeat_impl<T,2> {
static std::tuple<T,T> invoke(T&&x) { return {x,x};}
};
template <typename T> struct make_tuple_repeat_impl<T,3> {
static std::tuple<T,T,T> invoke(T&&x) { return {x,x,x};}
};
template <typename T> struct make_tuple_repeat_impl<T, 4> {
static std::tuple<T,T,T,T> invoke(T &&x) { return {x, x, x, x}; }
};
template <int R, typename T> auto make_tuple_repeat(T &&x) { return make_tuple_repeat_impl<T, R>::invoke(std::forward<T>(x)); }
/// _get_seq<T>() : from a tuple T, return the index sequence of the tuple length
template <typename T> std14::make_index_sequence<std::tuple_size<std14::decay_t<T>>::value> _get_seq() {
return {};
@ -120,6 +139,18 @@ namespace triqs { namespace tuple {
template <typename C, typename T>
AUTO_DECL apply_construct(T &&t) RETURN(apply_construct_impl<C>(std::forward<T>(t), _get_seq<T>()));
/**
* apply_construct_parenthesis<C>(t)
* C : a class
* t a tuple
* Returns : C(t0, t1, ....)
*/
template <typename C, typename T, size_t... Is>
AUTO_DECL apply_construct_parenthesis_impl(T &&t, std14::index_sequence<Is...>) RETURN(C(std::get<Is>(std::forward<T>(t))...));
template <typename C, typename T>
AUTO_DECL apply_construct_parenthesis(T &&t) RETURN(apply_construct_parenthesis_impl<C>(std::forward<T>(t), _get_seq<T>()));
/**
* called_on_tuple(f)
* f : a callable object
@ -455,7 +486,7 @@ namespace triqs { namespace tuple {
/*
* print a tuple
*/
inline void _triqs_print_tuple_impl(std::ostream &os) {}
/*inline void _triqs_print_tuple_impl(std::ostream &os) {}
template <typename T0, typename... T> void _triqs_print_tuple_impl(std::ostream &os, T0 const &x0, T const &... x) {
os << x0;
if (sizeof...(T) > 0) os << ',';
@ -466,6 +497,8 @@ namespace triqs { namespace tuple {
void _triqs_print_tuple(std::ostream &os, std::tuple<T...> const &t, std14::index_sequence<Is...>) {
_triqs_print_tuple_impl(os, std::get<Is>(t)...);
}
*/
}
}
@ -473,7 +506,8 @@ namespace std {
template <typename... T> std::ostream &operator<<(std::ostream &os, std::tuple<T...> const &t) {
os << "(";
triqs::tuple::_triqs_print_tuple(os, t, std14::make_index_sequence<sizeof...(T)>());
triqs::tuple::for_each(t, [&os, c=0](auto & x) mutable {if (c++) os << ','; os << x;});
//triqs::tuple::_triqs_print_tuple(os, t, std14::make_index_sequence<sizeof...(T)>());
return os << ")";
}