work on gf: scalar_valued, ...

- introducing scalar_valued gf
- Change Fourier routines to run on scalar_valued,
and then use those routines to run on matrix_valued.
- Tools for slices of 2 variables functions

test/triqs/gf/fourier1.cpp

@@ -25,8 +25,8 @@ int main() {
  auto G3 = make_gf<imfreq> (beta, Fermion, make_shape(2,2));
  auto Gt = make_gf<imtime> (beta, Fermion, make_shape(2,2));
 
- //auto gt = inverse_fourier(G);
+ auto gt = inverse_fourier(G);
- //auto gw = fourier(gt);
+ auto gw = fourier(gt);
 
  //gw() = lazy_fourier(gt);
  G() = lazy_fourier(Gt);

test/triqs/gf/gf_2times_b.cpp

@@ -27,7 +27,7 @@ int main(){
 
  std::cout << rg (1,1)<< std::endl ;
  std::cout << R.on_mesh(1,1)<< std::endl ;
- std::cout << R(1,1)<< std::endl ;
+ std::cout << R(0.001,0.001)<< std::endl ;
 
  auto R2 = R;
 

test/triqs/gf/pb_affichage.cpp

@@ -0,0 +1,39 @@
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+
+#include <triqs/utility/first_include.hpp>
+#include <iostream>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <fstream>
+#include <triqs/gf/gf.hpp>
+#include <triqs/gf/two_real_times.hpp>
+#include <complex>
+using namespace triqs::gf;
+using namespace std;
+using triqs::arrays::make_shape;
+
+int main(){
+
+ double dt=0.001;
+ double tmax=0.005;
+ int nt=tmax/dt;
+ auto  R= make_gf<two_real_times> (tmax,nt,make_shape(1,1));//results
+
+ for(auto point:R.mesh()) R(point)=0;
+
+ const auto rg = on_mesh(R);
+ R.on_mesh(1,1) = 10;
+
+ std::cout << rg (1,1)<< std::endl ;
+ std::cout << R.on_mesh(1,1)<< std::endl ;
+ std::cout << R(0.001,0.001)<< std::endl ;
+
+ auto R2 = R;
+
+ //on_mesh(R2)(1,1) =  on_mesh(R)(1,1) * on_mesh(R)(1,1);
+ on_mesh(R2)(1,1)() =  on_mesh(R)(1,1) * on_mesh(R)(1,1);
+
+ std::cout  << on_mesh(R2)(1,1)<< std::endl; 
+ return 0;
+};

test/triqs/gf/test_fourier_matsubara.cpp

@@ -0,0 +1,64 @@
+//#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+
+#include <triqs/gf/imfreq.hpp> 
+#include <triqs/gf/imtime.hpp> 
+#include <triqs/gf/local/fourier_matsubara.hpp> 
+
+namespace tql= triqs::clef;
+// namespace tqa= triqs::arrays;
+// using tqa::range;
+using triqs::arrays::make_shape;
+using triqs::gf::Fermion;
+using triqs::gf::imfreq;
+using triqs::gf::imtime;
+using triqs::gf::make_gf;
+using triqs::arrays::range;
+#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl;
+
+int main() {
+ 
+ double precision=10e-9;
+ H5::H5File file("test_fourier_matsubara.h5",H5F_ACC_TRUNC);
+ triqs::clef::placeholder<0> om_;
+ double beta =1;
+ int N=10000;
+ double E=1;
+ 
+ auto Gw1 = make_gf<imfreq> (beta, Fermion, make_shape(1,1), N);
+ Gw1(om_) << 1/(om_-E);
+//  for(auto const& w:Gw1.mesh()){
+//   std::cout<<"w="<<std::complex<double>(w)<<", Gw1=" << Gw1(w)(0,0)<<std::endl;
+//  }
+ h5_write(file, "Gw1", Gw1);   // the original lorentzian
+ 
+ auto Gt1 = make_gf<imtime> (beta, Fermion, make_shape(1,1), N);
+ inverse_fourier_impl( Gt1, Gw1, triqs::gf::matrix_valued() );
+//  for(auto const& t:Gt1.mesh()){
+//   std::cout<<"t="<<t<<",  expected="<<exp(-E*t) * ( (t>0?-1:0)+1/(1+exp(E*beta)) )<<std::endl;
+//  }
+ h5_write(file, "Gt1", Gt1);   // the lorentzian TF : lorentzian_inverse
+ 
+ ///verification that TF(TF^-1)=Id
+ auto Gw1b = make_gf<imfreq> (beta, Fermion, make_shape(1,1), N);
+ fourier_impl(Gw1b, Gt1, triqs::gf::matrix_valued());
+ for(auto const& w:Gw1.mesh()){
+//   std::cout<<"w="<<std::complex<double>(w)<<",Gw1b=" << Gw1b(w)(0,0)<<std::endl;
+//   std::cout<<"w="<<std::complex<double>(w)<<",Delta Gw1b=" << Gw1b(w)(0,0)-Gw1(w)(0,0)<<std::endl;
+  if ( std::abs(Gw1b(w)(0,0)-Gw1(w)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : w="<<std::complex<double>(w)<<" ,Gw1b="<<std::abs(Gw1b(w)(0,0))<<"\n";
+ }
+ h5_write(file,"Gw1b",Gw1b); // must be 0
+ 
+ ///verification that the TF is OK
+ for(auto const & t:Gt1.mesh()){
+  Gt1(t)-= exp(-E*t) * ( (t>0?-1:0)+1/(1+exp(E*beta)) );
+  if ( std::abs(Gt1(t)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_matsubara error : t="<<t<<" ,G1="<<std::abs(Gt1(t)(0,0))<<"\n";
+ }
+ h5_write(file,"Gt1b",Gt1); // must be 0
+ 
+ ///to verify that lazy_fourier computes
+ auto Gw2 = make_gf<imfreq> (beta, Fermion, make_shape(1,1));
+ Gw2() = lazy_fourier(Gt1);
+ 
+}
+
+

test/triqs/gf/test_fourier_real_time.cpp

@@ -0,0 +1,100 @@
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+
+#include <triqs/gf/refreq.hpp> 
+#include <triqs/gf/retime.hpp> 
+#include <triqs/gf/local/fourier_real.hpp> 
+#include <triqs/arrays.hpp>
+
+using triqs::arrays::make_shape;
+using triqs::gf::refreq;
+using triqs::gf::retime;
+using triqs::gf::make_gf;
+
+double lorentzian(double w, double a){
+  return 2*a / (w*w + a*a) ; 
+};
+std::complex<double> lorentzian_inverse(double t, double a){
+  return std::exp(-a*std::abs(t)) ; 
+};
+double theta(double x){
+  return x>0 ? 1.0 : ( x<0 ? 0.0 : 0.5 ) ; 
+};
+
+int main() {
+ 
+
+  double precision=10e-10;
+  H5::H5File file("fourier_real_time.h5",H5F_ACC_TRUNC);
+
+  std::complex<double> I(0,1);  
+  
+  //Test on the tail: GF in frequency that is a lorentzian, with its singularity, TF and TF^-1. 
+  
+  double wmax=10;
+  int Nw=1001;
+  
+  auto Gw1 = make_gf<refreq> (-wmax, wmax, Nw, make_shape(1,1),triqs::gf::full_bins);
+  double a = Gw1.mesh().delta() * sqrt( Gw1.mesh().size() );
+  for(auto const & w:Gw1.mesh()) Gw1(w)=lorentzian(w,a);
+  Gw1.singularity()(2)=triqs::arrays::matrix<double>{{2.0*a}};
+  h5_write(file,"Gw1",Gw1);   // the original lorentzian
+  
+  auto Gt1 = inverse_fourier(Gw1);  
+  h5_write(file,"Gt1",Gt1);   // the lorentzian TF : lorentzian_inverse
+  
+  // verification that TF(TF^-1)=Id
+  auto Gw1b = 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="<<w<<" ,G1="<<std::abs(Gw1b(w)(0,0))<<"\n";
+  }
+  h5_write(file,"Gw1b",Gw1b); // must be 0
+  
+  // verification that TF is the lorentzian_inverse function
+  for(auto const & t:Gt1.mesh()){
+    Gt1(t)-=lorentzian_inverse(t,a);
+    if ( std::abs(Gt1(t)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : t="<<t<<" ,G1="<<std::abs(Gt1(t)(0,0))<<"\n";
+  }
+  h5_write(file,"Gt1b",Gt1); // must be 0
+  
+  
+  
+  //Test on the tail: GF in time that is a decreasing exponential 
+  
+  double tmax=10.;
+  int Nt=501;
+  
+  auto Gt2 = make_gf<retime> (-tmax, tmax, Nt, make_shape(1,1));
+  a = 2*acos(-1.) / ( Gt2.mesh().delta() *sqrt( Gt2.mesh().size() ) );
+  for(auto const & t:Gt2.mesh()) Gt2(t) = 0.5 *I * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
+  //for(auto const & t:Gt2.mesh()) Gt2(t) = 0.5_j * ( lorentzian_inverse(-t,a)*theta(-t)-lorentzian_inverse(t,a)*theta(t) );
+  Gt2.singularity()(1)=triqs::arrays::matrix<double>{{1.0}};
+  h5_write(file,"Gt2",Gt2);
+  
+  auto Gw2 = fourier(Gt2);
+  h5_write(file,"Gw2",Gw2);
+  
+  for(auto const & w:Gw2.mesh()){
+    Gw2(w)-= 0.5/(w+a*I)+0.5/(w-a*I);
+    //Gw2(w)-= 0.5/(w+a*1_j)+0.5/(w-a*1_j);
+    if ( std::abs(Gw2(w)(0,0)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_real_time error : w="<<w<<" ,G2="<<std::abs(Gw2(w)(0,0))<<"\n";
+  }
+  h5_write(file,"Gw2b",Gw2);
+
+  
+  
+  //Test : GF in time is a simple trigonometric function, the result is a sum of Dirac functions 
+  
+  tmax=4*acos(-1.);
+  
+  auto Gt3 = make_gf<retime> (-tmax, tmax, Nt, make_shape(1,1));
+  for(auto const & t:Gt3.mesh()) Gt3(t) = 1.0 * std::cos(10*t) + 0.25*std::sin(4*t) + 0.5 * I*std::sin(8*t+0.3*acos(-1.)) ;
+  //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);
+  h5_write(file,"Gw3",Gw3);
+  
+}
+
+

test/triqs/gf/test_gf_triqs.cpp

@@ -35,7 +35,6 @@ void print_to_file(std::string const s, gf<imtime> const & gt){
 
 void test_0(){
  
- triqs::gf::freq_infty inf;
  int Ntau = 10001;
  double beta =1;
 

triqs/gf/gf.hpp

@@ -385,20 +385,33 @@ namespace triqs { namespace gf {
   }
 
  // tool for lazy transformation
- template<typename Tag, typename D> struct gf_keeper{ gf_view<D> g; gf_keeper (gf_view<D> const & g_) : g(g_) {} };
+ template<typename Tag, typename D, typename Target = matrix_valued> struct gf_keeper{ gf_view<D,Target> g; gf_keeper (gf_view<D,Target> const & g_) : g(g_) {} };
 
  // ---------------------------------- slicing ------------------------------------
  
+ //slice
  template<typename Variable, typename Target, typename Opt, bool V, typename... Args>
-  gf_view<Variable,Target,Opt> slice_target (gf_impl<Variable,Target,Opt,V> const & g, Args... args) {
+  gf_view<Variable,matrix_valued,Opt> slice_target (gf_impl<Variable,Target,Opt,V> const & g, Args... args) {
-   return gf_view<Variable,Target,Opt>(g.mesh(), g.data()(tqa::range(), args... ), slice_target (g.singularity(),args...), g.symmetry());
+   static_assert(std::is_same<Target,matrix_valued>::value, "slice_target only for matrix_valued GF's");
+   using arrays::range;
+   //auto sg=slice_target (g.singularity(),range(args,args+1)...);
+   return gf_view<Variable,matrix_valued,Opt>(g.mesh(), g.data()(tqa::range(), args... ),  slice_target (g.singularity(),args...) , g.symmetry());
   }
 
  template<typename Variable, typename Target, typename Opt, bool V, typename... Args>
-  gf_view<Variable,Target,Opt> slice_mesh (gf_impl<Variable,Target,Opt,V> const & g, Args... args) {
+  gf_view<Variable,scalar_valued,Opt> slice_target_to_scalar (gf_impl<Variable,Target,Opt,V> const & g, Args... args) {
-   return gf_view<Variable,Target,Opt>(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), g.singularity(), g.symmetry());
+   static_assert(std::is_same<Target,matrix_valued>::value, "slice_target only for matrix_valued GF's");
+   using arrays::range;
+   auto sg=slice_target (g.singularity(),range(args,args+1)...);
+   return gf_view<Variable,scalar_valued,Opt>(g.mesh(), g.data()(tqa::range(), args... ), sg , g.symmetry());
   }
 
+/*
+  template<typename Variable1,typename Variable2, typename Target, typename Opt, bool V, typename... Args>
+  gf_view<Variable2,Target,Opt> slice_mesh (gf_impl<Variable1,Target,Opt,V> const & g, Args... args) {
+   return gf_view<Variable2,Target,Opt>(g.mesh().slice(args...), g.data()(g.mesh().slice_get_range(args...),arrays::ellipsis()), g.singularity(), g.symmetry());
+  }*/
+
 }}
 
 #include "./gf_expr.hpp"

triqs/gf/imfreq.hpp

@@ -44,29 +44,33 @@ namespace triqs { namespace gf {
 
   //singularity
   template<typename Opt> struct singularity<imfreq,matrix_valued,Opt>  { typedef local::tail type;};
+  template<typename Opt> struct singularity<imfreq,scalar_valued,Opt>  { typedef local::tail type;};
   
   //h5 name
   template<typename Opt> struct h5_name<imfreq,matrix_valued,Opt>      { static std::string invoke(){ return "GfImFreq";}};
+  template<typename Opt> struct h5_name<imfreq,scalar_valued,Opt>      { static std::string invoke(){ return "GfImFreq_s";}};
 
   /// ---------------------------  evaluator ---------------------------------
-
+  template<typename Opt, typename Target>
-  template<typename Opt>
+   struct evaluator<imfreq,Target,Opt> {
-   struct evaluator<imfreq,matrix_valued,Opt> {
     static constexpr int arity = 1;
+    typedef typename std::conditional < std::is_same<Target, matrix_valued>::value, arrays::matrix_view<std::complex<double>>, std::complex<double>>::type rtype; 
     template<typename G>
-     arrays::matrix_view<std::complex<double> >  operator() (G const * g, long n)  const {return g->data()(n, arrays::range(), arrays::range()); }
+     rtype operator() (G const * g, long n)  const {return g->data()(n, arrays::ellipsis()); }
     template<typename G>
      local::tail_view operator()(G const * g, freq_infty const &) const {return g->singularity();}
    };
 
   /// ---------------------------  data access  ---------------------------------
-
   template<typename Opt> struct data_proxy<imfreq,matrix_valued,Opt> : data_proxy_array<std::complex<double>,3> {};
+  template<typename Opt> struct data_proxy<imfreq,scalar_valued,Opt> : data_proxy_array<std::complex<double>,1> {};
 
   // -------------------------------   Factories  --------------------------------------------------
   
+  // matrix_valued
   template<typename Opt> struct factories<imfreq,matrix_valued,Opt> { 
     typedef gf<imfreq,matrix_valued,Opt> gf_t;
+    typedef gf_view<imfreq,matrix_valued,Opt> gf_view_t;
     
     template<typename MeshType>
     static gf_t make_gf(MeshType && m, tqa::mini_vector<size_t,2> shape, local::tail_view const & t) {
@@ -83,6 +87,28 @@ namespace triqs { namespace gf {
       return make_gf(mesh<imfreq,Opt>::make(beta,S,Nmax), shape, t);
     }
   };
+  
+  // scalar_valued
+  template<typename Opt> struct factories<imfreq,scalar_valued,Opt> { 
+   typedef gf<imfreq,scalar_valued,Opt> gf_t;
+   typedef gf_view<imfreq,scalar_valued,Opt> gf_view_t;
+
+   template<typename MeshType>
+   static gf_t make_gf(MeshType && m, local::tail_view const & t) {
+    typename gf_t::data_non_view_t A(m.size()); A() =0;
+    return gf_t ( std::forward<MeshType>(m), std::move(A), t, nothing() ) ;
+   }
+   static gf_t make_gf(double beta, statistic_enum S) {
+     return make_gf(mesh<imfreq,Opt>::make(beta,S), local::tail(tqa::mini_vector<size_t,2> (1,1)));
+   }
+   static gf_t make_gf(double beta, statistic_enum S, size_t Nmax) {
+    return make_gf(mesh<imfreq,Opt>::make(beta,S,Nmax), local::tail(tqa::mini_vector<size_t,2> (1,1)));
+   }
+   static gf_t make_gf(double beta, statistic_enum S, size_t Nmax, local::tail_view const & t) {
+    return make_gf(mesh<imfreq,Opt>::make(beta,S,Nmax), t);
+   }
+  };
  } // gf_implementation 
+
 }}
 #endif

triqs/gf/imtime.hpp

@@ -43,14 +43,21 @@ namespace triqs { namespace gf {
  
   // singularity 
   template<typename Opt> struct singularity<imtime,matrix_valued,Opt>  { typedef local::tail type;};
+  template<typename Opt> struct singularity<imtime,scalar_valued,Opt>  { typedef local::tail type;};
   
   // h5 name
   template<typename Opt> struct h5_name<imtime,matrix_valued,Opt>      { static std::string invoke(){ return  "GfImTime";}};
+  template<typename Opt> struct h5_name<imtime,scalar_valued,Opt>      { static std::string invoke(){ return  "GfImTime_s";}};
+
+  /// ---------------------------  data access  ---------------------------------
+
+  template<typename Opt> struct data_proxy<imtime,matrix_valued,Opt> : data_proxy_array<double,3> {};
+  template<typename Opt> struct data_proxy<imtime,scalar_valued,Opt> : data_proxy_array<double,1> {};
 
    /// ---------------------------  closest mesh point on the grid ---------------------------------
 
-  template<typename Opt>
+  template<typename Opt, typename Target>
-   struct get_closest_point <imtime,matrix_valued,Opt> {
+   struct get_closest_point <imtime,Target,Opt> {
     // index_t is size_t
     template<typename G, typename T>
      static size_t invoke(G const * g, closest_pt_wrap<T> const & p) {
@@ -63,21 +70,11 @@ namespace triqs { namespace gf {
 
   /// ---------------------------  evaluator ---------------------------------
 
-  template<typename Opt>
-   struct evaluator<imtime,matrix_valued,Opt> {
-    private:
-     mutable arrays::matrix<double> _tmp;
-    public :
-     static constexpr int arity = 1;
-     evaluator() = default;
-     evaluator(size_t n1, size_t n2) : _tmp(n1,n2) {}
-     // WHAT happen in resize ??
-
   // NOT TESTED
   // TEST THE SPPED when q_view are incorporated...
   // true evaluator with interpolation ...
-     template<typename G>
+  template<typename G, typename ReturnType>
-      arrays::matrix<double> const & operator()(G const * g, double tau) const {
+   ReturnType evaluator_imtime_impl (G const * g, double tau, ReturnType && _tmp) {
     // interpolate between n and n+1, with weight
     double beta = g->mesh().domain().beta;
     int p = std::floor(tau/beta);
@@ -86,10 +83,11 @@ namespace triqs { namespace gf {
     long n = std::floor(a);
     double w = a-n;
     assert(n < g->mesh().size()-1);
+    auto _ = arrays::ellipsis();
     if ((g->mesh().domain().statistic == Fermion) && (p%2==1))
-	_tmp = - w*g->data()(n, arrays::range(), arrays::range()) - (1-w)*g->data()(n+1, arrays::range(), arrays::range());
+     _tmp = - w*g->data()(n, _) - (1-w)*g->data()(n+1, _);
     else
-	_tmp =   w*g->data()(n, arrays::range(), arrays::range()) + (1-w)*g->data()(n+1, arrays::range(), arrays::range());
+     _tmp =   w*g->data()(n, _) + (1-w)*g->data()(n+1, _);
     //else { // Speed test to redo when incoparated qview in main branch
     // _tmp(0,0) =   w*g->data()(n, 0,0) + (1-w)*g->data()(n+1, 0,0);
     // _tmp(0,1) =   w*g->data()(n, 0,1) + (1-w)*g->data()(n+1, 0,1);
@@ -99,32 +97,44 @@ namespace triqs { namespace gf {
     return _tmp;
    }
 
+  template<typename Opt>
+   struct evaluator<imtime,matrix_valued,Opt> {
+    private:
+     mutable arrays::matrix<double> _tmp;
+    public :
+     static constexpr int arity = 1;
+     evaluator() = default;
+     evaluator(size_t n1, size_t n2) : _tmp(n1,n2) {} // WHAT happen in resize ??
+
+     template<typename G>
+      arrays::matrix<double> const & operator()(G const * g, double tau) const { return evaluator_imtime_impl(g, tau, _tmp);}
+
      template<typename G>
       typename G::singularity_t const & operator()(G const * g,freq_infty const &) const {return g->singularity();}
    };
 
-  /// ---------------------------  data access  ---------------------------------
+  template<typename Opt>
+   struct evaluator<imtime,scalar_valued,Opt> {
+    public :
+     static constexpr int arity = 1;
 
-  template<typename Opt> struct data_proxy<imtime,matrix_valued,Opt> : data_proxy_array<double,3> {};
+     template<typename G> double operator()(G const * g, double tau) const { return evaluator_imtime_impl(g, tau, 0.0);}
+
+     template<typename G>
+      typename G::singularity_t const & operator()(G const * g,freq_infty const &) const {return g->singularity();}
+   };
 
  // -------------------------------   Factories  --------------------------------------------------
 
+  // matrix_valued
   template<typename Opt> struct factories<imtime,matrix_valued,Opt> { 
    typedef gf<imtime,matrix_valued,Opt> gf_t;
-
    template<typename MeshType>
     static gf_t make_gf(MeshType && m, tqa::mini_vector<size_t,2> shape, local::tail_view const & t) {
      typename gf_t::data_non_view_t A(shape.front_append(m.size())); A() =0;
      //return gf_t ( m, std::move(A), t, nothing() ) ;
      return gf_t (std::forward<MeshType>(m), std::move(A), t, nothing(), evaluator<imtime,matrix_valued,Opt>(shape[0],shape[1]) ) ;
     }
-   /*static gf_t make_gf(double beta, statistic_enum S, tqa::mini_vector<size_t,2> shape) {
-     return make_gf(make_mesh(beta,S,1025,half_bins), shape, local::tail(shape));
-     }
-     static gf_t make_gf(double beta, statistic_enum S, tqa::mini_vector<size_t,2> shape, size_t Nmax) {
-     return make_gf(make_mesh(beta,S,Nmax,half_bins), shape, local::tail(shape));
-     }
-     */
    static gf_t make_gf(double beta, statistic_enum S,  tqa::mini_vector<size_t,2> shape, size_t Nmax=1025, mesh_kind mk= half_bins) {
     return make_gf(mesh<imtime,Opt>::make(beta,S,Nmax,mk), shape, local::tail(shape));
    }
@@ -132,7 +142,24 @@ namespace triqs { namespace gf {
     return make_gf(mesh<imtime,Opt>::make(beta,S,Nmax,mk), shape, t);
    }
   };
- } // gf_implementation
+
+  // scalar_valued 
+  template<typename Opt> struct factories<imtime,scalar_valued,Opt> { 
+   typedef gf<imtime,scalar_valued,Opt> gf_t;
+   template<typename MeshType>
+    static gf_t make_gf(MeshType && m, local::tail_view const & t) {
+     typename gf_t::data_non_view_t A(m.size()); A() =0;
+     return gf_t (std::forward<MeshType>(m), std::move(A), t, nothing());
+    }
+   static gf_t make_gf(double beta, statistic_enum S, size_t Nmax=1025, mesh_kind mk= half_bins) {
+    return make_gf(mesh<imtime,Opt>::make(beta,S,Nmax,mk), local::tail(tqa::mini_vector<size_t,2> (1,1)));
+   }
+   static gf_t make_gf(double beta, statistic_enum S, size_t Nmax, mesh_kind mk, local::tail_view const & t) {
+    return make_gf(mesh<imtime,Opt>::make(beta,S,Nmax,mk), t);
+   }
+  };
+ } // gf_implementation.
+ 
 }}
 #endif
 

triqs/gf/local/fourier_base.cpp

@@ -40,13 +40,11 @@ namespace triqs { namespace gf { namespace details {
   const size_t imax = std::min(L,in.size());
   for (size_t i =0; i<imax; ++i) { inFFT[i][0] = real(in_ptr[i]); inFFT[i][1] = imag(in_ptr[i]);}
   for (size_t i =imax; i<L; ++i) { inFFT[i][0] = 0; inFFT[i][1] = 0;}
-
+//   for (size_t i =0; i<L; ++i) { inFFT[i][0] = real(in_ptr[i]); inFFT[i][1] = imag(in_ptr[i]);}
   p = fftw_plan_dft_1d(L, inFFT, outFFT, (direct ? FFTW_BACKWARD : FFTW_FORWARD), FFTW_ESTIMATE);
   fftw_execute(p); 
-
   const size_t jmax = std::min(L,out.size());
-  for (size_t j =0; j<jmax; ++j)  {out_ptr[j] = dcomplex(outFFT[j][0] ,outFFT[j][1]);}
+  for (size_t j =0; j<jmax; ++j) out_ptr[j] = dcomplex(outFFT[j][0] ,outFFT[j][1]);
-
   fftw_destroy_plan(p); 
   fftw_free(inFFT); fftw_free(outFFT);
  }

triqs/gf/local/fourier_base.hpp

@@ -30,6 +30,7 @@ namespace triqs { namespace gf {
    typedef std::complex<double> dcomplex;
    
    void fourier_base(const tqa::vector<dcomplex> &in, tqa::vector<dcomplex> &out, size_t L, bool direct);
+   void fourier_base(const tqa::vector<dcomplex> &in, tqa::vector<dcomplex> &out, size_t L1, size_t L2, bool direct);
 
  }
 

triqs/gf/local/fourier_matsubara.cpp

@@ -26,148 +26,162 @@ namespace triqs { namespace gf {
  
  namespace impl_local_matsubara {
   
-  inline dcomplex oneFermion(dcomplex a,double b,double tau,double Beta) {
+  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)) );
+   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) {
+  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)) );
+   return a*( b >=0 ? exp(-b*tau)/(exp(-beta*b)-1) : exp(b*(beta-tau))/(1-exp(b*beta)) );
   }
  }
  
- //--------------------------------------------------------------------------------------
-
- void fourier_impl  (gf_view<imfreq> &gw , gf_view<imtime> const & gt){
-
-  // set behavior according to mesh kind
-  double shift;
-  size_t L;
-  switch(gt.mesh().kind()) {
-    case half_bins: shift = 0.5; L = gt.mesh().size(); break;
-    case full_bins: shift = 0.0; L = gt.mesh().size()-1; break;
-    case without_last: shift = 0.0; L = gt.mesh().size(); break;
-  }
-
-  auto ta = gt(freq_infty());
-  long numberTimeSlices = gt.mesh().size();
-  double Beta = gt.domain().beta, Pi = std::acos(-1);
-  dcomplex I(0,1);
-  tqa::vector<dcomplex> g_in(gt.mesh().size()), g_out (gw.mesh().size()); 
-
-  using namespace impl_local_matsubara;
-  for (size_t n1=0; n1<gw.data().shape()[1];n1++)
-   for (size_t n2=0; n2<gw.data().shape()[2];n2++) {
-    dcomplex d= ta(1)(n1,n2), A= ta.get_or_zero(2)(n1,n2),B = ta.get_or_zero(3)(n1,n2);
-    //dcomplex d= ta(1)(n1,n2), A= ta(2)(n1,n2),B = ta(3)(n1,n2);
-    double b1, b2, b3;
-    dcomplex a1, a2, a3;
-    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;
-
-    if (gw.domain().statistic == Fermion){
-     for (auto & t : gt.mesh())  
-      g_in(t.index()) = exp(I*Pi*t/Beta)*( gt(t)(n1,n2) - ( 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()) =  gt(t)(n1,n2) -  ( oneBoson(a1,b1,t,Beta) + oneBoson(a2,b2,t,Beta) + oneBoson(a3,b3,t,Beta) );  
-    }
-
-    g_in *= Beta/numberTimeSlices;
-
-    details::fourier_base(g_in, g_out, L, true);
-
-    for (auto & w : gw.mesh()) {
-     gw(w)(n1,n2) = g_out(w.index())*exp(2*I*w.index()*shift*Pi/Beta*gt.mesh().delta()) + a1/(w-b1) + a2/(w-b2) + a3/(w-b3); 
-    }
-
-    // set tail
-    gw.singularity() = gt.singularity();
-
-   }
- }
- //---------------------------------------------------------------------------
-
  template<typename GfElementType> GfElementType convert_green         ( dcomplex const& x) { return x;}
  template<>                       double        convert_green<double> ( dcomplex const& x) { return real(x);}
  
- //---------------------------------------------------------------------------
+ //--------------------------------------------------------------------------------------
  
- void inverse_fourier_impl (gf_view<imtime> &gt,  gf_view<imfreq> const & gw) { 
+ struct impl_worker {
   
-  // set behavior according to mesh kind
+  tqa::vector<dcomplex> g_in, g_out;
-  double shift;
+  
-  size_t L;
+  void direct ( gf_view<imfreq,scalar_valued> &gw , gf_view<imtime,scalar_valued> const& gt) {
-  switch(gt.mesh().kind()) {
+   using namespace impl_local_matsubara;
-    case half_bins: shift = 0.5; L = gt.mesh().size(); break;
+   auto ta = gt(freq_infty());
-    case full_bins: shift = 0.0; L = gt.mesh().size()-1; break;
+   //TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO
-    case without_last: shift = 0.0; L = gt.mesh().size(); break;
+   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);
+   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; 
+   }
+   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) ) );  
+   }
+   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); 
+   }
+   gw.singularity() = gt.singularity();// set tail
   }
   
-  static bool Green_Function_Are_Complex_in_time = false;
-  auto ta = gw(freq_infty());
   
-  double Beta = gt.domain().beta, Pi = std::acos(-1);
-  dcomplex I(0,1);
-  tqa::vector<dcomplex> g_in(gw.mesh().size()), g_out (gt.mesh().size());
   
+  void inverse(gf_view<imtime,scalar_valued> gt, gf_view<imfreq,scalar_valued> const gw){
    using namespace impl_local_matsubara;
-  for (size_t n1=0; n1<gt.data().shape()[1];n1++)
+   static bool Green_Function_Are_Complex_in_time = false;
-   for (size_t n2=0; n2<gt.data().shape()[2];n2++) {
+   // If the Green function are NOT complex, then one use the symmetry property
-    dcomplex d= ta(1)(n1,n2), A= ta.get_or_zero(2)(n1,n2),B = ta.get_or_zero(3)(n1,n2);
+   // fold the sum and get a factor 2
-    //dcomplex d= ta(1)(n1,n2), A= ta(2)(n1,n2),B = ta(3)(n1,n2);
+   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);
    double b1, b2, b3;
    dcomplex a1, a2, a3;
-    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;
    
+   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;
+   }
+   g_in() = 0;
    for (auto & w: gw.mesh()) {
-     g_in(w.index()) =  exp(-I*2*w.index()*shift*Pi/Beta*gt.mesh().delta()) * ( gw(w)(n1,n2) - (a1/(w-b1) + a2/(w-b2) + a3/(w-b3)) );
+    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;
    
    details::fourier_base(g_in, g_out, L, false);
    
-    // If the Green function are NOT complex, then one use the symmetry property
-    // fold the sum and get a factor 2
-    double fact = (Green_Function_Are_Complex_in_time ? 1 : 2);
-    g_out = g_out*(fact/Beta);
-
    // CORRECT FOR COMPLEX G(tau) !!!
    typedef double gt_result_type;
-    //typedef boost::mpl::if_<gt_result_type;
    //typedef typename gf<imtime>::mesh_type::gf_result_type gt_result_type;
-
    if (gw.domain().statistic == Fermion){
-     for (auto & t : gt.mesh()) 
+    for (auto & t : gt.mesh()){
-       gt(t)(n1,n2) = convert_green<gt_result_type> (g_out(t.index())*exp(-I*Pi*t/Beta)
+     gt(t) = convert_green<gt_result_type> ( 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) );
+     + oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) );
+    }
    }
    else {
     for (auto & t : gt.mesh()) 
-       gt(t)(n1,n2) = convert_green<gt_result_type> (g_out(t.index())
+     gt(t) = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() )
-                         + oneBoson(a1,b1,t,Beta) + oneBoson(a2,b2,t,Beta) + oneBoson(a3,b3,t,Beta) );
+     + 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)(n1,n2) = -gt.on_mesh(0)(n1,n2)-convert_green<gt_result_type>(ta(1)(n1,n2));
+    gt.on_mesh(L) = -gt.on_mesh(0)-convert_green<gt_result_type>(ta(1)(0,0));
-
    // set tail
     gt.singularity() = gw.singularity();
+  }
   
+ };
+ 
+ void fourier_impl  (gf_view<imfreq,scalar_valued> gw , gf_view<imtime,scalar_valued> const gt, scalar_valued){
+  impl_worker w;
+  w.direct(gw,gt);
+ }
+ 
+ void fourier_impl  (gf_view<imfreq,matrix_valued> gw , gf_view<imtime,matrix_valued> const gt, matrix_valued){
+  impl_worker w;
+  for (size_t n1=0; n1<gt.data().shape()[1];n1++)
+   for (size_t n2=0; n2<gt.data().shape()[2];n2++){
+    auto gw_sl=slice_target_to_scalar(gw,n1,n2);
+    auto gt_sl=slice_target_to_scalar(gt,n1,n2);
+    w.direct(gw_sl, gt_sl);
    }
  }
  
  
- gf_keeper<tags::fourier,imtime> lazy_fourier         (gf_view<imtime> const & g) { return g;}
- gf_keeper<tags::fourier,imfreq> lazy_inverse_fourier (gf_view<imfreq> const & g) { return g;}
  
- void triqs_gf_view_assign_delegation( gf_view<imfreq> &g, gf_keeper<tags::fourier,imtime> const & L) { fourier_impl (g,L.g);}
+ //---------------------------------------------------------------------------
- void triqs_gf_view_assign_delegation( gf_view<imtime> &g, gf_keeper<tags::fourier,imfreq> const & L) { inverse_fourier_impl(g,L.g);}
+ 
+ void inverse_fourier_impl  (gf_view<imtime,scalar_valued> gt , gf_view<imfreq,scalar_valued> const gw, scalar_valued){ 
+  impl_worker w;
+  w.inverse(gt,gw);
+ }
+ 
+ void inverse_fourier_impl  (gf_view<imtime,matrix_valued> gt , gf_view<imfreq,matrix_valued> const gw, matrix_valued){
+  impl_worker w;
+  for (size_t n1=0; n1<gw.data().shape()[1];n1++)
+   for (size_t n2=0; n2<gw.data().shape()[2];n2++){
+    auto gt_sl=slice_target_to_scalar(gt, n1, n2);
+    auto gw_sl=slice_target_to_scalar(gw, n1, n2);
+    w.inverse( gt_sl, gw_sl);
+   }
+ }
+ 
+ 
+ 
+ //---------------------------------------------------------------------------
+ void triqs_gf_view_assign_delegation( gf_view<imfreq,scalar_valued> g, gf_keeper<tags::fourier,imtime,scalar_valued> const & L) { fourier_impl ( g ,L.g, scalar_valued() ); }
+ void triqs_gf_view_assign_delegation( gf_view<imfreq,matrix_valued> g, gf_keeper<tags::fourier,imtime,matrix_valued> const & L) { fourier_impl ( g, L.g, matrix_valued() ); }
+ void triqs_gf_view_assign_delegation( gf_view<imtime,scalar_valued> g, gf_keeper<tags::fourier,imfreq,scalar_valued> const & L) { inverse_fourier_impl( g, L.g, scalar_valued() ); }
+ void triqs_gf_view_assign_delegation( gf_view<imtime,matrix_valued> g, gf_keeper<tags::fourier,imfreq,matrix_valued> const & L) { inverse_fourier_impl( g, L.g, matrix_valued() ); }
  
 }}
 

triqs/gf/local/fourier_matsubara.hpp

@@ -28,14 +28,53 @@
 namespace triqs { namespace gf {
  
  // First the implementation of the fourier transform
- void fourier_impl         (gf_view<imfreq> &gw , gf_view<imtime> const & gt);
+ void fourier_impl         (gf_view<imfreq,scalar_valued> gw , gf_view<imtime,scalar_valued> const gt, scalar_valued);
- void inverse_fourier_impl (gf_view<imtime> &gt,  gf_view<imfreq> const & gw);
+ void fourier_impl         (gf_view<imfreq,matrix_valued> gw , gf_view<imtime,matrix_valued> const gt, matrix_valued);
+ void inverse_fourier_impl (gf_view<imtime,scalar_valued> gt,  gf_view<imfreq,scalar_valued> const gw, scalar_valued);
+ void inverse_fourier_impl (gf_view<imtime,matrix_valued> gt,  gf_view<imfreq,matrix_valued> const gw, matrix_valued);
  
- gf_keeper<tags::fourier,imtime> lazy_fourier         (gf_view<imtime> const & g);
+ inline gf_view<imfreq,matrix_valued> fourier (gf_view<imtime,matrix_valued> const & gt) {
- gf_keeper<tags::fourier,imfreq> lazy_inverse_fourier (gf_view<imfreq> const & g);
+  size_t L = (gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size() );
+  auto gw = make_gf<imfreq,matrix_valued>(gt.domain().beta, gt.domain().statistic , gt.data().shape().front_pop(), L);
+  auto V = gw();
+  fourier_impl(V, gt, matrix_valued());
+  return gw;
+ }
+ inline gf_view<imfreq,scalar_valued> fourier (gf_view<imtime,scalar_valued> const & gt) {
+  size_t L = (gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size() );
+  auto gw = make_gf<imfreq,scalar_valued>(gt.domain().beta, gt.domain().statistic, L);
+  auto V = gw();
+  fourier_impl(V, gt, scalar_valued());
+  return gw;
+ }
+ 
+ inline gf_view<imtime, matrix_valued> inverse_fourier (gf_view<imfreq, matrix_valued> const & gw, mesh_kind mk = half_bins) {
+  double pi = std::acos(-1);
+  size_t L = (mk == full_bins ? gw.mesh().size()+1 : gw.mesh().size() );
+  auto gt = make_gf<imtime,matrix_valued>(gw.domain().beta, gw.domain().statistic, gw.data().shape().front_pop(), L);
+  auto V = gt();
+  inverse_fourier_impl(V, gw, matrix_valued());
+  return gt;
+ }
+ inline gf_view<imtime,scalar_valued> inverse_fourier (gf_view<imfreq,scalar_valued> const & gw, mesh_kind mk = half_bins) {
+  double pi = std::acos(-1);
+  size_t L = (mk == full_bins ? gw.mesh().size()+1 : gw.mesh().size() );
+  auto gt = make_gf<imtime,scalar_valued>(gw.domain().beta, gw.domain().statistic, L);
+  auto V = gt();
+  inverse_fourier_impl(V, gw,scalar_valued());
+  return gt;
+ }
+ 
+ inline gf_keeper<tags::fourier,imtime,scalar_valued> lazy_fourier         (gf_view<imtime,scalar_valued> const & g) { return g;}
+ inline gf_keeper<tags::fourier,imfreq,scalar_valued> lazy_inverse_fourier (gf_view<imfreq,scalar_valued> const & g) { return g;}
+ inline gf_keeper<tags::fourier,imtime,matrix_valued> lazy_fourier         (gf_view<imtime,matrix_valued> const & g) { return g;}
+ inline gf_keeper<tags::fourier,imfreq,matrix_valued> lazy_inverse_fourier (gf_view<imfreq,matrix_valued> const & g) { return g;}
+ 
+ void triqs_gf_view_assign_delegation( gf_view<imfreq,scalar_valued> g, gf_keeper<tags::fourier,imtime,scalar_valued> const & L);
+ void triqs_gf_view_assign_delegation( gf_view<imfreq,matrix_valued> g, gf_keeper<tags::fourier,imtime,matrix_valued> const & L);
+ void triqs_gf_view_assign_delegation( gf_view<imtime,scalar_valued> g, gf_keeper<tags::fourier,imfreq,scalar_valued> const & L);
+ void triqs_gf_view_assign_delegation( gf_view<imtime,matrix_valued> g, gf_keeper<tags::fourier,imfreq,matrix_valued> const & L);
  
- void triqs_gf_view_assign_delegation( gf_view<imfreq> &g, gf_keeper<tags::fourier,imtime> const & L);
- void triqs_gf_view_assign_delegation( gf_view<imtime> &g, gf_keeper<tags::fourier,imfreq> const & L);
  
 }}
 #endif

triqs/gf/local/fourier_real.cpp

@@ -24,24 +24,28 @@
 
 namespace triqs { namespace gf { 
 
-  namespace impl_local_real {
+  namespace {
-    inline dcomplex th_expo(double t, double a ) { return (t < 0 ? 0 : -I * exp(-a*t)); }
+    double pi = std::acos(-1);
-    inline dcomplex th_expo_neg(double t, double a ) { return (t > 0 ? 0 : I * exp(a*t)); }
+    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) ; }
   }
 
   //--------------------------------------------------------------------------------------
 
-  void fourier_impl(gf_view<refreq> & gw, gf_view<retime> const & gt) {
+  void fourier_impl(gf_view<refreq> gw, gf_view<retime> const gt) {
-
-    using namespace impl_local_real;
 
     size_t L = gt.mesh().size();
-    if (gw.mesh().size() != gt.mesh().size()) TRIQS_RUNTIME_ERROR << "Meshes are different";
+    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() + (gt.mesh().kind() == half_bins ? 0.5 : 0.0) * gt.mesh().delta();
+    const double tmin = gt.mesh().x_min();
-    const double wmin = gw.mesh().x_min() + (gw.mesh().kind() == half_bins ? 0.5 : 0.0) * gw.mesh().delta();
+    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<dcomplex> g_in(L), g_out(L);
@@ -50,19 +54,19 @@ namespace triqs { namespace gf {
       for (size_t n2=0; n2<gw.data().shape()[2];n2++) {
 
         dcomplex t1 = ta(1)(n1,n2), t2= ta.get_or_zero(2)(n1,n2);
-        dcomplex a1 = (t1 - I * t2)/2, a2 = (t1 + I * t2)/2;
+        dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.;
 
         g_in() = 0;
 
-        for (auto & t : gt.mesh()) {
+        for (auto const & t : gt.mesh()) {
-          g_in(t.index()) = (gt(t)(n1,n2) - (a1*th_expo(t,1) + a2*th_expo_neg(t,1))) * std::exp(I*t*wmin);
+          g_in(t.index()) = (gt(t)(n1,n2) - (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 & w : gw.mesh()) {
+        for (auto const & w : gw.mesh()) {
-          gw(w)(n1,n2) = gt.mesh().delta() * std::exp(I*w*tmin) * std::exp(-I*wmin*tmin) * g_out(w.index())
+          gw(w)(n1,n2) = gt.mesh().delta() * std::exp(I*(w-wmin)*tmin) * g_out(w.index())
-                         + (a1/(w+I) + a2/(w-I));
+                         + a1*th_expo_inv(w,a) + a2*th_expo_neg_inv(w,a);
         }
       }
     }
@@ -74,17 +78,17 @@ namespace triqs { namespace gf {
 
   //---------------------------------------------------------------------------
 
-  void inverse_fourier_impl (gf_view<retime> & gt, gf_view<refreq> const & gw) { 
+  void inverse_fourier_impl (gf_view<retime> gt, gf_view<refreq> const gw) { 
-
-    using namespace impl_local_real;
 
     size_t L = gw.mesh().size();
-    if (gw.mesh().size() != gt.mesh().size()) TRIQS_RUNTIME_ERROR << "Meshes are different";
+    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() + (gt.mesh().kind() == half_bins ? 0.5 : 0.0) * gt.mesh().delta();
+    const double tmin = gt.mesh().x_min();
-    const double wmin = gw.mesh().x_min() + (gw.mesh().kind() == half_bins ? 0.5 : 0.0) * gw.mesh().delta();
+    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<dcomplex> g_in(L), g_out(L);
@@ -93,20 +97,18 @@ namespace triqs { namespace gf {
       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)/2, a2 = (t1 + I * t2)/2;
+        dcomplex a1 = (t1 + I * t2/a )/2., a2 = (t1 - I * t2/a )/2.;
         g_in() = 0;
 
-        for (auto & w: gw.mesh()) {
+        for (auto const & w: gw.mesh())
-          g_in(w.index()) = (gw(w)(n1,n2) - (a1/(w+I) + a2/(w-I))) * std::exp(-I*w*tmin);
+          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);
 
         const double corr = 1.0/(gt.mesh().delta()*L);
-        for (auto & t : gt.mesh()) {
+        for (auto const & t : gt.mesh())
-          gt(t)(n1,n2) = corr * std::exp(I*wmin*tmin) * std::exp(-I*wmin*t) *
+          gt(t)(n1,n2) = corr * std::exp(I*wmin*(tmin-t)) *
-                                g_out(t.index()) + a1*th_expo(t,1) + a2*th_expo_neg(t,1);
+                                 g_out( t.index() ) + a1 * th_expo(t,a) + a2 * th_expo_neg(t,a) ;
-        }
       }
     }
 
@@ -120,7 +122,7 @@ namespace triqs { namespace gf {
   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<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);}
+  void triqs_gf_view_assign_delegation( gf_view<retime> g, gf_keeper<tags::fourier,refreq> const & L) { inverse_fourier_impl(g,L.g);}
 
 }}

triqs/gf/local/fourier_real.hpp

@@ -27,29 +27,26 @@
 
 namespace triqs { namespace gf { 
 
- namespace impl_local_real {
-   dcomplex I(0,1);
-   double pi = std::acos(-1);
- }
-
  // 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> gw , gf_view<retime> const gt);
- void inverse_fourier_impl (gf_view<retime> &gt,  gf_view<refreq> const & gw);
+ void inverse_fourier_impl (gf_view<retime> gt,  gf_view<refreq> const gw);
 
- gf_view<refreq> fourier (gf_view<retime> const & gt) { 
+ inline gf_view<refreq> fourier (gf_view<retime> const & gt) { 
+   double pi = std::acos(-1);
    size_t L = gt.mesh().size();
-   double wmin = -impl_local_real::pi * (L-1) / (L*gt.mesh().delta());
+   double wmin = -pi * (L-1) / (L*gt.mesh().delta());
-   double wmax =  impl_local_real::pi * (L-1) / (L*gt.mesh().delta());
+   double wmax =  pi * (L-1) / (L*gt.mesh().delta());
    auto gw = make_gf<refreq>(wmin, wmax, L, gt.data().shape().front_pop());
    auto V = gw();
    fourier_impl(V, gt);
    return gw;
  }
 
- gf_view<retime> inverse_fourier (gf_view<refreq> const & gw) { 
+ inline gf_view<retime> inverse_fourier (gf_view<refreq> const & gw) { 
+   double pi = std::acos(-1);
    size_t L = gw.mesh().size();
-   double tmin = -impl_local_real::pi * (L-1) / (L*gw.mesh().delta());
+   double tmin = -pi * (L-1) / (L*gw.mesh().delta());
-   double tmax =  impl_local_real::pi * (L-1) / (L*gw.mesh().delta());
+   double tmax =  pi * (L-1) / (L*gw.mesh().delta());
    auto gt = make_gf<retime>(tmin, tmax, L, gw.data().shape().front_pop());
    auto V = gt();
    inverse_fourier_impl(V, gw);
@@ -59,8 +56,8 @@ namespace triqs { namespace gf {
  gf_keeper<tags::fourier,retime> lazy_fourier         (gf_view<retime> const & g);
  gf_keeper<tags::fourier,refreq> lazy_inverse_fourier (gf_view<refreq> const & 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);
 
 }}
 #endif

triqs/gf/meshes/product.hpp

@@ -23,6 +23,7 @@
 #include "./mesh_tools.hpp"
 #include "../domains/product.hpp"
 #include <triqs/utility/tuple_tools.hpp>
+#include <triqs/utility/mini_vector.hpp>
 namespace triqs { namespace gf {
 
  template<typename... Meshes> struct mesh_product : tag::composite {
@@ -57,6 +58,14 @@ namespace triqs { namespace gf {
    struct _aux3 { template<typename P, typename M> size_t operator()(M const & m, P const & p,size_t R) {return p.linear_index() + R * m.size();}};
    size_t mp_to_linear(m_pt_tuple_t const & mp) const { return triqs::tuple::fold_on_zip(_aux3(), m_tuple, mp, size_t(0)); }
 
+   //
+   struct _aux4 { template< typename M, typename V> V * operator()(M const & m, V * v) {*v = m.size(); return ++v;}};
+   utility::mini_vector<size_t,dim> all_size_as_mini_vector () const { 
+    utility::mini_vector<size_t,dim> res;
+    triqs::tuple::fold(_aux4(), m_tuple, &res[0] ); 
+    return res;
+   }
+
    // Same but a variadic list of mesh_point_t
    template<typename ... MP> size_t mesh_pt_components_to_linear(MP const &  ... mp) const {
     static_assert(std::is_same< std::tuple<MP...>, m_pt_tuple_t>::value, "Call incorrect ");
@@ -75,6 +84,7 @@ namespace triqs { namespace gf {
     mesh_point_t(mesh_product const & m_)                  : m(&m_), _c (triqs::tuple::apply(F1(), m_.m_tuple)), _atend(false)    {}
     m_pt_tuple_t const & components_tuple() const { return _c;}
     size_t linear_index() const { return m->mp_to_linear(_c);}
+    const mesh_product * mesh() const { return m;}
 
     typedef domain_pt_t cast_t;
     operator cast_t() const { return m->index_to_point(index);}
@@ -146,5 +156,35 @@ namespace triqs { namespace gf {
    domain_t _dom;
  };
 
+ //template<int pos, typename ... M> 
+ //typename std::tuple_element<pos,typename mesh_product<M...>::index_t>::type get_index1(typename mesh_product<M...>::mesh_point_t const & p) { return std::get<pos>(p.components_tuple());}
+
+ template<int pos, typename P> 
+  auto get_index(P const & p) DECL_AND_RETURN( std::get<pos>(p.components_tuple()).index());
+
+ template<int pos, typename P> 
+  auto get_point(P const & p) DECL_AND_RETURN( std::get<pos>( p.mesh()->components() ).index_to_point( std::get<pos>(p.components_tuple()).index() ) );
+ 
+ // C++14
+ //auto get_point(P const & p) { return std::get<pos> (p.mesh()->components()).index_to_point( std::get<pos>(p.components_tuple()));}
+
+ // Given a composite mesh m , and a linear array of storage A
+ // reinterpret_linear_array(m,A) returns a d-dimensionnal view of the array
+ // with indices egal to the indices of the components of the mesh.
+ // Very useful for slicing, currying functions.
+ template<typename ... Meshes, typename T, ull_t OptionsFlags >
+  arrays::array_view<T, sizeof...(Meshes),OptionsFlags, arrays::indexmaps::mem_layout::fortran_order(sizeof...(Meshes)) >
+  reinterpret_linear_array(mesh_product<Meshes...> const & m, arrays::array_view<T,1,OptionsFlags> const & A) {
+   static int constexpr rank = sizeof...(Meshes);
+   typedef arrays::array_view<T, sizeof...(Meshes),OptionsFlags, arrays::indexmaps::mem_layout::fortran_order(rank)> return_t;
+   typedef typename return_t::indexmap_type im_t;
+   auto l = m.all_size_as_mini_vector();
+   typename im_t::strides_type sv;
+   std::ptrdiff_t s= 1;
+   for (int u=0; u<rank; ++u) { sv[u] = s; s *= l[u];} // fortran type folding 
+   return return_t (im_t (l,sv,0) , A.storage());
+  }
+
+
 }}
 #endif

triqs/gf/refreq.hpp

@@ -35,41 +35,46 @@ namespace triqs { namespace gf {
   template<typename Opt> struct mesh<refreq,Opt>                { 
    typedef linear_mesh<R_domain> type;
    typedef typename type::domain_t domain_t;
-   static type make(double wmin, double wmax, size_t n_freq, mesh_kind mk) {
+   static type make(double wmin, double wmax, size_t n_freq, mesh_kind mk=full_bins) {
     return type(domain_t(), wmin, wmax, n_freq, mk);
    }
   };
   
+  // singularity 
   template<typename Opt> struct singularity<refreq,matrix_valued,Opt>  { typedef local::tail type;};
+  template<typename Opt> struct singularity<refreq,scalar_valued,Opt>  { typedef local::tail type;};
+  
+  // h5 name
   template<typename Opt> struct h5_name<refreq,matrix_valued,Opt>      { static std::string invoke(){ return "GfReFreq";}};
+  template<typename Opt> struct h5_name<refreq,scalar_valued,Opt>      { static std::string invoke(){ return "GfReFreq_s";}};
 
   /// ---------------------------  evaluator ---------------------------------
-
+  template<typename Opt, typename Target>
-  template<typename Opt>
+   struct evaluator<refreq,Target,Opt> {
-   struct evaluator<refreq,matrix_valued,Opt> {
     static constexpr int arity = 1;
+    typedef typename std::conditional < std::is_same<Target, matrix_valued>::value, arrays::matrix_view<std::complex<double>>, std::complex<double>>::type rtype; 
     template<typename G>
-     arrays::matrix_view<std::complex<double> >  operator() (G const * g,double w0)  const {
+     rtype operator() (G const * g,double w0)  const {
       auto & data = g->data();
       auto & mesh = g->mesh();
       size_t index; double w; bool in;
       std::tie(in, index, w) = windowing(mesh,w0);
       if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
-      arrays::matrix<std::complex<double> > res = w*data(mesh.index_to_linear(index), arrays::ellipsis()) + (1-w)*data(mesh.index_to_linear(index+1), arrays::ellipsis());
+      return w*data(mesh.index_to_linear(index), arrays::ellipsis()) + (1-w)*data(mesh.index_to_linear(index+1), arrays::ellipsis());
-      return res;
      }
     template<typename G>
      local::tail_view operator()(G const * g,freq_infty const &) const {return g->singularity();}
    };
 
   /// ---------------------------  data access  ---------------------------------
-
   template<typename Opt> struct data_proxy<refreq,matrix_valued,Opt> : data_proxy_array<std::complex<double>,3> {};
+  template<typename Opt> struct data_proxy<refreq,scalar_valued,Opt> : data_proxy_array<std::complex<double>,1> {};
    
   // -------------------------------   Factories  --------------------------------------------------
   
+  //matrix_valued
   template<typename Opt> struct factories<refreq, matrix_valued,Opt> {
-   typedef gf<refreq> gf_t;
+   typedef gf<refreq,matrix_valued> gf_t;
 
    template<typename MeshType>
     static gf_t make_gf(MeshType && m, tqa::mini_vector<size_t,2> shape, local::tail_view const & t) {
@@ -86,6 +91,27 @@ namespace triqs { namespace gf {
     typename gf_t::data_non_view_t A(shape.front_append(n_freq)); A() =0;
     return gf_t(mesh<refreq,Opt>::make(wmin, wmax, n_freq, mk), std::move(A), local::tail(shape), nothing());
    }
+  };
+   
+   //scalar_valued
+  template<typename Opt> struct factories<refreq,scalar_valued,Opt> {
+   typedef gf<refreq,scalar_valued> gf_t;
+
+   template<typename MeshType>
+    static gf_t make_gf(MeshType && m, local::tail_view const & t) {
+     typename gf_t::data_non_view_t A(m.size()); A() =0;
+     return gf_t ( std::forward<MeshType>(m), std::move(A), t, nothing() ) ;
+    }
+    
+    static gf_t make_gf(double wmin, double wmax, size_t n_freq) {
+     typename gf_t::data_non_view_t A(n_freq); A() =0;
+     return gf_t(mesh<refreq,Opt>::make(wmin, wmax, n_freq), std::move(A), local::tail(tqa::mini_vector<size_t,2>(1,1)), nothing());
+    }
+
+   static gf_t make_gf(double wmin, double wmax, size_t n_freq, mesh_kind mk) {
+    typename gf_t::data_non_view_t A(n_freq); A() =0;
+    return gf_t(mesh<refreq,Opt>::make(wmin, wmax, n_freq, mk), std::move(A), local::tail(tqa::mini_vector<size_t,2>(1,1)), nothing());
+   }
 
   };
  } // gf_implementation

triqs/gf/retime.hpp

@@ -36,53 +36,80 @@ namespace triqs { namespace gf {
    typedef linear_mesh<R_domain> type;
    typedef typename type::domain_t domain_t;
 
-   static type make(double tmin, double tmax, size_t n_points, mesh_kind mk) {
+   static type make(double tmin, double tmax, size_t n_points, mesh_kind mk=full_bins) {
      return type(domain_t(), tmin, tmax, n_points, mk);
    }
+   
   };
   
+  // singularity 
   template<typename Opt> struct singularity<retime,matrix_valued,Opt>  { typedef local::tail type;};
+  template<typename Opt> struct singularity<retime,scalar_valued,Opt>  { typedef local::tail type;};
+  
+  // h5 name
   template<typename Opt> struct h5_name<retime,matrix_valued,Opt>      { static std::string invoke(){ return  "GfReTime";}};
+  template<typename Opt> struct h5_name<retime,scalar_valued,Opt>      { static std::string invoke(){ return  "GfReTime_s";}};
 
   /// ---------------------------  evaluator ---------------------------------
-
+  template<typename Opt, typename Target>
-  template<typename Opt>
+   struct evaluator<retime,Target,Opt> {
-   struct evaluator<retime,matrix_valued,Opt> {
     static constexpr int arity = 1;
+    //typedef typename std::conditional < std::is_same<Target, matrix_valued>::value, arrays::matrix_view<std::complex<double>>, std::complex<double>>::type rtype; 
+    typedef typename std::conditional < std::is_same<Target, matrix_valued>::value, arrays::matrix<std::complex<double>>, std::complex<double>>::type rtype; 
     template<typename G>
-     arrays::matrix_view<std::complex<double> >  operator() (G const * g,double t0)  const {
+      rtype operator() (G const * g,double t0)  const {
       auto & data = g->data();
       auto & mesh = g->mesh();
       size_t index; double w; bool in;
       std::tie(in, index, w) = windowing(mesh,t0);
       if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
-      arrays::matrix<std::complex<double> > res = w*data(mesh.index_to_linear(index), arrays::ellipsis()) + (1-w)*data(mesh.index_to_linear(index+1), arrays::ellipsis());
+      return 
-      return res;
+        (1-w) * data(mesh.index_to_linear(index  ), arrays::ellipsis() ) 
+        + w *   data(mesh.index_to_linear(index+1), arrays::ellipsis() );
      }
     template<typename G>
      local::tail_view operator()(G const * g,freq_infty const &) const {return g->singularity();}
    };
 
   /// ---------------------------  data access  ---------------------------------
-
    template<typename Opt> struct data_proxy<retime,matrix_valued,Opt> : data_proxy_array<std::complex<double>,3> {};
+   template<typename Opt> struct data_proxy<retime,scalar_valued,Opt> : data_proxy_array<std::complex<double>,1> {};
 
   // -------------------------------   Factories  --------------------------------------------------
 
+  //matrix_valued
   template<typename Opt> struct factories<retime, matrix_valued,Opt> {
-   typedef gf<retime> gf_t;
+   typedef gf<retime,matrix_valued> gf_t;
-
-   static gf_t make_gf(double tmin, double tmax, size_t n_points, tqa::mini_vector<size_t,2> shape) {
-    typename gf_t::data_non_view_t A(shape.front_append(n_points)); A() =0;
-    return gf_t(mesh<retime,Opt>::make(tmin, tmax, n_points, full_bins), std::move(A), local::tail(shape), nothing());
-   }
    
    static gf_t make_gf(double tmin, double tmax, size_t n_points, tqa::mini_vector<size_t,2> shape, mesh_kind mk) {
     typename gf_t::data_non_view_t A(shape.front_append(n_points)); A() =0;
     return gf_t(mesh<retime,Opt>::make(tmin, tmax, n_points,mk), std::move(A), local::tail(shape), nothing());
    }
    
+   static gf_t make_gf(double tmin, double tmax, size_t n_points, tqa::mini_vector<size_t,2> shape) {
+    typename gf_t::data_non_view_t A(shape.front_append(n_points)); A() =0;
+    return gf_t(mesh<retime,Opt>::make(tmin, tmax, n_points), std::move(A), local::tail(shape), nothing());
+   }
+   
   };
+  
+  //scalar_valued
+  template<typename Opt> struct factories<retime, scalar_valued,Opt> {
+   typedef gf<retime,scalar_valued> gf_t;
+   
+   static gf_t make_gf(double tmin, double tmax, size_t n_points, mesh_kind mk) {
+    typename gf_t::data_non_view_t A(n_points); A() =0;
+    return gf_t(mesh<retime,Opt>::make(tmin, tmax, n_points,mk), std::move(A), local::tail(tqa::mini_vector<size_t,2>(1,1)), nothing());
+   }
+   
+   static gf_t make_gf(double tmin, double tmax, size_t n_points) {
+    typename gf_t::data_non_view_t A(n_points); A() =0;
+    return gf_t(mesh<retime,Opt>::make(tmin, tmax, n_points), std::move(A), local::tail(tqa::mini_vector<size_t,2>(1,1)), nothing());
+   }
+   
+  };
+  
+  
  } // gf_implementation
 }}
 #endif

triqs/gf/two_real_times.hpp

@@ -72,10 +72,16 @@ namespace triqs { namespace gf {
   struct evaluator<two_real_times,matrix_valued,Opt> {
     static constexpr int arity = 2;
     template<typename G>
-     arrays::matrix_view<std::complex<double> > operator() (G const * g, double t0, double t1)  const {
+    arrays::matrix<std::complex<double> > operator() (G const * g, double t0, double t1)  const {
-      auto & m0 = std::get<0>(g->mesh().components()); 
+      auto & data = g->data();
-      double s= m0.x_max()/m0.size();
+      auto & m = std::get<0>(g->mesh().components()); 
-      return g->data()(g->mesh().index_to_linear( typename G::mesh_t::index_t(t0*s, t1*s)), arrays::range(), arrays::range());//mesh.index_to_linear(mesh.point_to_index (t1,t2)));
+      size_t n0,n1; double w0,w1; bool in;
+      std::tie(in, n0, w0) = windowing(m,t0);
+      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+      std::tie(in, n1, w1) = windowing(m,t1);
+      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+      auto gg = [g,data]( size_t n0, size_t n1) {return data(g->mesh().index_to_linear(std::tuple<size_t,size_t>{n0,n1}), arrays::ellipsis());};
+      return  w0 * ( w1*gg(n0,n1) + (1-w1)*gg(n0,n1+1) ) + (1-w0) * ( w1*gg(n0+1,n1) + (1-w1)*gg(n0+1,n1+1)); 
     } 
   };
 
@@ -130,13 +136,14 @@ namespace triqs { namespace gf {
   // -------------------------------   Additionnal free function for this gf  --------------------------------------------------
   
   // from g(t,t') and t, return g(t-t') for any t'>t 
+  // 
  gf<retime> slice (gf_view<two_real_times> const & g, double t) { 
-  auto const & m = std::get<0> (g.mesh().components());
+  auto const & m = std::get<0> (g.mesh().components());   //one-time mesh
-  long it = get_closest_mesh_pt_index(m, t);
+  long it = get_closest_mesh_pt_index(m, t);              //index of t on this mesh
   long nt = m.size() - it;
-  if (it < nt) nt = it ;
+  if (it+1 < nt) nt = it+1 ;                              //nt=length of the resulting GF's mesh
   double dt = m.delta();
-  auto res = make_gf<retime>(0, nt*dt, nt, g(t,t).shape());
+  auto res = make_gf<retime>(0, 2*(nt-1)*dt, nt, g(t,t).shape());
   res() = 0;
   auto _ = arrays::range();// everyone
   for(long sh=0; sh<nt; sh++){