mirror of
https://github.com/triqs/dft_tools
synced 2024-12-25 13:53:40 +01:00
fourier_real: scalar_valued implementation
This commit is contained in:
parent
f4d42a4ec8
commit
a69e2f52aa
@ -1,5 +1,5 @@
|
|||||||
/*******************************************************************************
|
/*******************************************************************************
|
||||||
*
|
*
|
||||||
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
* TRIQS: a Toolbox for Research in Interacting Quantum Systems
|
||||||
*
|
*
|
||||||
* Copyright (C) 2011 by M. Ferrero, O. Parcollet
|
* Copyright (C) 2011 by M. Ferrero, O. Parcollet
|
||||||
@ -23,106 +23,128 @@
|
|||||||
#include <fftw3.h>
|
#include <fftw3.h>
|
||||||
|
|
||||||
namespace triqs { namespace gf {
|
namespace triqs { namespace gf {
|
||||||
|
|
||||||
namespace {
|
namespace {
|
||||||
double pi = std::acos(-1);
|
double pi = std::acos(-1);
|
||||||
dcomplex I(0,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(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_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_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) ; }
|
inline dcomplex th_expo_neg_inv(double w, double a ) { return 1./(w-I*a) ; }
|
||||||
|
}
|
||||||
|
|
||||||
|
struct impl_worker {
|
||||||
|
|
||||||
|
tqa::vector<dcomplex> g_in, g_out;
|
||||||
|
|
||||||
|
void direct (gf_view<refreq,scalar_valued> gw, gf_view<retime,scalar_valued> 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 inverse(gf_view<retime,scalar_valued> gt, gf_view<refreq,scalar_valued> const gw){
|
||||||
|
|
||||||
void fourier_impl(gf_view<refreq> gw, gf_view<retime> const gt) {
|
size_t L = gw.mesh().size();
|
||||||
|
if ( L != gt.mesh().size()) TRIQS_RUNTIME_ERROR << "Meshes are different";
|
||||||
size_t L = gt.mesh().size();
|
double test = std::abs(gt.mesh().delta() * gw.mesh().delta() * L / (2*pi) -1);
|
||||||
if (gw.mesh().size() != L) TRIQS_RUNTIME_ERROR << "Meshes are different";
|
if (test > 1.e-10) TRIQS_RUNTIME_ERROR << "Meshes are not compatible";
|
||||||
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();
|
||||||
const double tmin = gt.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 wmin = gw.mesh().x_min();
|
const double a = gw.mesh().delta() * sqrt( double(L) );
|
||||||
//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<dcomplex> g_in(L), g_out(L);
|
||||||
auto ta = gt(freq_infty());
|
|
||||||
tqa::vector<dcomplex> 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.;
|
||||||
for (size_t n1=0; n1<gw.data().shape()[1];n1++) {
|
g_in() = 0;
|
||||||
for (size_t n2=0; n2<gw.data().shape()[2];n2++) {
|
|
||||||
|
for (auto const & w: gw.mesh())
|
||||||
dcomplex t1 = ta(1)(n1,n2), t2= ta.get_or_zero(2)(n1,n2);
|
g_in(w.index()) = (gw(w) - a1*th_expo_inv(w,a) - a2*th_expo_neg_inv(w,a) ) * std::exp(-I*w*tmin);
|
||||||
dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.;
|
|
||||||
|
details::fourier_base(g_in, g_out, L, false);
|
||||||
g_in() = 0;
|
|
||||||
|
const double corr = 1.0/(gt.mesh().delta()*L);
|
||||||
for (auto const & t : gt.mesh()) {
|
for (auto const & t : gt.mesh())
|
||||||
g_in(t.index()) = (gt(t)(n1,n2) - (a1*th_expo(t,a) + a2*th_expo_neg(t,a))) * std::exp(I*t*wmin);
|
gt(t) = corr * std::exp(I*wmin*(tmin-t)) *
|
||||||
}
|
g_out( t.index() ) + a1 * th_expo(t,a) + a2 * th_expo_neg(t,a) ;
|
||||||
|
|
||||||
details::fourier_base(g_in, g_out, L, true);
|
// set tail
|
||||||
|
gt.singularity() = gw.singularity();
|
||||||
for (auto const & w : gw.mesh()) {
|
|
||||||
gw(w)(n1,n2) = 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);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
// set tail
|
|
||||||
gw.singularity() = gt.singularity();
|
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
//---------------------------------------------------------------------------
|
};
|
||||||
|
|
||||||
void inverse_fourier_impl (gf_view<retime> gt, gf_view<refreq> const gw) {
|
|
||||||
|
//--------------------------------------------------------------------------------------
|
||||||
size_t L = gw.mesh().size();
|
|
||||||
if ( L != gt.mesh().size()) TRIQS_RUNTIME_ERROR << "Meshes are different";
|
void fourier_impl(gf_view<refreq,scalar_valued> gw, gf_view<retime,scalar_valued> const gt, scalar_valued){
|
||||||
double test = std::abs(gt.mesh().delta() * gw.mesh().delta() * L / (2*pi) -1);
|
impl_worker w;
|
||||||
if (test > 1.e-10) TRIQS_RUNTIME_ERROR << "Meshes are not compatible";
|
w.direct(gw, gt);
|
||||||
|
}
|
||||||
const double tmin = gt.mesh().x_min();
|
|
||||||
const double wmin = gw.mesh().x_min();
|
void fourier_impl(gf_view<refreq,matrix_valued> gw, gf_view<retime,matrix_valued> const gt, matrix_valued){
|
||||||
//a is a number very larger than delta_w and very smaller than wmax-wmin, used in the tail computation
|
impl_worker w;
|
||||||
const double a = gw.mesh().delta() * sqrt( double(L) );
|
for (size_t n1=0; n1<gw.data().shape()[1];n1++)
|
||||||
|
for (size_t n2=0; n2<gw.data().shape()[2];n2++) {
|
||||||
auto ta = gw(freq_infty());
|
auto gw_sl=slice_target_to_scalar(gw, n1, n2);
|
||||||
tqa::vector<dcomplex> g_in(L), g_out(L);
|
auto gt_sl=slice_target_to_scalar(gt, n1, n2);
|
||||||
|
w.direct(gw_sl, gt_sl);
|
||||||
for (size_t n1=0; n1<gt.data().shape()[1];n1++) {
|
}
|
||||||
for (size_t n2=0; n2<gt.data().shape()[2];n2++) {
|
}
|
||||||
|
|
||||||
dcomplex t1 = ta(1)(n1,n2), t2 = ta.get_or_zero(2)(n1,n2);
|
//---------------------------------------------------------------------------
|
||||||
dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.;
|
|
||||||
g_in() = 0;
|
void inverse_fourier_impl (gf_view<retime,scalar_valued> gt, gf_view<refreq,scalar_valued> const gw, scalar_valued){
|
||||||
|
impl_worker w;
|
||||||
for (auto const & w: gw.mesh())
|
w.inverse(gt,gw);
|
||||||
g_in(w.index()) = (gw(w)(n1,n2) - 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);
|
void inverse_fourier_impl (gf_view<retime,matrix_valued> gt, gf_view<refreq,matrix_valued> const gw, matrix_valued){
|
||||||
|
impl_worker w;
|
||||||
const double corr = 1.0/(gt.mesh().delta()*L);
|
for (size_t n1=0; n1<gt.data().shape()[1];n1++)
|
||||||
for (auto const & t : gt.mesh())
|
for (size_t n2=0; n2<gt.data().shape()[2];n2++) {
|
||||||
gt(t)(n1,n2) = corr * std::exp(I*wmin*(tmin-t)) *
|
auto gt_sl=slice_target_to_scalar(gt, n1, n2);
|
||||||
g_out( t.index() ) + a1 * th_expo(t,a) + a2 * th_expo_neg(t,a) ;
|
auto gw_sl=slice_target_to_scalar(gw, n1, n2);
|
||||||
}
|
w.inverse(gt_sl, gw_sl);
|
||||||
}
|
}
|
||||||
|
}
|
||||||
// set tail
|
|
||||||
gt.singularity() = gw.singularity();
|
//---------------------------------------------------------------------------
|
||||||
|
|
||||||
}
|
void triqs_gf_view_assign_delegation( gf_view<refreq,matrix_valued> g, gf_keeper<tags::fourier,retime,matrix_valued> const & L) { fourier_impl (g,L.g, matrix_valued());}
|
||||||
|
void triqs_gf_view_assign_delegation( gf_view<refreq,scalar_valued> g, gf_keeper<tags::fourier,retime,scalar_valued> const & L) { fourier_impl (g,L.g, scalar_valued());}
|
||||||
//---------------------------------------------------------------------------
|
void triqs_gf_view_assign_delegation( gf_view<retime,matrix_valued> g, gf_keeper<tags::fourier,refreq,matrix_valued> const & L) { inverse_fourier_impl(g,L.g, matrix_valued());}
|
||||||
|
void triqs_gf_view_assign_delegation( gf_view<retime,scalar_valued> g, gf_keeper<tags::fourier,refreq,scalar_valued> const & L) { inverse_fourier_impl(g,L.g, scalar_valued());}
|
||||||
gf_keeper<tags::fourier,retime> lazy_fourier (gf_view<retime> const & g) { return g;}
|
|
||||||
gf_keeper<tags::fourier,refreq> lazy_inverse_fourier (gf_view<refreq> const & g) { return g;}
|
|
||||||
|
|
||||||
void triqs_gf_view_assign_delegation( gf_view<refreq> g, gf_keeper<tags::fourier,retime> const & L) { fourier_impl (g,L.g);}
|
|
||||||
void triqs_gf_view_assign_delegation( gf_view<retime> g, gf_keeper<tags::fourier,refreq> const & L) { inverse_fourier_impl(g,L.g);}
|
|
||||||
|
|
||||||
}}
|
}}
|
||||||
|
@ -28,33 +28,57 @@
|
|||||||
namespace triqs { namespace gf {
|
namespace triqs { namespace gf {
|
||||||
|
|
||||||
// First the implementation of the fourier transform
|
// First the implementation of the fourier transform
|
||||||
void fourier_impl (gf_view<refreq> gw , gf_view<retime> const gt);
|
void fourier_impl (gf_view<refreq,scalar_valued> gw , gf_view<retime,scalar_valued> const gt, scalar_valued);
|
||||||
void inverse_fourier_impl (gf_view<retime> gt, gf_view<refreq> const gw);
|
void fourier_impl (gf_view<refreq,matrix_valued> gw , gf_view<retime,matrix_valued> const gt, matrix_valued);
|
||||||
|
void inverse_fourier_impl (gf_view<retime,scalar_valued> gt, gf_view<refreq,scalar_valued> const gw, scalar_valued);
|
||||||
inline gf_view<refreq> fourier (gf_view<retime> const & gt) {
|
void inverse_fourier_impl (gf_view<retime,matrix_valued> gt, gf_view<refreq,matrix_valued> const gw, matrix_valued);
|
||||||
double pi = std::acos(-1);
|
|
||||||
size_t L = gt.mesh().size();
|
inline gf_view<refreq,matrix_valued> fourier (gf_view<retime, matrix_valued> const gt) {
|
||||||
double wmin = -pi * (L-1) / (L*gt.mesh().delta());
|
double pi = std::acos(-1);
|
||||||
double wmax = pi * (L-1) / (L*gt.mesh().delta());
|
size_t L = gt.mesh().size();
|
||||||
auto gw = make_gf<refreq>(wmin, wmax, L, gt.data().shape().front_pop());
|
double wmin = -pi * (L-1) / (L*gt.mesh().delta());
|
||||||
auto V = gw();
|
double wmax = pi * (L-1) / (L*gt.mesh().delta());
|
||||||
fourier_impl(V, gt);
|
auto gw = make_gf<refreq,matrix_valued>(wmin, wmax, L, gt.data().shape().front_pop());
|
||||||
return gw;
|
auto V = gw();
|
||||||
|
fourier_impl(V, gt, matrix_valued());
|
||||||
|
return gw;
|
||||||
}
|
}
|
||||||
|
inline gf_view<refreq,scalar_valued> fourier (gf_view<retime, scalar_valued> const gt) {
|
||||||
inline gf_view<retime> inverse_fourier (gf_view<refreq> const & gw) {
|
double pi = std::acos(-1);
|
||||||
double pi = std::acos(-1);
|
size_t L = gt.mesh().size();
|
||||||
size_t L = gw.mesh().size();
|
double wmin = -pi * (L-1) / (L*gt.mesh().delta());
|
||||||
double tmin = -pi * (L-1) / (L*gw.mesh().delta());
|
double wmax = pi * (L-1) / (L*gt.mesh().delta());
|
||||||
double tmax = pi * (L-1) / (L*gw.mesh().delta());
|
auto gw = make_gf<refreq,scalar_valued>(wmin, wmax, L);
|
||||||
auto gt = make_gf<retime>(tmin, tmax, L, gw.data().shape().front_pop());
|
auto V = gw();
|
||||||
auto V = gt();
|
fourier_impl(V, gt, scalar_valued());
|
||||||
inverse_fourier_impl(V, gw);
|
return gw;
|
||||||
return gt;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
gf_keeper<tags::fourier,retime> lazy_fourier (gf_view<retime> const & g);
|
inline gf_view<retime,matrix_valued> inverse_fourier (gf_view<refreq,matrix_valued> const gw) {
|
||||||
gf_keeper<tags::fourier,refreq> lazy_inverse_fourier (gf_view<refreq> const & g);
|
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<retime,matrix_valued>(tmin, tmax, L, gw.data().shape().front_pop());
|
||||||
|
auto V = gt();
|
||||||
|
inverse_fourier_impl(V, gw, matrix_valued());
|
||||||
|
return gt;
|
||||||
|
}
|
||||||
|
inline gf_view<retime,scalar_valued> inverse_fourier (gf_view<refreq,scalar_valued> 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<retime,scalar_valued>(tmin, tmax, L);
|
||||||
|
auto V = gt();
|
||||||
|
inverse_fourier_impl(V, gw, scalar_valued());
|
||||||
|
return gt;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline gf_keeper<tags::fourier,retime,scalar_valued> lazy_fourier (gf_view<retime,scalar_valued> const & g) { return g;}
|
||||||
|
inline gf_keeper<tags::fourier,refreq,scalar_valued> lazy_inverse_fourier (gf_view<refreq,scalar_valued> const & g) { return g;}
|
||||||
|
inline gf_keeper<tags::fourier,retime,matrix_valued> lazy_fourier (gf_view<retime,matrix_valued> const & g) { return g;}
|
||||||
|
inline gf_keeper<tags::fourier,refreq,matrix_valued> lazy_inverse_fourier (gf_view<refreq,matrix_valued> const & g) { return g;}
|
||||||
|
|
||||||
void triqs_gf_view_assign_delegation( gf_view<refreq> g, gf_keeper<tags::fourier,retime> const & L);
|
void triqs_gf_view_assign_delegation( gf_view<refreq> g, gf_keeper<tags::fourier,retime> const & L);
|
||||||
void triqs_gf_view_assign_delegation( gf_view<retime> g, gf_keeper<tags::fourier,refreq> const & L);
|
void triqs_gf_view_assign_delegation( gf_view<retime> g, gf_keeper<tags::fourier,refreq> const & L);
|
||||||
|
Loading…
Reference in New Issue
Block a user