From a69e2f52aa0cad57b3de1087c91fee9a4cb8d61c Mon Sep 17 00:00:00 2001 From: Laura Messio Date: Fri, 19 Jul 2013 10:55:37 +0200 Subject: [PATCH] fourier_real: scalar_valued implementation --- triqs/gf/local/fourier_real.cpp | 224 ++++++++++++++++++-------------- triqs/gf/local/fourier_real.hpp | 74 +++++++---- 2 files changed, 172 insertions(+), 126 deletions(-) diff --git a/triqs/gf/local/fourier_real.cpp b/triqs/gf/local/fourier_real.cpp index 664f0202..6c1ddc09 100644 --- a/triqs/gf/local/fourier_real.cpp +++ b/triqs/gf/local/fourier_real.cpp @@ -1,5 +1,5 @@ /******************************************************************************* - * + * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * * Copyright (C) 2011 by M. Ferrero, O. Parcollet @@ -23,106 +23,128 @@ #include namespace triqs { namespace gf { - - namespace { - double pi = std::acos(-1); - dcomplex I(0,1); - inline dcomplex th_expo(double t, double a ) { return (t > 0 ? -I * exp(-a*t) : ( t < 0 ? 0 : -0.5 * I * exp(-a*t) ) ) ; } - inline dcomplex th_expo_neg(double t, double a ) { return (t < 0 ? I * exp( a*t) : ( t > 0 ? 0 : 0.5 * I * exp( a*t) ) ) ; } - inline dcomplex th_expo_inv(double w, double a ) { return 1./(w+I*a) ; } - inline dcomplex th_expo_neg_inv(double w, double a ) { return 1./(w-I*a) ; } + + namespace { + double pi = std::acos(-1); + dcomplex I(0,1); + inline dcomplex th_expo(double t, double a ) { return (t > 0 ? -I * exp(-a*t) : ( t < 0 ? 0 : -0.5 * I * exp(-a*t) ) ) ; } + inline dcomplex th_expo_neg(double t, double a ) { return (t < 0 ? I * exp( a*t) : ( t > 0 ? 0 : 0.5 * I * exp( a*t) ) ) ; } + inline dcomplex th_expo_inv(double w, double a ) { return 1./(w+I*a) ; } + inline dcomplex th_expo_neg_inv(double w, double a ) { return 1./(w-I*a) ; } + } + + struct impl_worker { + + tqa::vector g_in, g_out; + + void direct (gf_view gw, gf_view const gt){ + + size_t L = gt.mesh().size(); + if (gw.mesh().size() != L) TRIQS_RUNTIME_ERROR << "Meshes are different"; + double test = std::abs(gt.mesh().delta() * gw.mesh().delta() * L / (2*pi) -1); + if (test > 1.e-10) TRIQS_RUNTIME_ERROR << "Meshes are not compatible"; + + const double tmin = gt.mesh().x_min(); + const double wmin = gw.mesh().x_min(); + //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()); + g_in.resize(L); + g_in() = 0; + g_out.resize(L); + + dcomplex t1 = ta(1)(0,0), t2= ta.get_or_zero(2)(0,0); + dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.; + + for (auto const & t : gt.mesh()) + g_in(t.index()) = (gt(t) - (a1*th_expo(t,a) + a2*th_expo_neg(t,a))) * std::exp(I*t*wmin); + + details::fourier_base(g_in, g_out, L, true); + + for (auto const & w : gw.mesh()) + gw(w) = gt.mesh().delta() * std::exp(I*(w-wmin)*tmin) * g_out(w.index()) + + a1*th_expo_inv(w,a) + a2*th_expo_neg_inv(w,a); + + gw.singularity() = gt.singularity();// set tail + } - - //-------------------------------------------------------------------------------------- - - void fourier_impl(gf_view gw, gf_view const gt) { - - size_t L = gt.mesh().size(); - if (gw.mesh().size() != L) TRIQS_RUNTIME_ERROR << "Meshes are different"; - double test = std::abs(gt.mesh().delta() * gw.mesh().delta() * L / (2*pi) -1); - if (test > 1.e-10) TRIQS_RUNTIME_ERROR << "Meshes are not compatible"; - - const double tmin = gt.mesh().x_min(); - const double wmin = gw.mesh().x_min(); - //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()); - tqa::vector g_in(L), g_out(L); - - for (size_t n1=0; n1 gt, gf_view const gw){ + + size_t L = gw.mesh().size(); + if ( L != gt.mesh().size()) TRIQS_RUNTIME_ERROR << "Meshes are different"; + double test = std::abs(gt.mesh().delta() * gw.mesh().delta() * L / (2*pi) -1); + if (test > 1.e-10) TRIQS_RUNTIME_ERROR << "Meshes are not compatible"; + + const double tmin = gt.mesh().x_min(); + const double wmin = gw.mesh().x_min(); + //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()); + tqa::vector g_in(L), g_out(L); + + dcomplex t1 = ta(1)(0,0), t2 = ta.get_or_zero(2)(0,0); + dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.; + g_in() = 0; + + for (auto const & w: gw.mesh()) + g_in(w.index()) = (gw(w) - a1*th_expo_inv(w,a) - a2*th_expo_neg_inv(w,a) ) * std::exp(-I*w*tmin); + + details::fourier_base(g_in, g_out, L, false); + + const double corr = 1.0/(gt.mesh().delta()*L); + for (auto const & t : gt.mesh()) + gt(t) = corr * std::exp(I*wmin*(tmin-t)) * + g_out( t.index() ) + a1 * th_expo(t,a) + a2 * th_expo_neg(t,a) ; + + // set tail + gt.singularity() = gw.singularity(); } - - //--------------------------------------------------------------------------- - - void inverse_fourier_impl (gf_view gt, gf_view const gw) { - - size_t L = gw.mesh().size(); - if ( L != gt.mesh().size()) TRIQS_RUNTIME_ERROR << "Meshes are different"; - double test = std::abs(gt.mesh().delta() * gw.mesh().delta() * L / (2*pi) -1); - if (test > 1.e-10) TRIQS_RUNTIME_ERROR << "Meshes are not compatible"; - - const double tmin = gt.mesh().x_min(); - const double wmin = gw.mesh().x_min(); - //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()); - tqa::vector g_in(L), g_out(L); - - for (size_t n1=0; n1 lazy_fourier (gf_view const & g) { return g;} - gf_keeper lazy_inverse_fourier (gf_view const & g) { return 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 fourier_impl(gf_view gw, gf_view const gt, scalar_valued){ + impl_worker w; + w.direct(gw, gt); + } + + void fourier_impl(gf_view gw, gf_view const gt, matrix_valued){ + impl_worker w; + for (size_t n1=0; n1 gt, gf_view const gw, scalar_valued){ + impl_worker w; + w.inverse(gt,gw); + } + + void inverse_fourier_impl (gf_view gt, gf_view const 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/gf/local/fourier_real.hpp b/triqs/gf/local/fourier_real.hpp index 81909ebc..f3bb4d47 100644 --- a/triqs/gf/local/fourier_real.hpp +++ b/triqs/gf/local/fourier_real.hpp @@ -28,33 +28,57 @@ namespace triqs { namespace gf { // First the implementation of the fourier transform - void fourier_impl (gf_view gw , gf_view const gt); - void inverse_fourier_impl (gf_view gt, gf_view const gw); - - inline gf_view fourier (gf_view const & gt) { - double pi = std::acos(-1); - size_t L = gt.mesh().size(); - double wmin = -pi * (L-1) / (L*gt.mesh().delta()); - double wmax = pi * (L-1) / (L*gt.mesh().delta()); - auto gw = make_gf(wmin, wmax, L, gt.data().shape().front_pop()); - auto V = gw(); - fourier_impl(V, gt); - return gw; + void fourier_impl (gf_view gw , gf_view const gt, scalar_valued); + void fourier_impl (gf_view gw , gf_view const gt, matrix_valued); + void inverse_fourier_impl (gf_view gt, gf_view const gw, scalar_valued); + void inverse_fourier_impl (gf_view gt, gf_view const gw, matrix_valued); + + inline gf_view fourier (gf_view const gt) { + double pi = std::acos(-1); + size_t L = gt.mesh().size(); + double wmin = -pi * (L-1) / (L*gt.mesh().delta()); + double wmax = pi * (L-1) / (L*gt.mesh().delta()); + auto gw = make_gf(wmin, wmax, L, gt.data().shape().front_pop()); + auto V = gw(); + fourier_impl(V, gt, matrix_valued()); + return gw; } - - inline gf_view inverse_fourier (gf_view const & gw) { - double pi = std::acos(-1); - size_t L = gw.mesh().size(); - double tmin = -pi * (L-1) / (L*gw.mesh().delta()); - double tmax = pi * (L-1) / (L*gw.mesh().delta()); - auto gt = make_gf(tmin, tmax, L, gw.data().shape().front_pop()); - auto V = gt(); - inverse_fourier_impl(V, gw); - return gt; + inline gf_view fourier (gf_view const gt) { + double pi = std::acos(-1); + size_t L = gt.mesh().size(); + double wmin = -pi * (L-1) / (L*gt.mesh().delta()); + double wmax = pi * (L-1) / (L*gt.mesh().delta()); + auto gw = make_gf(wmin, wmax, L); + auto V = gw(); + fourier_impl(V, gt, scalar_valued()); + return gw; } - - gf_keeper lazy_fourier (gf_view const & g); - gf_keeper lazy_inverse_fourier (gf_view const & g); + + inline gf_view inverse_fourier (gf_view const gw) { + double pi = std::acos(-1); + size_t L = gw.mesh().size(); + double tmin = -pi * (L-1) / (L*gw.mesh().delta()); + double tmax = pi * (L-1) / (L*gw.mesh().delta()); + auto gt = make_gf(tmin, tmax, L, gw.data().shape().front_pop()); + auto V = gt(); + inverse_fourier_impl(V, gw, matrix_valued()); + return gt; + } + inline gf_view inverse_fourier (gf_view const gw) { + double pi = std::acos(-1); + size_t L = gw.mesh().size(); + double tmin = -pi * (L-1) / (L*gw.mesh().delta()); + double tmax = pi * (L-1) / (L*gw.mesh().delta()); + auto gt = make_gf(tmin, tmax, L); + auto V = gt(); + inverse_fourier_impl(V, gw, scalar_valued()); + return gt; + } + + inline gf_keeper lazy_fourier (gf_view const & g) { return g;} + inline gf_keeper lazy_inverse_fourier (gf_view const & g) { return g;} + inline gf_keeper lazy_fourier (gf_view const & g) { return g;} + inline gf_keeper lazy_inverse_fourier (gf_view const & g) { return g;} 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);