diff --git a/doc/reference/c++/gf/clef.rst b/doc/reference/c++/gf/clef.rst new file mode 100644 index 00000000..f7ca3faa --- /dev/null +++ b/doc/reference/c++/gf/clef.rst @@ -0,0 +1,48 @@ +.. highlight:: c + +.. _gf_clef: + +Interaction with CLEF expressions +============================================ + +* The gf containers and their view classes can be used with the :doc:`../clef/contents` library : + + * They can be called with CLEF expressions. + * :doc:`Automatic assignment<../clef/assign>` has been set up. + +* Using the CLEF library offers a quick and efficient way to fill an array with multiple advantages : + + * It is simpler and more readeable than a series of for loops. + * It is usually more optimal since the for loops are automatically written in the TraversalOrder of the array. + + +* **Example** : + +.. compileblock:: + + #include + using namespace triqs::gfs; using triqs::clef::placeholder; + int main(){ + + // Cf gf specialisation page for the constructor + double beta=10; int Nfreq =100; + auto g = gf { {beta,Fermion,Nfreq}, {1,1} }; + + // Filling the gf with something... + placeholder<0> wn_; + g(wn_) << 1/ (wn_ + 2); + g(wn_) << 1/ (wn_ + 2 + g(wn_) ); + } + +.. note:: + + The syntax uses a <<, not = since the array is not assigned to an expression + but filled by the evaluation thereof. + + The LHS uses () and not brackets, even though it is on the mesh, because of the strange C++ limitation + that [] can not be overloaded for several variables... + + + + + diff --git a/doc/reference/c++/gf/fourier.rst b/doc/reference/c++/gf/fourier.rst index 0703acd3..aeaff1f4 100644 --- a/doc/reference/c++/gf/fourier.rst +++ b/doc/reference/c++/gf/fourier.rst @@ -3,152 +3,79 @@ Fourier transforms ################### +The Fourier transforms from real and imaginary frequencies to times, and inverse, are currently implemented, +along with the analogous transformation from the Legendre expansion to imaginary time and frequencies. + +Synopsis and example +====================== + +**Synopsis** :: + + lazy_object fourier (gf const &) + lazy_object fourier (gf_const_view const &) + + lazy_object inverse_fourier (gf const &) + lazy_object inverse_fourier (gf_const_view const &) + +The fourier/inverse_fourier functions do **not** perform the Fourier transformation, +but returns a small lazy object (basically saying "Fourier Transform of XXX"), +which is then used in an assignment of a *view* of a gf. + +Example + +.. compileblock:: + + #include + using namespace triqs::gfs; + int main() { + double beta =1, a=1; + int N=10000; + auto gw = gf {{beta, Fermion, N}, {1,1}}; + auto gt = gf {{beta, Fermion, N}, {1,1}}; + + triqs::clef::placeholder<0> om_; + gw (om_) << 1/(om_-a); + + gt() = inverse_fourier(gw); // fills the *View* with the contents of the FFT. + // NB : the mesh must have the same size. + + // make a new fresh gf, with the same size mesh, from the FFT of gt + auto gw2 = make_gf_from_fourier(gt); + } + +Note that : + +* the LHS of the = must be a view, since the RHS can not compute the domain of the function + (how many points to use ?). + +* In the make_gf_from_fourier function, choice is explicitly made to generate a new gf with the same number of points in the mesh. -Convention -========== +Convention +=========== For real time/frequency: - :label: _TF_R .. math:: \tilde G(\omega)=\int_{-\infty}^\infty dt G(t)e^{i\omega t} - :label: _inv_TF_R .. math:: G(t)=\int_{-\infty}^\infty \frac{d\omega}{2\pi} \tilde G(\omega)e^{-i\omega t} For Matsubara (imaginary) time/frequency: - :label: _TF_I .. math:: \tilde G(i\omega_n)=\int_{0}^\beta d\tau G(t)e^{i\omega_n \tau} - :label: _inv_TF_I .. math:: G(\tau)=\sum_{n=-\infty}^\infty \frac{1}{\beta} \tilde G(i\omega_n)e^{-i\omega_n \tau} The :math:`\omega_n`'s are :math:`\frac{(2n+1)\pi}{\beta}` for fermions, :math:`\frac{2n\pi}{\beta}` for bosons (as :math:`G(\tau+\beta)=-G(\tau)` for fermions, :math:`G(\tau)` for bosons). - -The FFTW library -================ - -Documentation on FFTW is on https://www.fftw.org. -FFTW is a C subroutine library for computing the discrete Fourier transform (DFT) in one or more dimensions, of arbitrary input size, and of both real and complex data. -It will be used to calculate the (inverse) Fourier transform, in real/imaginary time/frequency, using the fact that the GF values are stored for a finite amount of regularly spaced values. - -The DFT transforms of a sequence of :math:`N` complex numbers :math:`f_0...,f_{N-1}` into a sequence of :math:`N` complex numbers :math:`\tilde f_0...,\tilde f_{N-1}` according to the formula: - :label: _DFT - .. math:: \tilde f_k = \sum_{n=0}^{N-1} f_n e^{-i 2 \pi k n / N}. -The inverse DFT formula is - :label: _inv_DFT - .. math:: f_n = \frac{1}{N} \sum_{k=0}^{N-1} \tilde f_k e^{i 2 \pi k n / N}. - - - -Implementation in real time/frequency using FFTW -================================================ - -The real time mesh parameters are :math:`t_{min}`, :math:`\delta t` and :math:`N_t`. -For the real frequency mesh, they are :math:`\omega_{min}`, :math:`\delta \omega` and :math:`N_\omega`. -The Fourier transform requires :math:`N_\omega=N_t` and :math:`\delta t \delta \omega= \frac{2\pi}{N_t}`. -The times are :math:`t_k=t_{min}+k\delta t` and the frequencies :math:`\omega_m=\omega_{min}+m\delta \omega`. - -By approximating Eq. :ref:`TF_R` by - .. math:: \tilde G(\omega_m) = \delta t \sum_{k=0}^{N_t} G(t_k) e^{i\omega_m t_k}, -we recognize an inverse DFT (Eq. :ref:`inv_DFT`). To calculate it using FFTW, we first need to prepare the input :math:`\tilde f_k`, then to do the DFT and finally to modify the output to obtain :math:`\tilde G(\omega_m)` using the two formulas: - .. math:: \tilde f_k = G(t_k) e^{i \omega_{min}t_k}, - .. math:: \tilde G(\omega_m) = \delta t f_m e^{i t_{min}(\omega_m-\omega_{min})}. - -Similarly, the inverse transformation is obtained by approximating Eq. :ref:`eq_inv_TF_R` by - .. math:: G(t_k)=\frac{\delta\omega}{2\pi}\sum_{m=0}^{N_\omega} \tilde G(\omega_m)e^{-i\omega_m t_k}, -we recognize a DFT (Eq. :ref:`DFT`). To calculate it using FFTW, we first need to prepare the input :math:`f_m`, then to do the inverse DFT and finally to modify the output to obtain :math:`G(t_k)`: - .. math:: f_m = \tilde G(\omega_m) e^{-i t_{min}\omega_m}, - .. math:: G(t_k) = \frac{1}{N_t \delta t}\tilde f_k e^{-i \omega_{min}(t_k-t_{min})}. - - - -Implementation in imaginary time/frequency using FFTW -===================================================== - - -The imaginary time mesh parameters are :math:`\beta` and :math:`N_\tau`, plus a tag ``half_bins``, ``full_bins`` or ``without_last``. -In the ``full_bins`` case, one point of the time GF has to be removed for the fourier transform. -From these parameters, we deduce :math:`\delta\tau=\beta/N_\tau`. - -For the imaginary frequency mesh, the mesh parameters are :math:`\beta`, :math:`n_{min}` and :math:`N_{\omega_n}`. -From them, we deduce :math:`\delta\omega=\frac{2\pi}{\beta}`. - -The Fourier transform requires :math:`N_\omega=N_\tau`. -The times are :math:`\tau_k=\tau_{min}+k\delta\tau` and the frequencies :math:`\omega_n=\omega_{min}+n\delta \omega`. -:math:`\tau_{min}` is either 0 or :math:`\delta\tau/2` depending on the mesh kind. -:math:`\omega_{min}` is either :math:`\frac{2\pi(n_{min}+1)}{\beta}` or :math:`\frac{2\pi n_{min}}{\beta}` depending on the statistics. - -We approximate the TF and its inverse by - .. math:: \tilde G(i\omega_n) = \delta\tau \sum_{k=0}^{N_\tau} G(\tau_k)e^{i\omega_n \tau_k} - .. math:: G(\tau_k) = \sum_{n=n_{min}}^{N_\tau} \frac{1}{\beta} \tilde G(i\omega_n)e^{-i\omega_n \tau_k} - -We use for the TF: - .. math:: \tilde f_k = G(\tau_k) e^{i \omega_{min}\tau_k}, - .. math:: \tilde G(i\omega_n) = \frac{\beta}{N_\tau} f_n e^{i \tau_{min}(\omega_n-\omega_{min})}. -and for the inverse TF: - .. math:: f_m = \frac{1}{\beta}\tilde G(i\omega_n) e^{-i t_{min}\omega_n}, - .. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})}, - - -Special case of real functions in time for fermions ----------------------------------------------------- - -In this case, :math:`G(i\omega_n)=conj(G(i\omega_n))` and we only store the values of :math:`G(i\omega_n)` for :math:`\omega_n > 0`. -The Eq. :ref:`inv_DFT_I` becomes: - - :label: _inv_TF_I_real_fermion - .. math:: G(\tau)=\sum_{n=0}^\infty \frac{2}{\beta} \tilde G(i\omega_n)\cos(\omega_n \tau) - -The inverse TF formulas are in this case - .. math:: f_m = \frac{2}{\beta}\tilde G(i\omega_n) e^{-i t_{min}\omega_n}, - .. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})}, - -Special case of real functions in time for bosons --------------------------------------------------- - -In this case, :math:`G(i\omega_n)=conj(G(i\omega_n))` and we only store the values of :math:`G(i\omega_n)` for :math:`\omega_n \ge 0`. -The Eq. :ref:`inv_DFT_I` becomes: - - :label: _inv_TF_I_real_bosons - .. math:: G(\tau)=\frac{1}{\beta} \tilde G(0)+\sum_{n=1}^\infty \frac{2}{\beta} \tilde G(i\omega_n)\cos(\omega_n \tau) - -The inverse TF formulas are in this case - .. math:: f_0 = \frac{1}{\beta}\tilde G(0), - .. math:: f_m = \frac{2}{\beta}\tilde G(i\omega_n) \cos(t_{min}\omega_n), - .. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})}, - - -Effect of a TF on the tail -=========================== - -The tail is unchanged during a TF, but its value is used to limit the errors. - -The components 1 and 2 of the tail (:math:`t_1` and :math:`t_2`) are used to improve the computation of the GF in the following way: -in the large :math:`\omega` limit, - .. math:: G(\omega)\simeq \frac{t_1}{\omega}+\frac{t_2}{\omega^2}\simeq \frac{a_1}{\omega+i}+\frac{a_2}{\omega-i} -with :math:`a_1=\frac{t_1+it_2}{2}` and :math:`a_2=\frac{t_1-it_2}{2}`. - -As these large w terms are badly taken into account if we naively Fourier transform the function described by its values on the mesh in w, we substract them, do the Fourier transform and add their Fourier transform to the result. - -We use the following Fourier tranforms: - - .. math:: \frac{2a}{\omega^2+a^2} \leftrightarrow e^{-a|x|} - - .. math:: \frac{1}{\omega+i} \leftrightarrow -i e^{-x} \theta(x) - - .. math:: \frac{1}{\omega-i} \leftrightarrow i e^{x} \theta(-x) - -For the inverse Fourier transform, the inverse procedure is used. - -In the library, :math:`a` is optimized according to the mesh properties (its size :math:`L=G.mesh().size()` and its precision :math:`\delta = G.mesh().delta()`). -The requirements are :math:`a \gg \delta\omega` and :math:`a \ll L\delta\omega`, or equivalently :math:`a \gg \delta t` and :math:`a \ll L\delta t`. -Thus, we chose :math:`a=\sqrt{L}\delta\omega` - - - +* +.. toctree:: + :maxdepth: 1 + + fourier_impl_notes + + + diff --git a/doc/reference/c++/gf/fourier_impl_notes.rst b/doc/reference/c++/gf/fourier_impl_notes.rst new file mode 100644 index 00000000..39661457 --- /dev/null +++ b/doc/reference/c++/gf/fourier_impl_notes.rst @@ -0,0 +1,126 @@ +.. highlight:: c + +Fourier : implementation notes +####################################### + + +The FFTW library +--------------------- + +Documentation on FFTW is on https://www.fftw.org. +FFTW is a C subroutine library for computing the discrete Fourier transform (DFT) in one or more dimensions, of arbitrary input size, and of both real and complex data. +It will be used to calculate the (inverse) Fourier transform, in real/imaginary time/frequency, using the fact that the GF values are stored for a finite amount of regularly spaced values. + +The DFT transforms of a sequence of :math:`N` complex numbers :math:`f_0...,f_{N-1}` into a sequence of :math:`N` complex numbers :math:`\tilde f_0...,\tilde f_{N-1}` according to the formula: + .. math:: \tilde f_k = \sum_{n=0}^{N-1} f_n e^{-i 2 \pi k n / N}. +The inverse DFT formula is + .. math:: f_n = \frac{1}{N} \sum_{k=0}^{N-1} \tilde f_k e^{i 2 \pi k n / N}. + + + +Implementation in real time/frequency using FFTW +--------------------------------------------------------------- + +The real time mesh parameters are :math:`t_{min}`, :math:`\delta t` and :math:`N_t`. +For the real frequency mesh, they are :math:`\omega_{min}`, :math:`\delta \omega` and :math:`N_\omega`. +The Fourier transform requires :math:`N_\omega=N_t` and :math:`\delta t \delta \omega= \frac{2\pi}{N_t}`. +The times are :math:`t_k=t_{min}+k\delta t` and the frequencies :math:`\omega_m=\omega_{min}+m\delta \omega`. + +By approximating Eq. :ref:`TF_R` by + .. math:: \tilde G(\omega_m) = \delta t \sum_{k=0}^{N_t} G(t_k) e^{i\omega_m t_k}, +we recognize an inverse DFT (Eq. :ref:`inv_DFT`). To calculate it using FFTW, we first need to prepare the input :math:`\tilde f_k`, then to do the DFT and finally to modify the output to obtain :math:`\tilde G(\omega_m)` using the two formulas: + .. math:: \tilde f_k = G(t_k) e^{i \omega_{min}t_k}, + .. math:: \tilde G(\omega_m) = \delta t f_m e^{i t_{min}(\omega_m-\omega_{min})}. + +Similarly, the inverse transformation is obtained by approximating Eq. :ref:`eq_inv_TF_R` by + .. math:: G(t_k)=\frac{\delta\omega}{2\pi}\sum_{m=0}^{N_\omega} \tilde G(\omega_m)e^{-i\omega_m t_k}, +we recognize a DFT (Eq. :ref:`DFT`). To calculate it using FFTW, we first need to prepare the input :math:`f_m`, then to do the inverse DFT and finally to modify the output to obtain :math:`G(t_k)`: + .. math:: f_m = \tilde G(\omega_m) e^{-i t_{min}\omega_m}, + .. math:: G(t_k) = \frac{1}{N_t \delta t}\tilde f_k e^{-i \omega_{min}(t_k-t_{min})}. + + + +Implementation in imaginary time/frequency using FFTW +--------------------------------------------------------------- + + +The imaginary time mesh parameters are :math:`\beta` and :math:`N_\tau`, plus a tag ``half_bins``, ``full_bins`` or ``without_last``. +In the ``full_bins`` case, one point of the time GF has to be removed for the fourier transform. +From these parameters, we deduce :math:`\delta\tau=\beta/N_\tau`. + +For the imaginary frequency mesh, the mesh parameters are :math:`\beta`, :math:`n_{min}` and :math:`N_{\omega_n}`. +From them, we deduce :math:`\delta\omega=\frac{2\pi}{\beta}`. + +The Fourier transform requires :math:`N_\omega=N_\tau`. +The times are :math:`\tau_k=\tau_{min}+k\delta\tau` and the frequencies :math:`\omega_n=\omega_{min}+n\delta \omega`. +:math:`\tau_{min}` is either 0 or :math:`\delta\tau/2` depending on the mesh kind. +:math:`\omega_{min}` is either :math:`\frac{2\pi(n_{min}+1)}{\beta}` or :math:`\frac{2\pi n_{min}}{\beta}` depending on the statistics. + +We approximate the TF and its inverse by + .. math:: \tilde G(i\omega_n) = \delta\tau \sum_{k=0}^{N_\tau} G(\tau_k)e^{i\omega_n \tau_k} + .. math:: G(\tau_k) = \sum_{n=n_{min}}^{N_\tau} \frac{1}{\beta} \tilde G(i\omega_n)e^{-i\omega_n \tau_k} + +We use for the TF: + .. math:: \tilde f_k = G(\tau_k) e^{i \omega_{min}\tau_k}, + .. math:: \tilde G(i\omega_n) = \frac{\beta}{N_\tau} f_n e^{i \tau_{min}(\omega_n-\omega_{min})}. +and for the inverse TF: + .. math:: f_m = \frac{1}{\beta}\tilde G(i\omega_n) e^{-i t_{min}\omega_n}, + .. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})}, + + +Special case of real functions in time for fermions +........................................................... + +In this case, :math:`G(i\omega_n)=conj(G(i\omega_n))` and we only store the values of :math:`G(i\omega_n)` for :math:`\omega_n > 0`. +The Eq. :ref:`inv_DFT_I` becomes: + + .. math:: G(\tau)=\sum_{n=0}^\infty \frac{2}{\beta} \tilde G(i\omega_n)\cos(\omega_n \tau) + +The inverse TF formulas are in this case + .. math:: f_m = \frac{2}{\beta}\tilde G(i\omega_n) e^{-i t_{min}\omega_n}, + .. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})}, + +Special case of real functions in time for bosons +........................................................... + +In this case, :math:`G(i\omega_n)=conj(G(i\omega_n))` and we only store the values of :math:`G(i\omega_n)` for :math:`\omega_n \ge 0`. +The Eq. :ref:`inv_DFT_I` becomes: + + .. math:: G(\tau)=\frac{1}{\beta} \tilde G(0)+\sum_{n=1}^\infty \frac{2}{\beta} \tilde G(i\omega_n)\cos(\omega_n \tau) + +The inverse TF formulas are in this case + .. math:: f_0 = \frac{1}{\beta}\tilde G(0), + .. math:: f_m = \frac{2}{\beta}\tilde G(i\omega_n) \cos(t_{min}\omega_n), + .. math:: G(t_k) = \tilde f_k e^{-i \omega_{min}(\tau_k-\tau_{min})}, + + +Effect of a TF on the tail +------------------------------------------ + +The tail is unchanged during a TF, but its value is used to limit the errors. + +The components 1 and 2 of the tail (:math:`t_1` and :math:`t_2`) are used to improve the computation of the GF in the following way: +in the large :math:`\omega` limit, + .. math:: G(\omega)\simeq \frac{t_1}{\omega}+\frac{t_2}{\omega^2}\simeq \frac{a_1}{\omega+i}+\frac{a_2}{\omega-i} +with :math:`a_1=\frac{t_1+it_2}{2}` and :math:`a_2=\frac{t_1-it_2}{2}`. + +As these large w terms are badly taken into account if we naively Fourier transform the function described by its values on the mesh in w, we substract them, do the Fourier transform and add their Fourier transform to the result. + +We use the following Fourier tranforms: + + .. math:: \frac{2a}{\omega^2+a^2} \leftrightarrow e^{-a|x|} + + .. math:: \frac{1}{\omega+i} \leftrightarrow -i e^{-x} \theta(x) + + .. math:: \frac{1}{\omega-i} \leftrightarrow i e^{x} \theta(-x) + +For the inverse Fourier transform, the inverse procedure is used. + +In the library, :math:`a` is optimized according to the mesh properties (its size :math:`L=G.mesh().size()` and its precision :math:`\delta = G.mesh().delta()`). +The requirements are :math:`a \gg \delta\omega` and :math:`a \ll L\delta\omega`, or equivalently :math:`a \gg \delta t` and :math:`a \ll L\delta t`. +Thus, we chose :math:`a=\sqrt{L}\delta\omega` + + + + + diff --git a/pytriqs/gf/local/functions.pxd b/pytriqs/gf/local/functions.pxd index f9e1bac7..8e35d880 100644 --- a/pytriqs/gf/local/functions.pxd +++ b/pytriqs/gf/local/functions.pxd @@ -8,22 +8,22 @@ cdef extern from "triqs/gfs/local/functions.hpp": ############### Fourier ######################### cdef extern from "triqs/gfs/local/fourier_matsubara.hpp" : - gf_imfreq lazy_fourier (gf_imtime & ) - gf_imtime lazy_inverse_fourier (gf_imfreq & ) + gf_imfreq fourier (gf_imtime & ) + gf_imtime inverse_fourier (gf_imfreq & ) cdef extern from "triqs/gfs/local/fourier_real.hpp" : - gf_refreq lazy_fourier (gf_retime & ) - gf_retime lazy_inverse_fourier (gf_refreq & ) - gf_refreq fourier (gf_retime & ) except + - gf_retime inverse_fourier (gf_refreq & ) except + + gf_refreq fourier (gf_retime & ) + gf_retime inverse_fourier (gf_refreq & ) + #gf_refreq fourier (gf_retime & ) except + + #gf_retime inverse_fourier (gf_refreq & ) except + ############### Legendre ######################### cdef extern from "triqs/gfs/local/legendre_matsubara.hpp" : - gf_imfreq lazy_legendre_imfreq (gf_legendre &) - gf_imtime lazy_legendre_imtime (gf_legendre &) - gf_legendre lazy_imfreq_legendre (gf_imfreq &) - gf_legendre lazy_imtime_legendre (gf_imtime &) + gf_imfreq legendre_to_imfreq (gf_legendre &) + gf_imtime legendre_to_imtime (gf_legendre &) + gf_legendre imfreq_to_legendre (gf_imfreq &) + gf_legendre imtime_to_legendre (gf_imtime &) ############### Pade ######################### diff --git a/pytriqs/gf/local/imfreq.pyx b/pytriqs/gf/local/imfreq.pyx index 76c3d712..83f8cf9e 100644 --- a/pytriqs/gf/local/imfreq.pyx +++ b/pytriqs/gf/local/imfreq.pyx @@ -12,11 +12,11 @@ cdef class GfImFreq_cython: def set_from_fourier(self,GfImTime_cython gt) : """Fills self with the Fourier transform of gt""" - self._c << lazy_fourier( gt._c ) + self._c << fourier( gt._c ) def set_from_legendre(self, GfLegendre_cython gl) : """Fills self with the Legendre transform of gl""" - self._c << lazy_legendre_imfreq(gl._c) + self._c << legendre_to_imfreq(gl._c) def density(self): return density(self._c).to_python() diff --git a/pytriqs/gf/local/imtime.pyx b/pytriqs/gf/local/imtime.pyx index d8d54c5e..ff20fa29 100644 --- a/pytriqs/gf/local/imtime.pyx +++ b/pytriqs/gf/local/imtime.pyx @@ -12,11 +12,11 @@ cdef class GfImTime_cython: def set_from_inverse_fourier(self,GfImFreq_cython gw) : """Fills self with the Inverse Fourier transform of gw""" - self._c << lazy_inverse_fourier( gw._c) + self._c << inverse_fourier( gw._c) def set_from_legendre(self, GfLegendre_cython gl) : """Fills self with the Legendre transform of gl""" - self._c << lazy_legendre_imtime(gl._c) + self._c << legendre_to_imtime(gl._c) def __dealloc__ (self): pass diff --git a/pytriqs/gf/local/legendre.pyx b/pytriqs/gf/local/legendre.pyx index 65ae5a59..1697183e 100644 --- a/pytriqs/gf/local/legendre.pyx +++ b/pytriqs/gf/local/legendre.pyx @@ -10,11 +10,11 @@ cdef class GfLegendre_cython: def set_from_imtime(self, GfImTime_cython gt) : """Fills self with the Legendre transform of gt""" - self._c << lazy_imtime_legendre(gt._c) + self._c << imtime_to_legendre(gt._c) def set_from_imfreq(self, GfImFreq_cython gw) : """Fills self with the Legendre transform of gw""" - self._c << lazy_imfreq_legendre(gw._c) + self._c << imfreq_to_legendre(gw._c) def density(self): return density(self._c).to_python() diff --git a/pytriqs/gf/local/refreq.pyx b/pytriqs/gf/local/refreq.pyx index 117ffc08..7b03b040 100644 --- a/pytriqs/gf/local/refreq.pyx +++ b/pytriqs/gf/local/refreq.pyx @@ -10,10 +10,11 @@ cdef class GfReFreq_cython: def set_from_fourier(self, GfReTime_cython gt) : """Fills self with the Fourier transform of gt""" - self._c << lazy_fourier( gt._c ) + self._c << fourier( gt._c ) - def inverse_fourier(self): - return make_GfReTime(inverse_fourier(self._c)) + # put if back with make_gf_from_fourier when approved + #def inverse_fourier(self): + # return make_GfReTime(inverse_fourier(self._c)) def set_from_pade(self, GfImFreq_cython gw, n_points = 100, freq_offset = 0.0) : pade(self._c, gw._c, n_points, freq_offset) diff --git a/pytriqs/gf/local/retime.pyx b/pytriqs/gf/local/retime.pyx index c2cc690a..26212767 100644 --- a/pytriqs/gf/local/retime.pyx +++ b/pytriqs/gf/local/retime.pyx @@ -8,12 +8,13 @@ cdef class GfReTime_cython: def __write_hdf5_cython__ (self, gr , char * key) : h5_write (make_h5_group(gr), key, self._c) - def fourier(self): - return make_GfReFreq(fourier(self._c)) + # Put it back ? + #def fourier(self): + # return make_GfReFreq(fourier(self._c)) def set_from_inverse_fourier(self, GfReFreq_cython gw) : """Fills self with the Inverse Fourier transform of gw""" - self._c << lazy_inverse_fourier( gw._c) + self._c << inverse_fourier( gw._c) def __dealloc__ (self): pass diff --git a/test/triqs/gfs/curry_and_fourier.cpp b/test/triqs/gfs/curry_and_fourier.cpp index bb39bd11..6755e62c 100644 --- a/test/triqs/gfs/curry_and_fourier.cpp +++ b/test/triqs/gfs/curry_and_fourier.cpp @@ -35,11 +35,11 @@ try { auto G_w_wn_curry0 = curry<0>(G_w_wn); auto G_w_tau_curry0 = curry<0>(G_w_tau); - for (auto const & w : G_w_wn_curry0.mesh()) G_w_wn_curry0[w] = lazy_fourier(G_w_tau_curry0[w]); + for (auto const & w : G_w_wn_curry0.mesh()) G_w_wn_curry0[w] = fourier(G_w_tau_curry0[w]); - G_w_wn_curry0[w_] << lazy_fourier(G_w_tau_curry0[w_]); + G_w_wn_curry0[w_] << fourier(G_w_tau_curry0[w_]); - curry<0>(G_w_wn) [w_] << lazy_fourier(curry<0>(G_w_tau)[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; diff --git a/test/triqs/gfs/fourier1.cpp b/test/triqs/gfs/fourier1.cpp index ae0e85e5..9d61b17a 100644 --- a/test/triqs/gfs/fourier1.cpp +++ b/test/triqs/gfs/fourier1.cpp @@ -12,13 +12,13 @@ int main() { auto G = gf {{beta, Fermion}, {2,2}}; auto Gc = G; auto G3 = G; - auto Gt = gf {{beta, Fermion}, {2,2}}; + auto Gt = gf {{beta, Fermion,100}, {2,2}}; - auto gt = inverse_fourier(G); - auto gw = fourier(gt); + auto gt = make_gf_from_inverse_fourier(G); + auto gw = make_gf_from_fourier(gt); - //gw() = lazy_fourier(gt); - G() = lazy_fourier(Gt); + //gw() = fourier(gt); + G() = fourier(Gt); } diff --git a/test/triqs/gfs/gf_notail.cpp b/test/triqs/gfs/gf_notail.cpp new file mode 100644 index 00000000..46944b00 --- /dev/null +++ b/test/triqs/gfs/gf_notail.cpp @@ -0,0 +1,38 @@ +#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK +#include +using namespace triqs::gfs; +using namespace triqs::arrays; +#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl; +#include + +int main() { + try { + double beta = 1, a = 1; + + // construct + auto gw_n = gf{{beta, Fermion}, {2, 2}}; + auto gt_n = gf{{beta, Fermion, 100}, {2, 2}}; + + // make a view from a g with a tail + auto gw = gf{{beta, Fermion}, {2, 2}}; + auto gt = gf{gf_mesh{beta, Fermion, 100}, {2, 2}}; + + auto vw = make_gf_view_without_tail(gw); + auto vt = make_gf_view_without_tail(gt); + + triqs::clef::placeholder<0> tau_; + // should not and does not compile : improve error message ??? + // gt(tau_) << exp(-a * tau_) / (1 + exp(-beta * a)); + vt(tau_) << exp(-a * tau_) / (1 + exp(-beta * a)); + + // test hdf5 + H5::H5File file("ess_g_notail.h5", H5F_ACC_TRUNC); + h5_write(file, "g", vt); + + // rebuilding a new gf... + auto g3 = make_gf_from_g_and_tail(vw, gw.singularity()); + // need to test all this.... + + } + TRIQS_CATCH_AND_ABORT; +} diff --git a/test/triqs/gfs/gfv2.cpp b/test/triqs/gfs/gfv2.cpp index 733d5eb8..e55e68ca 100644 --- a/test/triqs/gfs/gfv2.cpp +++ b/test/triqs/gfs/gfv2.cpp @@ -20,7 +20,7 @@ int main() { auto G = gf{ {beta, Fermion}, {2,2} }; auto Gc = gf{ {beta, Fermion}, {2,2} }; auto G3 = gf{ {beta, Fermion}, {2,2} }; - auto Gt = gf{ {beta, Fermion}, {2,2} }; + auto Gt = gf{ {beta, Fermion,100}, {2,2} }; auto Gv = G(); TEST( G( 0) ) ; diff --git a/test/triqs/gfs/test_fourier_matsubara.cpp b/test/triqs/gfs/test_fourier_matsubara.cpp index 1afb0d68..2e98779c 100644 --- a/test/triqs/gfs/test_fourier_matsubara.cpp +++ b/test/triqs/gfs/test_fourier_matsubara.cpp @@ -5,6 +5,14 @@ using namespace triqs::arrays; #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) < +namespace triqs { namespace gfs { + // defined in the cpp file + void fourier_impl (gf_view gw , gf_const_view gt, scalar_valued); + void fourier_impl (gf_view gw , gf_const_view gt, matrix_valued); + void inverse_fourier_impl (gf_view gt, gf_const_view gw, scalar_valued); +void inverse_fourier_impl (gf_view gt, gf_const_view gw, matrix_valued); +}} + int main() { double precision=10e-9; @@ -45,9 +53,9 @@ int main() { } h5_write(file,"Gt1b",Gt1); // must be 0 - ///to verify that lazy_fourier computes + ///to verify that fourier computes auto Gw2 = gf {{beta, Fermion}, {1,1}}; - Gw2() = lazy_fourier(Gt1); + Gw2() = fourier(Gt1); } diff --git a/test/triqs/gfs/test_fourier_real_time.cpp b/test/triqs/gfs/test_fourier_real_time.cpp index 354a4c4b..3580ca84 100644 --- a/test/triqs/gfs/test_fourier_real_time.cpp +++ b/test/triqs/gfs/test_fourier_real_time.cpp @@ -33,11 +33,11 @@ int main() { Gw1.singularity()(2)=triqs::arrays::matrix{{2.0*a}}; h5_write(file,"Gw1",Gw1); // the original lorentzian - auto Gt1 = inverse_fourier(Gw1); + auto Gt1 = make_gf_from_inverse_fourier(Gw1); h5_write(file,"Gt1",Gt1); // the lorentzian TF : lorentzian_inverse // verification that TF(TF^-1)=Id - auto Gw1b = fourier(Gt1); + auto Gw1b = make_gf_from_fourier(Gt1); for(auto const & w:Gw1b.mesh()){ Gw1b[w]-=Gw1[w]; if ( std::abs(Gw1b[w](0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<{{1.0}}; h5_write(file,"Gt2",Gt2); - auto Gw2 = fourier(Gt2); + auto Gw2 = make_gf_from_fourier(Gt2); h5_write(file,"Gw2",Gw2); for(auto const & w:Gw2.mesh()){ @@ -86,7 +86,7 @@ int main() { //for(auto const & t:Gt3.mesh()) Gt3[t] = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5_j*std::sin(8*t+0.3*acos(-1.)) ; h5_write(file,"Gt3",Gt3); - auto Gw3 = fourier(Gt3); + auto Gw3 = make_gf_from_fourier(Gt3); h5_write(file,"Gw3",Gw3); } diff --git a/test/triqs/gfs/test_gf_triqs.cpp b/test/triqs/gfs/test_gf_triqs.cpp index fe982e17..5e57b1f6 100644 --- a/test/triqs/gfs/test_gf_triqs.cpp +++ b/test/triqs/gfs/test_gf_triqs.cpp @@ -33,7 +33,7 @@ void test_0(){ /* ---------- Fourier transform ---------------------*/ auto Gt = gf {{beta, Fermion, Ntau, full_bins}, {1,1}}; - Gt() = lazy_inverse_fourier(G); + Gt() = inverse_fourier(G); TEST(Gt(0.0)); @@ -52,7 +52,7 @@ void test_1(){ auto Gw = gf {{beta, Fermion, 100}, {1,1}}; Gw.singularity()(1) = 1; - Gt() = lazy_inverse_fourier(Gw); + Gt() = inverse_fourier(Gw); } int main() { diff --git a/triqs/gfs.hpp b/triqs/gfs.hpp index 8ac5a13f..e87be1a1 100644 --- a/triqs/gfs.hpp +++ b/triqs/gfs.hpp @@ -30,5 +30,8 @@ #include #include +#include +#include +#include #endif diff --git a/triqs/gfs/domains/matsubara.hpp b/triqs/gfs/domains/matsubara.hpp index de757ad4..5ef02afe 100644 --- a/triqs/gfs/domains/matsubara.hpp +++ b/triqs/gfs/domains/matsubara.hpp @@ -30,7 +30,8 @@ namespace triqs { namespace gfs { typedef typename mpl::if_c, double>::type point_t; double beta; statistic_enum statistic; - matsubara_domain (double Beta=1, statistic_enum s = Fermion): beta(Beta), statistic(s){ + matsubara_domain() = default; + matsubara_domain (double Beta, statistic_enum s): beta(Beta), statistic(s){ if(beta<0)TRIQS_RUNTIME_ERROR<<"Matsubara domain construction : beta <0 : beta ="<< beta <<"\n"; } matsubara_domain(matsubara_domain const &) = default; diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp index c4b7552e..5adc9bb5 100644 --- a/triqs/gfs/gf.hpp +++ b/triqs/gfs/gf.hpp @@ -94,6 +94,9 @@ namespace triqs { namespace gfs { template auto get_gf_data_shape(G const &g) DECL_AND_RETURN(g.get_data_shape()); + template + auto get_target_shape(gf_impl const &g) DECL_AND_RETURN(g.data().shape().front_pop()); + // ---------------------- implementation -------------------------------- // overload get_shape for a vector to simplify code below in gf block case. @@ -334,8 +337,10 @@ namespace triqs { namespace gfs { template friend void triqs_clef_auto_assign_subscript (gf_impl & g, RHS rhs) { triqs_clef_auto_assign(g,rhs);} private: - template void triqs_clef_auto_assign_impl (RHS const & rhs, std::integral_constant) { - for (auto const & w: this->mesh()) (*this)[w] = rhs(w); + template void triqs_clef_auto_assign_impl(RHS const &rhs, std::integral_constant) { + for (auto const &w : this->mesh()) { + (*this)[w] = rhs(w); + } //for (auto const & w: this->mesh()) (*this)[w] = rhs(typename B::mesh_t::mesh_point_t::cast_t(w)); } template void triqs_clef_auto_assign_impl(RHS const &rhs, std::integral_constant) { @@ -582,7 +587,11 @@ namespace triqs { namespace gfs { template struct factories { typedef gf gf_t; - struct target_shape_t {}; + struct target_shape_t { + target_shape_t front_pop() const { // this make the get_target_shape function works in this case... + return {}; + } + }; typedef typename gf_t::mesh_t mesh_t; static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) { @@ -620,6 +629,7 @@ namespace triqs { namespace gfs { }; } // gfs_implementation + }} // same as for arrays : views can not be swapped by the std::swap. Delete it diff --git a/triqs/gfs/imfreq.hpp b/triqs/gfs/imfreq.hpp index e68cb44b..8a16f731 100644 --- a/triqs/gfs/imfreq.hpp +++ b/triqs/gfs/imfreq.hpp @@ -23,6 +23,7 @@ #include "./tools.hpp" #include "./gf.hpp" #include "./local/tail.hpp" +#include "./local/no_tail.hpp" #include "./domains/matsubara.hpp" #include "./meshes/linear.hpp" #include "./evaluators.hpp" @@ -34,16 +35,17 @@ namespace triqs { namespace gfs { typedef linear_mesh> B; static double m1(double beta) { return std::acos(-1)/beta;} gf_mesh() = default; + gf_mesh(B const &x) : B(x) {} // enables also construction from another Opt gf_mesh (typename B::domain_t const & d, int Nmax = 1025) : B(d, d.statistic==Fermion?m1(d.beta):0, d.statistic==Fermion?(2*Nmax+1)*m1(d.beta): 2*Nmax*m1(d.beta), Nmax, without_last){} - gf_mesh (double beta, statistic_enum S, int Nmax = 1025) : gf_mesh({beta,S}, Nmax){} + gf_mesh (double beta, statistic_enum S, int Nmax = 1025) : gf_mesh({beta,S}, Nmax){} }; namespace gfs_implementation { //singularity - template struct singularity { typedef local::tail type;}; - template struct singularity { typedef local::tail type;}; + template<> struct singularity { typedef local::tail type;}; + template<> struct singularity { typedef local::tail type;}; //h5 name template struct h5_name { static std::string invoke(){ return "ImFreq";}}; diff --git a/triqs/gfs/imtime.hpp b/triqs/gfs/imtime.hpp index 828af163..a6447aa9 100644 --- a/triqs/gfs/imtime.hpp +++ b/triqs/gfs/imtime.hpp @@ -23,6 +23,7 @@ #include "./tools.hpp" #include "./gf.hpp" #include "./local/tail.hpp" +#include "./local/no_tail.hpp" #include "./domains/matsubara.hpp" #include "./meshes/linear.hpp" #include "./evaluators.hpp" @@ -35,15 +36,16 @@ namespace triqs { namespace gfs { template struct gf_mesh : linear_mesh> { typedef linear_mesh> B; gf_mesh() = default; + gf_mesh(B const &x) : B(x) {} // enables also construction from another Opt gf_mesh(typename B::domain_t d, int n_time_slices, mesh_kind mk = half_bins) : B(d, 0, d.beta, n_time_slices, mk) {} gf_mesh(double beta, statistic_enum S, int n_time_slices, mesh_kind mk = half_bins) : gf_mesh({beta, S}, n_time_slices, mk) {} }; namespace gfs_implementation { - // singularity - template struct singularity { typedef local::tail type;}; - template struct singularity { typedef local::tail type;}; + // singularity. If no_tail is given, then it is the default (nothing) + template<> struct singularity { typedef local::tail type;}; + template<> struct singularity { typedef local::tail type;}; // h5 name template struct h5_name { static std::string invoke(){ return "ImTime";}}; @@ -92,6 +94,7 @@ namespace triqs { namespace gfs { template struct evaluator : evaluator_one_var{}; } // gfs_implementation. + }} #endif diff --git a/triqs/gfs/local/fourier_matsubara.cpp b/triqs/gfs/local/fourier_matsubara.cpp index 78bc574c..a4aeea26 100644 --- a/triqs/gfs/local/fourier_matsubara.cpp +++ b/triqs/gfs/local/fourier_matsubara.cpp @@ -1,5 +1,5 @@ /******************************************************************************* - * + * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * * Copyright (C) 2011 by M. Ferrero, O. Parcollet @@ -18,168 +18,203 @@ * TRIQS. If not, see . * ******************************************************************************/ - +#include "fourier_base.hpp" #include "fourier_matsubara.hpp" #include -namespace triqs { namespace gfs { - - namespace impl_local_matsubara { - - inline dcomplex oneFermion(dcomplex a,double b,double tau,double beta) { - return -a*( b >=0 ? exp(-b*tau)/(1+exp(-beta*b)) : exp(b*(beta-tau))/(1+exp(beta*b)) ); - } - - inline dcomplex oneBoson(dcomplex a,double b,double tau,double beta) { - return a*( b >=0 ? exp(-b*tau)/(exp(-beta*b)-1) : exp(b*(beta-tau))/(1-exp(b*beta)) ); - } - } - - template GfElementType convert_green ( dcomplex const& x) { return x;} - template<> double convert_green ( dcomplex const& x) { return real(x);} - +namespace triqs { +namespace gfs { + + template GfElementType convert_green(dcomplex const& x) { return x; } + template <> double convert_green(dcomplex const& x) { return real(x); } + //-------------------------------------------------------------------------------------- - + struct impl_worker { - + tqa::vector g_in, g_out; - - void direct (gf_view gw, gf_const_view gt) { - using namespace impl_local_matsubara; + + dcomplex oneFermion(dcomplex a, double b, double tau, double beta) { + return -a * (b >= 0 ? exp(-b * tau) / (1 + exp(-beta * b)) : exp(b * (beta - tau)) / (1 + exp(beta * b))); + } + + dcomplex oneBoson(dcomplex a, double b, double tau, double beta) { + return a * (b >= 0 ? exp(-b * tau) / (exp(-beta * b) - 1) : exp(b * (beta - tau)) / (1 - exp(b * beta))); + } + + //------------------------------------- + + void direct(gf_view gw, gf_const_view gt) { auto ta = gt(freq_infty()); - //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; + 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) =delete; + + //------------------------------------- + + private: + void direct_impl(gf_view gw, gf_const_view gt, + local::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; dcomplex a1, a2, a3; - double beta=gt.mesh().domain().beta; - auto L = ( gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size()); - double fact= beta/ gt.mesh().size(); - dcomplex iomega = dcomplex(0.0,1.0) * std::acos(-1) / beta; - dcomplex iomega2 = iomega * 2 * gt.mesh().delta() * ( gt.mesh().kind() == half_bins ? 0.5 : 0.0); + double beta = gt.mesh().domain().beta; + auto L = (gt.mesh().kind() == full_bins ? gt.mesh().size() - 1 : gt.mesh().size()); + double fact = beta / gt.mesh().size(); + dcomplex iomega = dcomplex(0.0, 1.0) * std::acos(-1) / beta; + dcomplex iomega2 = iomega * 2 * gt.mesh().delta() * (gt.mesh().kind() == half_bins ? 0.5 : 0.0); g_in.resize(gt.mesh().size()); g_out.resize(gw.mesh().size()); - if (gw.domain().statistic == Fermion){ - b1 = 0; b2 =1; b3 =-1; - a1 = d-B; a2 = (A+B)/2; a3 = (B-A)/2; + if (gw.domain().statistic == Fermion) { + b1 = 0; + b2 = 1; + b3 = -1; + a1 = d - B; + a2 = (A + B) / 2; + a3 = (B - A) / 2; + } else { + b1 = -0.5; + b2 = -1; + b3 = 1; + a1 = 4 * (d - B) / 3; + a2 = B - (d + A) / 2; + a3 = d / 6 + A / 2 + B / 3; } - else { - b1 = -0.5; b2 =-1; b3 =1; - a1 = 4*(d-B)/3; a2 = B-(d+A)/2; a3 = d/6+A/2+B/3; - } - if (gw.domain().statistic == Fermion){ - for (auto & t : gt.mesh()) - g_in[t.index()] = fact * exp(iomega*t) * ( gt[t] - ( oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) ) ); - } - else { - for (auto & t : gt.mesh()) - g_in[t.index()] = fact * ( gt[t] - ( oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) ) ); + if (gw.domain().statistic == Fermion) { + for (auto& t : gt.mesh()) + g_in[t.index()] = fact * exp(iomega * t) * + (gt[t] - (oneFermion(a1, b1, t, beta) + oneFermion(a2, b2, t, beta) + oneFermion(a3, b3, t, beta))); + } else { + for (auto& t : gt.mesh()) + g_in[t.index()] = fact * (gt[t] - (oneBoson(a1, b1, t, beta) + oneBoson(a2, b2, t, beta) + oneBoson(a3, b3, t, beta))); } details::fourier_base(g_in, g_out, L, true); - for (auto & w : gw.mesh()) { - gw[w] = g_out( w.index() ) * exp(iomega2*w.index() ) + a1/(w-b1) + a2/(w-b2) + a3/(w-b3); + for (auto& w : gw.mesh()) { + gw[w] = g_out(w.index()) * exp(iomega2 * w.index()) + a1 / (w - b1) + a2 / (w - b2) + a3 / (w - b3); } - gw.singularity() = gt.singularity();// set tail } - - void inverse(gf_view gt, gf_const_view gw){ - using namespace impl_local_matsubara; + + public: + //------------------------------------- + + void inverse(gf_view gt, gf_const_view gw) { 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()); - //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); + // 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; dcomplex a1, a2, a3; - - double beta=gw.domain().beta; - size_t L= gt.mesh().size() - ( gt.mesh().kind() == full_bins ? 1 : 0); //L can be different from gt.mesh().size() (depending on the mesh kind) and is given to the FFT algorithm - dcomplex iomega = dcomplex(0.0,1.0) * std::acos(-1) / beta; - dcomplex iomega2 = -iomega * 2 * gt.mesh().delta() * (gt.mesh().kind() == half_bins ? 0.5 : 0.0) ; - double fact = (Green_Function_Are_Complex_in_time ? 1 : 2)/beta; - g_in.resize( gw.mesh().size()); + + double beta = gw.domain().beta; + size_t L = gt.mesh().size() - (gt.mesh().kind() == full_bins ? 1 : 0); // L can be different from gt.mesh().size() (depending + // on the mesh kind) and is given to the FFT algorithm + dcomplex iomega = dcomplex(0.0, 1.0) * std::acos(-1) / beta; + dcomplex iomega2 = -iomega * 2 * gt.mesh().delta() * (gt.mesh().kind() == half_bins ? 0.5 : 0.0); + double fact = (Green_Function_Are_Complex_in_time ? 1 : 2) / beta; + g_in.resize(gw.mesh().size()); g_out.resize(gt.mesh().size()); - - if (gw.domain().statistic == Fermion){ - b1 = 0; b2 =1; b3 =-1; - a1 = d-B; a2 = (A+B)/2; a3 = (B-A)/2; - } - else { - b1 = -0.5; b2 =-1; b3 =1; - a1=4*(d-B)/3; a2=B-(d+A)/2; a3=d/6+A/2+B/3; + + if (gw.domain().statistic == Fermion) { + b1 = 0; + b2 = 1; + b3 = -1; + a1 = d - B; + a2 = (A + B) / 2; + a3 = (B - A) / 2; + } else { + b1 = -0.5; + b2 = -1; + b3 = 1; + a1 = 4 * (d - B) / 3; + a2 = B - (d + A) / 2; + a3 = d / 6 + A / 2 + B / 3; } g_in() = 0; - for (auto & w: gw.mesh()) { - g_in[ w.index() ] = fact * exp(w.index()*iomega2) * ( gw[w] - (a1/(w-b1) + a2/(w-b2) + a3/(w-b3)) ); + for (auto& w : gw.mesh()) { + g_in[w.index()] = fact * exp(w.index() * iomega2) * (gw[w] - (a1 / (w - b1) + a2 / (w - b2) + a3 / (w - b3))); } // for bosons GF(w=0) is divided by 2 to avoid counting it twice - if (gw.domain().statistic == Boson && !Green_Function_Are_Complex_in_time ) g_in(0) *= 0.5; - + if (gw.domain().statistic == Boson && !Green_Function_Are_Complex_in_time) g_in(0) *= 0.5; + details::fourier_base(g_in, g_out, L, false); - + // CORRECT FOR COMPLEX G(tau) !!! typedef double gt_result_type; - //typedef typename gf::mesh_type::gf_result_type gt_result_type; - if (gw.domain().statistic == Fermion){ - for (auto & t : gt.mesh()){ - gt[t] = convert_green ( g_out( t.index() == L ? 0 : t.index() ) * exp(-iomega*t) - + oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) ); + // typedef typename gf::mesh_type::gf_result_type gt_result_type; + if (gw.domain().statistic == Fermion) { + for (auto& t : gt.mesh()) { + gt[t] = + convert_green(g_out(t.index() == L ? 0 : t.index()) * exp(-iomega * t) + oneFermion(a1, b1, t, beta) + + oneFermion(a2, b2, t, beta) + oneFermion(a3, b3, t, beta)); } + } else { + for (auto& t : gt.mesh()) + gt[t] = convert_green(g_out(t.index() == L ? 0 : t.index()) + oneBoson(a1, b1, t, beta) + + oneBoson(a2, b2, t, beta) + oneBoson(a3, b3, t, beta)); } - else { - for (auto & t : gt.mesh()) - gt[t] = convert_green ( g_out( t.index() == L ? 0 : t.index() ) - + oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) ); - } - if (gt.mesh().kind() == full_bins) - gt.on_mesh(L) = -gt.on_mesh(0)-convert_green(ta(1)(0,0)); + if (gt.mesh().kind() == full_bins) gt.on_mesh(L) = -gt.on_mesh(0) - convert_green(ta(1)(0, 0)); // set tail - gt.singularity() = gw.singularity(); + gt.singularity() = gw.singularity(); } - - }; - - void fourier_impl (gf_view gw , gf_const_view gt, scalar_valued){ + + }; // class worker + + //-------------------------------------------- + + void fourier_impl(gf_view gw, gf_const_view gt, scalar_valued) { impl_worker w; - w.direct(gw,gt); + w.direct(gw, gt); } - - void fourier_impl (gf_view gw , gf_const_view gt, matrix_valued){ + + void fourier_impl(gf_view gw, gf_const_view gt, matrix_valued) { impl_worker w; - for (size_t n1=0; n1 gt , gf_const_view gw, scalar_valued){ + + void inverse_fourier_impl(gf_view gt, gf_const_view gw, scalar_valued) { impl_worker w; - w.inverse(gt,gw); + w.inverse(gt, gw); } - - void inverse_fourier_impl (gf_view gt , gf_const_view gw, matrix_valued){ + + 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, scalar_valued() ); } - void triqs_gf_view_assign_delegation( gf_view 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) { inverse_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) { + fourier_impl(g, L.g, scalar_valued()); + } + void triqs_gf_view_assign_delegation(gf_view 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) { + inverse_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()); + } +} +} diff --git a/triqs/gfs/local/fourier_matsubara.hpp b/triqs/gfs/local/fourier_matsubara.hpp index 5ba7045a..a940a69c 100644 --- a/triqs/gfs/local/fourier_matsubara.hpp +++ b/triqs/gfs/local/fourier_matsubara.hpp @@ -1,5 +1,5 @@ /******************************************************************************* - * + * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * * Copyright (C) 2011 by M. Ferrero, O. Parcollet @@ -18,67 +18,62 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_LOCAL_FOURIER_MATSU_H +#ifndef TRIQS_GF_LOCAL_FOURIER_MATSU_H #define TRIQS_GF_LOCAL_FOURIER_MATSU_H #include "fourier_base.hpp" -#include -#include +#include +#include -namespace triqs { namespace gfs { - - // First the implementation of the fourier transform - void fourier_impl (gf_view gw , gf_const_view gt, scalar_valued); - void fourier_impl (gf_view gw , gf_const_view gt, matrix_valued); - void inverse_fourier_impl (gf_view gt, gf_const_view gw, scalar_valued); - void inverse_fourier_impl (gf_view gt, gf_const_view gw, matrix_valued); - - inline gf_view fourier (gf_const_view gt) { - int L = (gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size() ); - auto gw = gf{ {gt.domain(),L}, gt.data().shape().front_pop() }; - auto V = gw(); - fourier_impl(V, gt, matrix_valued()); +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}; + } + + 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); + + template gf_mesh make_mesh_fourier_compatible(gf_mesh const& m) { + int L = m.size() - (m.kind() == full_bins ? 1 : 0); + return {m.domain(), L}; + } + + template + gf_mesh make_mesh_fourier_compatible(gf_mesh const& m, mesh_kind mk = full_bins) { + int L = m.size() + (mk == full_bins ? 1 : 0); + return {m.domain(), L}; + } + + template + gf_view make_gf_from_fourier(gf_impl const& gt) { + auto gw = gf{make_mesh_fourier_compatible(gt.mesh()), get_target_shape(gt)}; + gw() = fourier(gt); return gw; } - inline gf_view fourier (gf_const_view gt) { - int L = (gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size() ); - auto gw = gf{ {gt.domain(),L} }; - auto V = gw(); - fourier_impl(V, gt, scalar_valued()); - return gw; - } - - inline gf_view inverse_fourier (gf_const_view gw, mesh_kind mk = half_bins) { - double pi = std::acos(-1); - int L = (mk == full_bins ? gw.mesh().size()+1 : gw.mesh().size() ); - auto gt = gf{ {gw.domain(),L}, gw.data().shape().front_pop()}; - auto V = gt(); - inverse_fourier_impl(V, gw, matrix_valued()); - return gt; - } - inline gf_view inverse_fourier (gf_const_view gw, mesh_kind mk = half_bins) { - double pi = std::acos(-1); - int L = (mk == full_bins ? gw.mesh().size()+1 : gw.mesh().size() ); - auto gt = gf{ {gw.domain(),L} }; - auto V = gt(); - inverse_fourier_impl(V, gw,scalar_valued()); - return gt; - } - - inline gf_keeper lazy_fourier (gf_const_view g) { return {g};} - inline gf_keeper lazy_inverse_fourier (gf_const_view g) { return {g};} - inline gf_keeper lazy_fourier (gf_const_view g) { return {g};} - inline gf_keeper lazy_inverse_fourier (gf_const_view 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); - 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); -}} -namespace triqs { namespace clef { -TRIQS_CLEF_MAKE_FNT_LAZY (lazy_fourier); -TRIQS_CLEF_MAKE_FNT_LAZY (lazy_inverse_fourier); -}} + template + gf_view make_gf_from_inverse_fourier(gf_impl const& gw, mesh_kind mk = full_bins) { + auto gt = gf{make_mesh_fourier_compatible(gw.mesh(), mk), get_target_shape(gw)}; + gt() = inverse_fourier(gw); + return gt; + } +} +} + +namespace triqs { +namespace clef { + TRIQS_CLEF_MAKE_FNT_LAZY(fourier); + TRIQS_CLEF_MAKE_FNT_LAZY(inverse_fourier); +} +} #endif diff --git a/triqs/gfs/local/fourier_real.hpp b/triqs/gfs/local/fourier_real.hpp index fd4fe649..63d1a66c 100644 --- a/triqs/gfs/local/fourier_real.hpp +++ b/triqs/gfs/local/fourier_real.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 @@ -27,63 +27,50 @@ namespace triqs { namespace gfs { - // First the implementation of the fourier transform - void fourier_impl (gf_view gw , gf_const_view gt, scalar_valued); - void fourier_impl (gf_view gw , gf_const_view gt, matrix_valued); - void inverse_fourier_impl (gf_view gt, gf_const_view gw, scalar_valued); - void inverse_fourier_impl (gf_view gt, gf_const_view gw, matrix_valued); - - inline gf_view fourier (gf_const_view gt) { - double pi = std::acos(-1); - int L = gt.mesh().size(); - double wmin = -pi * (L-1) / (L*gt.mesh().delta()); - double wmax = pi * (L-1) / (L*gt.mesh().delta()); - auto gw = gf{ {wmin, wmax, L}, gt.data().shape().front_pop()}; - auto V = gw(); - fourier_impl(V, gt, matrix_valued()); - return gw; + template + gf_keeper fourier(gf_impl const& g) { + return {g}; } - inline gf_view fourier (gf_const_view gt) { - double pi = std::acos(-1); - int L = gt.mesh().size(); - double wmin = -pi * (L-1) / (L*gt.mesh().delta()); - double wmax = pi * (L-1) / (L*gt.mesh().delta()); - auto gw = gf{ {wmin, wmax, L} }; - auto V = gw(); - fourier_impl(V, gt, scalar_valued()); - return gw; - } - - inline gf_view inverse_fourier (gf_const_view gw) { - double pi = std::acos(-1); - int L = gw.mesh().size(); - double tmin = -pi * (L-1) / (L*gw.mesh().delta()); - double tmax = pi * (L-1) / (L*gw.mesh().delta()); - auto gt = 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_const_view gw) { - double pi = std::acos(-1); - int L = gw.mesh().size(); - double tmin = -pi * (L-1) / (L*gw.mesh().delta()); - double tmax = pi * (L-1) / (L*gw.mesh().delta()); - auto gt = gf{ {tmin, tmax, L} }; - auto V = gt(); - inverse_fourier_impl(V, gw, scalar_valued()); - return gt; + template + gf_keeper inverse_fourier(gf_impl const& g) { + return {g}; } - inline gf_keeper lazy_fourier(gf_const_view g) { return {g}; } - inline gf_keeper lazy_inverse_fourier(gf_const_view g) { return {g}; } - inline gf_keeper lazy_fourier(gf_const_view g) { return {g}; } - inline gf_keeper lazy_inverse_fourier(gf_const_view 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); + 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 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); + template gf_mesh make_mesh_fourier_compatible(gf_mesh const& m) { + int L = m.size(); + double pi = std::acos(-1); + double wmin = -pi * (L - 1) / (L * m.delta()); + double wmax = pi * (L - 1) / (L * m.delta()); + return {wmin, wmax, L}; + } + + template + gf_mesh make_mesh_fourier_compatible(gf_mesh const& m, mesh_kind mk = full_bins) { + double pi = std::acos(-1); + int L = m.size(); + double tmin = -pi * (L-1) / (L*m.delta()); + double tmax = pi * (L-1) / (L*m.delta()); + return {tmin, tmax, L}; + } + + template + gf_view make_gf_from_fourier(gf_impl const& gt) { + auto gw = gf{make_mesh_fourier_compatible(gt.mesh()), get_target_shape(gt)}; + gw() = fourier(gt); + return gw; + } + + template + gf_view make_gf_from_inverse_fourier(gf_impl const& gw) { + auto gt = gf{make_mesh_fourier_compatible(gw.mesh()), get_target_shape(gw)}; + gt() = inverse_fourier(gw); + return gt; + } }} #endif diff --git a/triqs/gfs/local/legendre_matsubara.cpp b/triqs/gfs/local/legendre_matsubara.cpp index 5fd95c1c..8583c64a 100644 --- a/triqs/gfs/local/legendre_matsubara.cpp +++ b/triqs/gfs/local/legendre_matsubara.cpp @@ -1,9 +1,8 @@ - /******************************************************************************* * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by L. Boehnke, M. Ferrero, O. Parcollet + * Copyright (C) 2011-2014 by L. Boehnke, 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 @@ -27,25 +26,63 @@ using namespace triqs::utility; -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { -void legendre_matsubara_direct(gf_view & gw, gf_const_view gl) { + // ---------------------------- + + void legendre_matsubara_direct(gf_view gw, gf_const_view gl) { gw() = 0.0; triqs::arrays::range R; // Use the transformation matrix - for (auto om: gw.mesh()) { - for (auto l: gl.mesh()) { - gw[om] += legendre_T(om.index(),l.index()) * gl[l]; - } + for (auto om : gw.mesh()) { + for (auto l : gl.mesh()) { + gw[om] += legendre_T(om.index(), l.index()) * gl[l]; + } } gw.singularity() = get_tail(gl, gw.singularity().size(), gw.singularity().order_min()); + } -} + // ---------------------------- -void legendre_matsubara_inverse (gf_view & gl, gf_const_view gw) { + void legendre_matsubara_direct(gf_view gt, gf_const_view gl) { + + gt() = 0.0; + legendre_generator L; + + for (auto t : gt.mesh()) { + L.reset(2 * t / gt.domain().beta - 1); + for (auto l : gl.mesh()) { + gt[t] += sqrt(2 * l.index() + 1) / gt.domain().beta * gl[l] * L.next(); + } + } + + gt.singularity() = get_tail(gl, gt.singularity().size(), gt.singularity().order_min()); + } + + // ---------------------------- + + void legendre_matsubara_inverse(gf_view gl, gf_const_view gt) { + + gl() = 0.0; + legendre_generator L; + + // Do the integral over imaginary time + for (auto t : gt.mesh()) { + L.reset(2 * t / gt.domain().beta - 1); + for (auto l : gl.mesh()) { + gl[l] += sqrt(2 * l.index() + 1) * L.next() * gt[t]; + } + } + gl.data() *= gt.mesh().delta(); + } + + // ---------------------------- + + void legendre_matsubara_inverse(gf_view gl, gf_const_view gw) { gl() = 0.0; @@ -53,64 +90,25 @@ void legendre_matsubara_inverse (gf_view & gl, gf_const_view // I set Nt time bins. This is ugly, one day we must code the direct // transformation without going through imaginary time int Nt = 50000; - auto gt = gf{ {gw.domain(),Nt, half_bins}, gw.data().shape().front_pop() }; + auto gt = gf{{gw.domain(), Nt, half_bins}, gw.data().shape().front_pop()}; // We first transform to imaginary time because it's been coded with the knowledge of the tails - gt() = lazy_inverse_fourier(gw); + gt() = inverse_fourier(gw); legendre_matsubara_inverse(gl, gt()); - -} + } -void legendre_matsubara_direct (gf_view & gt, gf_const_view gl) { - - gt() = 0.0; - legendre_generator L; - - for (auto t : gt.mesh()) { - L.reset( 2*t/gt.domain().beta-1 ); - for (auto l : gl.mesh()) { - gt[t] += sqrt(2*l.index()+1) / gt.domain().beta * gl[l] * L.next(); - } - } - - gt.singularity() = get_tail(gl, gt.singularity().size(), gt.singularity().order_min()); - -} - -void legendre_matsubara_inverse (gf_view & gl, gf_const_view gt) { - - gl() = 0.0; - legendre_generator L; - - // Do the integral over imaginary time - for (auto t : gt.mesh()) { - L.reset( 2*t/gt.domain().beta-1 ); - for (auto l : gl.mesh()) { - gl[l] += sqrt(2*l.index()+1) * L.next() * gt[t]; - } - } - gl.data() *= gt.mesh().delta(); - -} - - -gf_keeper lazy_legendre_imfreq (gf_const_view gl) { return {gl}; } -gf_keeper lazy_legendre_imtime (gf_const_view gl) { return {gl}; } -gf_keeper lazy_imfreq_legendre (gf_const_view gw) { return {gw}; } -gf_keeper lazy_imtime_legendre (gf_const_view gt) { return {gt}; } - -void triqs_gf_view_assign_delegation( gf_view &gw, gf_keeper const & L) { + void triqs_gf_view_assign_delegation(gf_view gw, gf_keeper const& L) { legendre_matsubara_direct(gw, L.g); -} -void triqs_gf_view_assign_delegation( gf_view >, gf_keeper const & L) { + } + void triqs_gf_view_assign_delegation(gf_view gt, gf_keeper const& L) { legendre_matsubara_direct(gt, L.g); -} -void triqs_gf_view_assign_delegation( gf_view &gl, gf_keeper const & L) { + } + void triqs_gf_view_assign_delegation(gf_view gl, gf_keeper const& L) { legendre_matsubara_inverse(gl, L.g); -} -void triqs_gf_view_assign_delegation( gf_view &gl, gf_keeper const & L) { + } + void triqs_gf_view_assign_delegation(gf_view gl, gf_keeper const& L) { legendre_matsubara_inverse(gl, L.g); + } +} } - -}} diff --git a/triqs/gfs/local/legendre_matsubara.hpp b/triqs/gfs/local/legendre_matsubara.hpp index 83a6de5f..38ffb95e 100644 --- a/triqs/gfs/local/legendre_matsubara.hpp +++ b/triqs/gfs/local/legendre_matsubara.hpp @@ -1,9 +1,8 @@ - /******************************************************************************* * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011 by L. Boehnke, M. Ferrero, O. Parcollet + * Copyright (C) 2011-2014 by L. Boehnke, 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 @@ -19,33 +18,37 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_GF_LOCAL_LEGENDRE_MATSU_H +#ifndef TRIQS_GF_LOCAL_LEGENDRE_MATSU_H #define TRIQS_GF_LOCAL_LEGENDRE_MATSU_H #include #include #include -namespace triqs { namespace gfs { +namespace triqs { +namespace gfs { - void legendre_matsubara_direct (gf_view &gw, gf_const_view gl); - void legendre_matsubara_inverse (gf_view &gl, gf_const_view gw); + namespace tags { + struct legendre {}; + } - void legendre_matsubara_direct (gf_view >, gf_const_view gl); - void legendre_matsubara_inverse (gf_view &gl, gf_const_view gt); + inline gf_keeper legendre_to_imfreq(gf_const_view gl) { + return {gl}; + } + inline gf_keeper legendre_to_imtime(gf_const_view gl) { + return {gl}; + } + inline gf_keeper imfreq_to_legendre(gf_const_view gw) { + return {gw}; + } + inline gf_keeper imtime_to_legendre(gf_const_view gt) { + return {gt}; + } - namespace tags { struct legendre{}; } - - gf_keeper lazy_legendre_imfreq (gf_const_view gl); - gf_keeper lazy_legendre_imtime (gf_const_view gl); - gf_keeper lazy_imfreq_legendre (gf_const_view gw); - gf_keeper lazy_imtime_legendre (gf_const_view gt); - - void triqs_gf_view_assign_delegation( gf_view &gw, gf_keeper const & L); - void triqs_gf_view_assign_delegation( gf_view >, gf_keeper const & L); - void triqs_gf_view_assign_delegation( gf_view &gl, gf_keeper const & L); - void triqs_gf_view_assign_delegation( gf_view &gl, gf_keeper const & L); - - -}} + void triqs_gf_view_assign_delegation(gf_view gw, gf_keeper const &L); + void triqs_gf_view_assign_delegation(gf_view gt, gf_keeper const &L); + void triqs_gf_view_assign_delegation(gf_view gl, gf_keeper const &L); + void triqs_gf_view_assign_delegation(gf_view gl, gf_keeper const &L); +} +} #endif diff --git a/triqs/gfs/local/no_tail.hpp b/triqs/gfs/local/no_tail.hpp new file mode 100644 index 00000000..b3ffa0e1 --- /dev/null +++ b/triqs/gfs/local/no_tail.hpp @@ -0,0 +1,48 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2013-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 "./tail.hpp" + +namespace triqs { +namespace gfs { + + struct no_tail {}; + + template + gf_view make_gf_view_without_tail(gf_impl const &g) { + return {g.mesh(), g.data(), {}, g.symmetry()}; + } + + template + gf_view make_gf_from_g_and_tail(gf_impl const &g, local::tail t) { + 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) { + return {g.mesh(), g.data(), t, g.symmetry()}; + } +} +} diff --git a/triqs/gfs/meshes/linear.hpp b/triqs/gfs/meshes/linear.hpp index 107d372e..c632d4da 100644 --- a/triqs/gfs/meshes/linear.hpp +++ b/triqs/gfs/meshes/linear.hpp @@ -43,7 +43,7 @@ namespace gfs { linear_mesh() : _dom(), L(0), a_pt(0), b_pt(0), xmin(0), xmax(0), del(0), meshk(half_bins) {} - linear_mesh(domain_t dom, double a, double b, size_t n_pts, mesh_kind mk) + explicit linear_mesh(domain_t dom, double a, double b, size_t n_pts, mesh_kind mk) : _dom(std::move(dom)), L(n_pts), a_pt(a), b_pt(b), meshk(mk) { switch (mk) { case half_bins: diff --git a/triqs/gfs/tools.hpp b/triqs/gfs/tools.hpp index 08e40be8..31175102 100644 --- a/triqs/gfs/tools.hpp +++ b/triqs/gfs/tools.hpp @@ -50,6 +50,8 @@ namespace triqs { namespace gfs { struct freq_infty{}; // the point at infinity + //------------------------------------------------------ + inline std::vector split(const std::string &s, char delim){ std::vector elems; std::stringstream ss(s);