From 17ab555213acbae2ae77a1d387030a72a1e515b7 Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Wed, 20 Nov 2013 15:11:29 +0100 Subject: [PATCH] gf: fix matsubara_freq and mesh_pt - cleaner implementation, solve some bugs... --- doc/reference/c++/parameters/parameters.rst | 10 +-- test/triqs/arrays/mapped_functions.cpp | 1 - test/triqs/arrays/mapped_functions.output | 6 +- triqs/arrays/mapped_functions.hpp | 7 ++- triqs/gfs/domains/matsubara.hpp | 36 ++--------- triqs/gfs/imfreq.hpp | 21 ++++++- triqs/gfs/meshes/discrete.hpp | 2 +- triqs/gfs/meshes/linear.hpp | 2 +- triqs/gfs/meshes/matsubara_freq.hpp | 23 +++---- triqs/gfs/meshes/mesh_tools.hpp | 27 +------- triqs/utility/arithmetic_ops_by_cast.hpp | 69 +++++++++++++++++++++ 11 files changed, 114 insertions(+), 90 deletions(-) create mode 100644 triqs/utility/arithmetic_ops_by_cast.hpp diff --git a/doc/reference/c++/parameters/parameters.rst b/doc/reference/c++/parameters/parameters.rst index 50ce6cd1..17a43b98 100644 --- a/doc/reference/c++/parameters/parameters.rst +++ b/doc/reference/c++/parameters/parameters.rst @@ -60,13 +60,9 @@ required and optional parameters. Default values are provided for optional param .. code-block:: c parameter_defaults pdef; - pdef.required - ( "Beta", double(), "Inverse temperature") - ; - - pdef.optional - ( "N_time", int(100), "Number of time points") - ( "N_freq", int(200), "Number of frequencies") + pdef.required("Beta", double(), "Inverse temperature") + .optional("N_time", int(100), "Number of time points") + .optional("N_freq", int(200), "Number of frequencies") ; The default_parameters object serves two main purposes: Firstly, the input parameters can diff --git a/test/triqs/arrays/mapped_functions.cpp b/test/triqs/arrays/mapped_functions.cpp index a506b563..77d95f58 100644 --- a/test/triqs/arrays/mapped_functions.cpp +++ b/test/triqs/arrays/mapped_functions.cpp @@ -62,7 +62,6 @@ template void test( T val=1 ) { int main(int argc, char **argv) { - conj (8); test(); test(); test(); diff --git a/test/triqs/arrays/mapped_functions.output b/test/triqs/arrays/mapped_functions.output index 5843c2db..39c78292 100644 --- a/test/triqs/arrays/mapped_functions.output +++ b/test/triqs/arrays/mapped_functions.output @@ -139,9 +139,9 @@ [[2.71828,7.38906] [20.0855,54.5982]] (make_matrix(conj(A))) ---> -[[(1,0),(3,0),(5,0)] - [(2,0),(4,0),(6,0)] - [(3,0),(5,0),(7,0)]] +[[1,3,5] + [2,4,6] + [3,5,7]] (A) ---> [[(1,0),(3,0),(5,0)] [(2,0),(4,0),(6,0)] diff --git a/triqs/arrays/mapped_functions.hpp b/triqs/arrays/mapped_functions.hpp index 32b3af73..51206c17 100644 --- a/triqs/arrays/mapped_functions.hpp +++ b/triqs/arrays/mapped_functions.hpp @@ -25,9 +25,10 @@ namespace triqs { namespace arrays { // complex conjugation for integers - inline int conj(int const& x) { return x; } - inline long conj(long const& x) { return x; } - inline long long conj(long long const& x) { return x; } + inline int conj(int x) { return x; } + inline long conj(long x) { return x; } + inline long long conj(long long x) { return x; } + inline double conj(double x) { return x; } //C++14 will simply be ... //template decltype(auto) abs(A && a) { return map( [](auto const &x) { using std::abs; return abs(a);}, std::forward(a));} diff --git a/triqs/gfs/domains/matsubara.hpp b/triqs/gfs/domains/matsubara.hpp index 4d6b2bf8..9157f78e 100644 --- a/triqs/gfs/domains/matsubara.hpp +++ b/triqs/gfs/domains/matsubara.hpp @@ -21,7 +21,7 @@ #ifndef TRIQS_GF_DOMAIN_MATSUBARA_H #define TRIQS_GF_DOMAIN_MATSUBARA_H #include "../tools.hpp" -#include "../meshes/mesh_tools.hpp" +#include namespace triqs { namespace gfs { @@ -29,10 +29,11 @@ namespace gfs { // --------------- One matsubara frequency, with its arithmetics ------------------------- // all operations are done by casting to complex, except addition and substraction of matsubara_freq - struct matsubara_freq : public arith_ops_by_cast> { + struct matsubara_freq : public utility::arithmetic_ops_by_cast_disable_same_type> { int n; double beta; statistic_enum statistic; + matsubara_freq() : n(0), beta(1), statistic(Fermion) {} matsubara_freq(int const &n_, double beta_, statistic_enum stat_) : n(n_), beta(beta_), statistic(stat_) {} using cast_t = std::complex; operator cast_t() const { @@ -40,42 +41,17 @@ namespace gfs { } }; - inline matsubara_freq operator+(matsubara_freq const &x, matsubara_freq const &y) { + inline matsubara_freq operator+(matsubara_freq const &x, matsubara_freq const &y) { return {x.n + y.n + (x.statistic & y.statistic), x.beta, ((x.statistic ^ y.statistic) == 1 ? Fermion : Boson)}; } inline matsubara_freq operator-(matsubara_freq const &x, matsubara_freq const &y) { - // std::cout << x.n - y.n - (~x.statistic & y.statistic)<< std::endl; return {x.n - y.n - (~x.statistic & y.statistic), x.beta, ((x.statistic ^ y.statistic) == 1 ? Fermion : Boson)}; } - /* -#define IMPL_OP(OP) \ - template auto operator OP(matsubara_freq const &x, Y &&y)->decltype(std::complex{}) { \ - return std::complex(x) OP std::forward(y); \ - } \ - template \ - auto operator OP(Y &&y, matsubara_freq const &x)->decltype(std::forward(y) OP std::complex{}) { \ - return std::forward(y) OP std::complex(x); \ + inline matsubara_freq operator-(matsubara_freq const &mp) { + return {-(mp.n + mp.statistic==Fermion ? 1: 0), mp.beta, mp.statistic}; } - IMPL_OP(+); - IMPL_OP(-); - IMPL_OP(*); - IMPL_OP(/ ); -#undef IMPL_OP -*/ - -/* - * template auto operator OP(matsubara_freq const &x, Y &&y)->decltype(std::complex{}) { \ - return std::complex(x) OP std::forward(y); \ - } \ - template \ - auto operator OP(Y &&y, matsubara_freq const &x) \ - ->TYPE_DISABLE_IF(decltype(std::forward(y) OP std::complex{}), \ - std::is_same::type>::type, matsubara_freq>) { \ - return std::forward(y) OP std::complex(x); \ - } -*/ //--------------------------------------------------------------------------------------------------------- /// The domain diff --git a/triqs/gfs/imfreq.hpp b/triqs/gfs/imfreq.hpp index e12f6bdd..eb755e8c 100644 --- a/triqs/gfs/imfreq.hpp +++ b/triqs/gfs/imfreq.hpp @@ -109,9 +109,26 @@ namespace gfs { return r; } + // does not work on gcc 4.8.1 ??? + /* template + auto operator()(G const *g, matsubara_freq const &f) const + DECL_AND_RETURN(_call_impl(g, f, Target{}, std::integral_constant::value>{})); + */ + template - auto operator()(G const *g, matsubara_freq const &f) const - DECL_AND_RETURN(_call_impl(g, f, Target{}, std::integral_constant::value>{})); + typename std::conditional::value, arrays::matrix_const_view>, + std::complex>::type + operator()(G const *g, matsubara_freq const &f) const { + return _call_impl(g, f, Target{}, std::integral_constant::value>{}); + } + +#ifdef __clang__ + // to generate a clearer error message ? . Only ok on clang ? + template struct error { static_assert(n>0, "Green function can not be evaluated on a complex number !");}; + + template + error<0> operator()(G const *g, std::complex) const { return {};} +#endif template typename G::singularity_t const &operator()(G const *g, freq_infty const &) const { return g->singularity(); diff --git a/triqs/gfs/meshes/discrete.hpp b/triqs/gfs/meshes/discrete.hpp index ef00f653..e71e31c3 100644 --- a/triqs/gfs/meshes/discrete.hpp +++ b/triqs/gfs/meshes/discrete.hpp @@ -42,7 +42,7 @@ namespace triqs { namespace gfs { size_t index_to_linear(index_t ind) const {return ind;} /// The wrapper for the mesh point - class mesh_point_t : tag::mesh_point,public arith_ops_by_cast { + class mesh_point_t : tag::mesh_point,public utility::arithmetic_ops_by_cast { discrete_mesh const * m; index_t _index; public: diff --git a/triqs/gfs/meshes/linear.hpp b/triqs/gfs/meshes/linear.hpp index 2128ede0..8fa1919f 100644 --- a/triqs/gfs/meshes/linear.hpp +++ b/triqs/gfs/meshes/linear.hpp @@ -82,7 +82,7 @@ namespace gfs { size_t index_to_linear(index_t ind) const { return ind; } /// The wrapper for the mesh point - class mesh_point_t : tag::mesh_point, public arith_ops_by_cast { + class mesh_point_t : tag::mesh_point, public utility::arithmetic_ops_by_cast { linear_mesh const *m; index_t _index; diff --git a/triqs/gfs/meshes/matsubara_freq.hpp b/triqs/gfs/meshes/matsubara_freq.hpp index 2778bc8a..b8803b3d 100644 --- a/triqs/gfs/meshes/matsubara_freq.hpp +++ b/triqs/gfs/meshes/matsubara_freq.hpp @@ -49,22 +49,14 @@ namespace gfs { size_t index_to_linear(index_t ind) const { return ind; } /// The wrapper for the mesh point - struct mesh_point_t : tag::mesh_point, public arith_ops_by_cast { - index_t n; - double beta; - statistic_enum statistic; + struct mesh_point_t : tag::mesh_point, matsubara_freq { index_t n_max; - - mesh_point_t()= default; + mesh_point_t() = default; mesh_point_t(matsubara_freq_mesh const &mesh, index_t const &index_) - : n(index_), beta(mesh.domain().beta), statistic(mesh.domain().statistic), n_max(mesh.size()) {} + : matsubara_freq(index_, mesh.domain().beta, mesh.domain().statistic), n_max(mesh.size()) {} mesh_point_t(index_t const &index_, double beta_, statistic_enum stat_, index_t n_max_) - : n(index_), beta(beta_), statistic(stat_), n_max(n_max_) {} + : matsubara_freq(index_, beta_, stat_), n_max(n_max_) {} void advance() { ++n; } - using cast_t = domain_pt_t; - operator cast_t() const { return 1_j * M_PI * (2 * n + (statistic == Fermion ? 1 : 0)) / beta; } - using cast1_t = matsubara_freq; - operator cast1_t() const { return {n,beta,statistic};} size_t linear_index() const { return n; } size_t index() const { return n; } bool at_end() const { return (n == n_max); } @@ -92,7 +84,7 @@ namespace gfs { h5::group gr = fg.create_group(subgroup_name); h5_write(gr, "domain", m.domain()); h5_write(gr, "size", m.size()); - //h5_write(gr, "start_at_0", m._positive_only); + // h5_write(gr, "start_at_0", m._positive_only); // kept for backward compatibility auto beta = m.domain().beta; h5_write(gr, "min", Fermion ? M_PI / beta : 0); @@ -108,7 +100,7 @@ namespace gfs { bool s = true; h5_read(gr, "domain", dom); h5_read(gr, "size", L); - //h5_read(gr, "start_at_0", s); + // h5_read(gr, "start_at_0", s); m = matsubara_freq_mesh{std::move(dom), L, s}; } @@ -129,7 +121,6 @@ namespace gfs { int n_max; bool _positive_only; }; - - } +} } diff --git a/triqs/gfs/meshes/mesh_tools.hpp b/triqs/gfs/meshes/mesh_tools.hpp index 7f0fe6f7..d7b1c75c 100644 --- a/triqs/gfs/meshes/mesh_tools.hpp +++ b/triqs/gfs/meshes/mesh_tools.hpp @@ -21,35 +21,11 @@ #ifndef TRIQS_GF_MESHTOOLS_H #define TRIQS_GF_MESHTOOLS_H #include "../tools.hpp" +#include namespace triqs { namespace gfs { - // Derive from this object using CRTP to provide arithmetic operation by casting the final object to C - // by not forwarding x, I assume the cast is a simple value, not a matrix, but this is ok - template struct arith_ops_by_cast {}; - -#define IMPL_OP(OP) \ - template \ - auto operator OP(arith_ops_by_cast const& x, Y&& y)->decltype(std::declval() OP std::forward(y)) { \ - return C(static_cast(x)) OP std::forward(y); \ - } \ - template \ - auto operator OP(Y&& y, arith_ops_by_cast const& x) \ - ->TYPE_DISABLE_IF(decltype(std::forward(y) OP std::declval()), \ - std::is_same::type>::type, D>) { \ - return std::forward(y) OP C(static_cast(x)); \ - } - - IMPL_OP(+); - IMPL_OP(-); - IMPL_OP(*); - IMPL_OP(/ ); -#undef IMPL_OP - - - //------------------------------------------------------ - template class mesh_pt_generator : public boost::iterator_facade< mesh_pt_generator, typename MeshType::mesh_point_t , boost::forward_traversal_tag, @@ -70,6 +46,5 @@ namespace gfs { typename MeshType::domain_t::point_t to_point() const { return pt;} }; - }} #endif diff --git a/triqs/utility/arithmetic_ops_by_cast.hpp b/triqs/utility/arithmetic_ops_by_cast.hpp new file mode 100644 index 00000000..0f18c18d --- /dev/null +++ b/triqs/utility/arithmetic_ops_by_cast.hpp @@ -0,0 +1,69 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 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 +namespace triqs { +namespace utility { + + // Derive from this object using CRTP to provide arithmetic operation by casting the final object to C + // by not forwarding x, I assume the cast is a simple value, not a matrix, but this is ok + template struct arithmetic_ops_by_cast {}; + +#define IMPL_OP(OP) \ + template \ + auto operator OP(arithmetic_ops_by_cast const &x, Y &&y)->decltype(std::declval() OP std::forward(y)) { \ + return C(static_cast(x)) OP std::forward(y); \ + } \ + template \ + auto operator OP(Y &&y, arithmetic_ops_by_cast const &x) \ + ->TYPE_DISABLE_IF(decltype(std::forward(y) OP std::declval()), \ + std::is_base_of::type>::type>) { \ + return std::forward(y) OP C(static_cast(x)); \ + } + + IMPL_OP(+); + IMPL_OP(-); + IMPL_OP(*); + IMPL_OP(/ ); +#undef IMPL_OP + + // Same thing, but disallow the operations with the type itself + template struct arithmetic_ops_by_cast_disable_same_type {}; +#define IMPL_OP(OP) \ + template \ + auto operator OP(arithmetic_ops_by_cast_disable_same_type const &x, Y &&y) \ + ->TYPE_DISABLE_IF(decltype(std::declval() OP std::forward(y)), \ + std::is_base_of::type>::type>) { \ + return C(static_cast(x)) OP std::forward(y); \ + } \ + template \ + auto operator OP(Y &&y, arithmetic_ops_by_cast_disable_same_type const &x) \ + ->TYPE_DISABLE_IF(decltype(std::forward(y) OP std::declval()), \ + std::is_base_of::type>::type>) { \ + return std::forward(y) OP C(static_cast(x)); \ + } + IMPL_OP(+); + IMPL_OP(-); + IMPL_OP(*); + IMPL_OP(/ ); +#undef IMPL_OP +} +} +