diff --git a/CMakeLists.txt b/CMakeLists.txt index 3208e16b..bc0dd43a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -174,11 +174,7 @@ if(CMAKE_COMPILER_IS_ICC) else() set(BOOST_VERSION 1.46) endif() -if (PythonSupport) - find_package(Boost ${BOOST_VERSION} REQUIRED COMPONENTS mpi python serialization system ) -else() - find_package(Boost ${BOOST_VERSION} REQUIRED COMPONENTS mpi serialization system ) -endif() +find_package(Boost ${BOOST_VERSION} REQUIRED COMPONENTS mpi serialization system ) if (NOT Boost_FOUND) message(FATAL_ERROR "Boost not found") endif() diff --git a/pytriqs/gf/local/gf_desc.py b/pytriqs/gf/local/gf_desc.py index 570c6b60..6af43830 100644 --- a/pytriqs/gf/local/gf_desc.py +++ b/pytriqs/gf/local/gf_desc.py @@ -8,7 +8,6 @@ module.add_include("") module.add_include("") module.add_using("namespace triqs::arrays") module.add_using("namespace triqs::gfs") -module.add_using("namespace triqs::gfs::local") module.add_using("triqs::utility::mini_vector") ######################## @@ -16,8 +15,8 @@ module.add_using("triqs::utility::mini_vector") ######################## t = class_( py_type = "TailGf", - c_type = "local::tail_view", - c_type_absolute = "triqs::gfs::local::tail_view", + c_type = "tail_view", + c_type_absolute = "triqs::gfs::tail_view", serializable= "tuple", is_printable= True, arithmetic = ("algebra","double") @@ -41,7 +40,7 @@ t.add_property(getter = cfunction("int order_max()"), doc = "Max order of the expansion") t.add_property(name = "mask", - getter = cfunction("array_view mask_view()"), + getter = cfunction("array_view mask()"), doc = "Access to the mask") t.add_method(name = "has_coef", @@ -51,7 +50,7 @@ t.add_method(name = "has_coef", # strange, I should not quality : ADL ?? t.add_method(name = "invert", - calling_pattern = "self_c = local::inverse(self_c)", + calling_pattern = "self_c = inverse(self_c)", signature = "void()", doc = "Invert") @@ -63,7 +62,7 @@ t.add_method(name = "zero", t.add_method_copy() t.add_method_copy_from() -t.add_call(calling_pattern = "auto result = self_c.evaluate(u)", +t.add_call(calling_pattern = "auto result = evaluate(self_c,u)", signature = "dcomplex(dcomplex u)", doc = "") @@ -252,7 +251,7 @@ def make_gf( py_type, c_tag, is_complex_data = True, is_im = False, has_tail = T if has_tail: g.add_property(name = "tail", - getter = cfunction(c_name="singularity", signature = "local::tail_view()"), + getter = cfunction(c_name="singularity", signature = "tail_view()"), doc ="The high frequency tail") g.add_property(name = "indices", @@ -332,7 +331,7 @@ def make_gf( py_type, c_tag, is_complex_data = True, is_im = False, has_tail = T g.number_protocol['multiply'].add_overload(calling_pattern = "*", signature = "gf<%s>(gf<%s> x,matrix<%s> y)"%(c_tag,c_tag,data_type)) g.add_method(name = "from_L_G_R", - calling_pattern = "self_c = L_G_R(l,g,r)", + calling_pattern = "self_c = L_G_R(l,g(),r)", signature = "void(matrix<%s> l,gf<%s> g,matrix<%s> r)"%(data_type,c_tag,data_type), doc = "self <<= l * g * r") diff --git a/test/triqs/gfs/g_r_k.cpp b/test/triqs/gfs/g_r_k.cpp new file mode 100644 index 00000000..38a1d48c --- /dev/null +++ b/test/triqs/gfs/g_r_k.cpp @@ -0,0 +1,57 @@ +#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK +#include +#include +#include +#include + +namespace h5 = triqs::h5; +using namespace triqs::gfs; +using namespace triqs::clef; +using namespace triqs::arrays; +using namespace triqs::lattice; + +#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl; + +int main() { + try { + double beta = 1; + int N = 10; + int S = 2; + placeholder<0> k_; + placeholder<1> r_; + + auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix(2)}}; + auto gk = gf{{bz_, N}, {S, S}}; + + gk(k_) << -2 * (cos(k_(0)) + cos(k_(1))); + + auto gr = gf{{N, N}, {S, S}}; + + gr() = inverse_fourier(gk); // ADD Security if gf = + + // for (auto k : gk.mesh()) if (max_element(abs(gk[k])) > 1.e-10) std::cout << k << gk[k] << std::endl; + + for (auto r : gr.mesh()) + if (max_element(abs(gr[r])) > 1.e-10) std::cout << r << gr[r] << std::endl; + + // reverse transfo + auto gk2 = gk; + gk2() = fourier(gr); + if (max_element(abs(gk.data() - gk2.data())) > 1.e-13) TRIQS_RUNTIME_ERROR << "fourier pb"; + + // check gr + auto gr_test = gf{{N, N}, {S, S}}; + gr_test[{1,0,0}] = -1; + gr_test[{-1,0,0}] = -1; + gr_test[{0,1,0}] = -1; + gr_test[{0,-1,0}] = -1; + + if (max_element(abs(gr.data() - gr_test.data())) > 1.e-13) TRIQS_RUNTIME_ERROR << "fourier inverse pb"; // << gr.data() << gr_test.data(); + + // hdf5 + h5::file file("ess_g_r_k.h5", H5F_ACC_TRUNC); + h5_write(file, "gr", gr); + h5_write(file, "gk", gk); + } + TRIQS_CATCH_AND_ABORT; +} diff --git a/test/triqs/gfs/gf.output b/test/triqs/gfs/gf.output deleted file mode 100644 index 8f13d685..00000000 --- a/test/triqs/gfs/gf.output +++ /dev/null @@ -1,156 +0,0 @@ -(G1( 0)) ---> empty - -(G( 0)) ---> -[[(0,0),(0,0)] - [(0,0),(0,0)]] - -(Gv2( 0)) ---> -[[(0,0)]] - -(Gv2( 0)) ---> -[[(10,0)]] - -(G( 0)) ---> -[[(10,0),(0,0)] - [(0,0),(0,0)]] - -(G(om_)) ---> gf_view(x0,) - -(tql::eval(G(om_), om_=0)) ---> -[[(10,0),(0,0)] - [(0,0),(0,0)]] - -(Gv(om_)) ---> gf_view(x0,) - -(tql::eval(Gv(om_), om_=0)) ---> -[[(10,0),(0,0)] - [(0,0),(0,0)]] - --------------lazy assign ------------------ -(G(0)) ---> -[[(2.3,3.14159),(0,0)] - [(0,0),(2.3,3.14159)]] - -(G(inf)) ---> tail/tail_view : Ordermin/max = -1 3 - ... Order -1 = -[[(1,0),(0,0)] - [(0,0),(1,0)]] - ... Order 0 = -[[(2.3,0),(0,0)] - [(0,0),(2.3,0)]] - ... Order 1 = -[[(0,0),(0,0)] - [(0,0),(0,0)]] - ... Order 2 = -[[(0,0),(0,0)] - [(0,0),(0,0)]] - ... Order 3 = -[[(0,0),(0,0)] - [(0,0),(0,0)]] - --------------lazy assign ------------------ -(G(0)) ---> -[[(0.151719,-0.207234),(0,0)] - [(0,0),(0.151719,-0.207234)]] - -(G(inf)) ---> tail/tail_view : Ordermin/max = -1 3 - ... Order -1 = -[[(1,0),(0,0)] - [(0,0),(1,0)]] - ... Order 0 = -[[(-2.3,0),(0,0)] - [(0,0),(-2.3,0)]] - ... Order 1 = -[[(5.29,0),(0,0)] - [(0,0),(5.29,0)]] - ... Order 2 = -[[(-12.167,0),(0,0)] - [(0,0),(-12.167,0)]] - ... Order 3 = -[[(27.9841,0),(0,0)] - [(0,0),(27.9841,0)]] - -(inverse(G(inf))) ---> tail_inv_lazy - -------------------------------------- -(Gv(om_)) ---> gf_view(x0,) - -(tql::eval(Gv(om_), om_=0)) ---> -[[(0.151719,-0.207234),(0,0)] - [(0,0),(0.151719,-0.207234)]] - -(t.order_min()) ---> -1 - -(t( 0)) ---> -[[(-2.3,0),(0,0)] - [(0,0),(-2.3,0)]] - -(Gv2(inf)( 0)) ---> -[[(-2.3,0)]] - -(G( 0)) ---> -[[(0.151719,-0.207234),(0,0)] - [(0,0),(0.151719,-0.207234)]] - -(Gc( 0)) ---> -[[(0.151719,-0.207234),(0,0)] - [(0,0),(0.151719,-0.207234)]] - -(G( 0)) ---> -[[(0.151719,-0.207234),(0,0)] - [(0,0),(0.151719,-0.207234)]] - -(G(inf)(0)) ---> -[[(-2.3,0),(0,0)] - [(0,0),(-2.3,0)]] - -(( G(inf) + G(inf) ) (0)) ---> ( -[[(-2.3,0),(0,0)] - [(0,0),(-2.3,0)]] + -[[(-2.3,0),(0,0)] - [(0,0),(-2.3,0)]]) - -(( G(inf) * G(inf) ) (0)) ---> -[[(15.87,0),(0,0)] - [(0,0),(15.87,0)]] - -((G + Gc)( inf)) ---> ( tail/tail_view : Ordermin/max = -1 3 - ... Order -1 = -[[(1,0),(0,0)] - [(0,0),(1,0)]] - ... Order 0 = -[[(-2.3,0),(0,0)] - [(0,0),(-2.3,0)]] - ... Order 1 = -[[(5.29,0),(0,0)] - [(0,0),(5.29,0)]] - ... Order 2 = -[[(-12.167,0),(0,0)] - [(0,0),(-12.167,0)]] - ... Order 3 = -[[(27.9841,0),(0,0)] - [(0,0),(27.9841,0)]] + tail/tail_view : Ordermin/max = -1 3 - ... Order -1 = -[[(1,0),(0,0)] - [(0,0),(1,0)]] - ... Order 0 = -[[(-2.3,0),(0,0)] - [(0,0),(-2.3,0)]] - ... Order 1 = -[[(5.29,0),(0,0)] - [(0,0),(5.29,0)]] - ... Order 2 = -[[(-12.167,0),(0,0)] - [(0,0),(-12.167,0)]] - ... Order 3 = -[[(27.9841,0),(0,0)] - [(0,0),(27.9841,0)]]) - -((t + 2.3)(0)) ---> ( -[[(-2.3,0),(0,0)] - [(0,0),(-2.3,0)]] + (2.3,0)) - -(t(-1)) ---> -[[(1,0),(0,0)] - [(0,0),(1,0)]] - diff --git a/test/triqs/gfs/gf_notail.cpp b/test/triqs/gfs/gf_notail.cpp index d1270970..9adae332 100644 --- a/test/triqs/gfs/gf_notail.cpp +++ b/test/triqs/gfs/gf_notail.cpp @@ -48,7 +48,7 @@ int main() { gw_n (tau_) << 1/(tau_-1.); auto gt_with_full_tail = make_gf_from_inverse_fourier(make_gf_from_g_and_tail(gw_n, gw.singularity())); TEST(gt_with_full_tail(.5)); - triqs::gfs::local::tail t(2,2); + triqs::gfs::tail t(2,2); t(1)=1; TEST(t); auto gt_tail_with_one_term = make_gf_from_inverse_fourier(make_gf_from_g_and_tail(gw_n, t)); diff --git a/test/triqs/gfs/gfv2.cpp b/test/triqs/gfs/gfv2.cpp index e8464fca..b2f663eb 100644 --- a/test/triqs/gfs/gfv2.cpp +++ b/test/triqs/gfs/gfv2.cpp @@ -14,7 +14,6 @@ namespace h5 = triqs::h5; int main() { try { - triqs::gfs::freq_infty inf; double beta =1; @@ -48,15 +47,15 @@ int main() { Gv(om_) << (0.2 + om_ + 2.1); TEST(G(0)); - TEST(G(inf)); + TEST(G.singularity()); std::cout <<"-------------lazy assign 2 ------------------"< tail/tail_view: min/smallest/max = -1 -1 8 +(G.singularity()) ---> tail/tail_view: min/smallest/max = -1 -1 8 ... Order -1 = [[(1,0),(0,0)] [(0,0),(1,0)]] @@ -74,7 +74,7 @@ [[(0.151719,-0.207234),(0,0)] [(0,0),(0.151719,-0.207234)]] -(G(inf)) ---> tail/tail_view: min/smallest/max = -1 1 8 +(G.singularity()) ---> tail/tail_view: min/smallest/max = -1 1 8 ... Order -1 = [[(0,0),(0,0)] [(0,0),(0,0)]] @@ -106,7 +106,7 @@ [[(-340.483,0),(0,0)] [(0,0),(-340.483,0)]] -(inverse(G(inf))) ---> tail/tail_view: min/smallest/max = -1 -1 6 +(inverse(G.singularity())) ---> tail/tail_view: min/smallest/max = -1 -1 6 ... Order -1 = [[(1,0),(0,0)] [(0,0),(1,0)]] @@ -145,7 +145,7 @@ [[(-2.3,0),(0,0)] [(0,0),(-2.3,0)]] -(Gv2(inf)( 2)) ---> +(Gv2.singularity()( 2)) ---> [[(-2.3,0)]] (G( 0)) ---> @@ -165,15 +165,15 @@ [[(0.151719,-0.207234),(0,0)] [(0,0),(0.151719,-0.207234)]] -(G(inf)(2)) ---> +(G.singularity()(2)) ---> [[(-2.3,0),(0,0)] [(0,0),(-2.3,0)]] -(( G(inf) + G(inf) ) (2)) ---> +(( G.singularity() + G.singularity() ) (2)) ---> [[(-4.6,0),(0,0)] [(0,0),(-4.6,0)]] -(( G(inf) * G(inf) ) (4)) ---> +(( G.singularity() * G.singularity() ) (4)) ---> [[(15.87,0),(0,0)] [(0,0),(15.87,0)]] diff --git a/test/triqs/gfs/multivar/curry_and_fourier.cpp b/test/triqs/gfs/multivar/curry_and_fourier.cpp index 8a156839..7870ec9a 100644 --- a/test/triqs/gfs/multivar/curry_and_fourier.cpp +++ b/test/triqs/gfs/multivar/curry_and_fourier.cpp @@ -1,43 +1,71 @@ #define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK +#include +#include +#include -#include -#include -#include - -#include -#include -#include - -using namespace triqs::gfs; namespace h5 = triqs::h5; +using namespace triqs::gfs; +using namespace triqs::clef; +using namespace triqs::arrays; +using namespace triqs::lattice; + +template +void assert_equal(T const& x, T const& y, std::string mess) { + if (std::abs(x - y) > 1.e-13) TRIQS_RUNTIME_ERROR << mess; +} + +template +void assert_equal_array(T1 const& x, T2 const& y, std::string mess) { + if (max_element(abs(x - y)) > 1.e-13) TRIQS_RUNTIME_ERROR << mess << "\n" << x << "\n" << y << "\n" << max_element(abs(x - y)); +} + +#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl; int main() { + try { + double beta = 1; + auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix(2)}}; -try { - double tmax=10; - double wmin=0.; - double wmax=1.0; - int n_re_freq=100; - int Nt=100; + int n_freq = 100; + int n_times = n_freq * 2; + int n_bz = 50; + auto gkw = gf>{{{bz_, n_bz}, {beta, Fermion, n_freq}}, {1, 1}}; + auto gkt = gf>{{{bz_, n_bz}, {beta, Fermion, n_times}}, {1, 1}}; + + placeholder<0> k_; + placeholder<1> w_; + placeholder<2> wn_; + placeholder<3> tau_; + + auto eps_k = -2 * (cos(k_(0)) + cos(k_(1))); + gkw(k_, w_) << 1 / (w_ - eps_k - 1 / (w_ + 2)); + + auto gk_w = curry<0>(gkw); + auto gk_t = curry<0>(gkt); + + gk_t[k_] << inverse_fourier(gk_w[k_]); + + // works also, but uses the evaluator which return to the same point + // gk_t(k_) << inverse_fourier(gk_w(k_)); + // check last assertion + for (auto k : gk_t.mesh()) assert_equal(k.linear_index(), gk_t.mesh().locate_neighbours(k).linear_index(), "k location point"); - triqs::clef::placeholder<0> w_; - triqs::clef::placeholder<1> wn_; - triqs::clef::placeholder<2> tau_; - - auto G_w_wn = gf, scalar_valued>( {gf_mesh(wmin, wmax, n_re_freq), gf_mesh(wmin, wmax, n_re_freq)}); - auto G_w_tau = gf, scalar_valued>( {gf_mesh(wmin, wmax, n_re_freq), gf_mesh(-tmax, tmax, Nt)}); + /// Testing the result + auto gk_w_test = gf{{beta, Fermion, n_freq}, {1, 1}}; + auto gk_t_test = gf{{beta, Fermion, n_times}, {1, 1}}; + assert_equal_array(gkt.singularity().data().data, gkw.singularity().data().data, "Error 05"); + for (auto & k : std::get<0>(gkw.mesh().components())) { + gk_w_test(w_) << 1 / (w_ - eval(eps_k, k_ = k) - 1 / (w_ + 2)); + gk_t_test() = inverse_fourier(gk_w_test); + assert_equal_array(gk_w_test.singularity().data(), gk_w[k].singularity().data(), "Error 0"); + assert_equal_array(gk_t_test.singularity().data(), gk_t[k].singularity().data(), "Error 0s"); + assert_equal_array(gk_w_test.data(), gk_w[k].data(), "Error 1"); + assert_equal_array(gk_t_test.data(), gk_t[k].data(), "Error 2"); + } - G_w_wn(w_,wn_)<<1/(wn_-1)/( pow(w_,3) ); - G_w_tau(w_,tau_)<< exp( -2*tau_ ) / (w_*w_ + 1 ); - - 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_]); -} -// temp fix : THE TEST DOES NOT RUN !! -//TRIQS_CATCH_AND_ABORT; -catch(std::exception const & e ) { std::cout << "error "<< e.what()<< std::endl;} + // hdf5 + h5::file file("ess_g_k_om.h5", H5F_ACC_TRUNC); + h5_write(file, "g", gkw); + } + TRIQS_CATCH_AND_ABORT; } diff --git a/test/triqs/gfs/multivar/g_k_om.cpp b/test/triqs/gfs/multivar/g_k_om.cpp index db2da644..3549b407 100644 --- a/test/triqs/gfs/multivar/g_k_om.cpp +++ b/test/triqs/gfs/multivar/g_k_om.cpp @@ -1,6 +1,7 @@ #define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK #include #include +#include namespace h5 = triqs::h5; using namespace triqs::gfs; @@ -8,34 +9,79 @@ using namespace triqs::clef; using namespace triqs::arrays; using namespace triqs::lattice; +template +void assert_equal(T1 const& x, T2 const& y, std::string mess) { + if (std::abs(x - y) > 1.e-13) TRIQS_RUNTIME_ERROR << mess; +} + #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl; int main() { try { double beta = 1; - - auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix(2)}}; + auto bz_ = brillouin_zone{bravais_lattice{{{1, 0}, {0, 1}}}}; auto g_eps = gf{{bz_, 100}, {1, 1}}; + // Default, with a tail. auto G = gf>{{{bz_, 100}, {beta, Fermion, 100}}, {1, 1}}; - // try to assign to expression + // Same : the gf_default_singularity is the same + static_assert(std::is_same, gf>::value, "!!"); + + // auto G = gf, matrix_valued, gf>{{{bz_, 100}, {beta, Fermion, 100}}, {1, 1}}; + + // The variant without tail + auto G_no_tail = gf, matrix_valued, no_tail>{{{bz_, 100}, {beta, Fermion, 100}}, {1, 1}}; + + // try to assign to expression placeholder<0> k_; placeholder<1> w_; + /// The calculation is 3x longer with the tail !! auto eps_k = -2 * (cos(k_(0)) + cos(k_(1))); G(k_, w_) << 1 / (w_ - eps_k - 1 / (w_ + 2)); - // same - // auto eps = make_expr( [](k_t const& k) { return -2 * (cos(k(0)) + cos(k(1))); }) ; - // G(k_, w_) << 1 / (w_ - eps(k_) - 1 / (w_ + 2)); + //G_no_tail(k_, w_) << 1 / (w_ - eps_k - 1 / (w_ + 2)); + //return 0; + if (0) { // to check speed of computation, which is long due to tail (inversion ?) + auto t = tail_omega(1, 1); + auto t2 = t; + for (int u = 0; u < 10000; ++u) { + t2 = (t + u - inverse(t - u)); + //t2 = inverse(t + u - inverse(t - u)); + } + } + + for (auto k : std::get<0>(G.mesh().components())) { + assert_equal(G.singularity()[k](1)(0, 0), 1, "should be 1/omega"); + assert_equal(G.singularity()[k](2)(0, 0), eval(eps_k, k_ = k), "should be eps_k/omega^2"); + } + + auto gk = curry<0>(G); + + std::cerr << G.singularity()[{0, 0, 0}] << std::endl; + std::cerr << gk[{0, 0, 0}].singularity() << std::endl; + + if (1) { + + for (auto k : gk.mesh()) { + // std::cout << k << std::endl; + auto p = gk[k]; + assert_equal(gk[k].singularity()(1)(0, 0), 1, "should be 1/omega"); + assert_equal(gk[k].singularity()(2)(0, 0), eval(eps_k, k_ = k), "should be eps_k/omega^2"); + } + + for (auto k : gk.mesh()) { + auto p = partial_eval<0>(G, 0); + assert_equal(gk[k].singularity()(1)(0, 0), 1, "should be 1/omega"); + assert_equal(gk[k].singularity()(2)(0, 0), eval(eps_k, k_ = k), "should be eps_k/omega^2"); + } + } // hdf5 h5::file file("ess_g_k_om.h5", H5F_ACC_TRUNC); h5_write(file, "g", G); - - } TRIQS_CATCH_AND_ABORT; } diff --git a/test/triqs/gfs/multivar/g_k_om2.cpp b/test/triqs/gfs/multivar/g_k_om2.cpp index f02b112c..642873b4 100644 --- a/test/triqs/gfs/multivar/g_k_om2.cpp +++ b/test/triqs/gfs/multivar/g_k_om2.cpp @@ -44,10 +44,9 @@ int main() { G_k_iom(k_, w_) << 1 / (w_ - eps(k_)); - auto G_loc = gf{{beta, Fermion, 100}, {1, 1}}; - auto r = G_k_iom(k_t{0, 0}, matsubara_freq{0, beta, Fermion}); + auto r = G_k_iom(k_t{0, 0, 0}, matsubara_freq{0, beta, Fermion}); auto r5 = sum_gf(k_ >> G_k_iom(k_,0), g_eps.mesh()); G_loc(w_) << sum_gf(k_ >> G_k_iom(k_,w_), g_eps.mesh()); diff --git a/test/triqs/gfs/multivar/g_k_tail.cpp b/test/triqs/gfs/multivar/g_k_tail.cpp index 3d46bb6a..004e46fa 100644 --- a/test/triqs/gfs/multivar/g_k_tail.cpp +++ b/test/triqs/gfs/multivar/g_k_tail.cpp @@ -8,7 +8,6 @@ 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; @@ -17,23 +16,16 @@ using local::tail; int main() { try { double beta = 1; - auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix(2)}}; auto t = tail{1,1}; - auto G = gf{{bz_, 100}}; + auto G = gf{{bz_, 100}, {1,1}}; - // 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_) << -2 * (cos(k_(0)) + cos(k_(1))) * t; - 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; diff --git a/test/triqs/gfs/test_fit_tail.cpp b/test/triqs/gfs/test_fit_tail.cpp index 0b9d7472..983ecf02 100644 --- a/test/triqs/gfs/test_fit_tail.cpp +++ b/test/triqs/gfs/test_fit_tail.cpp @@ -5,7 +5,7 @@ using triqs::arrays::make_shape; using namespace triqs::gfs; -using triqs::gfs::local::tail; +using triqs::gfs::tail; #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) < - -struct infty {}; - -template struct wrap { - T x; - wrap(T const & x_): x(x_) {} - operator infty() const { return infty();} - operator T () const { return x;} -}; - - -double g(infty) { return 3;} -double f(double) { return 30;} - -int main() { - - wrap w(2); - std::cout << w +1 <. * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_ALL_H -#define TRIQS_ARRAYS_ALL_H +#pragma once +// for python code generator, we need to know what to include... #define TRIQS_ARRAYS_INCLUDED // The basic classes @@ -31,6 +31,7 @@ // #include #include +#include #include #include @@ -42,7 +43,6 @@ #include // Linear algebra ?? Keep here ? -//#include +#include -#endif diff --git a/triqs/gfs/block.hpp b/triqs/gfs/block.hpp index 9bbf5167..c101a9c1 100644 --- a/triqs/gfs/block.hpp +++ b/triqs/gfs/block.hpp @@ -69,18 +69,14 @@ namespace gfs { // ------------------------------- Factories -------------------------------------------------- - template struct factories { + template struct data_factory { using mesh_t = gf_mesh; using gf_t = gf; using gf_view_t = gf_view; using aux_t = nothing; - struct target_shape_t {}; - static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t, aux_t) { return std::vector(m.size()); } - static typename gf_t::singularity_t make_singularity(mesh_t const &m, target_shape_t) { - return {}; - } + static typename gf_t::data_t make(mesh_t const &m, target_shape_t, aux_t) { return std::vector(m.size()); } }; } // gfs_implementation diff --git a/triqs/gfs/bz.hpp b/triqs/gfs/bz.hpp index fb420c69..7229ab39 100644 --- a/triqs/gfs/bz.hpp +++ b/triqs/gfs/bz.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2012-2013 by O. Parcollet + * 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 @@ -22,7 +22,8 @@ #include "./tools.hpp" #include "./gf.hpp" #include "./local/tail.hpp" -#include "../lattice/bz_mesh.hpp" +#include "./domains/R.hpp" +#include "../lattice/regular_bz_mesh.hpp" #include "./evaluators.hpp" namespace triqs { @@ -30,8 +31,8 @@ namespace gfs { struct bz {}; - template struct gf_mesh : lattice::bz_mesh { - template gf_mesh(T &&... x) : lattice::bz_mesh(std::forward(x)...) {} + template struct gf_mesh : lattice::regular_bz_mesh { + template gf_mesh(T &&... x) : lattice::regular_bz_mesh(std::forward(x)...) {} }; namespace gfs_implementation { @@ -49,25 +50,23 @@ namespace gfs { // simple evaluation : take the point on the grid... template <> struct evaluator_fnt_on_mesh { + // lattice::regular_bz_mesh::index_t n; size_t n; evaluator_fnt_on_mesh() = default; - template evaluator_fnt_on_mesh(MeshType const &m, lattice::k_t const &k) { - n = m.locate_neighbours(k); // TO BE IMPROVED + template evaluator_fnt_on_mesh(MeshType const &m, lattice::k_t const &k) { + // n = m.locate_neighbours(k).index(); + n = m.locate_neighbours(k).linear_index(); } - template auto operator()(F const &f) const DECL_AND_RETURN(f(n)); - //template decltype(auto) operator()(F const &f) const { return f(n); } + template AUTO_DECL operator()(F const &f) const RETURN(f(n)); + // template decltype(auto) operator()(F const &f) const { return f(n); } }; - // ------------- evaluator ------------------- - // handle the case where the matsu. freq is out of grid... + // -------------------------------------------------------------- template struct evaluator { static constexpr int arity = 1; - template - std::c14::conditional_t::value, arrays::matrix_const_view>, - std::complex> - operator()(G const *g, lattice::k_t const &k) const { - auto n = g->mesh().locate_neighbours(k); // TO BE IMPROVED + template auto operator()(G const *g, lattice::k_t const &k) const { + auto n = g->mesh().locate_neighbours(k); return (*g)[n]; } }; diff --git a/triqs/gfs/bz1.hpp b/triqs/gfs/bz1.hpp new file mode 100644 index 00000000..bab93bd6 --- /dev/null +++ b/triqs/gfs/bz1.hpp @@ -0,0 +1,78 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2012-2013 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 . + * + ******************************************************************************/ +#pragma once +#include "./tools.hpp" +#include "./gf.hpp" +#include "./local/tail.hpp" +#include "./domains/R.hpp" +#include "../lattice/bz_mesh.hpp" +#include "./evaluators.hpp" + +namespace triqs { +namespace gfs { + + struct bz {}; + + template struct gf_mesh : lattice::bz_mesh { + template gf_mesh(T &&... x) : lattice::bz_mesh(std::forward(x)...) {} + }; + + namespace gfs_implementation { + + // h5 name + template struct h5_name { + static std::string invoke() { return "BZ"; } + }; + + /// --------------------------- data access --------------------------------- + template struct data_proxy : data_proxy_array, 3> {}; + template struct data_proxy : data_proxy_array, 1> {}; + + /// --------------------------- evaluator --------------------------------- + + // simple evaluation : take the point on the grid... + template <> struct evaluator_fnt_on_mesh { + size_t n; + evaluator_fnt_on_mesh() = default; + template evaluator_fnt_on_mesh(MeshType const &m, lattice::k_t const &k) { + n = m.locate_neighbours(k); // TO BE IMPROVED + } + template auto operator()(F const &f) const DECL_AND_RETURN(f(n)); + //template decltype(auto) operator()(F const &f) const { return f(n); } + }; + + // ------------- evaluator ------------------- + // handle the case where the matsu. freq is out of grid... + template struct evaluator { + static constexpr int arity = 1; + + template + std::c14::conditional_t::value, arrays::matrix_const_view>, + std::complex> + operator()(G const *g, lattice::k_t const &k) const { + auto n = g->mesh().locate_neighbours(k); // TO BE IMPROVED + return (*g)[n]; + } + }; + } +} +} + diff --git a/triqs/gfs/curry.hpp b/triqs/gfs/curry.hpp index d589e7cd..5c9176c0 100644 --- a/triqs/gfs/curry.hpp +++ b/triqs/gfs/curry.hpp @@ -37,7 +37,7 @@ namespace triqs { namespace gfs { /// --------------------------- Factories --------------------------------- template - struct factories, lambda_valued, nothing, Opt> {}; + struct data_factory, lambda_valued, nothing, Opt> {}; /// --------------------------- partial_eval --------------------------------- // partial_eval<0> (g, 1) returns : x -> g(1,x) diff --git a/triqs/gfs/cyclic_lattice.hpp b/triqs/gfs/cyclic_lattice.hpp new file mode 100644 index 00000000..a258fe7a --- /dev/null +++ b/triqs/gfs/cyclic_lattice.hpp @@ -0,0 +1,74 @@ +/******************************************************************************* + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include "./tools.hpp" +#include "./gf.hpp" +#include "./local/tail.hpp" +#include "./domains/R.hpp" +#include "../lattice/cyclic_lattice.hpp" +#include "./evaluators.hpp" + +namespace triqs { +namespace gfs { + + struct cyclic_lattice {}; + + template struct gf_mesh : lattice::cyclic_lattice_mesh { + template gf_mesh(T &&... x) : lattice::cyclic_lattice_mesh(std::forward(x)...) {} + }; + + namespace gfs_implementation { + + // h5 name + template struct h5_name { + static std::string invoke() { return "R"; } + }; + + /// --------------------------- data access --------------------------------- + template struct data_proxy : data_proxy_array, 3> {}; + template struct data_proxy : data_proxy_array, 1> {}; + + /// --------------------------- evaluator --------------------------------- + + // simple evaluation : take the point on the grid... + template <> struct evaluator_fnt_on_mesh { + size_t n; + evaluator_fnt_on_mesh() = default; + template evaluator_fnt_on_mesh(MeshType const &m, R const &r) { + n = m.modulo_reduce(r).linear_index(); + } + template AUTO_DECL operator()(F const &f) const RETURN(f(n)); + //template decltype(auto) operator()(F const &f) const { return f(n); } + }; + + // -------------------------------------------------------------- + template struct evaluator { + static constexpr int arity = 1; + + template auto operator()(G const *g, R const &r) const { + auto n = g->mesh().modulo_reduce(r); + return (*g)[n]; + } + }; + } +} +} + diff --git a/triqs/gfs/evaluators.hpp b/triqs/gfs/evaluators.hpp index d834f029..92c946b4 100644 --- a/triqs/gfs/evaluators.hpp +++ b/triqs/gfs/evaluators.hpp @@ -70,8 +70,9 @@ namespace gfs { evaluator_one_var() = default; template auto operator()(G const *g, double x) const DECL_AND_RETURN(evaluator_fnt_on_mesh(g -> mesh(), x)(on_mesh(*g))); - template typename G::singularity_t const &operator()(G const *g, freq_infty const &) const { - return g->singularity(); + + template typename G::singularity_t operator()(G const *g, tail_view t) const { + return compose(g->singularity(), t); } }; } diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp index 7dcd0d04..436ad09a 100644 --- a/triqs/gfs/gf.hpp +++ b/triqs/gfs/gf.hpp @@ -39,6 +39,7 @@ namespace gfs { using arrays::array_view; using arrays::matrix; using arrays::matrix_view; + using arrays::matrix_const_view; using triqs::make_clone; // The default target, for each Variable. @@ -111,8 +112,9 @@ namespace gfs { template struct h5_name; // value is a const char template struct h5_rw; - // factories regroup all factories (constructors..) for all types of gf. Defaults implemented below. - template struct factories; + // factories (constructors..) for all types of gf. Defaults implemented below. + template struct data_factory; + template struct singularity_factory; } // gfs_implementation @@ -156,6 +158,7 @@ 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 linear_mesh_index_t = typename mesh_t::linear_index_t; using symmetry_t = typename gfs_implementation::symmetry::type; using indices_t = typename gfs_implementation::indices::type; using evaluator_t = gfs_implementation::evaluator; @@ -166,12 +169,11 @@ namespace gfs { using data_const_view_t = typename data_proxy_t::storage_const_view_t; using data_t = std14::conditional_t, data_regular_t>; - //using singularity_non_view_t = typename gfs_implementation::singularity::type; - using singularity_non_view_t = Singularity; - using singularity_view_t = typename view_type_if_exists_else_type::type; - using singularity_t = std14::conditional_t; - - //using is_container_t = void; // a tag to recognize the container + using singularity_regular_t = typename Singularity::regular_type; + using singularity_view_t = typename Singularity::view_type; + using singularity_const_view_t = typename Singularity::const_view_type; + using singularity_t = std14::conditional_t, + singularity_regular_t>; mesh_t const &mesh() const { return _mesh; } domain_t const &domain() const { return _mesh.domain(); } @@ -314,6 +316,15 @@ namespace gfs { } /// --------------------- A direct access to the grid point -------------------------- + + template r_type get_from_linear_index(Args &&... args) { + return _data_proxy(_data, linear_mesh_index_t(std::forward(args)...)); + } + + template cr_type get_from_linear_index(Args &&... args) const { + return _data_proxy(_data, linear_mesh_index_t(std::forward(args)...)); + } + template r_type on_mesh(Args &&... args) { return _data_proxy(_data, _mesh.index_to_linear(mesh_index_t(std::forward(args)...))); } @@ -391,10 +402,9 @@ namespace gfs { template void triqs_clef_auto_assign(gf_impl &g, RHS const &rhs) { triqs_clef_auto_assign_impl(g, rhs, typename std::is_base_of>::type()); - assign_from_expression(g.singularity(), rhs); + assign_singularity_from_function(g.singularity(), rhs); // access to the data . Beware, we view it as a *matrix* NOT an array... (crucial for assignment to scalars !) - // if f is an expression, replace the placeholder with a simple tail. If f is a function callable on freq_infty, - // it uses the fact that tail_non_view_t can be casted into freq_infty + // if f is an expression, replace the placeholder with a simple tail. } // enable the writing g[om_] << .... also @@ -436,7 +446,8 @@ namespace gfs { template class gf : public gf_impl { using B = gf_impl; - using factory = gfs_implementation::factories; + using data_factory = gfs_implementation::data_factory; + using singularity_factory = gfs_implementation::singularity_factory; public: gf() : B() {} @@ -458,12 +469,12 @@ namespace gfs { typename B::indices_t const &ind, std::string name = "") : B(std::move(m), std::move(dat), si, s, ind, name, typename B::evaluator_t{}) {} - using target_shape_t = typename factory::target_shape_t; + using target_shape_t = typename data_factory::target_shape_t; // with aux, and indices - gf(typename B::mesh_t m, target_shape_t shape, typename factory::aux_t aux, + gf(typename B::mesh_t m, target_shape_t shape, typename data_factory::aux_t aux, typename B::indices_t const &ind = typename B::indices_t{}, std::string name = "") - : B(std::move(m), factory::make_data(m, shape, aux), factory::make_singularity(m, shape), typename B::symmetry_t{}, ind, + : B(std::move(m), data_factory::make(m, shape, aux), singularity_factory::make(m, shape), typename B::symmetry_t{}, ind, name, // clean unncessary types typename B::evaluator_t{}) { if (this->_indices.is_empty()) this->_indices = typename B::indices_t(shape); @@ -472,7 +483,7 @@ namespace gfs { // without the aux (defaulted) gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{}, typename B::indices_t const &ind = typename B::indices_t{}, std::string name = "") - : B(std::move(m), factory::make_data(m, shape, typename factory::aux_t{}), factory::make_singularity(m, shape), + : B(std::move(m), data_factory::make(m, shape, typename data_factory::aux_t{}), singularity_factory::make(m, shape), typename B::symmetry_t{}, ind, name, // clean unncessary types typename B::evaluator_t{}) { if (this->_indices.is_empty()) this->_indices = typename B::indices_t(shape); @@ -800,7 +811,7 @@ namespace gfs { } template - gf L_G_R(matrix l, gf g, matrix r) { + gf L_G_R(matrix l, gf_view g, matrix r) { auto res = gf{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; @@ -811,31 +822,35 @@ namespace gfs { // ------------------------- default factories --------------------- - template struct factories_default_impl { + // Factory for the data + template struct data_factory_default_impl { using gf_t = gf; using target_shape_t = arrays::mini_vector; using mesh_t = typename gf_t::mesh_t; using aux_t = arrays::memory_layout< get_n_variables(Var{}) + VarDim>; // - static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape, aux_t ml) { + static typename gf_t::data_t make(mesh_t const &m, target_shape_t shape, aux_t ml) { typename gf_t::data_t A(gf_t::data_proxy_t::join_shape(m.size_of_components(), shape), ml); A() = 0; return A; } - // - static typename gf_t::singularity_t make_singularity(mesh_t const &m, target_shape_t shape) { - return typename gf_t::singularity_t{shape}; - } }; template - struct factories, Singularity, Opt> : factories_default_impl, R, Singularity, Opt> {}; + struct data_factory, Singularity, Opt> : data_factory_default_impl, R, Singularity, Opt> {}; template - struct factories : factories_default_impl {}; + struct data_factory : data_factory_default_impl {}; template - struct factories : factories_default_impl {}; + struct data_factory : data_factory_default_impl {}; + + // Factory for the singularity + template struct singularity_factory { + template static Singularity make(gf_mesh const &m, TargetShape shape) { + return Singularity{shape}; + } + }; // --------------------- hdf5 --------------------------------------- // scalar function : just add a _s diff --git a/triqs/gfs/imfreq.hpp b/triqs/gfs/imfreq.hpp index 63480ec9..801efade 100644 --- a/triqs/gfs/imfreq.hpp +++ b/triqs/gfs/imfreq.hpp @@ -37,10 +37,10 @@ namespace gfs { // singularity template <> struct gf_default_singularity { - using type = local::tail; + using type = tail; }; template <> struct gf_default_singularity { - using type = local::tail; + using type = tail; }; namespace gfs_implementation { @@ -68,7 +68,7 @@ namespace gfs { if ((p.n >= m.first_index()) && (p.n < m.size()+m.first_index())) {w=1; n =p.n;} else {w=0; n=0;} } - template auto operator()(F const &f) const DECL_AND_RETURN(w*f(n)); + template AUTO_DECL operator()(F const &f) const RETURN(w*f(n)); }; // ------------- evaluator ------------------- @@ -83,14 +83,15 @@ namespace gfs { 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::singularity_t const &operator()(G const *g, freq_infty const &) const { - return g->singularity(); + template typename G::singularity_t operator()(G const *g, tail_view t) const { + return compose(g->singularity(),t); + //return g->singularity(); } }; // --- various 4 specializations // scalar_valued, tail - template struct evaluator : _eval_imfreq_base_impl { + template struct evaluator : _eval_imfreq_base_impl { using _eval_imfreq_base_impl::operator(); @@ -101,7 +102,7 @@ namespace gfs { } 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); + return evaluate(g->singularity(),f)(0, 0); } }; @@ -122,7 +123,7 @@ namespace gfs { }; // matrix_valued, tail - template struct evaluator : _eval_imfreq_base_impl { + template struct evaluator : _eval_imfreq_base_impl { using _eval_imfreq_base_impl::operator(); @@ -134,7 +135,7 @@ namespace gfs { } 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); + return evaluate(g->singularity(), f); } }; diff --git a/triqs/gfs/imtime.hpp b/triqs/gfs/imtime.hpp index d5a126ce..0ed88bb6 100644 --- a/triqs/gfs/imtime.hpp +++ b/triqs/gfs/imtime.hpp @@ -39,10 +39,10 @@ namespace gfs { // singularity template <> struct gf_default_singularity { - using type = local::tail; + using type = tail; }; template <> struct gf_default_singularity { - using type = local::tail; + using type = tail; }; namespace gfs_implementation { diff --git a/triqs/gfs/local/fit_tail.cpp b/triqs/gfs/local/fit_tail.cpp index 250c4a17..57f7c764 100644 --- a/triqs/gfs/local/fit_tail.cpp +++ b/triqs/gfs/local/fit_tail.cpp @@ -1,5 +1,5 @@ #include "./fit_tail.hpp" -namespace triqs { namespace gfs { namespace local { +namespace triqs { namespace gfs { tail fit_tail_impl(gf_view gf, const tail_view known_moments, int n_moments, int n_min, int n_max) { @@ -44,7 +44,7 @@ namespace triqs { namespace gfs { namespace local { B(k, 0) = imag(gf.data()(gf.mesh().index_to_linear(n), i, j)); // subtract known tail if present if (known_moments.size() > 0) - B(k, 0) -= imag(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)).evaluate(iw)(0, 0)); + B(k, 0) -= imag(evaluate(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)), iw)(0, 0)); for (int l = 0; l < size_odd; l++) { int order = omin_odd + 2 * l; @@ -67,7 +67,7 @@ namespace triqs { namespace gfs { namespace local { B(k, 0) = real(gf.data()(gf.mesh().index_to_linear(n), i, j)); // subtract known tail if present if (known_moments.size() > 0) - B(k, 0) -= real(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)).evaluate(iw)(0, 0)); + B(k, 0) -= real(evaluate(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)), iw)(0, 0)); for (int l = 0; l < size_even; l++) { int order = omin_even + 2 * l; @@ -81,7 +81,7 @@ namespace triqs { namespace gfs { namespace local { } } } - res.mask_view()=n_moments; + res.mask()()=n_moments; return res; // return tail } @@ -92,7 +92,7 @@ namespace triqs { namespace gfs { namespace local { if (replace_by_fit) { // replace data in the fitting range by the values from the fitted tail int i = 0; for (auto iw : gf.mesh()) { // (arrays::range(n_min,n_max+1)) { - if (i >= n_min) gf[iw] = gf.singularity().evaluate(iw); + if (i >= n_min) gf[iw] = evaluate(gf.singularity(),iw); i++; } } @@ -109,4 +109,4 @@ namespace triqs { namespace gfs { namespace local { fit_tail(reinterpret_scalar_valued_gf_as_matrix_valued(gf), known_moments, n_moments, n_min, n_max, replace_by_fit ); } - }}} // namespace + }} // namespace diff --git a/triqs/gfs/local/fit_tail.hpp b/triqs/gfs/local/fit_tail.hpp index 69e0d0db..5b54dde6 100644 --- a/triqs/gfs/local/fit_tail.hpp +++ b/triqs/gfs/local/fit_tail.hpp @@ -24,14 +24,12 @@ #include #include -namespace triqs { namespace gfs { namespace local { +namespace triqs { namespace gfs { using triqs::gfs::imfreq; using triqs::gfs::block_index; using triqs::gfs::block_index; - namespace tgl = triqs::gfs::local; - // routine for fitting the tail (singularity) of a Matsubara Green's function // this is a *linear* least squares problem (with non-linear basis functions) // which is solved by singular value decomposition of the design matrix diff --git a/triqs/gfs/local/fourier_base.hpp b/triqs/gfs/local/fourier_base.hpp index df749126..97addb0b 100644 --- a/triqs/gfs/local/fourier_base.hpp +++ b/triqs/gfs/local/fourier_base.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by M. Ferrero, O. Parcollet + * Copyright (C) 2011-2014 by M. Ferrero, O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -18,28 +18,42 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_LOCAL_FOURIER_BASE_H -#define TRIQS_GF_LOCAL_FOURIER_BASE_H +#pragma once #include +#include namespace triqs { namespace gfs { namespace details { namespace tqa = triqs::arrays; - typedef std::complex dcomplex; - + using dcomplex = std::complex; + void fourier_base(const tqa::vector &in, tqa::vector &out, size_t L, bool direct); void fourier_base(const tqa::vector &in, tqa::vector &out, size_t L1, size_t L2, bool direct); - } namespace tags { struct fourier{}; } + // ------------------------------------------------------------------- + + // TO BE REPLACED BY A DIRECT CALL to many_fft in fftw, like for lattice case. + // The implementation of the Fourier transformation + // Reduce Matrix case to the scalar case. + template + void _fourier_impl(gf_view gw, gf_const_view gt) { + if (gt.data().shape().front_pop() != gw.data().shape().front_pop()) + TRIQS_RUNTIME_ERROR << "Fourier : matrix size of target mismatch"; + for (size_t n1 = 0; n1 < gt.data().shape()[1]; n1++) + for (size_t n2 = 0; n2 < gt.data().shape()[2]; n2++) { + auto gw_sl = slice_target_to_scalar(gw, n1, n2); + auto gt_sl = slice_target_to_scalar(gt, n1, n2); + _fourier_impl(gw_sl, gt_sl); + } + } + + template + void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const &L) { + _fourier_impl(g, L.g); + } }} - - - -#endif - - diff --git a/triqs/gfs/local/fourier_lattice.cpp b/triqs/gfs/local/fourier_lattice.cpp new file mode 100644 index 00000000..cab3cb6f --- /dev/null +++ b/triqs/gfs/local/fourier_lattice.cpp @@ -0,0 +1,89 @@ +/******************************************************************************* + * + * 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 . + * + ******************************************************************************/ +#include "./fourier_lattice.hpp" +#include + +#define ASSERT_EQUAL(X,Y,MESS) if (X!=Y) TRIQS_RUNTIME_ERROR << MESS; + +namespace triqs { +namespace gfs { + + // We rewrite the scalar in term of the matrix + // We will change this for other FFT later when this FFT on lattice is checked. + void _fourier_impl(gf_view gk, gf_const_view gr) { + _fourier_impl(reinterpret_scalar_valued_gf_as_matrix_valued(gk), reinterpret_scalar_valued_gf_as_matrix_valued(gr)); + } + //-------------------------------------------------------------------------------------- + + void _fourier_impl(gf_view gr, gf_const_view gk) { + _fourier_impl(reinterpret_scalar_valued_gf_as_matrix_valued(gr), reinterpret_scalar_valued_gf_as_matrix_valued(gk)); + } + + //-------------------------------------------------------------------------------------- + + // The implementation is almost the same in both cases... + template + void __impl(int sign, gf_view g_out, gf_const_view g_in) { + + ASSERT_EQUAL(g_out.data().shape(), g_in.data().shape(), "Meshes are different"); + ASSERT_EQUAL(g_out.data().indexmap().strides()[1], g_out.data().shape()[1], "Unexpected strides in fourier implementation"); + ASSERT_EQUAL(g_out.data().indexmap().strides()[2], 1, "Unexpected strides in fourier implementation"); + ASSERT_EQUAL(g_in.data().indexmap().strides()[1], g_in.data().shape()[1], "Unexpected strides in fourier implementation"); + ASSERT_EQUAL(g_in.data().indexmap().strides()[2], 1, "Unexpected strides in fourier implementation"); + + auto L = g_in.mesh().get_dimensions(); + auto rank = g_in.mesh().rank(); + auto in_FFT = reinterpret_cast(g_in.data().data_start()); + auto outFFT = reinterpret_cast(g_out.data().data_start()); + // auto p = fftw_plan_dft(rank, L.ptr(), in_FFT, outFFT, FFTW_BACKWARD, FFTW_ESTIMATE); + + // use the general routine that can do all the matrices at once. + auto p = fftw_plan_many_dft(rank, // rank + L.ptr(), // the dimension + g_in.data().shape()[1] * g_in.data().shape()[2], // how many FFT : here 1 + in_FFT, // in data + NULL, // embed : unused. Doc unclear ? + g_in.data().indexmap().strides()[0], // stride of the in data + 1, // in : shift for multi fft. + outFFT, // out data + NULL, // embed : unused. Doc unclear ? + g_out.data().indexmap().strides()[0], // stride of the out data + 1, // out : shift for multi fft. + sign, FFTW_ESTIMATE); + + + fftw_execute(p); + fftw_destroy_plan(p); + } + + //-------------------------------------------------------------------------------------- + + void _fourier_impl(gf_view gk, gf_const_view gr) { + __impl(FFTW_BACKWARD, gk, gr); + } + + //-------------------------------------------------------------------------------------- + void _fourier_impl(gf_view gr, gf_const_view gk) { + __impl(FFTW_FORWARD, gr, gk); + gr.data() /= gk.mesh().size(); + } +} +} diff --git a/triqs/gfs/local/fourier_lattice.hpp b/triqs/gfs/local/fourier_lattice.hpp new file mode 100644 index 00000000..ddb7fc68 --- /dev/null +++ b/triqs/gfs/local/fourier_lattice.hpp @@ -0,0 +1,48 @@ +/******************************************************************************* + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include "fourier_base.hpp" +#include +#include + +namespace triqs { +namespace gfs { + + template + gf_keeper + fourier(gf_impl const& g) { + return {g}; + } + + template + gf_keeper inverse_fourier(gf_impl const& g) { + return {g}; + } + + // implementation + void _fourier_impl(gf_view gk, gf_const_view gr); + void _fourier_impl(gf_view gk, gf_const_view gr); + + void _fourier_impl(gf_view gr, gf_const_view gk); + void _fourier_impl(gf_view gr, gf_const_view gk); + +}} + diff --git a/triqs/gfs/local/fourier_matsubara.cpp b/triqs/gfs/local/fourier_matsubara.cpp index 3560ba60..de776de7 100644 --- a/triqs/gfs/local/fourier_matsubara.cpp +++ b/triqs/gfs/local/fourier_matsubara.cpp @@ -18,7 +18,6 @@ * TRIQS. If not, see . * ******************************************************************************/ -#include "fourier_base.hpp" #include "fourier_matsubara.hpp" #include @@ -45,13 +44,13 @@ namespace gfs { //------------------------------------- void direct(gf_view gw, gf_const_view gt) { - auto ta = gt(freq_infty()); + auto ta = gt.singularity(); direct_impl(make_gf_view_without_tail(gw), make_gf_view_without_tail(gt), ta); gw.singularity() = gt.singularity(); // set tail } void direct(gf_view gw, gf_const_view gt) { - auto ta = local::tail{1,1}; + auto ta = tail{1,1}; direct_impl(gw, gt, ta); } @@ -59,7 +58,7 @@ namespace gfs { private: void direct_impl(gf_view gw, gf_const_view gt, - local::tail const& ta) { + tail const& ta) { // TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO dcomplex d = ta(1)(0, 0), A = ta.get_or_zero(2)(0, 0), B = ta.get_or_zero(3)(0, 0); double b1 = 0, b2 = 0, b3 = 0; @@ -112,7 +111,7 @@ namespace gfs { static bool Green_Function_Are_Complex_in_time = false; // If the Green function are NOT complex, then one use the symmetry property // fold the sum and get a factor 2 - auto ta = gw(freq_infty()); + auto ta = gw.singularity(); // TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO dcomplex d = ta(1)(0, 0), A = ta.get_or_zero(2)(0, 0), B = ta.get_or_zero(3)(0, 0); double b1, b2, b3; @@ -180,65 +179,22 @@ namespace gfs { //-------------------------------------------- - template - void fourier_impl(gf_view gw, gf_const_view gt) { + // Direct transformation imtime -> imfreq, with a tail + void _fourier_impl(gf_view gw, gf_const_view gt) { impl_worker w; w.direct(gw, gt); } - template - void fourier_impl(gf_view gw, gf_const_view gt) { + void _fourier_impl(gf_view gw, gf_const_view gt) { impl_worker w; - for (size_t n1 = 0; n1 < gt.data().shape()[1]; n1++) - for (size_t n2 = 0; n2 < gt.data().shape()[2]; n2++) { - auto gw_sl = slice_target_to_scalar(gw, n1, n2); - auto gt_sl = slice_target_to_scalar(gt, n1, n2); - w.direct(gw_sl, gt_sl); - } + w.direct(gw, gt); } - //--------------------------------------------------------------------------- - - void inverse_fourier_impl(gf_view gt, gf_const_view gw) { + // Inverse transformation imfreq -> imtime: tail is mandatory + void _fourier_impl(gf_view gt, gf_const_view gw) { impl_worker w; w.inverse(gt, gw); } - - void inverse_fourier_impl(gf_view gt, gf_const_view gw) { - impl_worker w; - for (size_t n1 = 0; n1 < gw.data().shape()[1]; n1++) - for (size_t n2 = 0; n2 < gw.data().shape()[2]; n2++) { - auto gt_sl = slice_target_to_scalar(gt, n1, n2); - auto gw_sl = slice_target_to_scalar(gw, n1, n2); - w.inverse(gt_sl, gw_sl); - } - } - - //--------------------------------------------------------------------------- - void triqs_gf_view_assign_delegation(gf_view g, - gf_keeper const& L) { - fourier_impl(g, L.g); - } - void triqs_gf_view_assign_delegation(gf_view g, - gf_keeper const& L) { - fourier_impl(g, L.g); - } - void triqs_gf_view_assign_delegation(gf_view g, - gf_keeper const& L) { - inverse_fourier_impl(g, L.g); - } - void triqs_gf_view_assign_delegation(gf_view g, - gf_keeper const& L) { - inverse_fourier_impl(g, L.g); - } - void triqs_gf_view_assign_delegation(gf_view g, - gf_keeper const& L) { - fourier_impl(g, L.g); - } - void triqs_gf_view_assign_delegation(gf_view g, - gf_keeper const& L) { - fourier_impl(g, L.g); - } } } diff --git a/triqs/gfs/local/fourier_matsubara.hpp b/triqs/gfs/local/fourier_matsubara.hpp index de3e49ca..62575b69 100644 --- a/triqs/gfs/local/fourier_matsubara.hpp +++ b/triqs/gfs/local/fourier_matsubara.hpp @@ -27,6 +27,8 @@ namespace triqs { namespace gfs { + // only a few functions allowed: + template gf_keeper fourier(gf_impl const& g) { return {g}; @@ -36,18 +38,12 @@ namespace gfs { return {g}; } - // The fourier transform with the tail - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); - - // The version without tail : only possible in one direction - void triqs_gf_view_assign_delegation(gf_view g, - gf_keeper const& L); - void triqs_gf_view_assign_delegation(gf_view g, - gf_keeper const& L); + /// + void _fourier_impl(gf_view gw, gf_const_view gt); + void _fourier_impl(gf_view gw, gf_const_view gt); + void _fourier_impl(gf_view gt, gf_const_view gw); + /// A few helper functions template gf make_gf_from_fourier(gf_impl const& gt, int n_iw) { auto m = gf_mesh{gt.mesh().domain(), n_iw}; diff --git a/triqs/gfs/local/fourier_real.cpp b/triqs/gfs/local/fourier_real.cpp index aa600eb0..14bd2613 100644 --- a/triqs/gfs/local/fourier_real.cpp +++ b/triqs/gfs/local/fourier_real.cpp @@ -49,7 +49,7 @@ namespace triqs { namespace gfs { //a is a number very larger than delta_w and very smaller than wmax-wmin, used in the tail computation const double a = gw.mesh().delta() * sqrt( double(L) ); - auto ta = gt(freq_infty()); + auto ta = gt.singularity(); g_in.resize(L); g_in() = 0; g_out.resize(L); @@ -82,7 +82,7 @@ namespace triqs { namespace gfs { //a is a number very larger than delta_w and very smaller than wmax-wmin, used in the tail computation const double a = gw.mesh().delta() * sqrt( double(L) ); - auto ta = gw(freq_infty()); + auto ta = gw.singularity(); arrays::vector g_in(L), g_out(L); dcomplex t1 = ta(1)(0,0), t2 = ta.get_or_zero(2)(0,0); @@ -104,47 +104,17 @@ namespace triqs { namespace gfs { } }; - - + //-------------------------------------------------------------------------------------- - - void fourier_impl(gf_view gw, gf_const_view gt, scalar_valued){ + + void _fourier_impl(gf_view gw, gf_const_view gt) { impl_worker w; w.direct(gw, gt); } - - void fourier_impl(gf_view gw, gf_const_view gt, matrix_valued){ + + void _fourier_impl(gf_view gt, gf_const_view gw) { impl_worker w; - for (size_t n1=0; n1 gt, gf_const_view gw, scalar_valued){ - impl_worker w; - w.inverse(gt,gw); - } - - void inverse_fourier_impl (gf_view gt, gf_const_view gw, matrix_valued){ - impl_worker w; - for (size_t n1=0; n1 g, gf_keeper const & L) { fourier_impl (g,L.g, matrix_valued());} - void triqs_gf_view_assign_delegation( gf_view g, gf_keeper const & L) { fourier_impl (g,L.g, scalar_valued());} - void triqs_gf_view_assign_delegation( gf_view g, gf_keeper const & L) { inverse_fourier_impl(g,L.g, matrix_valued());} - void triqs_gf_view_assign_delegation( gf_view g, gf_keeper const & L) { inverse_fourier_impl(g,L.g, scalar_valued());} - -}} +} +} diff --git a/triqs/gfs/local/fourier_real.hpp b/triqs/gfs/local/fourier_real.hpp index 9ddf4166..f450075d 100644 --- a/triqs/gfs/local/fourier_real.hpp +++ b/triqs/gfs/local/fourier_real.hpp @@ -34,11 +34,10 @@ namespace triqs { namespace gfs { return {g}; } - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); - void triqs_gf_view_assign_delegation(gf_view g, gf_keeper const& L); + void _fourier_impl(gf_view gw, gf_const_view gt); + void _fourier_impl(gf_view gt, gf_const_view gw); + // helper functions template gf_mesh make_mesh_fourier_compatible(gf_mesh const& m) { int L = m.size(); double pi = std::acos(-1); diff --git a/triqs/gfs/local/functions.cpp b/triqs/gfs/local/functions.cpp index cf9ff46e..132fa3dd 100644 --- a/triqs/gfs/local/functions.cpp +++ b/triqs/gfs/local/functions.cpp @@ -33,7 +33,7 @@ namespace triqs { namespace gfs { dcomplex I(0,1); auto sh = G.data().shape().front_pop(); auto Beta = G.domain().beta; - local::tail_view t = G(freq_infty()); + tail_view t = G.singularity(); if (!t.is_decreasing_at_infinity()) TRIQS_RUNTIME_ERROR<<" density computation : Green Function is not as 1/omega or less !!!"; const size_t N1=sh[0], N2 = sh[1]; arrays::array dens_part(sh), dens_tail(sh), dens(sh); @@ -89,11 +89,11 @@ namespace triqs { namespace gfs { // compute a tail from the Legendre GF // this is Eq. 8 of our paper - local::tail_view get_tail(gf_const_view gl, int size = 10, int omin = -1) { + tail_view get_tail(gf_const_view gl, int size = 10, int omin = -1) { auto sh = gl.data().shape().front_pop(); - local::tail t(sh, size, omin); - t.data() = 0.0; + tail t(sh, size, omin); + t.data()() = 0.0; for (int p=1; p<=t.order_max(); p++) for (auto l : gl.mesh()) diff --git a/triqs/gfs/local/functions.hpp b/triqs/gfs/local/functions.hpp index e68016ae..c6d62999 100644 --- a/triqs/gfs/local/functions.hpp +++ b/triqs/gfs/local/functions.hpp @@ -37,7 +37,7 @@ namespace triqs { arrays::matrix density(gf_view const & g); - local::tail_view get_tail(gf_const_view gl, int size, int omin); + tail_view get_tail(gf_const_view gl, int size, int omin); void enforce_discontinuity(gf_view & gl, triqs::arrays::array_view disc); diff --git a/triqs/gfs/local/no_tail.hpp b/triqs/gfs/local/no_tail.hpp index f23f7251..b6bc4a1e 100644 --- a/triqs/gfs/local/no_tail.hpp +++ b/triqs/gfs/local/no_tail.hpp @@ -48,15 +48,15 @@ namespace gfs { } template - gf_view make_gf_from_g_and_tail(gf_impl const &g, local::tail t) { + gf_view make_gf_from_g_and_tail(gf_impl const &g, tail t) { details::_equal_or_throw(t.shape(), get_target_shape(g)); auto g2 = gf{g}; // copy the function without tail return {std::move(g2.mesh()), std::move(g2.data()), std::move(t), g2.symmetry()}; } template - gf_view make_gf_view_from_g_and_tail(gf_impl const &g, - local::tail_view t) { + gf_view make_gf_view_from_g_and_tail(gf_impl const &g, + tail_view t) { details::_equal_or_throw(t.shape(), get_target_shape(g)); return {g.mesh(), g.data(), t, g.symmetry()}; } diff --git a/triqs/gfs/local/tail.cpp b/triqs/gfs/local/tail.cpp new file mode 100644 index 00000000..f5775643 --- /dev/null +++ b/triqs/gfs/local/tail.cpp @@ -0,0 +1,214 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2012 by M. Ferrero, O. Parcollet + * + * TRIQS is free software: you can redistribute it and/or modify it under the + * terms of the GNU General Public License as published by the Free Software + * Foundation, either version 3 of the License, or (at your option) any later + * version. + * + * TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY + * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more + * details. + * + * You should have received a copy of the GNU General Public License along with + * TRIQS. If not, see . + * + ******************************************************************************/ +#include "./tail.hpp" + +namespace triqs { +namespace gfs { + + // tail_const_view class + void tail_const_view::rebind(tail_const_view const &X) { + omin = X.omin; + _mask.rebind(X._mask); + _data.rebind(X._data); + } + + // tail_view class + void tail_view::rebind(tail_view const &X) { + omin = X.omin; + _mask.rebind(X._mask); + _data.rebind(X._data); + } + + tail_view &tail_view::operator=(std::complex x) { + _data() = 0.0; + mv_type(_data(-omin, ellipsis())) = x; + mask() = omin + size() - 1; + return *this; + } + + /// --- functions + + std::ostream &operator<<(std::ostream &out, tail_const_view x) { + if (x.data().is_empty()) return out << "empty tail" << std::endl; + out << "tail/tail_view: min/smallest/max = " << x.order_min() << " " << x.smallest_nonzero() << " " << x.order_max(); + for (int u = x.order_min(); u <= x.order_max(); ++u) out << "\n ... Order " << u << " = " << x(u); + return out; + } + + /// Evaluate the tail to sum_{n=order_min}^ordermax M_n/omega^n + matrix evaluate(tail_const_view t, dcomplex const &omega) { + auto r = arrays::matrix{t.shape()}; + r() = 0; + auto omin = t.order_min(); + auto omax = t.order_max(); // precompute since long to do... + auto _ = arrays::range{}; + for (int u = omax; u >= omin; --u) + r = r / omega + matrix_view{t.data()(u - omin, ellipsis())}; // need to make a matrix view because otherwise + is not defined + r /= pow(omega, omin); + return r; + } + + /// Evaluate the tail to sum_{n=order_min}^ordermax M_n/omega^n + tail compose(tail_const_view x, tail_const_view t) { + auto r = tail(x.shape(), x.size(), x.order_min()); // a new tail of same size, init to 0 + auto t_inv = inverse(t); + auto omin = x.order_min(); + auto omax = x.order_max(); // precompute since long to do... + if ((omin<=0) && (omax>=0)) r = t * x(0); + tail z = t; + for (int u = -1; u >= omin; --u) { + r += x(u) * z; + z = z * t; + } + z = t_inv; + for (int u = 1; u <= omax; ++u) { + r += x(u) * z; + z = z * t_inv; + } + return r; + } + + /// Factories + tail tail_omega(int N1, int N2, int size_, int order_min) { + tail t(N1, N2, size_, order_min); + t(-1) = 1; + return t; + } + tail_view tail_omega(tail::shape_type const &sh, int size_, int order_min) { return tail_omega(sh[0], sh[1], size_, order_min); } + tail_view tail_omega(tail_view t) { return tail_omega(t.shape(), t.size(), t.order_min()); } + + // Ops + tail conj(tail_const_view const &t) { + return {conj(t.data()), t.mask(), t.order_min()}; + } + + tail transpose(tail_const_view const &t) { + return {transposed_view(t.data(), 0, 2, 1), transposed_view(t.mask(), 1, 0), t.order_min()}; + } + + /// Slice in orbital space + tail_view slice_target(tail_view t, range R1, range R2) { + return {t.data()(range(), R1, R2), t.mask()(R1, R2), t.order_min()}; + } + + tail_const_view slice_target(tail_const_view const &t, range R1, range R2) { + return {t.data()(range(), R1, R2), t.mask()(R1, R2), t.order_min()}; + } + + // inverse + tail inverse(tail_const_view const &t) { + int omin1 = -t.smallest_nonzero(); + int omax1 = std::min(t.order_max() + 2 * omin1, t.order_min() + int(t.size()) - 1); + int si = omax1 - omin1 + 1; + + tail res(t.shape(), t.size(), t.order_min()); + res.mask()() = omax1; + + res(omin1) = inverse(t(-omin1)); + + // optimize for the case 1x1 + if (1 && (t.shape()==make_shape(1,1))) { + for (int n = 1; n < si; n++) { + for (int p = 0; p < n; p++) { + res(omin1 + n)(0,0) -= t(n - omin1 - p)(0,0) * res(omin1 + p)(0,0); + } + res(omin1 + n)(0,0) = res(omin1)(0,0) * res(omin1 + n)(0,0); + } + } + else { + for (int n = 1; n < si; n++) { + for (int p = 0; p < n; p++) { + res(omin1 + n) -= t(n - omin1 - p) * res(omin1 + p); + } + res(omin1 + n) = res(omin1) * res(omin1 + n); + //res(omin1 + n) = res(omin1) * make_clone(res(omin1 + n)); + } + } + return res; + } + + // Arithmetic ops + + tail operator*(matrix const &a, tail_const_view const &b) { + auto res = tail(first_dim(a), b.shape()[1], b.size(), b.order_min()); // no {} constructor to avoid narrowing:... + for (int n = res.order_min(); n <= res.order_max(); ++n) res(n) = a * b(n); + return res; + } + + tail operator*(tail_const_view const &a, matrix const &b) { + auto res = tail(a.shape()[0], second_dim(b), a.size(), a.order_min()); + for (int n = res.order_min(); n <= res.order_max(); ++n) res(n) = a(n) * b; + return res; + } + + tail operator*(tail_const_view const &l, tail_const_view const &r) { + if (l.shape()[1] != r.shape()[0] || l.order_min() != r.order_min() || l.size() != r.size()) + TRIQS_RUNTIME_ERROR << "tail multiplication: shape mismatch"; + // int omin1 = l.smallest_nonzero() + r.smallest_nonzero(); + int omax1 = std::min(std::min(r.order_max() + l.smallest_nonzero(), l.order_max() + r.smallest_nonzero()), + r.order_min() + int(r.size()) - 1); + // int si = omax1-omin1+1; + tail res(l.shape(), l.size(), l.order_min()); + res.mask()() = omax1; + for (int n = res.order_min(); n <= res.order_max(); ++n) { + // sum_{p}^n a_p b_{n-p}. p <= a.n_max, p >= a.n_min and n-p <=b.n_max and n-p >= b.n_min + // hence p <= min ( a.n_max, n-b.n_min ) and p >= max ( a.n_min, n- b.n_max) + const int pmin = std::max(l.smallest_nonzero(), n - r.order_max()); + const int pmax = std::min(l.order_max(), n - r.smallest_nonzero()); + for (int p = pmin; p <= pmax; ++p) { + res(n) += l(p) * r(n - p); + } + } + return res; + } + + tail operator*(dcomplex a, tail_const_view const &r) { + tail res(r); + res.data() *= a; + return res; + } + + tail operator/(tail_const_view const &r, dcomplex a) { + tail res(r); + res.data() /= a; + return res; + } + + tail operator+(tail_const_view const &l, tail_const_view const &r) { + if (l.shape() != r.shape() || l.order_min() != r.order_min() || (l.size() != r.size())) + TRIQS_RUNTIME_ERROR << "tail addition: shape mismatch" << l << r; + tail res(l.shape(), l.size(), l.order_min()); + res.mask()() = std::min(l.order_max(), r.order_max()); + for (int i = res.order_min(); i <= res.order_max(); ++i) res(i) = l(i) + r(i); + return res; + } + + tail operator-(tail_const_view const &l, tail_const_view const &r) { + if (l.shape() != r.shape() || l.order_min() != r.order_min() || (l.size() != r.size())) + TRIQS_RUNTIME_ERROR << "tail addition: shape mismatch" << l << r; + tail res(l.shape(), l.size(), l.order_min()); + res.mask()() = std::min(l.order_max(), r.order_max()); + for (int i = res.order_min(); i <= res.order_max(); ++i) res(i) = l(i) - r(i); + return res; + } + +} +} diff --git a/triqs/gfs/local/tail.hpp b/triqs/gfs/local/tail.hpp index ddd2ebbd..80bba127 100644 --- a/triqs/gfs/local/tail.hpp +++ b/triqs/gfs/local/tail.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2012 by M. Ferrero, O. Parcollet + * Copyright (C) 2012-2014 by M. Ferrero, O. Parcollet * * TRIQS is free software: you can redistribute it and/or modify it under the * terms of the GNU General Public License as published by the Free Software @@ -25,380 +25,302 @@ #include #include -namespace triqs { namespace gfs { namespace local { +namespace triqs { +namespace gfs { - namespace details { - static constexpr double small = 1.e-10; + using dcomplex = std::complex; + using arrays::array; + using arrays::array_view; + using arrays::matrix_view; + using arrays::matrix; + using arrays::make_shape; + using arrays::range; + using arrays::ellipsis; + using arrays::mini_vector; + + class tail; // the value class + class tail_view; // the view class + class tail_const_view; // the const view class + + // ---------------------- implementation -------------------------------- + + enum rvc_enum { + Regular, + View, + ConstView + }; + + template struct __get_rvc; + template using get_rvc = typename __get_rvc::type; + template struct __get_rvc { + using type = T; + }; + template struct __get_rvc { + using type = typename T::view_type; + }; + template struct __get_rvc { + using type = typename T::const_view_type; + }; + + /// A common implementation class. + template class tail_impl { + public: + TRIQS_MPI_IMPLEMENTED_VIA_BOOST; + using view_type = tail_view; + using const_view_type = tail_view; + using regular_type = tail; + + using data_type = get_rvc>; + using mask_type = get_rvc>; + + using mv_type = matrix_view; + using const_mv_type = matrix_view; + + data_type & data() { return _data; } + data_type const & data() const { return _data; } + mask_type & mask() { return _mask; } + mask_type const & mask() const { return _mask; } + + int order_min() const { return omin; } + int order_max() const { return min_element(_mask); } + size_t size() const { return _data.shape()[0]; } + int smallest_nonzero() const { + int om = omin; + while ((om < order_max()) && (max_element(abs(_data(om - omin, ellipsis()))) < 1.e-10)) om++; + return om; } - using arrays::matrix_view; - using arrays::matrix; - namespace tqa= triqs::arrays; //namespace tql= triqs::clef; namespace mpl= boost::mpl; - typedef std::complex dcomplex; + using shape_type = mini_vector; + shape_type shape() const { return shape_type(_data.shape()[1], _data.shape()[2]); } + bool is_decreasing_at_infinity() const { return (smallest_nonzero() >= 1); } - class tail; // the value class - class tail_view; // the view class + protected: + int omin = -1; + mask_type _mask; + data_type _data; - template struct LocalTail : mpl::false_{}; // a boolean trait to identify the objects modelling the concept LocalTail - template<> struct LocalTail : mpl::true_{}; - template<> struct LocalTail: mpl::true_{}; + tail_impl() = default; + tail_impl(int N1, int N2, int size_, int order_min) + : omin(order_min), _mask(arrays::make_shape(N1, N2)), _data(make_shape(size_, N1, N2)) { + _mask() = order_min + size_ - 1; + _data() = 0; + } + tail_impl(data_type const &d, mask_type const &m, int omin_) : omin(omin_), _mask(m), _data(d) {} + template tail_impl(tail_impl const &x) : omin(x.order_min()), _mask(x.mask()), _data(x.data()) {} - // a trait to find the scalar of the algebra i.e. the true scalar and the matrix ... - template struct is_scalar_or_element : mpl::or_< tqa::ImmutableMatrix, utility::is_in_ZRC > {}; + friend class tail; + friend class tail_view; + friend class tail_const_view; + friend tail compose(tail_const_view x, tail_const_view omega); - // ---------------------- implementation -------------------------------- + // operator = for views + void operator=(tail_const_view const &rhs); // implemented later + + public: - /// A common implementation class. Idiom: ValueView - template class tail_impl { - 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 data_regular_type; - typedef arrays::array_view data_view_type; - typedef typename mpl::if_c::type data_type; - - typedef arrays::array mask_regular_type; - typedef arrays::array_view mask_view_type; - typedef typename mpl::if_c::type mask_type; - - typedef arrays::matrix_view mv_type; - typedef arrays::matrix_view const_mv_type; - - data_view_type data() { return _data;} - const data_view_type data() const { return _data;} - mask_view_type mask_view() { return mask;} - const mask_view_type mask_view() const { return mask;} - - long order_min() const {return omin;} - long order_max() const {return min_element(mask);} - size_t size() const {return _data.shape()[0];} - long smallest_nonzero() const { - long om = omin; - while ((om < this->order_max()) && (max_element(abs(_data(om-omin,tqa::range(),tqa::range()))) < details::small)) om++; - return om; - } - - typedef tqa::mini_vector shape_type; - shape_type shape() const { return shape_type(_data.shape()[1], _data.shape()[2]);} - bool is_decreasing_at_infinity() const { return (smallest_nonzero() >=1);} - - protected: - - long omin; - mask_type mask; - data_type _data; - - // All constructors - tail_impl(): omin(0), mask(), _data() {} // all arrays of zero size (empty) - tail_impl(size_t N1, size_t N2, size_t size_, long order_min): - omin(order_min), mask(tqa::make_shape(N1,N2)), _data(tqa::make_shape(size_,N1,N2)) { - mask() = order_min+size_-1; - _data() = 0; - } - tail_impl(data_type const &d, mask_type const &m, long omin_) : omin(omin_), mask(m), _data(d) {} - tail_impl(tail_impl const & x): omin(x.omin), mask(x.mask), _data(x._data) {} - tail_impl(tail_impl const &) = default; - tail_impl(tail_impl &&) = default; - - friend class tail_impl; - - public: - - mv_type operator() (int n) { - if (n>this->order_max()) TRIQS_RUNTIME_ERROR<<" n > Max Order. n= "< Max Order. n= "< - void serialize(Archive & ar, const unsigned int version) { - ar & TRIQS_MAKE_NVP("omin",omin); - ar & TRIQS_MAKE_NVP("mask",mask); - ar & TRIQS_MAKE_NVP("data",_data); - } - - friend std::ostream & operator << (std::ostream & out, tail_impl const & x) { - if (x.data().is_empty()) return out << "empty tail"< const & x) { - _data() = 0.0; - mv_type(_data(-omin, tqa::range(), tqa::range())) = x; - mask() = omin+size()-1; - return *this; - } - - using B::operator(); // import all previously defined operator() for overloading - friend std::ostream & triqs_nvl_formal_print(std::ostream & out, tail_view const & x) { return out<<"tail_view"; } - - }; - - // ----------------------------- - // the regular class - class tail: public tail_impl { - typedef tail_impl B; - friend class tail_view; - public: - tail():B() {} - typedef tqa::mini_vector 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) : 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) {} - tail(tail &&) = default; - - // operator = for values - tail & operator = (tail_view const & rhs) { - omin = rhs.omin; - mask = rhs.mask; - _data = rhs._data; - return *this; - } - tail & operator = (tail const & rhs) { - omin = rhs.omin; - mask = rhs.mask; - _data = rhs._data; - return *this; - } - - using B::operator(); - - /// The simplest tail corresponding to omega - static tail_view omega(size_t N1, size_t N2, size_t size_, long order_min) { - tail t(N1, N2, size_, order_min); - t(-1) = 1; - return t; - } - - /// The simplest tail corresponding to omega, constructed from a shape for convenience - static tail_view omega(tail::shape_type const & sh, size_t size_, long order_min) { return omega(sh[0],sh[1],size_,order_min); } - - }; - - template void assign_from_expression(tail_view t,RHS const & rhs) { t = rhs( tail::omega(t.shape(),t.size(),t.order_min()) ); } - - inline void tail_view::rebind(tail const &X) { - omin = X.omin; - mask.rebind(X.mask); - _data.rebind(X._data); + mv_type operator()(int n) { + if (n > order_max()) TRIQS_RUNTIME_ERROR << " n > Max Order. n= " << n << ", Max Order = " << order_max(); + if (n < order_min()) TRIQS_RUNTIME_ERROR << " n < Min Order. n= " << n << ", Min Order = " << order_min(); + return _data(n - omin, ellipsis()); } - inline tail_view & tail_view::operator = (const tail & rhs) { - if ((_data.shape()[1] != rhs._data.shape()[1]) || (_data.shape()[2] != rhs._data.shape()[2]) || (omin != rhs.omin)) - TRIQS_RUNTIME_ERROR<<"tails are incompatible"; - mask = rhs.mask; - _data = rhs._data; - return *this; + const_mv_type operator()(int n) const { + if (n > order_max()) TRIQS_RUNTIME_ERROR << " n > Max Order. n= " << n << ", Max Order = " << order_max(); + if (n < order_min()) { + mv_type::regular_type r(shape()); + r() = 0; + return r; + } + return _data(n - omin, ellipsis()); } - inline tail conj(tail_view t) { return {conj(t.data()), t.mask_view(),t.order_min()};} - - inline tail transpose(tail_view t) { return {transposed_view(t.data(),0,2,1), transposed_view(t.mask_view(),1,0),t.order_min()};} - - /// Slice in orbital space - //template tail_view slice_target(tail_impl const & t, tqa::range R1, tqa::range R2) - inline tail_view slice_target(tail_view t, tqa::range R1, tqa::range R2) { - return tail_view(t.data()(tqa::range(),R1,R2), t.mask_view()(R1,R2), t.order_min()); + /// same as (), but if n is too large, then returns 0 instead of raising an exception + const_mv_type get_or_zero(int n) const { + if ((n > order_max()) || (n < order_min())) { + mv_type::regular_type r(shape()); + r() = 0; + return r; + } + return _data(n - omin, ellipsis()); } - inline tail inverse(tail_view const & t) { - long omin1 = - t.smallest_nonzero(); - long omax1 = std::min(t.order_max() + 2*omin1, t.order_min()+long(t.size())-1); - size_t si = omax1-omin1+1; + friend std::string get_triqs_hdf5_data_scheme(tail_impl const &g) { return "TailGf"; } - tail res(t.shape(), t.size(), t.order_min()); - res.mask_view() = omax1; - - res(omin1) = inverse(t(-omin1)); - - for (size_t n=1; n void serialize(Archive &ar, const unsigned int version) { + ar &TRIQS_MAKE_NVP("omin", omin); + ar &TRIQS_MAKE_NVP("mask", _mask); + ar &TRIQS_MAKE_NVP("data", _data); + } + }; + + // ----------------------------- + class tail_const_view : public tail_impl { + using B = tail_impl; + // friend class tail; + + public: + tail_const_view(B::data_type d, B::mask_type m, int order_min) : B(std::move(d), std::move(m), order_min) {} + tail_const_view(tail_const_view const &) = default; + tail_const_view(tail_const_view &&) = default; + tail_const_view(tail_impl const &t) : B(t) {} + tail_const_view(tail_impl const &t) : B(t) {} + + void rebind(tail_const_view const &X); + }; + + // ----------------------------- + class tail_view : public tail_impl { + using B = tail_impl; + // friend class tail; + + public: + tail_view(B::data_type d, B::mask_type m, int order_min) : B(std::move(d), std::move(m), order_min) {} + tail_view(tail_view const &) = default; + tail_view(tail_view &&) = default; + tail_view(tail_impl const &t) : B(t) {} + + void rebind(tail_view const &X); + + tail_view &operator=(tail_view const &x) { return B::operator=(tail_const_view(x)), *this; } + template tail_view &operator=(tail_impl const &x) { return B::operator=(tail_const_view(x)), *this; } + tail_view &operator=(dcomplex); + }; + + // ----------------------------- + // the regular class + class tail : public tail_impl { + using B = tail_impl; + // friend class tail_view; + + public: + using shape_type = mini_vector; - inline tail inverse(tail const & t) { return inverse(tail_view(t));} + tail() = default; + tail(int N1, int N2, int size_ = 10, int order_min = -1) : B(N1, N2, size_, order_min) {} + tail(shape_type const &sh, int size_ = 10, int order_min = -1) : B(sh[0], sh[1], size_, order_min) {} + tail(mini_vector) : tail(1, 1) {} + tail(B::data_type const &d, B::mask_type const &m, int order_min) : B(d, m, order_min) {} + tail(tail const &g) : B(g) {} + tail(tail_view const &g) : B(g) {} + tail(tail_const_view const &g) : B(g) {} + tail(tail &&) = default; - inline tail mult_impl(tail_view const & l, tail_view const& r) { - if (l.shape()[1] != r.shape()[0] || l.order_min() != r.order_min() || l.size() != r.size()) - TRIQS_RUNTIME_ERROR<< "tail multiplication: shape mismatch"; - //long omin1 = l.smallest_nonzero() + r.smallest_nonzero(); - long omax1 = std::min(std::min(r.order_max()+l.smallest_nonzero(), l.order_max()+r.smallest_nonzero()), r.order_min()+long(r.size())-1); - //size_t si = omax1-omin1+1; + tail &operator=(tail const &x) { return B::operator=(tail_const_view(x)), *this; } + template tail &operator=(tail_impl const &x) { return B::operator=(tail_const_view(x)), *this; } + }; - tail res(l.shape(), l.size(), l.order_min()); - res.mask_view() = omax1; - - for (long n=res.order_min(); n<=res.order_max(); ++n) { - // sum_{p}^n a_p b_{n-p}. p <= a.n_max, p >= a.n_min and n-p <=b.n_max and n-p >= b.n_min - // hence p <= min ( a.n_max, n-b.n_min ) and p >= max ( a.n_min, n- b.n_max) - const long pmin = std::max(l.smallest_nonzero(), n - r.order_max() ); - const long pmax = std::min(l.order_max(), n - r.smallest_nonzero() ); - for (long p = pmin; p <= pmax; ++p) { res(n) += l(p) * r(n-p);} - } - return res; + // ---- impl. of = for base class ------------- + template void tail_impl::operator=(tail_const_view const &rhs) { + if (RVC == Regular) { + omin = rhs.order_min(); + } else { + if ((_data.shape()[1] != rhs.data().shape()[1]) || (_data.shape()[2] != rhs.data().shape()[2]) || (omin != rhs.order_min())) + TRIQS_RUNTIME_ERROR << "tails are incompatible"; + } + _mask = rhs.mask(); + _data = rhs.data(); } - template - TYPE_ENABLE_IF(tail, mpl::and_, LocalTail>) operator*(T1 const &a, T2 const &b) { - return mult_impl(a, b); - } + // ---- Factories ------------- - template - TYPE_ENABLE_IF(tail, mpl::and_, LocalTail>) operator*(T1 const &a, T2 const &b) { - auto res = tail{first_dim(a), b.shape()[1], b.size(), b.order_min()}; - res.mask_view() = b.order_max(); - for (int n = res.order_min(); n <= res.order_max(); ++n) res(n) = a * b(n); - return res; - } + /// The simplest tail corresponding to omega + tail tail_omega(int N1, int N2, int size_ = 10, int order_min = -1); - template - TYPE_ENABLE_IF(tail, mpl::and_, tqa::ImmutableMatrix>) operator*(T1 const &a, T2 const &b) { - auto res = tail{a.shape()[0], second_dim(b), a.size(), a.order_min()}; - res.mask_view() = a.order_max(); - for (int n = res.order_min(); n <= res.order_max(); ++n) res(n) = a(n) * b; - return res; - } + /// The simplest tail corresponding to omega, constructed from a shape for convenience + tail_view tail_omega(tail::shape_type const &sh, int size_, int order_min); - inline tail operator * (dcomplex a, tail_view const & r) { tail res(r); res.data()*=a; return res;} - inline tail operator * (tail_view const & r, dcomplex a) { return a*r; } + /// The simplest tail corresponding to omega, built from the shape, size, ordermin of t + tail_view tail_omega(tail_view t); - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator/ (T1 const & a, T2 const & b) { return a * inverse(b); } + // ---- IO ------------- - inline tail operator / (tail_view const & r, dcomplex a) { tail res(r); res.data() /=a; return res;} - inline tail operator / (dcomplex a, tail_view const & r) { return a * inverse(r); } + std::ostream &operator<<(std::ostream &out, tail_const_view); - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator + (T1 const & l, T2 const& r) { - if (l.shape() != r.shape() || l.order_min() != r.order_min() || (l.size() != r.size())) - TRIQS_RUNTIME_ERROR<< "tail addition: shape mismatch"; - tail res(l.shape(), l.size(), l.order_min()); - res.mask_view() = std::min(l.order_max(), r.order_max()); - for (long i = res.order_min(); i<=res.order_max(); ++i) res(i) = l(i) + r(i); - return res; - } + // ---- ------------- - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator - (T1 const & l, T2 const& r) { - if (l.shape() != r.shape() || l.order_min() != r.order_min() || (l.size() != r.size())) - TRIQS_RUNTIME_ERROR<< "tail addition: shape mismatch"; - tail res(l.shape(), l.size(), l.order_min()); - res.mask_view() = std::min(l.order_max(), r.order_max()); - for (long i = res.order_min(); i<=res.order_max(); ++i) res(i) = l(i) - r(i); - return res; - } + // + template void assign_singularity_from_function(tail_view t, RHS const &rhs) { + t = rhs(tail_omega(t.shape(), t.size(), t.order_min())); + } - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator + (T1 const & a, T2 const & t) { - tail res(t); - res(0) += a; - return res; - } + // ---- Evaluate and compose ------------- - template TYPE_ENABLE_IF(tail,mpl::and_, is_scalar_or_element>) - operator + (T1 const & t, T2 const & a) { return a+t;} + /// Evaluate the tail to sum_{n=order_min}^ordermax M_n/omega^n + matrix evaluate(tail_const_view t, dcomplex const &omega); - template TYPE_ENABLE_IF(tail,mpl::and_, LocalTail>) - operator - (T1 const & a, T2 const & t) { return (-a) + t;} + tail compose(tail_const_view x, tail_const_view t); + + // ---- Ops ------------- + + tail conj(tail_const_view const &t); + tail transpose(tail_const_view const &t); + + /// Slice in orbital space + tail_const_view slice_target(tail_const_view const &t, range R1, range R2); + tail_view slice_target(tail_view t, range R1, range R2); + inline tail_view slice_target(tail &t, range R1, range R2) { return slice_target(tail_view{t}, R1, R2);} + inline tail_const_view slice_target(tail const &t, range R1, range R2) { return slice_target(tail_const_view{t}, R1, R2);} + + tail inverse(tail_const_view const &t); + + /// ------------------- Arithmetic operations ------------------------------ + + tail operator+(tail_const_view const &l, tail_const_view const &r); + tail operator-(tail_const_view const &l, tail_const_view const &r); + + // +/- with scalar or matrix. + template + struct is_scalar_or_matrix : std::integral_constant>::value || utility::is_in_ZRC::value> {}; + + template std14::enable_if_t::value, tail> operator+(T const &a, tail_const_view const &t) { + tail res(t); + res(0) += a; + return res; + } + + template std14::enable_if_t::value, tail> operator+(tail_const_view const &t, T const &a) { + return a + t; + } + + template std14::enable_if_t::value, tail> operator-(T const &a, tail_const_view const &t) { + return (-a) + t; + } + + template std14::enable_if_t::value, tail> operator-(tail_const_view const &t, T const &a) { + return (-a) + t; + } + + tail operator*(tail_const_view const &l, tail_const_view const &r); + + tail operator*(matrix const &a, tail_const_view const &b); + tail operator*(tail_const_view const &a, matrix const &b); + template tail operator*(matrix_view const &a, tail_const_view const &b) { return matrix(a) * b; } + template tail operator*(tail_const_view const &a, matrix_view const &b) { return a * matrix(b); } + + tail operator*(dcomplex a, tail_const_view const &r); + inline tail operator*(tail_const_view const &r, dcomplex a) { return a * r; } + + inline tail operator/(tail_const_view const &a, tail_const_view const &b) { return a * inverse(b); } + tail operator/(tail_const_view const &r, dcomplex a); + inline tail operator/(dcomplex a, tail_view r) { return a * inverse(r); } - template TYPE_ENABLE_IF(tail,mpl::and_, is_scalar_or_element>) - operator - (T1 const & t, T2 const & a) { return (-a) + t;} // inplace operators #define DEFINE_OPERATOR(OP1, OP2) \ @@ -411,6 +333,6 @@ namespace triqs { namespace gfs { namespace local { DEFINE_OPERATOR(/=, / ); #undef DEFINE_OPERATOR - -}} } +} + diff --git a/triqs/gfs/m_tail.hpp b/triqs/gfs/m_tail.hpp index b894ceb3..ba5bbfa6 100644 --- a/triqs/gfs/m_tail.hpp +++ b/triqs/gfs/m_tail.hpp @@ -22,83 +22,141 @@ #include "./tools.hpp" #include "./gf.hpp" #include "./local/tail.hpp" -#include "./meshes/discrete.hpp" namespace triqs { namespace gfs { + template + void assign_singularity_from_function(gf_impl &s, RHS const &rhs) { + auto t = tail_omega(s.get_from_linear_index(0)); + // a bit faster to first replace (some part of expression are precomputed). + clef::placeholder<0> x_; + //auto expr = rhs(x_, t); + //for (auto x : s.mesh()) s[x] = eval(expr, x_ = x); + for (auto x : s.mesh()) s[x] = rhs(x, t); + } + + /// --------------------------- singularity --------------------------------- + + template struct gf_default_singularity, Target> { + using S1 = typename gf_default_singularity::type; + using S2 = typename gf_default_singularity::type; + // type is nothing unless S1 is nothing and S2 is not, or 1 <-> 2 + using type = std14::conditional_t < is_nothing() && (!is_nothing()), gf, + std14::conditional_t < is_nothing() && (!is_nothing()), gf, nothing >> + ; + }; + namespace gfs_implementation { - /*template constexpr bool has_a_singularity() { - return std::is_same::type, local::tail>::value; - } - - /// --------------------------- singularity --------------------------------- - - template struct singularity, Target, Opt> { - using type = std14::conditional_t < (!has_a_singularity()) && has_a_singularity(), gf, nothing > ; - }; - */ /// --------------------------- hdf5 --------------------------------- - - template struct h5_name { + + template struct h5_name { static std::string invoke() { return "xxxxx"; } }; - template struct h5_rw { + template struct h5_rw { - static void write(h5::group gr, gf_const_view g) { - for (size_t i = 0; i < g.mesh().size(); ++i) h5_write(gr, std::to_string(i), g._data[i]); + static void write(h5::group gr, gf_const_view 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 static void read(h5::group gr, gf_impl &g) { + template static void read(h5::group gr, gf_impl &g) { // does not work : need to read the block name and remake the mesh... - //g._mesh = gf_mesh(gr.get_all_subgroup_names()); - //g._data.resize(g._mesh.size()); + // g._mesh = gf_mesh(gr.get_all_subgroup_names()); + // g._data.resize(g._mesh.size()); // if (g._data.size() != g._mesh.size()) TRIQS_RUNTIME_ERROR << "h5 read block gf : number of block mismatch"; - //for (size_t i = 0; i < g.mesh().size(); ++i) h5_read(gr, g.mesh().domain().names()[i], g._data[i]); + // 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 --------------------------------- + // --------------------------- data access --------------------------------- - template struct data_proxy : data_proxy_vector {}; + template struct data_proxy { - // ------------------------------- Factories -------------------------------------------------- + struct storage_t { + arrays::array data; + arrays::array mask; + int omin; + }; - template struct factories { - using mesh_t=gf_mesh ; - using gf_t=gf ; - //using gf_view_t=gf_view ; - struct target_shape_t {}; - using aux_t = nothing; + struct storage_view_t { + arrays::array_view data; + arrays::array_view mask; + int omin; + template storage_view_t(S &s): data(s.data), mask(s.mask), omin(s.omin){} + }; - static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t, aux_t) { return std::vector(m.size()); } - static nothing make_singularity(mesh_t const &m, target_shape_t) { - return {}; + struct storage_const_view_t { + arrays::array_const_view data; + arrays::array_const_view mask; + int omin; + template storage_const_view_t(S &s): data(s.data), mask(s.mask), omin(s.omin){} + }; + + // from the shape of the mesh and the target, make the shape of the array. default is to glue them + // template static auto join_shape(S1 const &s1, S2 const &s2) { return make_shape(join(s1, s2); } + + template static void assign_to_scalar(S &data, RHS &&rhs) { data() = std::forward(rhs); } + template static void rebind(ST &data, RHS &&rhs) { data.rebind(rhs.data()); } + + template tail_view operator()(S &data, int i) const { + return {data.data(i, arrays::ellipsis()), data.mask, data.omin}; + } + + template tail_const_view operator()(S const &data, int i) const { + return {data.data(i, arrays::ellipsis()), data.mask, data.omin}; } }; - // partial_eval - template struct partial_eval_impl { + // ------------------------------- Factory for data -------------------------------------------------- - using gv_t = gf_view; + template struct data_factory { + using mesh_t = gf_mesh; + using gf_t = gf; + using target_shape_t = arrays::mini_vector; + // struct target_shape_t {}; + using aux_t = nothing; + + static typename gf_t::data_t make(mesh_t const &m, target_shape_t sh, aux_t) { + auto t = tail(sh); // build a defaut tail + // and duplicate it over the mesh size + return {arrays::array{t.data().shape().front_append(m.size())}, arrays::array{t.shape()}, + t.order_min()}; + } + }; + + // ------------------------------- Factory for singularity -------------------------------------------------- + + template + struct singularity_factory, Target, gf, Opt> { + template + static gf make(gf_mesh, Opt> const &m, TargetShape shape) { + return {std::get<0>(m.components()), shape}; + } + }; + + // ------------------------------- partial_eval -------------------------------------------------- + + template + struct partial_eval_impl { + + using gv_t = gf_view; template static auto invoke(gv_t g, T const &... x) { return invoke_impl(g, std14::index_sequence(), x...); } - template 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 static tail_view invoke_impl(gv_t g, std14::index_sequence<0>, T const &x) { + return g.get_from_linear_index(x); } - template static nothing invoke_impl(gv_t g, std14::index_sequence<1>, T const &x) { - return {}; + template static gv_t invoke_impl(gv_t g, std14::index_sequence<1>, T const &x) { + return g; } }; } // gfs_implementation - }} diff --git a/triqs/gfs/meshes/linear.hpp b/triqs/gfs/meshes/linear.hpp index 2d8646b0..4a5bec3a 100644 --- a/triqs/gfs/meshes/linear.hpp +++ b/triqs/gfs/meshes/linear.hpp @@ -65,10 +65,10 @@ namespace gfs { } domain_t const &domain() const { return _dom; } - size_t size() const { return L; } + long size() const { return L; } utility::mini_vector size_of_components() const { - return {size()}; + return {size_t(size())}; } double delta() const { return del; } @@ -79,7 +79,6 @@ namespace gfs { /// Conversions point <-> index <-> linear_index domain_pt_t index_to_point(index_t ind) const { return xmin + ind * del; } - public: long index_to_linear(index_t ind) const { return ind; } /// The wrapper for the mesh point diff --git a/triqs/gfs/meshes/matsubara_freq.hpp b/triqs/gfs/meshes/matsubara_freq.hpp index 04ea9cbd..715f9256 100644 --- a/triqs/gfs/meshes/matsubara_freq.hpp +++ b/triqs/gfs/meshes/matsubara_freq.hpp @@ -150,7 +150,7 @@ namespace gfs { const_iterator cend() const { return const_iterator(this, true); } bool operator==(matsubara_freq_mesh const &M) const { - return (std::make_tuple(_dom, _n_pts, _positive_only) == std::make_tuple(M._dom, M._n_pts, M._positive_only)); + return (std::tie(_dom, _n_pts, _positive_only) == std::tie(M._dom, M._n_pts, M._positive_only)); } bool operator!=(matsubara_freq_mesh const &M) const { return !(operator==(M)); } diff --git a/triqs/gfs/product.hpp b/triqs/gfs/product.hpp index b9b0d308..e042a906 100644 --- a/triqs/gfs/product.hpp +++ b/triqs/gfs/product.hpp @@ -70,14 +70,14 @@ namespace gfs { /// --------------------------- hdf5 --------------------------------- // h5 name : name1_x_name2_..... - template struct h5_name, matrix_valued, nothing, Opt> { + template struct h5_name, matrix_valued, S, 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::invoke()...), std::string()); + reverse(std::make_tuple(h5_name::invoke()...)), std::string()); } }; - template - struct h5_name, tensor_valued, nothing, Opt> : h5_name, matrix_valued, nothing, Opt> {}; + template + struct h5_name, tensor_valued, S, Opt> : h5_name, matrix_valued, S, Opt> {}; /// --------------------------- evaluator --------------------------------- @@ -110,22 +110,27 @@ namespace gfs { template struct binder; template binder make_binder(G const *g, Tn tn, Ev const &ev) { - return binder{g, std::move(tn), ev}; + return {g, std::move(tn), ev}; } template struct binder { G const *g; Tn tn; Ev const &evals; - auto operator()(size_t p) const - DECL_AND_RETURN(std::get(evals)(make_binder(g, triqs::tuple::push_front(tn, p), evals))); + template decltype(auto) impl(size_t p, std14::index_sequence) const { + return std::get(evals)(make_binder(g, std::make_tuple(p, std::get(tn)...), evals)); + } + decltype(auto) operator()(size_t p) const { return impl(p, std14::make_index_sequence::value>()); } }; template struct binder { G const *g; Tn tn; Ev const &evals; - auto operator()(size_t p) const DECL_AND_RETURN(triqs::tuple::apply(on_mesh(*g), triqs::tuple::push_front(tn, p))); + template decltype(auto) impl(size_t p, std14::index_sequence) const { + return g->get_from_linear_index(p, std::get(tn)...); + } + decltype(auto) operator()(size_t p) const { return impl(p, std14::make_index_sequence::value>()); } }; // now the multi d evaluator itself. @@ -133,21 +138,11 @@ namespace gfs { static constexpr int arity = sizeof...(Ms); mutable std::tuple...> evals; - struct _poly_lambda { // replace by a polymorphic lambda in C++14 - template void operator()(A &a, B const &b, C const &c) const { - a = A{b, c}; - } - }; - - template - // std::complex operator() (G const * g, Args && ... args) const { - auto operator()(G const *g, Args &&... args) - const -> decltype(std::get(evals)(make_binder(g, std::make_tuple(), evals))) - // when do we get C++14 decltype(auto) ...!? - { + template decltype(auto) operator()(G const *g, Args &&... args) const { static constexpr int R = sizeof...(Args); // build the evaluators, as a tuple of ( evaluator ( mesh_component, args)) - triqs::tuple::for_each_zip(_poly_lambda(), evals, g->mesh().components(), std::make_tuple(args...)); + auto l = [](auto &a, auto &b, auto &c) { a = std14::decay_t{b, c}; }; + triqs::tuple::for_each_zip(l, evals, g->mesh().components(), std::make_tuple(args...)); return std::get(evals)(make_binder(g, std::make_tuple(), evals)); } }; diff --git a/triqs/gfs/refreq.hpp b/triqs/gfs/refreq.hpp index c11a96ee..fc8c00e1 100644 --- a/triqs/gfs/refreq.hpp +++ b/triqs/gfs/refreq.hpp @@ -37,10 +37,10 @@ namespace gfs { // singularity template <> struct gf_default_singularity { - using type = local::tail; + using type = tail; }; template <> struct gf_default_singularity { - using type = local::tail; + using type = tail; }; namespace gfs_implementation { diff --git a/triqs/gfs/retime.hpp b/triqs/gfs/retime.hpp index d0a0d264..d9c7fc43 100644 --- a/triqs/gfs/retime.hpp +++ b/triqs/gfs/retime.hpp @@ -37,10 +37,10 @@ namespace gfs { // singularity template <> struct gf_default_singularity { - using type = local::tail; + using type = tail; }; template <> struct gf_default_singularity { - using type = local::tail; + using type = tail; }; namespace gfs_implementation { diff --git a/triqs/gfs/tools.hpp b/triqs/gfs/tools.hpp index f9a58e63..fb491eac 100644 --- a/triqs/gfs/tools.hpp +++ b/triqs/gfs/tools.hpp @@ -55,8 +55,6 @@ namespace gfs { Fermion }; - struct freq_infty {}; // the point at infinity - //------------------------------------------------------ template struct closest_pt_wrap; @@ -171,6 +169,7 @@ namespace gfs { struct nothing { template explicit nothing(Args &&...) {} // takes anything, do nothing.. nothing() {} + using const_view_type = nothing; using view_type = nothing; using regular_type = nothing; void rebind(nothing) {} @@ -181,15 +180,19 @@ namespace gfs { friend class boost::serialization::access; template void serialize(Archive &ar, const unsigned int version) {} friend nothing operator+(nothing, nothing) { return nothing(); } - template friend void assign_from_expression(nothing &, RHS) {} + template friend void assign_singularity_from_function(nothing &, RHS) {} template bool check_size(A) {return true;} bool is_empty() const { return false;} }; + // Check if T is nothing + template constexpr bool is_nothing() { return std::is_same::value; } + template nothing partial_eval(nothing, T&&...) { return {};} inline nothing transpose(nothing) { return {};} inline nothing inverse(nothing) { return {};} inline nothing conj(nothing) { return {};} + template nothing compose(nothing,T&) { return {};} template nothing slice_target(nothing, T...) { return nothing(); } template nothing operator+(nothing, T const &) { return nothing(); } template nothing operator-(nothing, T const &) { return nothing(); } diff --git a/triqs/h5.hpp b/triqs/h5.hpp index 75e06005..78d486e0 100644 --- a/triqs/h5.hpp +++ b/triqs/h5.hpp @@ -40,6 +40,13 @@ namespace h5 { template struct has_h5_write(), std::string(), std::declval()))> : std::true_type { }; + +// A generic read + template T h5_read(group gr, std::string const& name) { + T x; + h5_read(gr, name, x); + return x; + } } } diff --git a/triqs/h5/vector.cpp b/triqs/h5/vector.cpp index 94c21628..b6117a20 100644 --- a/triqs/h5/vector.cpp +++ b/triqs/h5/vector.cpp @@ -131,8 +131,11 @@ namespace h5 { } // impl namespace + void h5_write(group f, std::string const &name, std::vector const &V) { h5_write_vector_impl(f, name, V); } void h5_write(group f, std::string const &name, std::vector const &V) { h5_write_vector_impl(f, name, V); } void h5_write(group f, std::string const &name, std::vector> const &V) { h5_write_vector_impl(f, name, V); } + + void h5_read(group f, std::string const &name, std::vector &V) { h5_read_impl(f, name, V); } void h5_read(group f, std::string const &name, std::vector &V) { h5_read_impl(f, name, V); } void h5_read(group f, std::string const &name, std::vector> &V) { h5_read_impl(f, name, V); } } diff --git a/triqs/h5/vector.hpp b/triqs/h5/vector.hpp index e1759fc7..413ecddb 100644 --- a/triqs/h5/vector.hpp +++ b/triqs/h5/vector.hpp @@ -38,9 +38,11 @@ namespace triqs { void h5_write (group f, std::string const & name, std::vector const & V); void h5_read (group f, std::string const & name, std::vector & V); + void h5_write (group f, std::string const & name, std::vector const & V); void h5_write (group f, std::string const & name, std::vector const & V); void h5_write (group f, std::string const & name, std::vector> const & V); + void h5_read (group f, std::string const & name, std::vector & V); void h5_read (group f, std::string const & name, std::vector & V); void h5_read (group f, std::string const & name, std::vector> & V); diff --git a/triqs/lattice/bz_mesh.hpp b/triqs/lattice/bz_mesh.hpp index bca6c980..058e7871 100644 --- a/triqs/lattice/bz_mesh.hpp +++ b/triqs/lattice/bz_mesh.hpp @@ -60,6 +60,7 @@ namespace lattice { /// The wrapper for the mesh point class mesh_point_t : gfs::tag::mesh_point, public utility::arithmetic_ops_by_cast { + public : bz_mesh const *m; index_t _index; diff --git a/triqs/lattice/cyclic_lattice.hpp b/triqs/lattice/cyclic_lattice.hpp new file mode 100644 index 00000000..0a57b446 --- /dev/null +++ b/triqs/lattice/cyclic_lattice.hpp @@ -0,0 +1,135 @@ +/******************************************************************************* + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include +#include +//#include +//#include + +namespace triqs { +namespace lattice { + + class cyclic_lattice_mesh { + + utility::mini_vector dims; // the size in each dimension + size_t _size = dims[0] * dims[1] * dims[2]; // total size + long s1 = dims[0]; // stride + long s2 = dims[0] * dims[1]; // stride + + long _modulo(long r, int i) const { + long res = r % dims[i]; + return (res >= 0 ? res : res + dims[i]); + } + + public: + cyclic_lattice_mesh(int L1 = 1, int L2 = 1, int L3 = 1) : dims{L1, L2, L3} {} + + int rank() const { return (dims[2] > 1 ? 3 : (dims[1] > 1 ? 2 : 1)); } + + utility::mini_vector get_dimensions() const { return dims; } + + /// ---------- Model the domain concept --------------------- + using point_t = arrays::vector; // domain concept. PUT on STACK + + /// ----------- Model the mesh concept ---------------------- + using domain_t = cyclic_lattice_mesh; + domain_t const& domain() const { return *this; } + + using index_t = utility::mini_vector; + using linear_index_t = long; + + size_t size() const { return _size; } + + utility::mini_vector size_of_components() const { + return {size()}; + } + + point_t index_to_point(index_t const& i) const { + return {i[0], i[1], i[2]}; // not very good. + } + + /// flatten the index + linear_index_t index_to_linear(index_t const& i) const { + return _modulo(i[0], 0) + _modulo(i[1], 1) * s1 + _modulo(i[2], 2) * s2; + } + // linear_index_t index_to_linear(index_t const& i) const { return i[0] + i[1] * s1 + i[2] * s2; } + + /// The wrapper for the mesh point + class mesh_point_t : public utility::index3_generator, public utility::arithmetic_ops_by_cast { + cyclic_lattice_mesh const* m = nullptr; + + public: + mesh_point_t() = default; + // explicit is important for g[ {1,2}] to disambiguate the mesh_point_t and the mesh_index_t + explicit mesh_point_t(cyclic_lattice_mesh const& mesh, index_t const& index) + : index3_generator(mesh.get_dimensions(), index), m(&mesh) {} + mesh_point_t(cyclic_lattice_mesh const& mesh) : mesh_point_t(mesh, {0, 0, 0}) {} + operator point_t() const { return m->index_to_point(index()); } + // linear_index_t linear_index() const { return m->index_to_linear(index()); } + // The mesh point behaves like a vector + long operator()(int i) const { return index()[i]; } + long operator[](int i) const { return operator()(i); } + friend std::ostream& operator<<(std::ostream& out, mesh_point_t const& x) { return out << x.index(); } + }; + + /// Accessing a point of the mesh + mesh_point_t operator[](index_t i) const { + return mesh_point_t{*this, i}; + } + + /// Iterating on all the points... + using const_iterator = gfs::mesh_pt_generator; + const_iterator begin() const { return const_iterator(this); } + const_iterator end() const { return const_iterator(this, true); } + const_iterator cbegin() const { return const_iterator(this); } + const_iterator cend() const { return const_iterator(this, true); } + + /// ----------- End mesh concept ---------------------- + + /// Reduce point modulo to the lattice. + mesh_point_t modulo_reduce(index_t const& r) const { + return mesh_point_t{*this, {_modulo(r[0], 0), _modulo(r[1], 1), _modulo(r[2], 2)}}; + } + + /// Write into HDF5 + friend void h5_write(h5::group fg, std::string subgroup_name, cyclic_lattice_mesh const& m) { + h5::group gr = fg.create_group(subgroup_name); + h5_write(gr, "dims", m.dims.to_vector()); + } + + /// Read from HDF5 + friend void h5_read(h5::group fg, std::string subgroup_name, cyclic_lattice_mesh& m) { + h5::group gr = fg.open_group(subgroup_name); + auto dims = h5::h5_read>(gr, "dims"); + m = cyclic_lattice_mesh(dims[0], dims[1], dims[2]); + } + + // BOOST Serialization + friend class boost::serialization::access; + template void serialize(Archive& ar, const unsigned int version) { + ar& TRIQS_MAKE_NVP("dims", dims); + ar& TRIQS_MAKE_NVP("_size", _size); + ar& TRIQS_MAKE_NVP("s2", s2); + ar& TRIQS_MAKE_NVP("s1", s1); + } + }; +} +} diff --git a/triqs/lattice/regular_bz_mesh.hpp b/triqs/lattice/regular_bz_mesh.hpp new file mode 100644 index 00000000..f75f32da --- /dev/null +++ b/triqs/lattice/regular_bz_mesh.hpp @@ -0,0 +1,150 @@ +/******************************************************************************* + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include +#include +#include "./brillouin_zone.hpp" +#include "../gfs/tools.hpp" +#include "../gfs/meshes/mesh_tools.hpp" + +namespace triqs { +namespace lattice { + + class regular_bz_mesh { + + brillouin_zone bz; // + utility::mini_vector dims = {1, 1, 1}; // the size in each dimension + size_t _size = dims[0] * dims[1] * dims[2]; // total size + long s1 = dims[0]; // stride + long s2 = dims[0] * dims[1]; // stride + utility::mini_vector step = {2 * M_PI / dims[0], 2 * M_PI / dims[1], 2 * M_PI / dims[2]}; + + public: + regular_bz_mesh() = default; + + regular_bz_mesh(brillouin_zone const& bz, int n_l) + : bz(bz), dims{n_l, (bz.lattice().dim() >= 2 ? n_l : 1), (bz.lattice().dim() >= 3 ? n_l : 1)} {} + + int rank() const { return (dims[2] > 1 ? 3 : (dims[1] > 1 ? 2 : 1)); } + + utility::mini_vector get_dimensions() const { return dims; } + + /// ----------- Model the mesh concept ---------------------- + using domain_t = brillouin_zone; + using domain_pt_t = typename domain_t::point_t; + domain_t const& domain() const { return bz; } + + using index_t = utility::mini_vector; + using linear_index_t = long; + + size_t size() const { return _size; } + + utility::mini_vector size_of_components() const { + return {size()}; + } + + k_t index_to_point(index_t const& i) const { + return {i[0] * step[0], i[1] * step[1], i[2] * step[2]}; + //return {(i[0] + 0.5) * step[0], (i[1] + 0.5) * step[1], (i[2] + 0.5) * step[2]}; + } + + // flatten the index + linear_index_t index_to_linear(index_t const& i) const { return i[0] + i[1] * s1 + i[2] * s2; } + + // f (k) -> void where k is a k_t, a point in the BZ + template friend void foreach(regular_bz_mesh const& m, F f) { + k_t k = {0,0,0}; //{0.5 * m.step[0], 0.5 * m.step[1], 0.5 * m.step[2]}; + for (long i2 = 0; i2 < m.dims[2]; ++i2, k(2) += m.step[2]) + for (long i1 = 0; i1 < m.dims[1]; ++i1, k(1) += m.step[1]) + for (long i0 = 0; i0 < m.dims[0]; ++i0, k(0) += m.step[0]) f(k); + } + + /// The wrapper for the mesh point + class mesh_point_t : public utility::index3_generator, public utility::arithmetic_ops_by_cast { + regular_bz_mesh const* m = nullptr; + + public: + mesh_point_t() = default; + mesh_point_t(regular_bz_mesh const& mesh, index_t const& index) : index3_generator(mesh.get_dimensions(), index), m(&mesh) {} + mesh_point_t(regular_bz_mesh const& mesh) : mesh_point_t(mesh, {0,0,0}) {} + operator domain_pt_t() const { return m->index_to_point(index()); } + linear_index_t linear_index() const { return m->index_to_linear(index()); } + // The mesh point behaves like a vector // NOT GOOD : take the ith componenet, this is SLOW + double operator()(int i) const { return index()[i] * m->step[i]; } + //double operator()(int i) const { return (index()[i] + 0.5) * m->step[i]; } + double operator[](int i) const { return operator()(i);} + friend std::ostream& operator<<(std::ostream& out, mesh_point_t const& x) { return out << (domain_pt_t)x; } + }; + + /// Accessing a point of the mesh + mesh_point_t operator[](index_t i) const { + return {*this, i}; + } + + /// Iterating on all the points... + using const_iterator = gfs::mesh_pt_generator; + const_iterator begin() const { return const_iterator(this); } + const_iterator end() const { return const_iterator(this, true); } + const_iterator cbegin() const { return const_iterator(this); } + const_iterator cend() const { return const_iterator(this, true); } + + /// ----------- End mesh concept ---------------------- + + /// locate the closest point + mesh_point_t locate_neighbours(k_t const& k) const { + auto l = [&](int i) { + long r = std::lround(k(i) / step[i]) % dims[i]; + return (r >= 0 ? r : r + dims[i]); + }; + return {*this, {l(0), l(1), l(2)}}; + } + + /// Write into HDF5 + friend void h5_write(h5::group fg, std::string subgroup_name, regular_bz_mesh const& m) { + h5::group gr = fg.create_group(subgroup_name); + h5_write(gr, "domain", m.domain()); + h5_write(gr, "n_pts", m.dims[2]); + //h5_write(gr, "dims", m.dims.to_vector()); + } + + /// Read from HDF5 + friend void h5_read(h5::group fg, std::string subgroup_name, regular_bz_mesh& m) { + h5::group gr = fg.open_group(subgroup_name); + auto bz = h5::h5_read(gr, "domain"); + auto dims = h5::h5_read(gr, "n_pts"); + m = regular_bz_mesh(bz, dims); + //auto dims = h5::h5_read>(gr, "dims"); + //m = regular_bz_mesh(bz, {dims[0], dims[1], dims[2]}); // NOT CORRECT IN GENERAL + } + + // BOOST Serialization + friend class boost::serialization::access; + template void serialize(Archive& ar, const unsigned int version) { + ar& TRIQS_MAKE_NVP("dims", dims); + ar& TRIQS_MAKE_NVP("_size", _size); + ar& TRIQS_MAKE_NVP("s2", s2); + ar& TRIQS_MAKE_NVP("s1", s1); + ar& TRIQS_MAKE_NVP("step", step); + } + }; +} +} + diff --git a/triqs/utility/index_generator.hpp b/triqs/utility/index_generator.hpp new file mode 100644 index 00000000..decf97cf --- /dev/null +++ b/triqs/utility/index_generator.hpp @@ -0,0 +1,61 @@ +/******************************************************************************* + * + * 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 . + * + ******************************************************************************/ +#pragma once +#include "./mini_vector.hpp" + +namespace triqs { +namespace utility { + + class index3_generator { + mini_vector i, d; // rely on mini_vector initialization + long i_flat =0; + bool _at_end = false; + + public: + index3_generator() = default; + index3_generator(mini_vector const& dims, mini_vector const &i) : d(dims),i(i) {} + void advance() { + ++i_flat; + ++i[2]; + if (i[2] < d[2]) return; + i[2] = 0; + ++i[1]; + if (i[1] < d[1]) return; + i[1] = 0; + ++i[0]; + if (i[0] < d[0]) return; + // i[0]=0; + _at_end = true; + } + AUTO_DECL index() const { return i; } + long linear_index() const { return i_flat;} + bool at_end() const { return _at_end; } + void reset() { + _at_end = false; + i_flat =0; + i[0] = 0; + i[1] = 0; + i[2] = 0; + } + }; +} +} +