diff --git a/test/speed/gf_s.cpp b/test/speed/gf_s.cpp
index 7fa2f483..7915d14d 100644
--- a/test/speed/gf_s.cpp
+++ b/test/speed/gf_s.cpp
@@ -1,49 +1,45 @@
 #define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
-#include <triqs/gfs.hpp> 
+
+//#define TRIQS_ARRAYS_CHECK_WEAK_REFS
+
+#include <triqs/gfs.hpp>
 using namespace triqs::gfs;
 using namespace triqs::arrays;
-#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl;
-#include <triqs/gfs/local/functions.hpp> 
+#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
+#include <triqs/gfs/local/functions.hpp>
 
 const int nl_interne = 1000;
 const int N = 1000;
 
 typedef double VALUE_TYPE;
-//typedef int VALUE_TYPE;
-inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8.0);}
-//inline VALUE_TYPE fnt(size_t i) { return i;} //*(i+2.0)*(i-8);}
-//inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8);}
+// typedef int VALUE_TYPE;
+// inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8.0);}
+inline VALUE_TYPE fnt(size_t i) { return i; } //*(i+2.0)*(i-8);}
+// inline VALUE_TYPE fnt(size_t i) { return i*(i+2.0)*(i-8);}
 
-struct with_sliding_view {
+struct with_g {
  void operator()() {
 
- double beta =1;
- auto G =  gf<imfreq> { {beta, Fermion,N}, {2,2}};
- G() =0;
-
- //auto slv = G.data_getter.slv;
-
- for (int u =0; u<nl_interne; ++u)
-   //for (int i =0; i<N-1; ++i) { G.data_getter.slv[i](0,0) = fnt(i);}
-   //for (int i =0; i<N-1; ++i) { slv.set(i); slv(0,0) = fnt(i);}
-   //for (int i =0; i<N-1; ++i) G.data_getter(i)(0,0) = fnt(i);
-  //for (int i =0; i<N-1; ++i) G.on_mesh(i)(0,0) = fnt(i);
-  for (int i =0; i<N-1; ++i) G[i](0,0) = fnt(i);
+  double beta = 1;
+  auto G = gf<imfreq>{{beta, Fermion, N}, {2, 2}};
+  G() = 0;
 
+  for (int u = 0; u < nl_interne; ++u)
+   // for (int i =0; i<N-1; ++i) G.on_mesh(i)(0,0) = fnt(i);
+   for (int i = 0; i < N - 1; ++i) G[i](0, 0) = fnt(i);
  }
 };
 
 struct array_code {
  void operator()() {
 
- double beta =1;
- auto G =  gf<imfreq> { {beta, Fermion,N}, {2,2}};
- G() =0;
- auto V = G.data();
-
- for (int u =0; u<nl_interne; ++u)
-   for (int i =0; i<N-1; ++i) V(i,0,0) = fnt(i);
+  double beta = 1;
+  auto G = gf<imfreq>{{beta, Fermion, N}, {2, 2}};
+  G() = 0;
+  auto V = G.data();
 
+  for (int u = 0; u < nl_interne; ++u)
+   for (int i = 0; i < N - 1; ++i) V(i, 0, 0) = fnt(i);
  }
 };
 
@@ -51,10 +47,9 @@ struct array_code {
 
 #include "./speed_tester.hpp"
 int main() {
- try { 
- speed_tester<array_code> (500);
- speed_tester<with_sliding_view> (500);
- //speed_tester<with_slices> (5000);
+ try {
+  speed_tester<array_code>(100);
+  speed_tester<with_g>(100);
  }
  TRIQS_CATCH_AND_ABORT;
 }
diff --git a/test/triqs/gfs/curry_and_fourier.cpp b/test/triqs/gfs/curry_and_fourier.cpp
index d22c0ec6..8a156839 100644
--- a/test/triqs/gfs/curry_and_fourier.cpp
+++ b/test/triqs/gfs/curry_and_fourier.cpp
@@ -33,11 +33,9 @@ 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] =  fourier(G_w_tau_curry0[w]);
- 
- G_w_wn_curry0[w_] << fourier(G_w_tau_curry0[w_]);
- 
- curry<0>(G_w_wn) [w_] << fourier(curry<0>(G_w_tau)[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_] << fourier(G_w_tau_curry0[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/g_k_om.cpp b/test/triqs/gfs/g_k_om.cpp
index 7554f8d7..db2da644 100644
--- a/test/triqs/gfs/g_k_om.cpp
+++ b/test/triqs/gfs/g_k_om.cpp
@@ -16,7 +16,7 @@ int main() {
 
   auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix<double>(2)}};
 
-  auto g_eps = gf<bz>{{bz_, 1000}, {1, 1}};
+  auto g_eps = gf<bz>{{bz_, 100}, {1, 1}};
 
   auto G = gf<cartesian_product<bz, imfreq>>{{{bz_, 100}, {beta, Fermion, 100}}, {1, 1}};
 
diff --git a/test/triqs/gfs/g_k_tail.cpp b/test/triqs/gfs/g_k_tail.cpp
new file mode 100644
index 00000000..3d46bb6a
--- /dev/null
+++ b/test/triqs/gfs/g_k_tail.cpp
@@ -0,0 +1,40 @@
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+#include <triqs/gfs.hpp>
+#include <triqs/gfs/bz.hpp>
+#include <triqs/gfs/m_tail.hpp>
+
+namespace h5 = triqs::h5;
+using namespace triqs::gfs;
+using namespace triqs::clef;
+using namespace triqs::arrays;
+using namespace triqs::lattice;
+using local::tail;
+
+#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
+
+// THE NAME bz is TOO SHORT
+
+int main() {
+ try {
+  double beta = 1;
+
+  auto bz_ = brillouin_zone{bravais_lattice{make_unit_matrix<double>(2)}};
+
+  auto t = tail{1,1};
+  auto G = gf<bz, tail>{{bz_, 100}};
+
+  // try to assign to expression 
+  placeholder<0> k_;
+  placeholder<1> w_;
+  auto eps = make_expr( [](k_t const& k) { return -2 * (cos(k(0)) + cos(k(1))); }) ;
+
+  G(k_) << eps(k_) * t;
+
+  // hdf5 
+  h5::file file("ess_g_k_tail.h5", H5F_ACC_TRUNC );
+  h5_write(file, "g", G);
+
+ 
+ }
+ TRIQS_CATCH_AND_ABORT;
+}
diff --git a/test/triqs/gfs/g_nu_nup.cpp b/test/triqs/gfs/g_nu_nup.cpp
new file mode 100644
index 00000000..1fb9217b
--- /dev/null
+++ b/test/triqs/gfs/g_nu_nup.cpp
@@ -0,0 +1,40 @@
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+#include <triqs/gfs.hpp>
+
+namespace h5 = triqs::h5;
+using namespace triqs::gfs;
+using namespace triqs::clef;
+using namespace triqs::arrays;
+
+#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
+
+int main() {
+ try {
+  double beta = 1;
+
+  //  auto G = gf<cartesian_product<imfreq, imfreq>, tensor_valued<4>, matrix_operations>{
+  //    {{beta, Fermion, 100}, {beta, Fermion, 100}}, {2, 2, 2, 2}};
+
+  /*
+  auto G = gf<imfreq_x2>{{beta, Fermion, 100}, {2, 2}};
+
+  // try to assign to expression 
+  placeholder<0> nu_;
+  placeholder<1> nup_;
+
+  G(nu_, nup_) << 1 / (nu_ + nup_ + 1);
+
+  for (auto w : G.mesh())
+   std::cerr  << G[w](0,0) <<std::endl;
+
+  // inverse 
+  auto Ginv = inverse(G);
+
+  // hdf5 
+  h5::file file("ess_g_nu_nup.h5", H5F_ACC_TRUNC );
+  h5_write(file, "g", G);
+ */
+
+ }
+ TRIQS_CATCH_AND_ABORT;
+}
diff --git a/test/triqs/gfs/g_om_nu_nup.cpp b/test/triqs/gfs/g_om_nu_nup.cpp
new file mode 100644
index 00000000..53cb8a94
--- /dev/null
+++ b/test/triqs/gfs/g_om_nu_nup.cpp
@@ -0,0 +1,44 @@
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+#include <triqs/gfs.hpp>
+
+namespace h5 = triqs::h5;
+using namespace triqs::gfs;
+using namespace triqs::clef;
+using namespace triqs::arrays;
+
+#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
+
+int main() {
+ try {
+  double beta = 1;
+
+  /*
+  auto G = gf<imfreq_bff, tensor_valued<4>>{{{beta, Boson, 3}, {beta, Fermion, 10}, {beta, Fermion, 10}}, {2, 2, 2, 2}};
+
+  // try to assign to expression 
+  placeholder<2> om_;
+  placeholder<0> nu_;
+  placeholder<1> nup_;
+
+  G(om_, nu_, nup_) << 1 / (nu_ + nup_ + om_ + 1);
+
+  for (auto w : G.mesh())
+   std::cerr  << G[w] <<std::endl;
+
+  // hdf5 
+  h5::file file("ess_g_om_nu_nup.h5", H5F_ACC_TRUNC );
+  h5_write(file, "g", G);
+
+  // Curry
+  //auto g_cur = curry<0>(G);
+
+  auto g0 = partial_eval<0>(G,0);
+  //auto gg = g_cur[0];
+
+  //h5_write(file, "g[0]c", gg);
+  h5_write(file, "g[0]", g0);
+
+ */
+  }
+ TRIQS_CATCH_AND_ABORT;
+}
diff --git a/test/triqs/gfs/gf_scalar.cpp b/test/triqs/gfs/gf_scalar.cpp
index 92063764..2880f8c3 100644
--- a/test/triqs/gfs/gf_scalar.cpp
+++ b/test/triqs/gfs/gf_scalar.cpp
@@ -1,4 +1,4 @@
-//#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
 #include <triqs/gfs.hpp> 
 #include <triqs/gfs/local/functions.hpp> 
 using namespace triqs::gfs;
diff --git a/test/triqs/gfs/gfv2.cpp b/test/triqs/gfs/gfv2.cpp
index bf831532..e8464fca 100644
--- a/test/triqs/gfs/gfv2.cpp
+++ b/test/triqs/gfs/gfv2.cpp
@@ -76,6 +76,8 @@ int main() {
   TEST( G( 0) ) ;
   TEST( Gc( 0) ) ;
 
+  // error
+  //TEST (G(1.0_j));
 
   // operations on gf
   G3 = G +2* Gc;
diff --git a/triqs/arrays/linalg/det_and_inverse.hpp b/triqs/arrays/linalg/det_and_inverse.hpp
index 8407a36b..6f544b23 100644
--- a/triqs/arrays/linalg/det_and_inverse.hpp
+++ b/triqs/arrays/linalg/det_and_inverse.hpp
@@ -41,10 +41,16 @@ namespace arrays {
   */
  template <typename A> struct inverse_lazy;
 
- template <class A> inverse_lazy<typename utility::remove_rvalue_ref<A>::type> inverse(A &&a) {
+ template <class A>
+ std14::enable_if_t<ImmutableMatrix<std14::remove_reference_t<A>>::value,
+                    inverse_lazy<typename utility::remove_rvalue_ref<A>::type>>
+ inverse(A &&a) {
   return {std::forward<A>(a)};
  }
 
+ template <class A>
+ std14::enable_if_t<ImmutableArray<std14::remove_reference_t<A>>::value> inverse(A &&a) = delete; // arrays can not be inverted.
+
  // ----------------- implementation -----------------------------------------
 
  // worker takes a contiguous view and compute the det and inverse in two steps.
diff --git a/triqs/arrays/matrix.hpp b/triqs/arrays/matrix.hpp
index b8f4e0e5..37d9a4c1 100644
--- a/triqs/arrays/matrix.hpp
+++ b/triqs/arrays/matrix.hpp
@@ -217,7 +217,7 @@ namespace triqs { namespace arrays {
   }
 
  template<typename ArrayType>
-  matrix_view<typename ArrayType::value_type, ArrayType::opt_flags, ArrayType::traversal_order>
+  matrix_view<typename ArrayType::value_type, ArrayType::opt_flags, ArrayType::traversal_order, true>
   make_matrix_view(ArrayType const & a) { 
    static_assert(ArrayType::rank ==2, "make_matrix_view only works for array of rank 2");
    return a;
diff --git a/triqs/gfs/block.hpp b/triqs/gfs/block.hpp
index 415983c3..21869af2 100644
--- a/triqs/gfs/block.hpp
+++ b/triqs/gfs/block.hpp
@@ -18,8 +18,7 @@
  * TRIQS. If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef TRIQS_GF_BLOCK_H
-#define TRIQS_GF_BLOCK_H
+#pragma once
 #include "./tools.hpp"
 #include "./gf.hpp"
 #include "./local/tail.hpp"
@@ -31,7 +30,7 @@ namespace gfs {
  struct block_index {};
 
  template <typename Opt> struct gf_mesh<block_index, Opt> : discrete_mesh<discrete_domain> {
-  typedef discrete_mesh<discrete_domain> B;
+  using B = discrete_mesh<discrete_domain>;
   gf_mesh() = default;
   gf_mesh(int s) : B(s) {}
   gf_mesh(discrete_domain const &d) : B(d) {}
@@ -42,20 +41,20 @@ namespace gfs {
 
   /// ---------------------------  hdf5 ---------------------------------
 
-  template <typename Target, typename Opt> struct h5_name<block_index, Target, Opt> {
+  template <typename Target> struct h5_name<block_index, Target, nothing, void> {
    static std::string invoke() { return "BlockGf"; }
   };
 
-  template <typename Target, typename Opt> struct h5_rw<block_index, Target, Opt> {
+  template <typename Target> struct h5_rw<block_index, Target, nothing, void> {
 
-   static void write(h5::group gr, gf_const_view<block_index, Target, Opt> g) {
+   static void write(h5::group gr, gf_const_view<block_index, Target, nothing, void> g) {
     for (size_t i = 0; i < g.mesh().size(); ++i) h5_write(gr, g.mesh().domain().names()[i], g._data[i]);
     // h5_write(gr,"symmetry",g._symmetry);
    }
 
-   template <bool IsView> static void read(h5::group gr, gf_impl<block_index, Target, Opt, IsView, false> &g) {
+   template <bool IsView> static void read(h5::group gr, gf_impl<block_index, Target, nothing, void, IsView, false> &g) {
     // does not work : need to read the block name and remake the mesh...
-    g._mesh = gf_mesh<block_index, Opt>(gr.get_all_subgroup_names());
+    g._mesh = gf_mesh<block_index>(gr.get_all_subgroup_names());
     g._data.resize(g._mesh.size());
     // if (g._data.size() != g._mesh.size()) TRIQS_RUNTIME_ERROR << "h5 read block gf : number of block mismatch";
     for (size_t i = 0; i < g.mesh().size(); ++i) h5_read(gr, g.mesh().domain().names()[i], g._data[i]);
@@ -65,15 +64,15 @@ namespace gfs {
 
   /// ---------------------------  data access  ---------------------------------
 
-  template <typename Target, typename Opt>
-  struct data_proxy<block_index, Target, Opt> : data_proxy_vector<typename regular_type_if_exists_else_type<Target>::type> {};
+  template <typename Target>
+  struct data_proxy<block_index, Target, void> : data_proxy_vector<typename regular_type_if_exists_else_type<Target>::type> {};
 
   // -------------------------------   Factories  --------------------------------------------------
 
-  template <typename Target, typename Opt> struct factories<block_index, Target, Opt> {
-   typedef gf_mesh<block_index, Opt> mesh_t;
-   typedef gf<block_index, Target> gf_t;
-   typedef gf_view<block_index, Target> gf_view_t;
+  template <typename Target> struct factories<block_index, Target, nothing, void> {
+   using mesh_t = gf_mesh<block_index>;
+   using gf_t = gf<block_index, Target>;
+   using gf_view_t = gf_view<block_index, Target>;
    struct target_shape_t {};
 
    static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t) { return std::vector<Target>(m.size()); }
@@ -93,54 +92,50 @@ namespace gfs {
  // -------------------------------   Free Factories for regular type  --------------------------------------------------
 
  // from a number and a gf to be copied
- template <typename Variable, typename Target, typename Opt>
- block_gf<Variable, Target, Opt> make_block_gf(int n, gf<Variable, Target, Opt> const &g) {
-  auto V = std::vector<gf<Variable, Target, Opt>>{};
+ template <typename... A> block_gf<A...> make_block_gf(int n, gf<A...> const &g) {
+  auto V = std::vector<gf<A...>>{};
   for (int i = 0; i < n; ++i) V.push_back(g);
   return {{n}, std::move(V), nothing{}, nothing{}, nothing{}};
  }
 
  // from a vector of gf (moving directly)
- template <typename Variable, typename Target, typename Opt>
- block_gf<Variable, Target, Opt> make_block_gf(std::vector<gf<Variable, Target, Opt>> V) {
+ template <typename... A> block_gf<A...> make_block_gf(std::vector<gf<A...>> V) {
   int s = V.size(); // DO NOT use V.size in next statement, the V is moved and the order of arg. evaluation is undefined.
   return {{s}, std::move(V), nothing{}, nothing{}, nothing{}};
  }
 
+ /*
  // from a vector of gf : generalized to have a different type of gf in the vector (e.g. views...)
- template <typename Variable, typename Target, typename Opt, typename GF2>
- block_gf<Variable, Target, Opt> make_block_gf(std::vector<GF2> const &V) {
-  auto V2 = std::vector<gf<Variable, Target, Opt>>{};
+ template <typename Variable, typename Target, typename Singularity, typename Opt, typename GF2>
+ block_gf<Variable, Target, Singularity, Opt> make_block_gf(std::vector<GF2> const &V) {
+  auto V2 = std::vector<gf<Variable, Target, Singularity, Opt>>{};
   for (auto const &g : V) V2.push_back(g);
   return {{int(V.size())}, std::move(V2), nothing{}, nothing{}, nothing{}};
  }
+*/
 
  // from a init list of GF with the correct type
- template <typename Variable, typename Target, typename Opt>
- block_gf<Variable, Target, Opt> make_block_gf(std::initializer_list<gf<Variable, Target, Opt>> const &V) {
+ template <typename... A> block_gf<A...> make_block_gf(std::initializer_list<gf<A...>> const &V) {
   return {{int(V.size())}, V, nothing{}, nothing{}, nothing{}};
  }
 
  // from vector<string> and a gf to be copied
- template <typename Variable, typename Target, typename Opt>
- block_gf<Variable, Target, Opt> make_block_gf(std::vector<std::string> block_names, gf<Variable, Target, Opt> const &g) {
-  auto V = std::vector<gf<Variable, Target, Opt>>{};
+ template <typename... A> block_gf<A...> make_block_gf(std::vector<std::string> block_names, gf<A...> const &g) {
+  auto V = std::vector<gf<A...>>{};
   for (int i = 0; i < block_names.size(); ++i) V.push_back(g);
   return {{block_names}, std::move(V), nothing{}, nothing{}, nothing{}};
  }
 
  // from vector<string>, vector<gf>
- template <typename Variable, typename Target, typename Opt>
- block_gf<Variable, Target, Opt> make_block_gf(std::vector<std::string> block_names, std::vector<gf<Variable, Target, Opt>> V) {
+ template <typename... A> block_gf<A...> make_block_gf(std::vector<std::string> block_names, std::vector<gf<A...>> V) {
   if (block_names.size() != V.size())
    TRIQS_RUNTIME_ERROR << "make_block_gf(vector<string>, vector<gf>) : the two vectors do not have the same size !";
   return {{block_names}, std::move(V), nothing{}, nothing{}, nothing{}};
  }
 
  // from vector<string>, init_list<GF>
- template <typename Variable, typename Target, typename Opt>
- block_gf<Variable, Target, Opt> make_block_gf(std::vector<std::string> block_names,
-                                               std::initializer_list<gf<Variable, Target, Opt>> const &V) {
+ template <typename... A>
+ block_gf<A...> make_block_gf(std::vector<std::string> block_names, std::initializer_list<gf<A...>> const &V) {
   if (block_names.size() != V.size()) TRIQS_RUNTIME_ERROR << "make_block_gf(vector<string>, init_list) : size mismatch !";
   return {{block_names}, V, nothing{}, nothing{}, nothing{}};
  }
@@ -148,10 +143,9 @@ namespace gfs {
  // -------------------------------   Free Factories for view type  --------------------------------------------------
 
  template <typename G0, typename... G>
- gf_view<block_index, typename std::remove_reference<G0>::type::view_type> make_block_gf_view(G0 &&g0, G &&... g) {
-  auto V = std::vector<typename std::remove_reference<G0>::type::view_type>{std::forward<G0>(g0), std::forward<G>(g)...};
+ gf_view<block_index, typename std14::decay_t<G0>::view_type> make_block_gf_view(G0 &&g0, G &&... g) {
+  auto V = std::vector<typename std14::decay_t<G0>::view_type>{std::forward<G0>(g0), std::forward<G>(g)...};
   return {{int(V.size())}, std::move(V), nothing{}, nothing{}, nothing{}};
-  // return { gf_mesh<block_index, Opt> {int(V.size())}, std::move(V), nothing{}, nothing{} } ;
  }
 
  template <typename GF> gf_view<block_index, typename GF::regular_type> make_block_gf_view_from_vector(std::vector<GF> V) {
@@ -167,44 +161,46 @@ namespace gfs {
 
  // ------------------------------- Extend reinterpret_scalar_valued_gf_as_matrix_valued for block gf   ------
 
- template <typename Variable, typename Opt, typename Opt2, bool IsConst>
- gf_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2, IsConst>
- reinterpret_scalar_valued_gf_as_matrix_valued(gf_view<block_index, gf<Variable, scalar_valued, Opt>, Opt2, IsConst> bg) {
-  std::vector<gf_view<Variable, matrix_valued, Opt>> V;
+ // TODO simplify ?
+ template <typename Variable, typename Singularity, typename Opt, bool IsConst>
+ gf_view<block_index, gf<Variable, matrix_valued, Singularity, Opt>, nothing, void, IsConst>
+ reinterpret_scalar_valued_gf_as_matrix_valued(
+     gf_view<block_index, gf<Variable, scalar_valued, Singularity, Opt>, nothing, void, IsConst> bg) {
+  std::vector<gf_view<Variable, matrix_valued, Singularity, Opt>> V;
   for (auto &g : bg) V.push_back(reinterpret_scalar_valued_gf_as_matrix_valued(g));
   return make_block_gf_view_from_vector(std::move(V));
  }
 
- template <typename Variable, typename Opt, typename Opt2>
- gf_const_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2>
- reinterpret_scalar_valued_gf_as_matrix_valued(gf<block_index, gf<Variable, scalar_valued, Opt>, Opt2> const &bg) {
+ template <typename Variable, typename Singularity, typename Opt>
+ block_gf_const_view<Variable, matrix_valued, Singularity, Opt>
+ reinterpret_scalar_valued_gf_as_matrix_valued(block_gf<Variable, scalar_valued, Singularity, Opt> const &bg) {
   return reinterpret_scalar_valued_gf_as_matrix_valued(bg());
  }
 
- template <typename Variable, typename Opt, typename Opt2>
- gf_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2>
- reinterpret_scalar_valued_gf_as_matrix_valued(gf<block_index, gf<Variable, scalar_valued, Opt>, Opt2> &bg) {
+ template <typename Variable, typename Singularity, typename Opt>
+ block_gf_view<Variable, matrix_valued, Singularity, Opt>
+ reinterpret_scalar_valued_gf_as_matrix_valued(block_gf<Variable, scalar_valued, Singularity, Opt> &bg) {
   return reinterpret_scalar_valued_gf_as_matrix_valued(bg());
  }
 
  // -------------------------------   Free functions   --------------------------------------------------
 
  // a simple function to get the number of blocks
- template <typename T> size_t n_blocks(gf<block_index, T> const &g) { return g.mesh().size(); }
- template <typename T> size_t n_blocks(gf_view<block_index, T> const &g) { return g.mesh().size(); }
+ template <typename... T> size_t n_blocks(gf<block_index, T...> const &g) { return g.mesh().size(); }
+ template <typename... T> size_t n_blocks(gf_view<block_index, T...> const &g) { return g.mesh().size(); }
 
  // -------------------------------   an iterator over the blocks --------------------------------------------------
 
- template<typename T> using __get_target = typename std::remove_reference<decltype(std::declval<T>()[0])>::type;
+ template <typename T> using __get_target = std14::remove_reference_t<decltype(std::declval<T>()[0])>;
 
  // iterator
  template <typename G>
  class block_gf_iterator
      : public boost::iterator_facade<block_gf_iterator<G>, __get_target<G>, boost::forward_traversal_tag, __get_target<G> &> {
   friend class boost::iterator_core_access;
-  typedef typename std::remove_reference<G>::type big_gf_t;
+  using big_gf_t = typename std::remove_reference<G>::type;
   big_gf_t &big_gf;
-  typedef typename big_gf_t::mesh_t::const_iterator mesh_iterator_t;
+  using mesh_iterator_t = typename big_gf_t::mesh_t::const_iterator;
   mesh_iterator_t mesh_it;
 
   __get_target<G> &dereference() const { return big_gf[*mesh_it]; }
@@ -217,25 +213,27 @@ namespace gfs {
  };
 
  //------------
- template <typename Target, typename Opt, bool B, bool C>
- block_gf_iterator<gf_impl<block_index, Target, Opt, B, C>> begin(gf_impl<block_index, Target, Opt, B, C> &bgf) {
+ template <typename Target, typename Singularity, typename Opt, bool B, bool C>
+ block_gf_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
+ begin(gf_impl<block_index, Target, Singularity, Opt, B, C> &bgf) {
   return {bgf, false};
  }
 
  //------------
- template <typename Target, typename Opt, bool B, bool C>
- block_gf_iterator<gf_impl<block_index, Target, Opt, B, C>> end(gf_impl<block_index, Target, Opt, B, C> &bgf) {
+ template <typename Target, typename Singularity, typename Opt, bool B, bool C>
+ block_gf_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
+ end(gf_impl<block_index, Target, Singularity, Opt, B, C> &bgf) {
   return {bgf, true};
  }
 
  //----- const iterator
  template <typename G>
- class block_gf_const_iterator
-     : public boost::iterator_facade<block_gf_const_iterator<G>, __get_target<G>, boost::forward_traversal_tag, __get_target<G> const &> {
+ class block_gf_const_iterator : public boost::iterator_facade<block_gf_const_iterator<G>, __get_target<G>,
+                                                               boost::forward_traversal_tag, __get_target<G> const &> {
   friend class boost::iterator_core_access;
-  typedef typename std::remove_reference<G>::type big_gf_t;
+  using big_gf_t = std14::remove_reference_t<G>;
   big_gf_t const &big_gf;
-  typedef typename big_gf_t::mesh_t::const_iterator mesh_iterator_t;
+  using mesh_iterator_t = typename big_gf_t::mesh_t::const_iterator;
   mesh_iterator_t mesh_it;
 
   __get_target<G> const &dereference() const { return big_gf[*mesh_it]; }
@@ -247,26 +245,29 @@ namespace gfs {
   bool at_end() const { return mesh_it.at_end(); }
  };
 
- template <typename Target, typename Opt, bool B, bool C>
- block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> begin(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
+ template <typename Target, typename Singularity, typename Opt, bool B, bool C>
+ block_gf_const_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
+ begin(gf_impl<block_index, Target, Singularity, Opt, B, C> const &bgf) {
   return {bgf, false};
  }
 
- template <typename Target, typename Opt, bool B, bool C>
- block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> end(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
+ template <typename Target, typename Singularity, typename Opt, bool B, bool C>
+ block_gf_const_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
+ end(gf_impl<block_index, Target, Singularity, Opt, B, C> const &bgf) {
   return {bgf, true};
  }
 
- template <typename Target, typename Opt, bool B, bool C>
- block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> cbegin(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
+ template <typename Target, typename Singularity, typename Opt, bool B, bool C>
+ block_gf_const_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
+ cbegin(gf_impl<block_index, Target, Singularity, Opt, B, C> const &bgf) {
   return {bgf, false};
  }
 
- template <typename Target, typename Opt, bool B, bool C>
- block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> cend(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
+ template <typename Target, typename Singularity, typename Opt, bool B, bool C>
+ block_gf_const_iterator<gf_impl<block_index, Target, Singularity, Opt, B, C>>
+ cend(gf_impl<block_index, Target, Singularity, Opt, B, C> const &bgf) {
   return {bgf, true};
  }
 }
 }
-#endif
 
diff --git a/triqs/gfs/bz.hpp b/triqs/gfs/bz.hpp
index 2558d2e8..fb420c69 100644
--- a/triqs/gfs/bz.hpp
+++ b/triqs/gfs/bz.hpp
@@ -37,7 +37,7 @@ namespace gfs {
  namespace gfs_implementation {
 
   // h5 name
-  template <typename Opt> struct h5_name<bz, matrix_valued, Opt> {
+  template <typename Singularity, typename Opt> struct h5_name<bz, matrix_valued, Singularity, Opt> {
    static std::string invoke() { return "BZ"; }
   };
 
@@ -60,7 +60,7 @@ namespace gfs {
 
   // ------------- evaluator  -------------------
   // handle the case where the matsu. freq is out of grid...
-  template <typename Target, typename Opt> struct evaluator<bz, Target, Opt> {
+  template <typename Target, typename Singularity, typename Opt> struct evaluator<bz, Target, Singularity, Opt> {
    static constexpr int arity = 1;
 
    template <typename G>
diff --git a/triqs/gfs/curry.hpp b/triqs/gfs/curry.hpp
index 8384aae4..d589e7cd 100644
--- a/triqs/gfs/curry.hpp
+++ b/triqs/gfs/curry.hpp
@@ -18,13 +18,16 @@
  * TRIQS. If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef TRIQS_GF_CURRY_H
-#define TRIQS_GF_CURRY_H
+#pragma once
 #include "./product.hpp"
 namespace triqs { namespace gfs { 
 
  template<typename F> struct lambda_valued {};
 
+ template <typename Var, typename M, typename L> gf_view<Var, lambda_valued<L>> make_gf_view_lambda_valued(M m, L l) {
+  return {std::move(m), l, nothing(), nothing(), {}};
+ }
+
  namespace gfs_implementation { 
   
   /// ---------------------------  data access  ---------------------------------
@@ -34,7 +37,7 @@ namespace triqs { namespace gfs {
   /// ---------------------------  Factories ---------------------------------
 
   template<typename F, typename Opt, typename ... Ms>
-   struct factories<cartesian_product<Ms...>, lambda_valued<F>, Opt> {};
+   struct factories<cartesian_product<Ms...>, lambda_valued<F>, nothing, Opt> {};
 
   /// ---------------------------  partial_eval ---------------------------------
   // partial_eval<0> (g, 1) returns :  x -> g(1,x)
@@ -42,82 +45,85 @@ namespace triqs { namespace gfs {
 
   // a technical trait: from a tuple of mesh, return the mesh (either M if it is a tuple of size 1, or the corresponding cartesian_product<M..>).
   template<typename ... Ms> struct cart_prod_impl;
-  template<typename ... Ms> using cart_prod = typename cart_prod_impl<Ms...>::type;
   template<typename ... Ms> struct cart_prod_impl<std::tuple<Ms...>> { using type = cartesian_product<Ms...>;};
   template<typename M> struct cart_prod_impl<std::tuple<M>> { using type = M;};
+  template<typename ... Ms> using cart_prod = typename cart_prod_impl<Ms...>::type;
 
-  template<typename M0, typename M1, typename ...M> auto rm_tuple_of_size_one(std::tuple<M0,M1,M...> const & t) DECL_AND_RETURN(t);
-  template<typename M> auto rm_tuple_of_size_one(std::tuple<M> const & t) DECL_AND_RETURN(std::get<0>(t));
+  // The implementation (can be overloaded for some types), so put in a struct to have partial specialization
+  template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst> struct partial_eval_impl;
 
-  // as_tuple leaves a tuple intact and wrap everything else in a tuple...
-  template<typename T> std::tuple<T> as_tuple(T && x) { return std::tuple<T> {std::forward<T>(x)};}
-  template<typename ... T> std::tuple<T...> as_tuple(std::tuple<T...> && x) { return std::forward<T...>(x);}
-  template<typename ... T> std::tuple<T...> const &  as_tuple(std::tuple<T...> const & x) { return x;}
-  template<typename ... T> std::tuple<T...> & as_tuple(std::tuple<T...> & x) { return x;}
-
-  template <int... pos, typename Opt, typename Target, bool IsConst, typename IT, typename... Ms>
-  gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, IsConst>
-  partial_eval(gf_view<cartesian_product<Ms...>, Target, Opt, IsConst> g, IT index) {
-    // meshes of the returned gf_view : just drop the mesh of the evaluated variables
-    auto meshes_tuple_partial = triqs::tuple::filter_out<pos...>(g.mesh().components());
-    // a view of the array of g, with the dimension sizeof...(Ms)
-    auto arr = reinterpret_linear_array(g.mesh(),g.data()); // NO the second () forces a view
-    // now rebuild a tuple of the size sizeof...(Ms), containing the indices and range at the position of evaluated variables.
-    auto arr_args = triqs::tuple::inverse_filter<sizeof...(Ms),pos...>(as_tuple(index), arrays::range());
-    // from it, we make a slice of the array of g, corresponding to the data of the returned gf_view
-    auto arr2 =  triqs::tuple::apply(arr, std::tuple_cat(arr_args, std::make_tuple(arrays::ellipsis{})));
-    // finally, we build the view on this data. 
-    using r_t = gf_view< cart_prod< triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>> ,Target, Opt,IsConst>;
-    return r_t{ rm_tuple_of_size_one(meshes_tuple_partial), arr2,  typename r_t::singularity_non_view_t{}, typename r_t::symmetry_t{} };
-   }
-
-  template <int... pos, typename Opt, typename Target, typename IT, typename... Ms>
-  gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, false>
-  partial_eval(gf<cartesian_product<Ms...>, Target, Opt> & g, IT index) {
-   return partial_eval<pos...>(g(),index);
+  // The user function
+  template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt, bool C, typename... T>
+  auto partial_eval(gf_view<Variable, Target, Singularity, Opt, C> g, T&&... x) {
+   return partial_eval_impl<Variable, Target, Singularity, Opt, C>::template invoke<pos...>(g(), std::forward<T>(x)...);
   }
 
-  template <int... pos, typename Opt, typename Target, typename IT, typename... Ms>
-  gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, true>
-  partial_eval(gf<cartesian_product<Ms...>, Target, Opt> const& g, IT index) {
-   return partial_eval<pos...>(g(),index);
+  template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt, typename... T>
+  auto partial_eval(gf<Variable, Target, Singularity, Opt>& g, T&&... x) {
+   return partial_eval_impl<Variable, Target, Singularity, Opt, false>::template invoke<pos...>(g(), std::forward<T>(x)...);
   }
- 
+
+  template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt, typename... T>
+  auto partial_eval(gf<Variable, Target, Singularity, Opt> const& g, T&&... x) {
+   return partial_eval_impl<Variable, Target, Singularity, Opt, true>::template invoke<pos...>(g(), std::forward<T>(x)...);
+  }
+
   /// ---------------------------  curry  ---------------------------------
   // curry<0>(g) returns : x-> y... -> g(x,y...)
   // curry<1>(g) returns : y-> x,z... -> g(x,y,z...)
-  
-  // to adapt the partial_eval as a polymorphic lambda (replace by a lambda in c++14)
-  template<typename Gview, int ... pos> struct curry_polymorphic_lambda { 
-   Gview g;
-   template<typename ...I> auto operator()(I ... i) const DECL_AND_RETURN(partial_eval<pos...>(g,std::make_tuple(i...)));
-   friend int get_shape(curry_polymorphic_lambda const&) { return 0;}// no shape here, but needed for compilation
-   //void resize(int){}
+
+  // The implementation (can be overloaded for some types)
+  template <int... pos, typename Target, typename Singularity, typename Opt, bool IsConst, typename... Ms>
+  auto curry_impl(gf_view<cartesian_product<Ms...>, Target, Singularity, Opt, IsConst> g) {
+   // pick up the meshed corresponding to the curryed variables
+   auto meshes_tuple = triqs::tuple::filter<pos...>(g.mesh().components());
+   using var_t = cart_prod<triqs::tuple::filter_t<std::tuple<Ms...>, pos...>>;
+   auto m = triqs::tuple::apply_construct<gf_mesh<var_t>>(meshes_tuple);
+   auto l = [g](auto&&... x) { return partial_eval<pos...>(g, x...); };
+   return make_gf_view_lambda_valued<var_t>(m, l);
   };
 
-  // curry function ...
-  template <int... pos, typename Target, typename Opt, bool IsConst, typename... Ms>
-  gf_view<cart_prod<triqs::tuple::filter_t<std::tuple<Ms...>, pos...>>,
-                lambda_valued<curry_polymorphic_lambda<gf_view<cartesian_product<Ms...>, Target, Opt,IsConst>, pos...>>, Opt, IsConst>
-  curry(gf_view<cartesian_product<Ms...>, Target, Opt, IsConst> g) {
-    // pick up the meshed corresponding to the curryed variables
-    auto meshes_tuple = triqs::tuple::filter<pos...>(g.mesh().components());
-    // building the view
-    return {rm_tuple_of_size_one(meshes_tuple),curry_polymorphic_lambda<gf_view<cartesian_product<Ms...>, Target,Opt,IsConst>, pos ...>{g}, nothing(), nothing()}; 
-    //using m_t = gf_mesh< cart_prod< triqs::tuple::filter_t<std::tuple<Ms...>,pos...>>>; 
-    //return {triqs::tuple::apply_construct<m_t>(meshes_tuple),curry_polymorphic_lambda<gf_view<cartesian_product<Ms...>, Target,Opt>, pos ...>{g}, nothing(), nothing()}; 
-   };
-  
-  template <int... pos, typename Target, typename Opt, typename... Ms>
-  auto curry(gf<cartesian_product<Ms...>, Target, Opt> & g) DECL_AND_RETURN(curry<pos...>(g()));
-  
-  template <int... pos, typename Target, typename Opt, typename... Ms>
-  auto curry(gf<cartesian_product<Ms...>, Target, Opt> const & g) DECL_AND_RETURN(curry<pos...>(g()));
+  // The user function
+  template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst>
+  auto curry(gf_view<Variable, Target, Singularity, Opt, IsConst> g) {
+   return curry_impl<pos...>(g());
+  }
+  template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt>
+  auto curry(gf<Variable, Target, Singularity, Opt>& g) {
+   return curry_impl<pos...>(g());
+  }
+  template <int... pos, typename Variable, typename Target, typename Singularity, typename Opt>
+  auto curry(gf<Variable, Target, Singularity, Opt> const& g) {
+   return curry_impl<pos...>(g());
+  }
+
+ //---------------------------------------------
+
+  // A generic impl. for cartesian product
+  template <typename Target, typename Singularity, typename Opt, bool IsConst, typename... Ms>
+  struct partial_eval_impl<cartesian_product<Ms...>, Target, Singularity, Opt, IsConst> {
+
+   template <int... pos, typename... T>
+   static auto invoke(gf_view<cartesian_product<Ms...>, Target, Singularity, Opt, IsConst> g, T const&... x) {
+    using var_t = cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>;
+    // meshes of the returned gf_view : just drop the mesh of the evaluated variables
+    auto meshes_tuple_partial = triqs::tuple::filter_out<pos...>(g.mesh().components());
+    auto m = triqs::tuple::apply_construct<gf_mesh<var_t>>(meshes_tuple_partial);
+    // now rebuild a tuple of the size sizeof...(Ms), containing the indices and range at the position of evaluated variables.
+    auto arr_args = triqs::tuple::inverse_filter<sizeof...(Ms), pos...>(std::make_tuple(x...), arrays::range());
+    // from it, we make a slice of the array of g, corresponding to the data of the returned gf_view
+    auto arr2 = triqs::tuple::apply(g.data(), std::tuple_cat(arr_args, std::make_tuple(arrays::ellipsis{})));
+    auto singv = partial_eval<pos...>(g.singularity(), x...);
+    using r_sing_t = typename decltype(singv)::regular_type;
+    // finally, we build the view on this data.
+    using r_t = gf_view<var_t, Target, r_sing_t, Opt, IsConst>;
+    return r_t{m, arr2, singv, {}, {}};
+   }
+  };
 
  } // gf_implementation
  using gfs_implementation::partial_eval;
  using gfs_implementation::curry;
 }}
-#endif
 
 
diff --git a/triqs/gfs/data_proxies.hpp b/triqs/gfs/data_proxies.hpp
index c61ed7fe..35de0686 100644
--- a/triqs/gfs/data_proxies.hpp
+++ b/triqs/gfs/data_proxies.hpp
@@ -23,96 +23,98 @@
 #include <triqs/utility/first_include.hpp>
 #include <utility>
 #include <triqs/arrays.hpp>
-//#include "./matrix_view_proxy.hpp"
 #include "../arrays/matrix_tensor_proxy.hpp"
 
+#define TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS
+
 namespace triqs { namespace gfs {
 
- //---------------------------- generic case array of dim R----------------------------------
- 
- template<typename T, int R> struct data_proxy_array { 
-  /// The storage
-  typedef arrays::array<T, R> storage_t;
-  typedef typename storage_t::view_type storage_view_t;
-  typedef typename storage_t::const_view_type storage_const_view_t;
+ //---------------------------- common stuff for array proxies ----------------------------------
 
-  /// The data access 
-  auto operator()(storage_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i));
-  auto operator()(storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
-  auto operator()(storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i));
-  auto operator()(storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
-  auto operator()(storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
-  auto operator()(storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
+ template <typename T, int D> struct data_proxy_array_common {
+  using storage_t = arrays::array<T, D>;
+  using storage_view_t = typename storage_t::view_type;
+  using storage_const_view_t = typename storage_t::const_view_type;
 
-#ifdef TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS
-  auto operator()(storage_t            & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  auto operator()(storage_t const      & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  auto operator()(storage_view_t       & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  auto operator()(storage_view_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  auto operator()(storage_const_view_t       & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  auto operator()(storage_const_view_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-#endif
+  // from the shape of the mesh and the target, make the shape of the array. default is to glue them
+  template <typename S1, typename S2> static auto join_shape(S1 const& s1, S2 const& s2) RETURN(join(s1, s2));
 
   template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs)           { data() = std::forward<RHS>(rhs);}
   template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); }
  };
 
- //---------------------------- 3d array : returns matrices in this case ! ----------------------------------
- 
- template<typename T> struct data_proxy_array<T,3> { 
-  /// The storage
-  typedef arrays::array<T,3>             storage_t;
-  typedef typename storage_t::view_type  storage_view_t;
-  typedef typename storage_t::const_view_type storage_const_view_t;
-
-  /// The data access
-  auto operator()(storage_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i));
-  auto operator()(storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
-  auto operator()(storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i));
-  auto operator()(storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
-  auto operator()(storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
-  auto operator()(storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
-
-#ifdef TRIQS_DATA_PROXIES_OLD_MATRIX_VIEW_PROXY
-  arrays::matrix_view_proxy<storage_t,0>            operator()(storage_t       & data, size_t i) const { return arrays::matrix_view_proxy<storage_t,0>(data,i); } 
-  arrays::const_matrix_view_proxy<storage_t,0>      operator()(storage_t const & data, size_t i) const { return arrays::const_matrix_view_proxy<storage_t,0>(data,i); } 
-  arrays::matrix_view_proxy<storage_view_t,0>       operator()(storage_view_t       & data, size_t i) const { return arrays::matrix_view_proxy<storage_view_t,0>(data,i); } 
-  arrays::const_matrix_view_proxy<storage_view_t,0> operator()(storage_view_t const & data, size_t i) const { return arrays::const_matrix_view_proxy<storage_view_t,0>(data,i); } 
-#endif
+ //---------------------------- generic case array of dim R----------------------------------
 
+ template <typename T, int R> struct data_proxy_array : data_proxy_array_common<T, R> {
+  using B = data_proxy_array_common<T, R>;
+/// The data access
 #ifdef TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS
-  auto operator()(storage_t            & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  auto operator()(storage_t const      & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
+  template <typename S> auto operator()(S& data, long i) const DECL_AND_RETURN(data(i, arrays::ellipsis()));
+#else
+  auto operator()(B::storage_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i));
+  auto operator()(B::storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
+  auto operator()(B::storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_tensor_proxy(data, i));
+  auto operator()(B::storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
+  auto operator()(B::storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
+  auto operator()(B::storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_tensor_proxy(data, i));
+#endif
+ };
 
-  auto operator()(storage_view_t       & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  auto operator()(storage_view_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  
-  auto operator()(storage_const_view_t       & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
-  auto operator()(storage_const_view_t const & data, size_t i) const DECL_AND_RETURN(data(i,arrays::ellipsis()));
+ //---------------------------- 3d array : returns matrices in this case ! ----------------------------------
+
+ template <typename T> struct data_proxy_array<T, 3> : data_proxy_array_common<T, 3> {
+  using B = data_proxy_array_common<T, 3>;
+#ifdef TRIQS_GF_DATA_PROXIES_WITH_SIMPLE_VIEWS
+  template <typename S> auto operator()(S & data, long i) const RETURN(make_matrix_view(data(i, arrays::ellipsis())));
+#else
+  /// The data access
+  auto operator()(B::storage_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i));
+  auto operator()(B::storage_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
+  auto operator()(B::storage_view_t& data, long i) const DECL_AND_RETURN(arrays::make_matrix_proxy(data, i));
+  auto operator()(B::storage_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
+  auto operator()(B::storage_const_view_t& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
+  auto operator()(B::storage_const_view_t const& data, long i) const DECL_AND_RETURN(arrays::make_const_matrix_proxy(data, i));
 #endif
 
-  template <typename S, typename RHS> static void assign_to_scalar(S& data, RHS&& rhs) { data() = std::forward<RHS>(rhs); }
-  template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); }
  };
 
  //---------------------------- 1d array ----------------------------------
 
- template<typename T> struct data_proxy_array<T,1>{ 
-  /// The storage
-  typedef arrays::array<T,1>             storage_t;
-  typedef typename storage_t::view_type  storage_view_t;
-  typedef typename storage_t::const_view_type storage_const_view_t;
+ template <typename T> struct data_proxy_array<T, 1> : data_proxy_array_common<T, 1> {
+  template <typename S> AUTO_DECL operator()(S& data, long i) const RETURN(data(i));
+ };
 
-  /// The data access 
-  auto operator()(storage_t       & data,size_t i) const -> decltype(data(i)) { return data(i);}
-  auto operator()(storage_t const & data,size_t i) const -> decltype(data(i)) { return data(i);}
-  auto operator()(storage_view_t       & data,size_t i) const -> decltype(data(i)) { return data(i);}
-  auto operator()(storage_view_t const & data,size_t i) const -> decltype(data(i)) { return data(i);}
-  auto operator()(storage_const_view_t       & data,size_t i) const -> decltype(data(i)) { return data(i);}
-  auto operator()(storage_const_view_t const & data,size_t i) const -> decltype(data(i)) { return data(i);}
+ //---------------------------- multi variable ----------------------------------
 
-  template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs)           { data() = std::forward<RHS>(rhs);}
-  template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) { data.rebind(rhs.data()); }
+ template <typename T, int TotalDim> struct data_proxy_array_multivar : data_proxy_array_common<T, TotalDim> {
+  // using the standard technique from tuple::apply with a sequence
+  template <typename S, typename Tu, size_t... Is>
+  AUTO_DECL _impl(S& data, Tu const& tu, std14::index_sequence<Is...>) const RETURN(data(std::get<Is>(tu)..., arrays::ellipsis()));
+  template <typename S, typename Tu>
+  AUTO_DECL operator()(S& data, Tu const& tu) const RETURN(_impl(data, tu, triqs::tuple::_get_seq<Tu>()));
+ };
+
+ //---------------------------- multi variable ----------------------------------
+
+ template <typename T, int TotalDim> struct data_proxy_array_multivar_matrix_valued : data_proxy_array_common<T, TotalDim> {
+  // using the standard technique from tuple::apply with a sequence
+  template <typename S, typename Tu, size_t... Is>
+  AUTO_DECL _impl(S& data, Tu const& tu, std14::index_sequence<Is...>) const RETURN(make_matrix_view(data(std::get<Is>(tu)..., arrays::range(), arrays::range())));
+  template <typename S, typename Tu>
+  AUTO_DECL operator()(S& data, Tu const& tu) const RETURN(_impl(data, tu, triqs::tuple::_get_seq<Tu>()));
+ };
+
+ //---------------------------- multi variable with index mixer----------------------------------
+ template <typename T, int Nvar, int TargetDim, typename IndexMixer>
+ struct data_proxy_array_index_mixer : data_proxy_array_common<T, Nvar + TargetDim> {
+
+  template <typename S1, typename S2> static utility::mini_vector<int, Nvar + TargetDim> join_shape(S1 const& s1, S2 const& s2) {
+   return tuple::apply_construct_parenthesis<utility::mini_vector<int, Nvar + TargetDim>>(IndexMixer::invoke(s1, s2));
+  }
+
+  template <typename S, typename Tu> auto operator()(S& data, Tu const& tu) const {
+   return tuple::apply(data, IndexMixer::invoke(tu, tuple::make_tuple_repeat<TargetDim>(arrays::range())));
+  }
  };
 
  //---------------------------- vector ----------------------------------
@@ -128,26 +130,25 @@ namespace triqs { namespace gfs {
   //template<typename X> view_proxy & operator = (X && x) { V::operator=( std::forward<X>(x) ); return *this;} 
  };
 
- template<typename T> struct data_proxy_vector { 
-  typedef typename T::view_type Tv;
-  typedef typename T::const_view_type Tcv;
+ template <typename T> struct data_proxy_vector {
+  using Tv = typename T::view_type;
+  using Tcv = typename T::const_view_type;
 
   /// The storage
-  typedef std::vector<T> storage_t;
-  typedef std::vector<view_proxy<Tv>> storage_view_t;
-  typedef std::vector<view_proxy<Tcv>> storage_const_view_t;
+  using storage_t = std::vector<T>;
+  using storage_view_t = std::vector<view_proxy<Tv>>;
+  using storage_const_view_t = std::vector<view_proxy<Tcv>>;
 
-  /// The data access 
+  /// The data access
+  template <typename S> AUTO_DECL operator()(S& data, size_t i) const  RETURN(data[i]);
+
+  /*
   T        &  operator()(storage_t  &           data, size_t i) { return data[i];}
   T  const &  operator()(storage_t      const & data, size_t i) const { return data[i];}
   Tv       &  operator()(storage_view_t &       data, size_t i)       { return data[i];}
   Tv const &  operator()(storage_view_t const & data, size_t i) const { return data[i];}
   Tcv       &  operator()(storage_const_view_t &       data, size_t i)       { return data[i];}
   Tcv const &  operator()(storage_const_view_t const & data, size_t i) const { return data[i];}
-/*Tv operator()(storage_view_t &       data, size_t i) const      { return data[i];}
-  Tv operator()(storage_view_t const & data, size_t i) const { return data[i];}
-  Tcv operator()(storage_const_view_t &       data, size_t i) const      { return data[i];}
-  Tcv operator()(storage_const_view_t const & data, size_t i) const { return data[i];}
 */
 
   template<typename S, typename RHS> static void assign_to_scalar   (S & data, RHS && rhs) {for (size_t i =0; i<data.size(); ++i) data[i] = rhs;}
@@ -156,16 +157,15 @@ namespace triqs { namespace gfs {
 
  //---------------------------- lambda ----------------------------------
 
- template<typename F> struct data_proxy_lambda { 
+ template <typename F> struct data_proxy_lambda {
 
   /// The storage
-  typedef F storage_t;
-  typedef F storage_view_t;
-  typedef F storage_const_view_t;
+  using storage_t = F;
+  using storage_view_t = F;
+  using storage_const_view_t = F;
 
-  /// The data access 
-  auto operator()(storage_t  &           data, size_t i)       DECL_AND_RETURN( data(i));
-  auto operator()(storage_t      const & data, size_t i) const DECL_AND_RETURN( data(i));
+  /// The data access
+  template <typename S, typename ... I> AUTO_DECL operator()(S& data, I const& ...i) const RETURN(data(i...));
 
   template<typename S, typename RHS> static void assign_to_scalar   (S & data, RHS && rhs) = delete;
   template <typename ST, typename RHS> static void rebind(ST& data, RHS&& rhs) = delete;
diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp
index 8103a558..320a1245 100644
--- a/triqs/gfs/gf.hpp
+++ b/triqs/gfs/gf.hpp
@@ -2,7 +2,7 @@
  *
  * TRIQS: a Toolbox for Research in Interacting Quantum Systems
  *
- * Copyright (C) 2012-2013 by M. Ferrero, O. Parcollet
+ * Copyright (C) 2012-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
@@ -42,40 +42,58 @@ namespace gfs {
  using arrays::matrix_view;
  using triqs::make_clone;
 
+ // The default target, for each Variable.
+ template <typename Variable> struct gf_default_target {
+  using type = matrix_valued;
+ };
+ template <typename Variable> using gf_default_target_t = typename gf_default_target<Variable>::type;
+
+ // The default singularity, for each Variable.
+ template <typename Variable, typename Target> struct gf_default_singularity {
+  using type = nothing;
+ };
+ template <typename Variable, typename Target>
+ using gf_default_singularity_t = typename gf_default_singularity<Variable, Target>::type;
+
  // the gf mesh
  template <typename Variable, typename Opt = void> struct gf_mesh;
 
  // The regular type
- template <typename Variable, typename Target = matrix_valued, typename Opt = void> class gf;
+ template <typename Variable, typename Target = gf_default_target_t<Variable>,
+           typename Singularity = gf_default_singularity_t<Variable, Target>, typename Opt = void>
+ class gf;
 
  // The view type
- template <typename Variable, typename Target = matrix_valued, typename Opt = void, bool IsConst = false> class gf_view;
+ template <typename Variable, typename Target = gf_default_target_t<Variable>,
+           typename Singularity = gf_default_singularity_t<Variable, Target>, typename Opt = void, bool IsConst = false>
+ class gf_view;
 
  // The const view type
- template <typename Variable, typename Target = matrix_valued, typename Opt = void>
- using gf_const_view = gf_view<Variable, Target, Opt, true>;
+ template <typename Variable, typename Target = gf_default_target_t<Variable>,
+           typename Singularity = gf_default_singularity_t<Variable, Target>, typename Opt = void>
+ using gf_const_view = gf_view<Variable, Target, Singularity, Opt, true>;
 
  // the implementation class
- template <typename Variable, typename Target, typename Opt, bool IsView, bool IsConst> class gf_impl;
+ template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst> class gf_impl;
 
  // various implementation traits
  namespace gfs_implementation { // never use using of this...
 
   // evaluator regroup functions to evaluate the function.
-  template <typename Variable, typename Target, typename Opt> struct evaluator {
+  template <typename Variable, typename Target, typename Singularity, typename Opt> struct evaluator {
    static constexpr int arity = 0;
   };
 
   // closest_point mechanism
-  template <typename Variable, typename Target, typename Opt> struct get_closest_point;
+  template <typename Variable, typename Target, typename Singularity, typename Opt> struct get_closest_point;
 
   // singularity
-  template <typename Variable, typename Target, typename Opt> struct singularity {
-   using type = nothing;
-  };
+  //template <typename Variable, typename Target, typename Opt> struct singularity {
+  // using type = nothing;
+  //};
 
   // symmetry
-  template <typename Variable, typename Target, typename Opt> struct symmetry {
+  template <typename Variable, typename Target, typename Singularity, typename Opt> struct symmetry {
    using type = nothing;
   };
 
@@ -94,21 +112,29 @@ namespace gfs {
   template <typename Variable, typename Target, typename Opt, typename Enable = void> struct data_proxy;
 
   // Traits to read/write in hdf5 files. Can be specialized for some case (Cf block). Defined below
-  template <typename Variable, typename Target, typename Opt> struct h5_name; // value is a const char
-  template <typename Variable, typename Target, typename Opt> struct h5_rw;
+  template <typename Variable, typename Target, typename Singularity, typename Opt> struct h5_name; // value is a const char
+  template <typename Variable, typename Target, typename Singularity, typename Opt> struct h5_rw;
 
   // factories regroup all factories (constructors..) for all types of gf. Defaults implemented below.
-  template <typename Variable, typename Target, typename Opt> struct factories;
+  template <typename Variable, typename Target, typename Singularity, typename Opt> struct factories;
 
  } // gfs_implementation
 
  // The trait that "marks" the Green function
  TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT(ImmutableGreenFunction);
 
- template <typename G> auto get_gf_data_shape(G const &g) DECL_AND_RETURN(g.get_data_shape());
+ template <typename Variable, typename Target, typename Singularity, typename Opt>
+ auto get_gf_data_shape(gf<Variable, Target, Singularity, Opt> const &g) RETURN(get_shape(g.data()));
 
- template <typename Variable, typename Target, typename Opt, bool IsView, bool IsConst>
- auto get_target_shape(gf_impl<Variable, Target, Opt, IsView, IsConst> const &g) DECL_AND_RETURN(g.data().shape().front_pop());
+ template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst>
+ auto get_gf_data_shape(gf_view<Variable, Target, Singularity, Opt, IsConst> const &g) RETURN(get_shape(g.data()));
+
+ template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst>
+ auto get_gf_data_shape(gf_impl<Variable, Target, Singularity, Opt, IsView, IsConst> const &g) RETURN(get_shape(g.data()));
+
+ template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst>
+ auto get_target_shape(gf_impl<Variable, Target, Singularity, Opt, IsView, IsConst> const &g)
+     DECL_AND_RETURN(g.data().shape().front_pop());
 
  // ---------------------- implementation --------------------------------
 
@@ -116,15 +142,15 @@ namespace gfs {
  template <typename T> long get_shape(std::vector<T> const &x) { return x.size(); }
 
  /// A common implementation class for gf and gf_view. They will only redefine contructor and = ...
- template <typename Variable, typename Target, typename Opt, bool IsView, bool IsConst>
+ template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsView, bool IsConst>
  class gf_impl : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction) {
   static_assert(!(!IsView && IsConst), "Internal error");
 
   public:
-  using mutable_view_type = gf_view<Variable, Target, Opt>;
-  using const_view_type = gf_const_view<Variable, Target, Opt>;
+  using mutable_view_type = gf_view<Variable, Target, Singularity, Opt>;
+  using const_view_type = gf_const_view<Variable, Target, Singularity, Opt>;
   using view_type = typename std::conditional<IsConst, const_view_type, mutable_view_type>::type;
-  using regular_type = gf<Variable, Target, Opt>;
+  using regular_type = gf<Variable, Target, Singularity, Opt>;
 
   using variable_t = Variable;
   using target_t = Target;
@@ -134,9 +160,9 @@ namespace gfs {
   using domain_t = typename mesh_t::domain_t;
   using mesh_point_t = typename mesh_t::mesh_point_t;
   using mesh_index_t = typename mesh_t::index_t;
-  using symmetry_t = typename gfs_implementation::symmetry<Variable, Target, Opt>::type;
+  using symmetry_t = typename gfs_implementation::symmetry<Variable, Target, Singularity, Opt>::type;
   using indices_t = typename gfs_implementation::indices<Target, Opt>::type;
-  using evaluator_t = gfs_implementation::evaluator<Variable, Target, Opt>;
+  using evaluator_t = gfs_implementation::evaluator<Variable, Target, Singularity, Opt>;
 
   using data_proxy_t = gfs_implementation::data_proxy<Variable, Target, Opt>;
   using data_regular_t = typename data_proxy_t::storage_t;
@@ -144,10 +170,13 @@ namespace gfs {
   using data_const_view_t = typename data_proxy_t::storage_const_view_t;
   using data_t = std14::conditional_t<IsView, std14::conditional_t<IsConst, data_const_view_t, data_view_t>, data_regular_t>;
 
-  using singularity_non_view_t = typename gfs_implementation::singularity<Variable, Target, Opt>::type;
+  //using singularity_non_view_t = typename gfs_implementation::singularity<Variable, Target, Opt>::type;
+  using singularity_non_view_t = Singularity; 
   using singularity_view_t = typename view_type_if_exists_else_type<singularity_non_view_t>::type;
   using singularity_t = std14::conditional_t<IsView, singularity_view_t, singularity_non_view_t>;
 
+  //using is_container_t = void; // a tag to recognize the container
+
   mesh_t const &mesh() const { return _mesh; }
   domain_t const &domain() const { return _mesh.domain(); }
   data_t &data() { return _data; }
@@ -158,8 +187,6 @@ namespace gfs {
   indices_t const &indices() const { return _indices; }
   evaluator_t const &get_evaluator() const { return _evaluator; }
 
-  auto get_data_shape() const DECL_AND_RETURN(get_shape(this -> data()));
-
   protected:
   mesh_t _mesh;
   data_t _data;
@@ -259,8 +286,8 @@ namespace gfs {
 
   // ------------- All the [] operators -----------------------------
   //  [] and access to the grid point
-  using r_type = typename std::result_of<data_proxy_t(data_t &, size_t)>::type;
-  using cr_type = typename std::result_of<data_proxy_t(data_t const &, size_t)>::type;
+  using r_type = std14::result_of_t<data_proxy_t(data_t &, typename mesh_t::linear_index_t)>;
+  using cr_type = std14::result_of_t<data_proxy_t(data_t const &, typename mesh_t::linear_index_t)>;
 
   r_type operator[](mesh_index_t const &arg) { return _data_proxy(_data, _mesh.index_to_linear(arg)); }
   cr_type operator[](mesh_index_t const &arg) const { return _data_proxy(_data, _mesh.index_to_linear(arg)); }
@@ -269,12 +296,12 @@ namespace gfs {
   cr_type operator[](mesh_point_t const &x) const { return _data_proxy(_data, x.linear_index()); }
 
   template <typename... U> r_type operator[](closest_pt_wrap<U...> const &p) {
-   return _data_proxy(_data,
-                      _mesh.index_to_linear(gfs_implementation::get_closest_point<Variable, Target, Opt>::invoke(this, p)));
+   return _data_proxy(
+       _data, _mesh.index_to_linear(gfs_implementation::get_closest_point<Variable, Target, Singularity, Opt>::invoke(this, p)));
   }
   template <typename... U> cr_type operator[](closest_pt_wrap<U...> const &p) const {
-   return _data_proxy(_data,
-                      _mesh.index_to_linear(gfs_implementation::get_closest_point<Variable, Target, Opt>::invoke(this, p)));
+   return _data_proxy(
+       _data, _mesh.index_to_linear(gfs_implementation::get_closest_point<Variable, Target, Singularity, Opt>::invoke(this, p)));
   }
 
   template <typename Arg>
@@ -321,17 +348,17 @@ namespace gfs {
   //----------------------------- HDF5 -----------------------------
 
   friend std::string get_triqs_hdf5_data_scheme(gf_impl const &g) {
-   auto s = gfs_implementation::h5_name<Variable, Target, Opt>::invoke();
+   auto s = gfs_implementation::h5_name<Variable, Target, Singularity, Opt>::invoke();
    return (s == "BlockGf" ? s : "Gf" + s);
   }
 
-  friend struct gfs_implementation::h5_rw<Variable, Target, Opt>;
+  friend struct gfs_implementation::h5_rw<Variable, Target, Singularity, Opt>;
 
   /// Write into HDF5
   friend void h5_write(h5::group fg, std::string subgroup_name, gf_impl const &g) {
    auto gr = fg.create_group(subgroup_name);
    gr.write_triqs_hdf5_data_scheme(g);
-   gfs_implementation::h5_rw<Variable, Target, Opt>::write(gr, g);
+   gfs_implementation::h5_rw<Variable, Target, Singularity, Opt>::write(gr, g);
   }
 
   /// Read from HDF5
@@ -343,7 +370,7 @@ namespace gfs {
    if (tag_file != tag_expected)
     TRIQS_RUNTIME_ERROR << "h5_read : mismatch of the tag TRIQS_HDF5_data_scheme tag in the h5 group : found " << tag_file
                         << " while I expected " << tag_expected;
-   gfs_implementation::h5_rw<Variable, Target, Opt>::read(gr, g);
+   gfs_implementation::h5_rw<Variable, Target, Singularity, Opt>::read(gr, g);
   }
 
   //-----------------------------  BOOST Serialization -----------------------------
@@ -365,8 +392,8 @@ namespace gfs {
  // -------------------------Interaction with the CLEF library : auto assignement implementation -----------------
  // auto assignment of the gf (gf(om_) << expression fills the functions by evaluation of expression)
 
- template <typename RHS, typename Variable, typename Target, typename Opt, bool IsView>
- void triqs_clef_auto_assign(gf_impl<Variable, Target, Opt, IsView, false> &g, RHS const &rhs) {
+ template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt, bool IsView>
+ void triqs_clef_auto_assign(gf_impl<Variable, Target, Singularity, Opt, IsView, false> &g, RHS const &rhs) {
   triqs_clef_auto_assign_impl(g, rhs, typename std::is_base_of<tag::composite, gf_mesh<Variable, Opt>>::type());
   assign_from_expression(g.singularity(), rhs);
   // access to the data . Beware, we view it as a *matrix* NOT an array... (crucial for assignment to scalars !)
@@ -375,8 +402,8 @@ namespace gfs {
  }
 
  // enable the writing g[om_] << .... also
- template <typename RHS, typename Variable, typename Target, typename Opt, bool IsView>
- void triqs_clef_auto_assign_subscript(gf_impl<Variable, Target, Opt, IsView, false> &g, RHS const &rhs) {
+ template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt, bool IsView>
+ void triqs_clef_auto_assign_subscript(gf_impl<Variable, Target, Singularity, Opt, IsView, false> &g, RHS const &rhs) {
   triqs_clef_auto_assign(g, rhs);
  }
 
@@ -389,8 +416,8 @@ namespace gfs {
   triqs_clef_auto_assign(std::forward<G>(g), std::forward<clef::make_fun_impl<Expr, Is...>>(rhs));
  }
 
- template <typename RHS, typename Variable, typename Target, typename Opt, bool IsView>
- void triqs_clef_auto_assign_impl(gf_impl<Variable, Target, Opt, IsView, false> &g, RHS const &rhs,
+ template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt, bool IsView>
+ void triqs_clef_auto_assign_impl(gf_impl<Variable, Target, Singularity, Opt, IsView, false> &g, RHS const &rhs,
                                   std::integral_constant<bool, false>) {
   for (auto const &w : g.mesh()) {
    triqs_gf_clef_auto_assign_impl_aux_assign(g[w], rhs(w));
@@ -398,8 +425,8 @@ namespace gfs {
   }
  }
 
- template <typename RHS, typename Variable, typename Target, typename Opt, bool IsView>
- void triqs_clef_auto_assign_impl(gf_impl<Variable, Target, Opt, IsView, false> &g, RHS const &rhs,
+ template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt, bool IsView>
+ void triqs_clef_auto_assign_impl(gf_impl<Variable, Target, Singularity, Opt, IsView, false> &g, RHS const &rhs,
                                   std::integral_constant<bool, true>) {
   for (auto const &w : g.mesh()) {
    triqs_gf_clef_auto_assign_impl_aux_assign(g[w], triqs::tuple::apply(rhs, w.components_tuple()));
@@ -409,19 +436,19 @@ namespace gfs {
 
  // -------------------------The regular class of GF --------------------------------------------------------
 
- template <typename Variable, typename Target, typename Opt> class gf : public gf_impl<Variable, Target, Opt, false, false> {
-  using B = gf_impl<Variable, Target, Opt, false, false>;
-  using factory = gfs_implementation::factories<Variable, Target, Opt>;
+ template <typename Variable, typename Target, typename Singularity, typename Opt> class gf : public gf_impl<Variable, Target, Singularity, Opt, false, false> {
+  using B = gf_impl<Variable, Target, Singularity, Opt, false, false>;
+  using factory = gfs_implementation::factories<Variable, Target, Singularity, Opt>;
 
   public:
   gf() : B() {}
   gf(gf const &g) : B(g) {}
   gf(gf &&g) noexcept : B(std::move(g)) {}
-  gf(gf_view<Variable, Target, Opt> const &g) : B(g, bool {}) {}
-  gf(gf_const_view<Variable, Target, Opt> const &g) : B(g, bool {}) {}
+  gf(gf_view<Variable, Target, Singularity, Opt> const &g) : B(g, bool {}) {}
+  gf(gf_const_view<Variable, Target, Singularity, Opt> const &g) : B(g, bool {}) {}
 
   template <typename GfType>
-  gf(GfType const &x, typename std::enable_if<ImmutableGreenFunction<GfType>::value>::type *dummy = 0)
+  gf(GfType const &x, std14::enable_if_t<ImmutableGreenFunction<GfType>::value> *dummy = 0)
      : B() {
    *this = x;
   }
@@ -481,27 +508,27 @@ namespace gfs {
  // But in chained clef expression, A(i_)(om_) where A is an array of gf
  // we need to call it with the gf, not gf_impl (or the template resolution find the deleted funciton in clef).
  // Another fix is to make gf, gf_view in the expression tree, but this requires using CRPT in gf_impl...
- template <typename RHS, typename Variable, typename Target, typename Opt>
- void triqs_clef_auto_assign(gf<Variable, Target, Opt> &g, RHS const &rhs) {
-  triqs_clef_auto_assign( static_cast<gf_impl<Variable, Target, Opt, false, false>&>(g), rhs);
+ template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt>
+ void triqs_clef_auto_assign(gf<Variable, Target, Singularity, Opt> &g, RHS const &rhs) {
+  triqs_clef_auto_assign( static_cast<gf_impl<Variable, Target, Singularity, Opt, false, false>&>(g), rhs);
  }
  
  // --------------------------The const View class of GF -------------------------------------------------------
 
- template <typename Variable, typename Target, typename Opt>
- class gf_view<Variable, Target, Opt, true> : public gf_impl<Variable, Target, Opt, true, true> {
-  using B = gf_impl<Variable, Target, Opt, true, true>;
+ template <typename Variable, typename Target, typename Singularity, typename Opt>
+ class gf_view<Variable, Target, Singularity, Opt, true> : public gf_impl<Variable, Target, Singularity, Opt, true, true> {
+  using B = gf_impl<Variable, Target, Singularity, Opt, true, true>;
 
   public:
   gf_view() = delete;
   gf_view(gf_view const &g) : B(g) {}
   gf_view(gf_view &&g) noexcept : B(std::move(g)) {}
 
-  gf_view(gf_impl<Variable, Target, Opt, true, true> const &g) : B(g, bool {}) {}                  // from a const_view
-  gf_view(gf_impl<Variable, Target, Opt, true, false> const &g) : B(g, bool {}) {}                 // from a view
-  gf_view(gf_impl<Variable, Target, Opt, false, false> const &g) : B(g, bool {}) {}                // from a const gf
-  gf_view(gf_impl<Variable, Target, Opt, false, false> &g) : B(g, bool {}) {}                      // from a gf &
-  gf_view(gf_impl<Variable, Target, Opt, false, false> &&g) noexcept : B(std::move(g), bool {}) {} // from a gf &&
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, true, true> const &g) : B(g, bool {}) {}   // from a const_view
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, true, false> const &g) : B(g, bool {}) {}  // from a view
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> const &g) : B(g, bool {}) {} // from a const gf
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> &g) : B(g, bool {}) {}       // from a gf &
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> &&g) noexcept : B(std::move(g), bool {}) {} // from a gf &&
 
   template <typename D>
   gf_view(typename B::mesh_t const &m, D const &dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s,
@@ -517,7 +544,7 @@ namespace gfs {
    this->name = X.name;
   }
 
-  void rebind(gf_view<Variable, Target, Opt, false> const &X) noexcept {
+  void rebind(gf_view<Variable, Target, Singularity, Opt, false> const &X) noexcept {
    rebind(gf_view{X});
   }
 
@@ -526,25 +553,25 @@ namespace gfs {
 
  // -------------------------  The View class of GF  -------------------------------------------------------
 
- template <typename Variable, typename Target, typename Opt>
- class gf_view<Variable, Target, Opt, false> : public gf_impl<Variable, Target, Opt, true, false> {
-  using B = gf_impl<Variable, Target, Opt, true, false>;
+ template <typename Variable, typename Target, typename Singularity, typename Opt>
+ class gf_view<Variable, Target, Singularity, Opt, false> : public gf_impl<Variable, Target, Singularity, Opt, true, false> {
+  using B = gf_impl<Variable, Target, Singularity, Opt, true, false>;
 
   public:
   gf_view() = delete;
   gf_view(gf_view const &g) : B(g) {}
   gf_view(gf_view &&g) noexcept : B(std::move(g)) {}
 
-  gf_view(gf_impl<Variable, Target, Opt, true, true> const &g) = delete;           // from a const view : impossible
-  gf_view(gf_impl<Variable, Target, Opt, true, false> const &g) : B(g, bool {}) {} // from a view
-  gf_view(gf_impl<Variable, Target, Opt, false, false> const &g) = delete;         // from a const gf : impossible
-  gf_view(gf_impl<Variable, Target, Opt, false, false> &g) : B(g, bool {}) {}      // from a gf &
-  gf_view(gf_impl<Variable, Target, Opt, false, false> &&g) noexcept : B(std::move(g), bool {}) {} // from a gf &&
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, true, true> const &g) = delete;           // from a const view : impossible
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, true, false> const &g) : B(g, bool {}) {} // from a view
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> const &g) = delete;         // from a const gf : impossible
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> &g) : B(g, bool {}) {}      // from a gf &
+  gf_view(gf_impl<Variable, Target, Singularity, Opt, false, false> &&g) noexcept : B(std::move(g), bool {}) {} // from a gf &&
 
   template <typename D>
   gf_view(typename B::mesh_t const &m, D &&dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s,
           typename B::indices_t const &ind = typename B::indices_t{}, std::string name = "")
-    : B(m, factory<typename B::data_t>(std::forward<D>(dat)), t, s, ind, name, typename B::evaluator_t{}) {}
+     : B(m, factory<typename B::data_t>(std::forward<D>(dat)), t, s, ind, name, typename B::evaluator_t{}) {}
 
   friend void swap(gf_view &a, gf_view &b) noexcept { a.swap_impl(b); }
 
@@ -569,69 +596,81 @@ namespace gfs {
  }; // class gf_view
 
  // delegate = so that I can overload it for specific RHS...
- template <typename Variable, typename Target, typename Opt, typename RHS>
- DISABLE_IF(arrays::is_scalar<RHS>) triqs_gf_view_assign_delegation(gf_view<Variable, Target, Opt> g, RHS const &rhs) {
+ template <typename Variable, typename Target, typename Singularity, typename Opt, typename RHS>
+ DISABLE_IF(arrays::is_scalar<RHS>) triqs_gf_view_assign_delegation(gf_view<Variable, Target, Singularity, Opt> g, RHS const &rhs) {
   if (!(g.mesh() == rhs.mesh()))
    TRIQS_RUNTIME_ERROR << "Gf Assignment in View : incompatible mesh" << g.mesh() << " vs " << rhs.mesh();
   for (auto const &w : g.mesh()) g[w] = rhs[w];
   g.singularity() = rhs.singularity();
  }
 
- template <typename Variable, typename Target, typename Opt, typename T>
- ENABLE_IF(arrays::is_scalar<T>) triqs_gf_view_assign_delegation(gf_view<Variable, Target, Opt> g, T const &x) {
-  gf_view<Variable, Target, Opt>::data_proxy_t::assign_to_scalar(g.data(), x);
+ template <typename Variable, typename Target, typename Singularity, typename Opt, typename T>
+ ENABLE_IF(arrays::is_scalar<T>) triqs_gf_view_assign_delegation(gf_view<Variable, Target, Singularity, Opt> g, T const &x) {
+  gf_view<Variable, Target, Singularity, Opt>::data_proxy_t::assign_to_scalar(g.data(), x);
   g.singularity() = x;
  }
 
  // tool for lazy transformation
- template <typename Tag, typename D, typename Target = matrix_valued, typename Opt = void> struct gf_keeper {
-  gf_const_view<D, Target, Opt> g;
+ // TODO  : clean this : typename G
+ template <typename Tag, typename D, typename Target = matrix_valued,
+           typename Singularity = gf_default_singularity_t<D, Target>, typename Opt = void>
+ struct gf_keeper {
+  gf_const_view<D, Target, Singularity, Opt> g;
  };
 
  // Cf gf
- template <typename RHS, typename Variable, typename Target, typename Opt>
- void triqs_clef_auto_assign(gf_view<Variable, Target, Opt> &g, RHS const &rhs) {
-  triqs_clef_auto_assign( static_cast<gf_impl<Variable, Target, Opt, true, false>&>(g), rhs);
+ template <typename RHS, typename Variable, typename Target, typename Singularity, typename Opt>
+ void triqs_clef_auto_assign(gf_view<Variable, Target, Singularity, Opt> &g, RHS const &rhs) {
+  triqs_clef_auto_assign( static_cast<gf_impl<Variable, Target, Singularity, Opt, true, false>&>(g), rhs);
  }
 
  // ---------------------------------- slicing ------------------------------------
 
  // slice
- template <typename Variable, typename Target, typename Opt, bool IsConst, typename... Args>
- gf_view<Variable, matrix_valued, Opt, IsConst> slice_target(gf_view<Variable, Target, Opt, IsConst> g, Args &&... args) {
+ template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst, typename... Args>
+ gf_view<Variable, matrix_valued, Singularity, Opt, IsConst> slice_target(gf_view<Variable, Target, Singularity, Opt, IsConst> g,
+                                                                          Args &&... args) {
   static_assert(std::is_same<Target, matrix_valued>::value, "slice_target only for matrix_valued GF's");
   using arrays::range;
   return {g.mesh(),                                                   g.data()(range(), std::forward<Args>(args)...),
-          slice_target(g.singularity(), std::forward<Args>(args)...), g.symmetry(), slice(g.indices(),std::forward<Args>(args)...), g.name };
+          slice_target(g.singularity(), std::forward<Args>(args)...), g.symmetry(),
+          slice(g.indices(), std::forward<Args>(args)...),            g.name};
  }
 
- template <typename Variable, typename Target, typename Opt, typename... Args>
- gf_view<Variable, matrix_valued, Opt> slice_target(gf<Variable, Target, Opt> &g, Args &&... args) {
+ template <typename Variable, typename Target, typename Singularity, typename Opt, typename... Args>
+ gf_view<Variable, matrix_valued, Singularity, Opt> slice_target(gf<Variable, Target, Singularity, Opt> &g, Args &&... args) {
   return slice_target(g(), std::forward<Args>(args)...);
  }
 
- template <typename Variable, typename Target, typename Opt, typename... Args>
- gf_const_view<Variable, matrix_valued, Opt> slice_target(gf<Variable, Target, Opt> const &g, Args &&... args) {
+ template <typename Variable, typename Target, typename Singularity, typename Opt, typename... Args>
+ gf_const_view<Variable, matrix_valued, Singularity, Opt> slice_target(gf<Variable, Target, Singularity, Opt> const &g,
+                                                                       Args &&... args) {
   return slice_target(g(), std::forward<Args>(args)...);
  }
 
  // slice to scalar
- template <typename Variable, typename Target, typename Opt, bool IsConst, typename... Args>
- gf_view<Variable, scalar_valued, Opt, IsConst> slice_target_to_scalar(gf_view<Variable, Target, Opt, IsConst> g,
-                                                                       Args &&... args) {
+ template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst, typename... Args>
+ gf_view<Variable, scalar_valued, Singularity, Opt, IsConst>
+ slice_target_to_scalar(gf_view<Variable, Target, Singularity, Opt, IsConst> g, Args &&... args) {
   static_assert(std::is_same<Target, matrix_valued>::value, "slice_target only for matrix_valued GF's");
   using arrays::range;
-  return {g.mesh(),                                                g.data()(range(), std::forward<Args>(args)...),
-          slice_target(g.singularity(), range(args, args + 1)...), g.symmetry(), {}, g.name};
+  return {g.mesh(),
+          g.data()(range(), std::forward<Args>(args)...),
+          slice_target(g.singularity(), range(args, args + 1)...),
+          g.symmetry(),
+          {},
+          g.name};
  }
 
- template <typename Variable, typename Target, typename Opt, typename... Args>
- gf_view<Variable, scalar_valued, Opt> slice_target_to_scalar(gf<Variable, Target, Opt> &g, Args &&... args) {
+ template <typename Variable, typename Target, typename Singularity, typename Opt, typename... Args>
+ gf_view<Variable, scalar_valued, Singularity, Opt> slice_target_to_scalar(gf<Variable, Target, Singularity, Opt> &g,
+                                                                           Args &&... args) {
   return slice_target_to_scalar(g(), std::forward<Args>(args)...);
  }
 
- template <typename Variable, typename Target, typename Opt, typename... Args>
- gf_const_view<Variable, scalar_valued, Opt> slice_target_to_scalar(gf<Variable, Target, Opt> const &g, Args &&... args) {
+ template <typename Variable, typename Target, typename Singularity, typename Opt, typename... Args>
+ gf_const_view<Variable, scalar_valued, Singularity, Opt> slice_target_to_scalar(gf<Variable, Target, Singularity, Opt> const &g,
+                                                                                 Args &&... args) {
   return slice_target_to_scalar(g(), std::forward<Args>(args)...);
  }
 
@@ -645,72 +684,69 @@ namespace gfs {
 */
 
  // ---------------------------------- target reinterpretation ------------------------------------
- 
+
  // a scalar_valued gf can be viewed as a 1x1 matrix
- template <typename Variable, typename Opt, bool IsConst>
- gf_view<Variable, matrix_valued, Opt, IsConst>
- reinterpret_scalar_valued_gf_as_matrix_valued(gf_view<Variable, scalar_valued, Opt, IsConst> g) {
-  using a_t = typename gf_view<Variable, matrix_valued, Opt, IsConst>::data_view_t;
+ template <typename Variable, typename Singularity, typename Opt, bool IsConst>
+ gf_view<Variable, matrix_valued, Singularity, Opt, IsConst>
+ reinterpret_scalar_valued_gf_as_matrix_valued(gf_view<Variable, scalar_valued, Singularity, Opt, IsConst> g) {
+  using a_t = typename gf_view<Variable, matrix_valued, Singularity, Opt, IsConst>::data_view_t;
   auto a = a_t{typename a_t::indexmap_type(join(g.data().shape(), make_shape(1, 1))), g.data().storage()};
-  return {g.mesh(), a, g.singularity(), g.symmetry(), {} , g.name};
+  return {g.mesh(), a, g.singularity(), g.symmetry(), {}, g.name};
  }
 
- template <typename Variable, typename Opt>
- gf_view<Variable, matrix_valued, Opt> reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Opt> &g) {
+ template <typename Variable, typename Singularity, typename Opt>
+ gf_view<Variable, matrix_valued, Singularity, Opt>
+ reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Singularity, Opt> &g) {
   return reinterpret_scalar_valued_gf_as_matrix_valued(g());
  }
 
- template <typename Variable, typename Opt>
- gf_const_view<Variable, matrix_valued, Opt>
- reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Opt> const &g) {
+ template <typename Variable, typename Singularity, typename Opt>
+ gf_const_view<Variable, matrix_valued, Singularity, Opt>
+ reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Singularity, Opt> const &g) {
   return reinterpret_scalar_valued_gf_as_matrix_valued(g());
  }
 
- /*
-    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());
-    }*/
-
  // ---------------------------------- some functions : invert, conjugate, transpose, ... ------------------------------------
 
- //  ---- inversion 
+ //  ---- inversion
  // auxiliary function : invert the data : one function for all matrix valued gf (save code).
- template <typename A3> void _gf_invert_data_in_place(A3 && a) {
-  for (int i = 0; i < first_dim(a); ++i) {// Rely on the ordering 
+ template <typename A3> void _gf_invert_data_in_place(A3 &&a) {
+  for (int i = 0; i < first_dim(a); ++i) { // Rely on the ordering
    auto v = make_matrix_view(a(i, arrays::range(), arrays::range()));
    v = triqs::arrays::inverse(v);
   }
  }
 
- template <typename Variable, typename Opt>
- void invert_in_place(gf_view<Variable, matrix_valued, Opt> g) {
+ template <typename Variable, typename Singularity, typename Opt>
+ void invert_in_place(gf_view<Variable, matrix_valued, Singularity, Opt> g) {
   _gf_invert_data_in_place(g.data());
-  g.singularity() = inverse(g.singularity()); 
+  g.singularity() = inverse(g.singularity());
  }
 
- template <typename Variable, typename Opt> gf<Variable, matrix_valued, Opt> inverse(gf<Variable, matrix_valued, Opt> const & g) {
+ template <typename Variable, typename Singularity, typename Opt>
+ gf<Variable, matrix_valued, Singularity, Opt> inverse(gf<Variable, matrix_valued, Singularity, Opt> const &g) {
   auto res = g;
-  gf_view<Variable, matrix_valued, Opt> v = res;
+  gf_view<Variable, matrix_valued, Singularity, Opt> v = res;
   invert_in_place(v);
   return res;
  }
 
- template <typename Variable, typename Opt, bool B> gf<Variable, matrix_valued, Opt> inverse(gf_view<Variable, matrix_valued, Opt, B> g) {
-  return inverse(gf<Variable, matrix_valued, Opt>(g));
+ template <typename Variable, typename Singularity, typename Opt, bool B>
+ gf<Variable, matrix_valued, Singularity, Opt> inverse(gf_view<Variable, matrix_valued, Singularity, Opt, B> g) {
+  return inverse(gf<Variable, matrix_valued, Singularity, Opt>(g));
  }
 
  //  ---- transpose : a new gf
 
- template <typename Variable, typename Opt>
- gf<Variable, matrix_valued, Opt> transpose(gf_view<Variable, matrix_valued, Opt> g) {
+ template <typename Variable, typename Singularity, typename Opt>
+ gf<Variable, matrix_valued, Singularity, Opt> transpose(gf_view<Variable, matrix_valued, Singularity, Opt> g) {
   return {g.mesh(), transposed_view(g.data(), 0, 2, 1), transpose(g.singularity()), g.symmetry(), transpose(g.indices()), g.name};
  }
- 
- //  ---- conjugate : always a new function -> changelog 
 
- template <typename Variable, typename Opt> gf<Variable, matrix_valued, Opt> conj(gf_view<Variable, matrix_valued, Opt> g) {
+ //  ---- conjugate : always a new function -> changelog
+
+ template <typename Variable, typename Singularity, typename Opt>
+ gf<Variable, matrix_valued, Singularity, Opt> conj(gf_view<Variable, matrix_valued, Singularity, Opt> g) {
   return {g.mesh(), conj(g.data()), conj(g.singularity()), g.symmetry(), g.indices(), g.name};
  }
 
@@ -725,7 +761,7 @@ namespace gfs {
 
  template <typename A3, typename T> void _gf_data_mul_L(matrix<T> const &l, A3 &&a) {
   for (int i = 0; i < first_dim(a); ++i) { // Rely on the ordering
-   matrix_view<T>  v = a(i, arrays::range(), arrays::range());
+   matrix_view<T> v = a(i, arrays::range(), arrays::range());
    v = l * v;
   }
  }
@@ -734,140 +770,108 @@ namespace gfs {
   auto res = typename A3::regular_type(first_dim(a), first_dim(l), second_dim(r));
   for (int i = 0; i < first_dim(a); ++i) { // Rely on the ordering
    auto _ = arrays::range{};
-   res(i, _, _) = l * make_matrix_view(a(i, _, _))* r;
+   res(i, _, _) = l * make_matrix_view(a(i, _, _)) * r;
   }
   return res;
  }
 
- template <typename Variable, typename Opt, typename T>
- gf<Variable, matrix_valued, Opt> operator*(gf<Variable, matrix_valued, Opt> g, matrix<T> r) {
+ template <typename Variable, typename Singularity, typename Opt, typename T>
+ gf<Variable, matrix_valued, Singularity, Opt> operator*(gf<Variable, matrix_valued, Singularity, Opt> g, matrix<T> r) {
   _gf_data_mul_R(g.data(), r);
   g.singularity() = g.singularity() * r;
   return g;
  }
 
- template <typename Variable, typename Opt, typename T>
- gf<Variable, matrix_valued, Opt> operator*(matrix<T> l, gf<Variable, matrix_valued, Opt> g) {
+ template <typename Variable, typename Singularity, typename Opt, typename T>
+ gf<Variable, matrix_valued, Singularity, Opt> operator*(matrix<T> l, gf<Variable, matrix_valued, Singularity, Opt> g) {
   _gf_data_mul_L(l, g.data());
   g.singularity() = l * g.singularity();
   return g;
  }
 
- template <typename Variable, typename Opt, typename T>
- gf<Variable, matrix_valued, Opt> L_G_R(matrix<T> l, gf<Variable, matrix_valued, Opt> g, matrix<T> r) {
-  auto res = gf<Variable, matrix_valued, Opt>{g.mesh(), {int(first_dim(l)), int(second_dim(r))}};
+ template <typename Variable, typename Singularity, typename Opt, typename T>
+ gf<Variable, matrix_valued, Singularity, Opt> L_G_R(matrix<T> l, gf<Variable, matrix_valued, Singularity, Opt> g, matrix<T> r) {
+  auto res = gf<Variable, matrix_valued, Singularity, Opt>{g.mesh(), {int(first_dim(l)), int(second_dim(r))}};
   res.data() = _gf_data_mul_LR(l, g.data(), r);
   res.singularity() = l * g.singularity() * r;
   return res;
  }
 
-
  namespace gfs_implementation { // implement some default traits
 
   // -------------------------  default factories ---------------------
 
-  // ----- tensor_valued
-  template <int R, typename Var, typename Opt> struct factories<Var, tensor_valued<R>, Opt> {
-   using gf_t = gf<Var, tensor_valued<R>, Opt>;
-   using target_shape_t = arrays::mini_vector<int, R>;
+  template <typename Var, typename Target, int VarDim, typename Singularity, typename Opt> struct factories_default_impl {
+   using gf_t = gf<Var, Target, Singularity, Opt>;
+   using target_shape_t = arrays::mini_vector<int, VarDim>;
    using mesh_t = typename gf_t::mesh_t;
-
+   //
    static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) {
-    typename gf_t::data_t A(shape.front_append(m.size()));
+    typename gf_t::data_t A(gf_t::data_proxy_t::join_shape(m.size_of_components(), shape));
     A() = 0;
     return A;
    }
-
+   //
    static typename gf_t::singularity_t make_singularity(mesh_t const &m, target_shape_t shape) {
-    return typename gf_t::singularity_t(shape);
+    return typename gf_t::singularity_t{shape};
    }
   };
 
-  // ----- matrix_valued
-  template <typename Var, typename Opt> struct factories<Var, matrix_valued, Opt> {
-   using gf_t = gf<Var, matrix_valued, Opt>;
-   using target_shape_t = arrays::mini_vector<int, 2>;
-   using mesh_t = typename gf_t::mesh_t;
+  template <int R, typename Var, typename Singularity, typename Opt>
+  struct factories<Var, tensor_valued<R>, Singularity, Opt> : factories_default_impl<Var, tensor_valued<R>, R, Singularity, Opt> {};
 
-   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) {
-    typename gf_t::data_t A(shape.front_append(m.size()));
-    A() = 0;
-    return A;
-   }
+  template <typename Var, typename Singularity, typename Opt>
+  struct factories<Var, matrix_valued, Singularity, Opt> : factories_default_impl<Var, matrix_valued, 2, Singularity, Opt> {};
 
-   static typename gf_t::singularity_t make_singularity(mesh_t const &m, target_shape_t shape) {
-    return typename gf_t::singularity_t(shape);
-   }
-  };
-
-  // ----- scalar_valued
-  template <typename Var, typename Opt> struct factories<Var, scalar_valued, Opt> {
-   using gf_t = gf<Var, scalar_valued, Opt>;
-   struct target_shape_t {
-    target_shape_t front_pop() const { // this make the get_target_shape function works in this case...
-     return {};
-    }
-    target_shape_t() = default;
-    template <typename T> target_shape_t(utility::mini_vector<T, 0>) {}
-   };
-
-   using mesh_t = typename gf_t::mesh_t;
-
-   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t shape) {
-    typename gf_t::data_t A(m.size());
-    A() = 0;
-    return A;
-   }
-   static typename gf_t::singularity_t make_singularity(mesh_t const &m, target_shape_t shape) {
-    return typename gf_t::singularity_t{1, 1};
-   }
-  };
+  template <typename Var, typename Singularity, typename Opt>
+  struct factories<Var, scalar_valued, Singularity, Opt> : factories_default_impl<Var, scalar_valued, 0, Singularity, Opt> {};
 
   // ---------------------   hdf5 ---------------------------------------
   // scalar function : just add a _s
-  template <typename Variable, typename Opt> struct h5_name<Variable, scalar_valued, Opt> {
-   static std::string invoke() { return h5_name<Variable, matrix_valued, Opt>::invoke() + "_s"; }
+  template <typename Variable, typename Singularity, typename Opt> struct h5_name<Variable, scalar_valued, Singularity, Opt> {
+   static std::string invoke() { return h5_name<Variable, matrix_valued, Singularity, Opt>::invoke() + "_s"; }
   };
 
   // the h5 write and read of gf members, so that we can specialize it e.g. for block gf
-  template <typename Variable, typename Target, typename Opt> struct h5_rw {
+  template <typename Variable, typename Target, typename Singularity, typename Opt> struct h5_rw {
 
-   static void write(h5::group gr, gf_const_view<Variable, Target, Opt> g) {
+   static void write(h5::group gr, gf_const_view<Variable, Target, Singularity, Opt> g) {
     h5_write(gr, "data", g._data);
     h5_write(gr, "singularity", g._singularity);
     h5_write(gr, "mesh", g._mesh);
     h5_write(gr, "symmetry", g._symmetry);
     h5_write(gr, "indices", g._indices);
-    //h5_write(gr, "name", g.name);
+    // h5_write(gr, "name", g.name);
    }
 
-   template <bool B> static void read(h5::group gr, gf_impl<Variable, Target, Opt, B, false> &g) {
+   template <bool B> static void read(h5::group gr, gf_impl<Variable, Target, Singularity, Opt, B, false> &g) {
     h5_read(gr, "data", g._data);
     h5_read(gr, "singularity", g._singularity);
     h5_read(gr, "mesh", g._mesh);
     h5_read(gr, "symmetry", g._symmetry);
     h5_read(gr, "indices", g._indices);
-    //h5_read(gr, "name", g.name);
+    // h5_read(gr, "name", g.name);
    }
   };
  } // gfs_implementation
 }
 
-namespace mpi { 
+namespace mpi {
 
-  template <typename Variable, typename Target, typename Opt>
- struct mpi_impl<gfs::gf<Variable, Target, Opt>, void> : mpi_impl_triqs_gfs<gfs::gf<Variable, Target, Opt>> {};
-
- template <typename Variable, typename Target, typename Opt, bool IsConst>
- struct mpi_impl<gfs::gf_view<Variable, Target, Opt, IsConst>, void> : mpi_impl_triqs_gfs<gfs::gf_view<Variable, Target, Opt, IsConst>> {};
+ template <typename Variable, typename Target, typename Singularity, typename Opt>
+ struct mpi_impl<gfs::gf<Variable, Target, Singularity, Opt>,
+                 void> : mpi_impl_triqs_gfs<gfs::gf<Variable, Target, Singularity, Opt>> {};
 
+ template <typename Variable, typename Target, typename Singularity, typename Opt, bool IsConst>
+ struct mpi_impl<gfs::gf_view<Variable, Target, Singularity, Opt, IsConst>,
+                 void> : mpi_impl_triqs_gfs<gfs::gf_view<Variable, Target, Singularity, Opt, IsConst>> {};
 }
-
 }
 
 // same as for arrays : views cannot be swapped by the std::swap. Delete it
 namespace std {
-template <typename Variable, typename Target, typename Opt, bool C1, bool C2>
-void swap(triqs::gfs::gf_view<Variable, Target, Opt, C1> &a, triqs::gfs::gf_view<Variable, Target, Opt, C2> &b) = delete;
+template <typename Variable, typename Target, typename Singularity, typename Opt, bool C1, bool C2>
+void swap(triqs::gfs::gf_view<Variable, Target, Singularity, Opt, C1> &a,
+          triqs::gfs::gf_view<Variable, Target, Singularity, Opt, C2> &b) = delete;
 }
 #include "./gf_expr.hpp"
diff --git a/triqs/gfs/gf_expr.hpp b/triqs/gfs/gf_expr.hpp
index 5f659cfc..d1f194e0 100644
--- a/triqs/gfs/gf_expr.hpp
+++ b/triqs/gfs/gf_expr.hpp
@@ -28,100 +28,114 @@ namespace triqs { namespace gfs {
 
  namespace gfs_expr_tools {
 
-  // a wrapper for scalars 
-  template<typename S> struct scalar_wrap {
-   typedef void variable_t;
-   typedef void target_t;
-   typedef void option_t;
-   S s; 
-   template<typename T> scalar_wrap(T && x):s(std::forward<T>(x)){}  
-   S singularity() const { return s;}
-   template<typename KeyType> S operator[](KeyType && key) const { return s;}
-   template<typename ... Args> inline S operator()(Args && ... args) const { return s;}
-   friend std::ostream &operator <<(std::ostream &sout, scalar_wrap const &expr){return sout << expr.s; }
+  // a wrapper for scalars
+  template <typename S> struct scalar_wrap {
+   using variable_t = void;
+   using target_t = void;
+   using option_t = void;
+   S s;
+   template <typename T> scalar_wrap(T &&x) : s(std::forward<T>(x)) {}
+   S singularity() const { return s; }
+   template <typename KeyType> S operator[](KeyType &&key) const { return s; }
+   template <typename... Args> inline S operator()(Args &&... args) const { return s; }
+   friend std::ostream &operator<<(std::ostream &sout, scalar_wrap const &expr) { return sout << expr.s; }
   };
 
   // Combine the two meshes of LHS and RHS : need to specialize where there is a scalar
   struct combine_mesh {
-   template<typename L, typename R> 
-    auto operator() (L && l, R && r) const -> decltype(std::forward<L>(l).mesh()) { 
-     if (!(l.mesh() == r.mesh())) TRIQS_RUNTIME_ERROR << "Mesh mismatch : in Green Function Expression "<< l.mesh()<<" vs" <<r.mesh();
-     return std::forward<L>(l).mesh();
-    }
-   template<typename S, typename R> auto operator() (scalar_wrap<S> const &, R && r) const DECL_AND_RETURN(std::forward<R>(r).mesh());
-   template<typename S, typename L> auto operator() (L && l, scalar_wrap<S> const &) const DECL_AND_RETURN(std::forward<L>(l).mesh());
+   template <typename L, typename R> auto operator()(L &&l, R &&r) const -> decltype(std::forward<L>(l).mesh()) {
+    if (!(l.mesh() == r.mesh()))
+     TRIQS_RUNTIME_ERROR << "Mesh mismatch : in Green Function Expression " << l.mesh() << " vs" << r.mesh();
+    return std::forward<L>(l).mesh();
+   }
+   template <typename S, typename R>
+   auto operator()(scalar_wrap<S> const &, R &&r) const DECL_AND_RETURN(std::forward<R>(r).mesh());
+   template <typename S, typename L>
+   auto operator()(L &&l, scalar_wrap<S> const &) const DECL_AND_RETURN(std::forward<L>(l).mesh());
   };
 
   // Same thing to get the data shape
   // NB : could be unified to one combine<F>, where F is a functor, but an easy usage requires polymorphic lambda ...
   struct combine_shape {
-   template<typename L, typename R> 
-    auto operator() (L && l, R && r) const -> decltype(get_gf_data_shape(std::forward<L>(l))) { 
-     if (!(get_gf_data_shape(l) == get_gf_data_shape(r))) 
-      TRIQS_RUNTIME_ERROR << "Shape mismatch in Green Function Expression: " << get_gf_data_shape(l) << " vs "<<  get_gf_data_shape(r);
-     return get_gf_data_shape(std::forward<L>(l));
-    }
-   template<typename S, typename R> auto operator() (scalar_wrap<S> const &, R && r) const DECL_AND_RETURN(get_gf_data_shape(std::forward<R>(r)));
-   template<typename S, typename L> auto operator() (L && l, scalar_wrap<S> const &) const DECL_AND_RETURN(get_gf_data_shape(std::forward<L>(l)));
+   template <typename L, typename R> auto operator()(L &&l, R &&r) const -> decltype(get_gf_data_shape(std::forward<L>(l))) {
+    if (!(get_gf_data_shape(l) == get_gf_data_shape(r)))
+     TRIQS_RUNTIME_ERROR << "Shape mismatch in Green Function Expression: " << get_gf_data_shape(l) << " vs "
+                         << get_gf_data_shape(r);
+    return get_gf_data_shape(std::forward<L>(l));
+   }
+   template <typename S, typename R>
+   auto operator()(scalar_wrap<S> const &, R &&r) const DECL_AND_RETURN(get_gf_data_shape(std::forward<R>(r)));
+   template <typename S, typename L>
+   auto operator()(L &&l, scalar_wrap<S> const &) const DECL_AND_RETURN(get_gf_data_shape(std::forward<L>(l)));
   };
 
-  template<typename T> struct node_t : std::conditional<utility::is_in_ZRC<T>::value, scalar_wrap<T>, typename remove_rvalue_ref<T>::type> {};
+  template <typename T> using node_t = std14::conditional_t<utility::is_in_ZRC<T>::value, scalar_wrap<T>, typename remove_rvalue_ref<T>::type>;
 
-  template <typename A, typename B> struct _or_ {typedef void type;};
-  template <typename A> struct _or_<A,A>        {typedef A type;};
-  template <typename A> struct _or_<void,A>     {typedef A type;};
-  template <typename A> struct _or_<A,void>     {typedef A type;};
-  template <> struct _or_<void,void>            {typedef void type;};
+  template <typename A, typename B> struct _or_ {using type=void ;};
+  template <typename A> struct _or_<A,A>        {using type=A ;};
+  template <typename A> struct _or_<void,A>     {using type=A ;};
+  template <typename A> struct _or_<A,void>     {using type=A ;};
+  template <> struct _or_<void,void>            {using type=void ;};
 
  }// gfs_expr_tools
 
- template<typename Tag, typename L, typename R> struct gf_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction){
-  typedef typename std::remove_reference<L>::type L_t;
-  typedef typename std::remove_reference<R>::type R_t;
-  typedef typename gfs_expr_tools::_or_<typename L_t::variable_t,typename R_t::variable_t>::type variable_t;
-  typedef typename gfs_expr_tools::_or_<typename L_t::target_t,typename R_t::target_t>::type target_t;
-  typedef typename gfs_expr_tools::_or_<typename L_t::option_t,typename R_t::option_t>::type option_t;
-  static_assert(!std::is_same<variable_t,void>::value, "Cannot combine two gf expressions with different variables");
-  static_assert(!std::is_same<target_t,void>::value, "Cannot combine two gf expressions with different target");
- 
-  L l; R r;
-  template<typename LL, typename RR> gf_expr(LL && l_, RR && r_):l(std::forward<LL>(l_)), r(std::forward<RR>(r_)) {}
-  
-  auto mesh() const DECL_AND_RETURN(gfs_expr_tools::combine_mesh()(l,r)); 
-  auto singularity() const DECL_AND_RETURN (utility::operation<Tag>()(l.singularity() , r.singularity()));
-  auto get_data_shape() const DECL_AND_RETURN (gfs_expr_tools::combine_shape()(l,r)); 
+ template <typename Tag, typename L, typename R> struct gf_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction) {
+  using L_t = typename std::remove_reference<L>::type;
+  using R_t = typename std::remove_reference<R>::type;
+  using variable_t = typename gfs_expr_tools::_or_<typename L_t::variable_t, typename R_t::variable_t>::type;
+  using target_t = typename gfs_expr_tools::_or_<typename L_t::target_t, typename R_t::target_t>::type;
+  using option_t = typename gfs_expr_tools::_or_<typename L_t::option_t, typename R_t::option_t>::type;
+  static_assert(!std::is_same<variable_t, void>::value, "Cannot combine two gf expressions with different variables");
+  static_assert(!std::is_same<target_t, void>::value, "Cannot combine two gf expressions with different target");
 
-  template<typename KeyType>  auto operator[](KeyType && key)   const DECL_AND_RETURN(utility::operation<Tag>()(l[std::forward<KeyType>(key)] , r[std::forward<KeyType>(key)]));
-  template<typename ... Args> auto operator()(Args && ... args) const DECL_AND_RETURN(utility::operation<Tag>()(l(std::forward<Args>(args)...) , r(std::forward<Args>(args)...)));
-  friend std::ostream &operator <<(std::ostream &sout, gf_expr const &expr){return sout << "("<<expr.l << " "<<utility::operation<Tag>::name << " "<<expr.r<<")" ; }
+  L l;
+  R r;
+  template <typename LL, typename RR> gf_expr(LL &&l_, RR &&r_) : l(std::forward<LL>(l_)), r(std::forward<RR>(r_)) {}
+
+  auto mesh() const DECL_AND_RETURN(gfs_expr_tools::combine_mesh()(l, r));
+  auto singularity() const DECL_AND_RETURN(utility::operation<Tag>()(l.singularity(), r.singularity()));
+
+  template <typename KeyType>
+  auto operator[](KeyType &&key) const
+      DECL_AND_RETURN(utility::operation<Tag>()(l[std::forward<KeyType>(key)], r[std::forward<KeyType>(key)]));
+  template <typename... Args>
+  auto operator()(Args &&... args) const
+      DECL_AND_RETURN(utility::operation<Tag>()(l(std::forward<Args>(args)...), r(std::forward<Args>(args)...)));
+  friend std::ostream &operator<<(std::ostream &sout, gf_expr const &expr) {
+   return sout << "(" << expr.l << " " << utility::operation<Tag>::name << " " << expr.r << ")";
+  }
  };
 
- // -------------------------------------------------------------------
- //a special case : the unary operator !
- template<typename L> struct gf_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction){
-  typedef typename std::remove_reference<L>::type L_t;
-  typedef typename L_t::variable_t variable_t;
-  typedef typename L_t::target_t target_t;
-  typedef typename L_t::option_t option_t;
+ template <typename Tag, typename L, typename R>
+ auto get_gf_data_shape(gf_expr<Tag, L, R> const &g) RETURN(gfs_expr_tools::combine_shape()(g.l, g.r));
 
-  L l; 
+ // -------------------------------------------------------------------
+ // a special case : the unary operator !
+ template <typename L> struct gf_unary_m_expr : TRIQS_CONCEPT_TAG_NAME(ImmutableGreenFunction) {
+  using L_t = typename std::remove_reference<L>::type;
+  using variable_t = typename L_t::variable_t;
+  using target_t = typename L_t::target_t;
+  using option_t = typename L_t::option_t;
+
+  L l;
   template<typename LL> gf_unary_m_expr(LL && l_) : l(std::forward<LL>(l_)) {}
   
   auto mesh() const DECL_AND_RETURN(l.mesh()); 
   auto singularity() const DECL_AND_RETURN(l.singularity());
-  AUTO_DECL get_data_shape() const RETURN (get_gf_data_shape(l));
   
   template<typename KeyType> auto operator[](KeyType&& key) const DECL_AND_RETURN( -l[key]); 
   template<typename ... Args> auto operator()(Args && ... args) const DECL_AND_RETURN( -l(std::forward<Args>(args)...));
   friend std::ostream &operator <<(std::ostream &sout, gf_unary_m_expr const &expr){return sout << '-'<<expr.l; }
  };
 
+ template <typename L> AUTO_DECL get_gf_data_shape(gf_unary_m_expr<L> const &g) RETURN(get_gf_data_shape(g.l));
+
  // -------------------------------------------------------------------
  // Now we can define all the C++ operators ...
 #define DEFINE_OPERATOR(TAG, OP, TRAIT1, TRAIT2) \
  template<typename A1, typename A2>\
- typename std::enable_if<TRAIT1<A1>::value && TRAIT2 <A2>::value, \
- gf_expr<utility::tags::TAG, typename gfs_expr_tools::node_t<A1>::type, typename gfs_expr_tools::node_t<A2>::type>>::type\
+ std14::enable_if_t<TRAIT1<A1>::value && TRAIT2 <A2>::value, \
+ gf_expr<utility::tags::TAG, gfs_expr_tools::node_t<A1>, gfs_expr_tools::node_t<A2>>>\
  operator OP (A1 && a1, A2 && a2) { return {std::forward<A1>(a1),std::forward<A2>(a2)};} 
 
  DEFINE_OPERATOR(plus,       +, ImmutableGreenFunction,ImmutableGreenFunction);
@@ -135,12 +149,10 @@ namespace triqs { namespace gfs {
 #undef DEFINE_OPERATOR
 
  // the unary is special
- template<typename A1> 
-  typename std::enable_if< 
-  ImmutableGreenFunction<A1>::value, 
-  gf_unary_m_expr<typename gfs_expr_tools::node_t<A1>::type >
-   >::type
-   operator - (A1 && a1) { return {std::forward<A1>(a1)};} 
+ template <typename A1>
+ std14::enable_if_t<ImmutableGreenFunction<A1>::value, gf_unary_m_expr<gfs_expr_tools::node_t<A1>>> operator-(A1 &&a1) {
+  return {std::forward<A1>(a1)};
+ }
 
 // Now the inplace operator. Because of expression template, there are useless for speed
 // we implement them trivially.
diff --git a/triqs/gfs/imfreq.hpp b/triqs/gfs/imfreq.hpp
index ac175997..63480ec9 100644
--- a/triqs/gfs/imfreq.hpp
+++ b/triqs/gfs/imfreq.hpp
@@ -32,24 +32,30 @@ namespace gfs {
 
  template <typename Opt> struct gf_mesh<imfreq, Opt> : matsubara_freq_mesh {
   template <typename... T> gf_mesh(T &&... x) : matsubara_freq_mesh(std::forward<T>(x)...) {}
-  //using matsubara_freq_mesh::matsubara_freq_mesh;
+  // using matsubara_freq_mesh::matsubara_freq_mesh;
+ };
+
+ // singularity
+ template <> struct gf_default_singularity<imfreq, matrix_valued> {
+  using type = local::tail;
+ };
+ template <> struct gf_default_singularity<imfreq, scalar_valued> {
+  using type = local::tail;
  };
 
  namespace gfs_implementation {
 
-  // singularity
-  template <> struct singularity<imfreq, matrix_valued, void> {
-   using type = local::tail;
-  };
-  template <> struct singularity<imfreq, scalar_valued, void> {
-   using type = local::tail;
-  };
-
-  // h5 name
-  template <typename Opt> struct h5_name<imfreq, matrix_valued, Opt> {
+  /// ---------------------------  hdf5 ---------------------------------
+  
+  template <typename S, typename Opt> struct h5_name<imfreq, matrix_valued, S, Opt> {
    static std::string invoke() { return "ImFreq"; }
   };
 
+  /// ---------------------------  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> {};
+
   /// ---------------------------  evaluator ---------------------------------
 
   // simple evaluation : take the point on the grid...
@@ -67,103 +73,89 @@ namespace gfs {
 
   // ------------- evaluator  -------------------
   // handle the case where the matsu. freq is out of grid...
-  template <typename Target, typename Opt> struct evaluator<imfreq, Target, Opt> {
-   static constexpr int arity = 1;
 
-   private:
+  struct _eval_imfreq_base_impl {
+   static constexpr int arity = 1;
    template <typename G> int sh(G const * g) const { return (g->mesh().domain().statistic == Fermion ? 1 : 0);}
 
-   // dispatch for 2x2 cases : matrix/scalar and tail/no_tail ( true means no_tail)
-   template <typename G>
-   std::complex<double> _call_impl(G const *g, matsubara_freq const &f, scalar_valued, std::false_type) const {
-    if (g->mesh().positive_only()){//only positive Matsubara frequencies
-     if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n];
-     if ((f.n < 0) && ((-f.n-sh(g)) < g->mesh().size())) return conj((*g)[-f.n-sh(g)]);
-    }
-    else{
-     if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size()+g->mesh().first_index())) return (*g)[f.n];
-    }
-    return g->singularity().evaluate(f)(0, 0);
-   }
-
-   template <typename G>
-   std::complex<double> _call_impl(G const *g, matsubara_freq const &f, scalar_valued, std::true_type) const {
-    if (g->mesh().positive_only()){//only positive Matsubara frequencies
-     if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n];
-     if ((f.n < 0) && ((-f.n-sh(g)) < g->mesh().size())) return conj((*g)[-f.n-sh(g)]);
-    }
-    else{
-     if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size()+g->mesh().first_index())) return (*g)[f.n];
-    }
-    return 0;
-   }
-
-   template <typename G>
-   arrays::matrix_const_view<std::complex<double>> _call_impl(G const *g, matsubara_freq const &f, matrix_valued,
-                                                              std::false_type) const {
-    if (g->mesh().positive_only()){//only positive Matsubara frequencies
-     if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n]();
-     if ((f.n < 0) && ((-f.n-sh(g)) < g->mesh().size()))
-      return arrays::matrix<std::complex<double>>{conj((*g)[-f.n-sh(g)]())};
-    }
-    else{
-     if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size()+g->mesh().first_index())) return (*g)[f.n];
-    }
-    return g->singularity().evaluate(f);
-   }
-
-   template <typename G>
-   arrays::matrix_const_view<std::complex<double>> _call_impl(G const *g, matsubara_freq const &f, matrix_valued,
-                                                              std::true_type) const {
-    if (g->mesh().positive_only()){//only positive Matsubara frequencies
-     if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n]();
-     if ((f.n < 0) && ((-f.n-sh(g)) < g->mesh().size()))
-      return arrays::matrix<std::complex<double>>{conj((*g)[-f.n-sh(g)]())};
-    }
-    else{
-     if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size()+g->mesh().first_index())) return (*g)[f.n];
-    }
-    auto r = arrays::matrix<std::complex<double>>{get_target_shape(*g)};
-    r() = 0;
-    return r;
-   }
-
-   // does not work on gcc 4.8.1 ???
-   /* template <typename G>
-      auto operator()(G const *g, matsubara_freq const &f) const
-      DECL_AND_RETURN(_call_impl(g, f, Target{}, std::integral_constant<bool, std::is_same<Opt, no_tail>::value>{}));
-      */
-   public:
-
-   template <typename G>
-   typename std::conditional<std::is_same<Target, matrix_valued>::value, arrays::matrix_const_view<std::complex<double>>,
-                             std::complex<double>>::type
-   operator()(G const *g, matsubara_freq const &f) const {
-    return _call_impl(g, f, Target{}, std::integral_constant<bool, std::is_same<Opt, no_tail>::value>{});
-   }
-
    // int -> replace by matsubara_freq
-   template <typename G> auto operator()(G const *g, int n) const DECL_AND_RETURN((*g)(matsubara_freq(n,g->mesh().domain().beta,g->mesh().domain().statistic)));
-
-#ifdef __clang__
-   // to generate a clearer error message ? . Only ok on clang ?
-   template <int n> struct error {
-    static_assert(n > 0, "Green function cannot be evaluated on a complex number !");
-   };
-
-   template <typename G> error<0> operator()(G const *g, std::complex<double>) const {
-    return {};
-   }
-#endif
+   template <typename G>
+   AUTO_DECL operator()(G const *g, int n) const
+       RETURN((*g)(matsubara_freq(n, g->mesh().domain().beta, g->mesh().domain().statistic)));
 
    template <typename G> typename G::singularity_t const &operator()(G const *g, freq_infty const &) const {
     return g->singularity();
    }
   };
+  // --- various 4 specializations
 
-  /// ---------------------------  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> {};
+  // scalar_valued, tail
+  template <typename Opt> struct evaluator<imfreq, scalar_valued, local::tail, Opt> : _eval_imfreq_base_impl {
+ 
+   using _eval_imfreq_base_impl::operator();
+
+   template <typename G> std::complex<double> operator()(G const *g, matsubara_freq const &f) const {
+    if (g->mesh().positive_only()) { // only positive Matsubara frequencies
+     if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n];
+     if ((f.n < 0) && ((-f.n - this->sh(g)) < g->mesh().size())) return conj((*g)[-f.n - this->sh(g)]);
+    } else {
+     if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size() + g->mesh().first_index())) return (*g)[f.n];
+    }
+    return g->singularity().evaluate(f)(0, 0);
+   }
+  };
+
+  // scalar_valued, no tail
+  template <typename Opt> struct evaluator<imfreq, scalar_valued, nothing, Opt> : _eval_imfreq_base_impl {
+
+   using _eval_imfreq_base_impl::operator();
+
+   template <typename G> std::complex<double> operator()(G const *g, matsubara_freq const &f) const {
+    if (g->mesh().positive_only()) { // only positive Matsubara frequencies
+     if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n];
+     if ((f.n < 0) && ((-f.n - this->sh(g)) < g->mesh().size())) return conj((*g)[-f.n - this->sh(g)]);
+    } else {
+     if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size() + g->mesh().first_index())) return (*g)[f.n];
+    }
+    return 0;
+   }
+  };
+
+  // matrix_valued, tail
+  template <typename Opt> struct evaluator<imfreq, matrix_valued, local::tail, Opt> : _eval_imfreq_base_impl {
+
+   using _eval_imfreq_base_impl::operator();
+
+   template <typename G> arrays::matrix_const_view<std::complex<double>> operator()(G const *g, matsubara_freq const &f) const {
+    if (g->mesh().positive_only()) { // only positive Matsubara frequencies
+     if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n]();
+     if ((f.n < 0) && ((-f.n - this->sh(g)) < g->mesh().size()))
+      return arrays::matrix<std::complex<double>>{conj((*g)[-f.n - this->sh(g)]())};
+    } else {
+     if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size() + g->mesh().first_index())) return (*g)[f.n];
+    }
+    return g->singularity().evaluate(f);
+   }
+  };
+
+  // matrix_valued, no tail
+  template <typename Opt> struct evaluator<imfreq, matrix_valued, nothing, Opt> : _eval_imfreq_base_impl {
+
+   using _eval_imfreq_base_impl::operator();
+
+   template <typename G> arrays::matrix_const_view<std::complex<double>> operator()(G const *g, matsubara_freq const &f) const {
+    if (g->mesh().positive_only()) { // only positive Matsubara frequencies
+     if ((f.n >= 0) && (f.n < g->mesh().size())) return (*g)[f.n]();
+     if ((f.n < 0) && ((-f.n - this->sh(g)) < g->mesh().size()))
+      return arrays::matrix<std::complex<double>>{conj((*g)[-f.n - this->sh(g)]())};
+    } else {
+     if ((f.n >= g->mesh().first_index()) && (f.n < g->mesh().size() + g->mesh().first_index())) return (*g)[f.n];
+    }
+    auto r = arrays::matrix<std::complex<double>>{get_target_shape(*g)};
+    r() = 0;
+    return r;
+   }
+  };
 
  } // gfs_implementation
 
diff --git a/triqs/gfs/imfreq_bff.hpp b/triqs/gfs/imfreq_bff.hpp
new file mode 100644
index 00000000..5803edbd
--- /dev/null
+++ b/triqs/gfs/imfreq_bff.hpp
@@ -0,0 +1,122 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "./product.hpp"
+#include "./meshes/matsubara_freq.hpp"
+namespace triqs {
+namespace gfs {
+
+ // short cut. Here only the change compare to default multi var implementation
+ using imfreq_bff = cartesian_product<imfreq, imfreq, imfreq>;
+ //struct imfreq_bff {}; // = cartesian_product<imfreq, imfreq, imfreq>;
+ using imfreq_mesh_3 = mesh_product<matsubara_freq_mesh, matsubara_freq_mesh, matsubara_freq_mesh>;
+
+ // The default target for this mesh
+ template <> struct gf_default_target<imfreq_bff> {
+  using type = tensor_valued<4>;
+ };
+
+ // reimplement a simpler constructor, which enforces that the 2 meshes are equal.
+ template <typename Opt>
+ struct gf_mesh<imfreq_bff, Opt> : imfreq_mesh_3 {
+  gf_mesh() = default;
+  gf_mesh(matsubara_freq_mesh const& m1, matsubara_freq_mesh const& m2) : imfreq_mesh_3{m1, m2, m2} {}
+ };
+
+ namespace gfs_implementation {
+
+  /// ---------------------------  data access  ---------------------------------
+
+  struct imfreq_bff_indices_mixer {
+   template <typename MI, typename TI> static auto invoke(MI const& m, TI const& t) {
+    return std::make_tuple(std::get<0>(m), std::get<1>(m), std::get<0>(t), std::get<1>(t), std::get<2>(m), std::get<2>(t),
+                           std::get<3>(t));
+   }
+  };
+
+  template <>
+  struct data_proxy<imfreq_bff, tensor_valued<4>, void> : data_proxy_array_index_mixer<std::complex<double>, 3, 4,
+                                                                                       imfreq_bff_indices_mixer> {
+
+   //template <typename S, typename Tu> auto operator()(S& data, Tu const& tu) const {
+   // return data(std::get<0>(tu), arrays::range(), std::get<1>(tu), arrays::range(), std::get<2>(tu), arrays::range());
+  // }
+  };
+
+
+  // -------------------------------  evaluator  --------------------------------------------------
+
+  //template <typename Target, typename Opt>
+  //struct evaluator<imfreq_bff, Target, nothing, Opt> : evaluator<cartesian_product<imfreq, imfreq, imfreq>, Target, nothing, Opt> {};
+
+  /// ---------------------------  hdf5  ---------------------------------
+
+  template <typename Opt> struct h5_name<imfreq_bff, tensor_valued<4>, nothing, Opt> {
+   static std::string invoke() { return "imfreq_bff"; }
+  };
+
+ /// ---------------------------  partial eval ---------------------------------
+
+  template <bool IsConst> struct partial_eval_impl<imfreq_x2, matrix_valued, nothing, void, IsConst> {
+
+   template <int pos, typename Omega>
+   gf_view<imfreq_x2, matrix_valued, nothing, void, IsConst>
+   partial_eval(gf_view<imfreq_bff, tensor_valued<4>, nothing, void, IsConst> g, Omega const& omega) {
+    static_assert(pos == 0, "EE");
+    auto& m = g.mesh().components();
+    auto av = g.data()(std::get<0>(m).index_to_linear(omega), arrays::ellipsis{});
+    return {{std::get<1>(m)}, av, {}, {}, {}};
+   }
+  };
+
+ } // gfs_implementation
+
+ gf_mesh<imfreq> get_bosonic_mesh(gf_view<imfreq_bff, tensor_valued<4>> g) { return std::get<0>(g.mesh().components()); }
+
+
+/* template <typename Omega>
+ auto partial_eval(gf<imfreq_bff, tensor_valued<4>>& g, Omega const& omega) RETURN(partial_eval(g(), omega));
+
+ template <typename Omega>
+ auto partial_eval(gf<imfreq_bff, tensor_valued<4>> const& g, Omega const& omega) RETURN(partial_eval(g(), omega));
+
+ /// ---------------------------  curry ---------------------------------
+
+ template <bool IsConst> auto curry_on_bosonic_freq(gf_view<imfreq_bff, tensor_valued<4>, nothing, void, IsConst> g) {
+  return make_gf_view_lambda_valued<imfreq>(std::get<0>(g.mesh().components()),
+                                            [g](auto&& x) { return partial_eval(g, x); });
+ }
+
+ auto curry_on_bosonic_freq(gf<imfreq_bff, tensor_valued<4>>& g) RETURN(curry_on_bosonic_freq(g()));
+ auto curry_on_bosonic_freq(gf<imfreq_bff, tensor_valued<4>> const& g) RETURN(curry_on_bosonic_freq(g()));
+*/
+ 
+ }
+}
+
+
+/* template<typename G>
+ struct curry_polymorphic_lambda2 {
+  G g;
+  template <typename... I> AUTO_DECL operator()(I &&... i) const RETURN(partial_eval(g, i...));
+ };
+*/
+
diff --git a/triqs/gfs/imfreq_x2.hpp b/triqs/gfs/imfreq_x2.hpp
new file mode 100644
index 00000000..215419d1
--- /dev/null
+++ b/triqs/gfs/imfreq_x2.hpp
@@ -0,0 +1,97 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "./product.hpp"
+#include "./meshes/matsubara_freq.hpp"
+namespace triqs {
+namespace gfs {
+
+ struct matrix_operations {};
+
+ // in this version, we store all frequencies, to be able to write matrix operation simply.
+ // memory usage is x2 too large in principle.
+ 
+ // short cut. Here only the change compare to default multi var implementation
+ using imfreq_x2 = cartesian_product<imfreq, imfreq>;
+
+ // reimplement a simpler constructor, which enforces that the 2 meshes are equal.
+ template <typename Opt> struct gf_mesh<imfreq_x2, Opt> : mesh_product<matsubara_freq_mesh, matsubara_freq_mesh> {
+  using B = mesh_product<matsubara_freq_mesh, matsubara_freq_mesh>;
+  gf_mesh() = default;
+  gf_mesh(matsubara_freq_mesh const& m) : B{m, m} {}
+  //gf_mesh(matsubara_freq_mesh const & m1, matsubara_freq_mesh const& m2) : B{m1,m2} {} //needed for curry. DO NOT document.
+  template <typename... T>
+  gf_mesh(T&&... x)
+     : B{matsubara_freq_mesh(std::forward<T>(x)...), matsubara_freq_mesh(std::forward<T>(x)...)} {}
+ };
+
+ // The default target for this mesh
+ template <> struct gf_default_target<imfreq_x2> {
+  using type = tensor_valued<4>;
+ };
+ namespace gfs_implementation {
+
+  /// ---------------------------  data access  ---------------------------------
+
+
+   struct imfreq_x2_indices_mixer {
+   template <typename MI, typename TI> static auto invoke(MI const& m, TI const& t) {
+    return std::make_tuple(std::get<0>(m), std::get<0>(t), std::get<1>(t), std::get<1>(m), std::get<2>(t), std::get<3>(t));
+   }
+  };
+
+  template <>
+  struct data_proxy<imfreq_x2, tensor_valued<4>, void> : data_proxy_array_index_mixer<std::complex<double>, 2, 4,
+                                                                                       imfreq_x2_indices_mixer> {
+
+   //template <typename S, typename Tu> auto operator()(S& data, Tu const& tu) const {
+   // return data(std::get<0>(tu), arrays::range(), std::get<1>(tu), arrays::range(), std::get<2>(tu), arrays::range());
+  // }
+  };
+/*
+  // Change the ordering of the indices in the array to allow matrix operations
+  template <typename Opt> struct data_proxy<imfreq_x2, tensor_valued<4>, Opt> : data_proxy_array_common<std::complex<double>, 6> {
+
+   template <typename S1, typename S2> static utility::mini_vector<int, 4> join_shape(S1 const& s1, S2 const& s2) {
+    return {int(s1[0]), s2[0], int(s1[1]), s2[1]};
+    // TODO : clean this size_t....
+   }
+
+   template <typename S, typename Tu>
+   AUTO_DECL operator()(S& data, Tu const& tu) const
+       RETURN(make_matrix_view(data(std::get<0>(tu), arrays::range(), std::get<1>(tu), arrays::range())));
+  };
+
+  */
+ } // gfs_implementation
+
+ /// ---------------------------  inverse  ---------------------------------
+
+ // the generic inverse is fine. We only need to redo the invert_in_place.
+ template <typename Opt> void invert_in_place(gf_view<imfreq_x2, matrix_valued, nothing, Opt> g) {
+
+ auto arr = make_matrix_view(group_indices_view(g.data(), {0, 1, 2}, {3, 4, 5})); // a view of the array properly regrouped
+ arr = inverse(arr);                                      // inverse as a big matrix (nu,n) x (nu', n')
+
+  // no singularity
+ }
+}
+}
diff --git a/triqs/gfs/imtime.hpp b/triqs/gfs/imtime.hpp
index ed7c38a0..d5a126ce 100644
--- a/triqs/gfs/imtime.hpp
+++ b/triqs/gfs/imtime.hpp
@@ -37,18 +37,18 @@ namespace gfs {
   // using matsubara_time_mesh::matsubara_time_mesh;
  };
 
+ // singularity
+ template <> struct gf_default_singularity<imtime, matrix_valued> {
+  using type = local::tail;
+ };
+ template <> struct gf_default_singularity<imtime, scalar_valued> {
+  using type = local::tail;
+ };
+
  namespace gfs_implementation {
 
-  // singularity. If no_tail is given, then it is the default (nothing)
-  template <> struct singularity<imtime, matrix_valued, void> {
-   using type = local::tail;
-  };
-  template <> struct singularity<imtime, scalar_valued, void> {
-   using type = local::tail;
-  };
-
   // h5 name
-  template <typename Opt> struct h5_name<imtime, matrix_valued, Opt> {
+  template <typename Singularity, typename Opt> struct h5_name<imtime, matrix_valued, Singularity, Opt> {
    static std::string invoke() { return "ImTime"; }
   };
 
@@ -59,7 +59,7 @@ namespace gfs {
 
   /// ---------------------------  closest mesh point on the grid ---------------------------------
 
-  template <typename Opt, typename Target> struct get_closest_point<imtime, Target, Opt> {
+  template <typename Singularity, typename Opt, typename Target> struct get_closest_point<imtime, Target, Singularity, Opt> {
    // index_t is int
    template <typename G, typename T> static int invoke(G const *g, closest_pt_wrap<T> const &p) {
     double x = (g->mesh().kind() == half_bins ? double(p.value) : double(p.value) + 0.5 * g->mesh().delta());
@@ -98,7 +98,8 @@ namespace gfs {
   };
 
   // now evaluator
-  template <typename Opt, typename Target> struct evaluator<imtime, Target, Opt> : evaluator_one_var<imtime> {};
+  template <typename Singularity, typename Opt, typename Target>
+  struct evaluator<imtime, Target, Singularity, Opt> : evaluator_one_var<imtime> {};
 
  } // gfs_implementation.
 }
diff --git a/triqs/gfs/legendre.hpp b/triqs/gfs/legendre.hpp
index aa5e0608..d0eebf33 100644
--- a/triqs/gfs/legendre.hpp
+++ b/triqs/gfs/legendre.hpp
@@ -26,39 +26,43 @@
 #include "./domains/legendre.hpp"
 #include "./meshes/discrete.hpp"
 
-namespace triqs { namespace gfs {
+namespace triqs {
+namespace gfs {
 
  struct legendre {};
 
  // mesh type and its factories
- template<typename Opt> struct gf_mesh<legendre,Opt> :discrete_mesh<legendre_domain> { 
+ template <typename Opt> struct gf_mesh<legendre, Opt> : discrete_mesh<legendre_domain> {
   typedef discrete_mesh<legendre_domain> B;
   gf_mesh() = default;
-  gf_mesh(double beta, statistic_enum S, size_t n_leg) : B(typename B::domain_t(beta,S,n_leg)) {}
+  gf_mesh(double beta, statistic_enum S, size_t n_leg) : B(typename B::domain_t(beta, S, n_leg)) {}
  };
 
- namespace gfs_implementation { 
+ namespace gfs_implementation {
 
   // h5 name
-  template<typename Opt> struct h5_name<legendre,matrix_valued,Opt>      { static std::string invoke(){ return  "Legendre";}};
+  template <typename Opt> struct h5_name<legendre, matrix_valued, nothing, Opt> {
+   static std::string invoke() { return "Legendre"; }
+  };
 
   /// ---------------------------  evaluator ---------------------------------
 
   // Not finished, not tested
-  template<typename Opt>
-   struct evaluator<legendre,matrix_valued,Opt> {
-    static constexpr int arity = 1;
-    //ERROR : give a double and interpolate
-    template<typename G>
-     arrays::matrix_view<double >  operator() (G const * g,long n)  const {return g->data()(n, arrays::range(), arrays::range()); }
-   };
+  template <typename Opt> struct evaluator<legendre, matrix_valued, nothing, Opt> {
+   static constexpr int arity = 1;
+   // ERROR : give a double and interpolate
+   template <typename G> arrays::matrix_view<double> operator()(G const* g, long n) const {
+    return g->data()(n, arrays::range(), arrays::range());
+   }
+  };
 
   /// ---------------------------  data access  ---------------------------------
 
-  template<typename Opt> struct data_proxy<legendre,matrix_valued,Opt> : data_proxy_array<double,3> {};
-  template<typename Opt> struct data_proxy<legendre,scalar_valued,Opt> : data_proxy_array<double,1> {};
+  template <typename Opt> struct data_proxy<legendre, matrix_valued, Opt> : data_proxy_array<double, 3> {};
+  template <typename Opt> struct data_proxy<legendre, scalar_valued, Opt> : data_proxy_array<double, 1> {};
 
  } // gfs_implementation
-}}
+}
+}
 #endif
 
diff --git a/triqs/gfs/local/fit_tail.cpp b/triqs/gfs/local/fit_tail.cpp
index 5b484291..250c4a17 100644
--- a/triqs/gfs/local/fit_tail.cpp
+++ b/triqs/gfs/local/fit_tail.cpp
@@ -25,8 +25,8 @@ namespace triqs { namespace gfs {  namespace local {
   int size1 = n_max - n_min + 1;
   // size2 is the number of moments
 
-  arrays::matrix<double, 2> A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT);
-  arrays::matrix<double, 2> B(size1, 1, FORTRAN_LAYOUT);
+  arrays::matrix<double> A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT);
+  arrays::matrix<double> B(size1, 1, FORTRAN_LAYOUT);
   arrays::vector<double> S(std::max(size_even, size_odd));
   const double rcond = 0.0;
   int rank;
diff --git a/triqs/gfs/local/fourier_matsubara.hpp b/triqs/gfs/local/fourier_matsubara.hpp
index f4880aa1..de3e49ca 100644
--- a/triqs/gfs/local/fourier_matsubara.hpp
+++ b/triqs/gfs/local/fourier_matsubara.hpp
@@ -18,8 +18,7 @@
  * TRIQS. If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef TRIQS_GF_LOCAL_FOURIER_MATSU_H
-#define TRIQS_GF_LOCAL_FOURIER_MATSU_H
+#pragma once
 
 #include "fourier_base.hpp"
 #include <triqs/gfs/imfreq.hpp>
@@ -28,12 +27,12 @@
 namespace triqs {
 namespace gfs {
 
- template <typename Target, typename Opt, bool V, bool C>
- gf_keeper<tags::fourier, imtime, Target, Opt> fourier(gf_impl<imtime, Target, Opt, V, C> const& g) {
+ template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf_keeper<tags::fourier, imtime, Target, Singularity, Opt> fourier(gf_impl<imtime, Target, Singularity, Opt, V, C> const& g) {
   return {g};
  }
- template <typename Target, typename Opt, bool V, bool C>
- gf_keeper<tags::fourier, imfreq, Target, Opt> inverse_fourier(gf_impl<imfreq, Target, Opt, V, C> const& g) {
+ template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf_keeper<tags::fourier, imfreq, Target, Singularity, Opt> inverse_fourier(gf_impl<imfreq, Target, Singularity, Opt, V, C> const& g) {
   return {g};
  }
 
@@ -49,30 +48,30 @@ namespace gfs {
  void triqs_gf_view_assign_delegation(gf_view<imfreq, matrix_valued, no_tail> g,
                                       gf_keeper<tags::fourier, imtime, matrix_valued, no_tail> const& L);
 
- template <typename Target, typename Opt, bool V, bool C>
- gf<imfreq, Target, Opt> make_gf_from_fourier(gf_impl<imtime, Target, Opt, V, C> const& gt, int n_iw) {
+ template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf<imfreq, Target, Singularity, Opt> make_gf_from_fourier(gf_impl<imtime, Target, Singularity, Opt, V, C> const& gt, int n_iw) {
   auto m = gf_mesh<imfreq, Opt>{gt.mesh().domain(), n_iw};
-  auto gw = gf<imfreq, Target, Opt>{m, get_target_shape(gt)};
+  auto gw = gf<imfreq, Target, Singularity, Opt>{m, get_target_shape(gt)};
   gw() = fourier(gt);
   return gw;
  }
 
- template <typename Target, typename Opt, bool V, bool C>
- gf<imfreq, Target, Opt> make_gf_from_fourier(gf_impl<imtime, Target, Opt, V, C> const& gt) {
+ template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf<imfreq, Target, Singularity, Opt> make_gf_from_fourier(gf_impl<imtime, Target, Singularity, Opt, V, C> const& gt) {
   return make_gf_from_fourier(gt, (gt.mesh().size() - (gt.mesh().kind() == full_bins ? 1 : 0)) / 2);
  }
 
- template <typename Target, typename Opt, bool V, bool C>
- gf<imtime, Target, Opt> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Opt, V, C> const& gw, int n_tau,
+ template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf<imtime, Target, Singularity, Opt> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Singularity, Opt, V, C> const& gw, int n_tau,
                                                       mesh_kind mk = full_bins) {
   auto m = gf_mesh<imtime, Opt>{gw.mesh().domain(), n_tau};
-  auto gt = gf<imtime, Target, Opt>{m, get_target_shape(gw)};
+  auto gt = gf<imtime, Target, Singularity, Opt>{m, get_target_shape(gw)};
   gt() = inverse_fourier(gw);
   return gt;
  }
 
- template <typename Target, typename Opt, bool V, bool C>
- gf<imtime, Target, Opt> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Opt, V, C> const& gw, mesh_kind mk = full_bins) {
+ template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf<imtime, Target, Singularity, Opt> make_gf_from_inverse_fourier(gf_impl<imfreq, Target, Singularity, Opt, V, C> const& gw, mesh_kind mk = full_bins) {
   return make_gf_from_inverse_fourier(gw, 2 * gw.mesh().size() + (mk == full_bins ? 1 : 0), mk);
  }
 }
@@ -82,5 +81,4 @@ namespace clef {
  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 63d1a66c..9ddf4166 100644
--- a/triqs/gfs/local/fourier_real.hpp
+++ b/triqs/gfs/local/fourier_real.hpp
@@ -18,21 +18,19 @@
  * TRIQS. If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef TRIQS_GF_LOCAL_FOURIER_REAL_H 
-#define TRIQS_GF_LOCAL_FOURIER_REAL_H
-
+#pragma once
 #include "fourier_base.hpp"
 #include <triqs/gfs/refreq.hpp> 
 #include <triqs/gfs/retime.hpp> 
 
 namespace triqs { namespace gfs { 
 
-  template <typename Target, typename Opt, bool V, bool C>
- gf_keeper<tags::fourier, retime, Target> fourier(gf_impl<retime, Target, Opt, V, C> const& g) {
+  template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf_keeper<tags::fourier, retime, Target, Singularity> fourier(gf_impl<retime, Target, Singularity, Opt, V, C> const& g) {
   return {g};
  }
- template <typename Target, typename Opt, bool V, bool C>
- gf_keeper<tags::fourier, refreq, Target> inverse_fourier(gf_impl<refreq, Target, Opt, V, C> const& g) {
+ template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf_keeper<tags::fourier, refreq, Target, Singularity> inverse_fourier(gf_impl<refreq, Target, Singularity, Opt, V, C> const& g) {
   return {g};
  }
 
@@ -58,20 +56,19 @@ namespace triqs { namespace gfs {
   return {tmin, tmax, L};
  }
 
-  template <typename Target, typename Opt, bool V, bool C>
- gf_view<refreq, Target> make_gf_from_fourier(gf_impl<retime, Target, Opt, V, C> const& gt) {
+  template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf_view<refreq, Target> make_gf_from_fourier(gf_impl<retime, Target, Singularity, Opt, V, C> const& gt) {
   auto gw = gf<refreq, Target>{make_mesh_fourier_compatible(gt.mesh()), get_target_shape(gt)};
   gw() = fourier(gt);
   return gw;
  }
 
- template <typename Target, typename Opt, bool V, bool C>
- gf_view<retime, Target> make_gf_from_inverse_fourier(gf_impl<refreq, Target, Opt, V, C> const& gw) {
+ template <typename Target, typename Singularity, typename Opt, bool V, bool C>
+ gf_view<retime, Target> make_gf_from_inverse_fourier(gf_impl<refreq, Target, Singularity, Opt, V, C> const& gw) {
   auto gt = gf<retime, Target>{make_mesh_fourier_compatible(gw.mesh()), get_target_shape(gw)};
   gt() = inverse_fourier(gw);
   return gt;
  }
 
 }}
-#endif
 
diff --git a/triqs/gfs/local/no_tail.hpp b/triqs/gfs/local/no_tail.hpp
index 6b38cc05..f23f7251 100644
--- a/triqs/gfs/local/no_tail.hpp
+++ b/triqs/gfs/local/no_tail.hpp
@@ -26,10 +26,11 @@
 namespace triqs {
 namespace gfs {
 
- struct no_tail {};
+ //struct no_tail {};
+ using no_tail = nothing;
 
- template <typename Variable, typename Target, bool V, bool C>
- gf_view<Variable, Target, no_tail, C> make_gf_view_without_tail(gf_impl<Variable, Target, void, V, C> const &g) {
+ template <typename Variable, typename Target, typename S, typename Opt, bool V, bool C>
+ gf_view<Variable, Target, no_tail, Opt, C> make_gf_view_without_tail(gf_impl<Variable, Target, S, Opt, V, C> const &g) {
   return {g.mesh(), g.data(), {}, g.symmetry(), g.indices(), g.name};
  }
 
@@ -46,15 +47,15 @@ namespace gfs {
   }
  }
 
- template <typename Variable, typename Target, bool V, bool C>
- gf_view<Variable, Target> make_gf_from_g_and_tail(gf_impl<Variable, Target, no_tail, V, C> const &g, local::tail t) {
+ template <typename Variable, typename Target, typename S, typename Opt, bool V, bool C>
+ gf_view<Variable, Target> make_gf_from_g_and_tail(gf_impl<Variable, Target, S, Opt, V, C> const &g, local::tail t) {
   details::_equal_or_throw(t.shape(), get_target_shape(g));
   auto g2 = gf<Variable, Target, no_tail>{g}; // copy the function without tail
   return {std::move(g2.mesh()), std::move(g2.data()), std::move(t), g2.symmetry()};
  }
 
- template <typename Variable, typename Target, bool V, bool C>
- gf_view<Variable, Target, void, C> make_gf_view_from_g_and_tail(gf_impl<Variable, Target, no_tail, V, C> const &g,
+ template <typename Variable, typename Target, typename S, typename Opt, bool V, bool C>
+ gf_view<Variable, Target, local::tail, Opt, C> make_gf_view_from_g_and_tail(gf_impl<Variable, Target, no_tail, Opt, V, C> const &g,
                                                                  local::tail_view t) {
   details::_equal_or_throw(t.shape(), get_target_shape(g));
   return {g.mesh(), g.data(), t, g.symmetry()};
diff --git a/triqs/gfs/local/tail.hpp b/triqs/gfs/local/tail.hpp
index 5de8e8e4..ddd2ebbd 100644
--- a/triqs/gfs/local/tail.hpp
+++ b/triqs/gfs/local/tail.hpp
@@ -53,6 +53,7 @@ namespace triqs { namespace gfs { namespace local {
     public:
       TRIQS_MPI_IMPLEMENTED_VIA_BOOST;
       typedef tail_view view_type;
+      typedef tail_view const_view_type; // not nice
       typedef tail      regular_type;
 
       typedef arrays::array      <dcomplex,3>                         data_regular_type;
@@ -232,6 +233,7 @@ namespace triqs { namespace gfs { namespace local {
     typedef tqa::mini_vector<size_t,2> shape_type;
     tail(size_t N1, size_t N2, size_t size_=10, long order_min=-1): B(N1,N2,size_,order_min) {}
     tail(shape_type const & sh, size_t size_=10, long order_min=-1): B(sh[0],sh[1],size_,order_min) {}
+    tail(tqa::mini_vector<int,0>) : tail(1,1) {}
     tail(B::data_type const &d, B::mask_type const &m, long order_min): B(d, m, order_min) {}
     tail(tail const & g): B(g) {}
     tail(tail_view const & g): B(g) {}
diff --git a/triqs/gfs/m_tail.hpp b/triqs/gfs/m_tail.hpp
new file mode 100644
index 00000000..1d5ef4e3
--- /dev/null
+++ b/triqs/gfs/m_tail.hpp
@@ -0,0 +1,103 @@
+/*******************************************************************************
+ *
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2012-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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#pragma once
+#include "./tools.hpp"
+#include "./gf.hpp"
+#include "./local/tail.hpp"
+#include "./meshes/discrete.hpp"
+
+namespace triqs {
+namespace gfs {
+
+ namespace gfs_implementation {
+
+  /*template <typename T> constexpr bool has_a_singularity() {
+   return std::is_same<typename singularity<V1, Target>::type, local::tail>::value;
+  }
+
+  /// ---------------------------  singularity ---------------------------------
+
+  template <typename V1, typename V2, typename Target, typename Opt> struct singularity<cartesian_product<V1,V2>, Target, Opt> {
+   using type = std14::conditional_t < (!has_a_singularity<V1>()) && has_a_singularity<V2>(), gf<V1, local::tail>, nothing > ;
+  };
+  */
+  /// ---------------------------  hdf5 ---------------------------------
+ 
+  template <typename Variable, typename Opt> struct h5_name<Variable, local::tail, nothing, Opt> {
+   static std::string invoke() { return "xxxxx"; }
+  };
+
+  template <typename Variable, typename Opt> struct h5_rw<Variable, local::tail, nothing, Opt> {
+
+   static void write(h5::group gr, gf_const_view<Variable, local::tail, nothing, Opt> g) {
+    for (size_t i = 0; i < g.mesh().size(); ++i) h5_write(gr, std::to_string(i), g._data[i]);
+    // h5_write(gr,"symmetry",g._symmetry);
+   }
+
+   template <bool IsView> static void read(h5::group gr, gf_impl<Variable, local::tail, nothing, Opt, IsView, false> &g) {
+    // does not work : need to read the block name and remake the mesh...
+    //g._mesh = gf_mesh<block_index, Opt>(gr.get_all_subgroup_names());
+    //g._data.resize(g._mesh.size());
+    // if (g._data.size() != g._mesh.size()) TRIQS_RUNTIME_ERROR << "h5 read block gf : number of block mismatch";
+    //for (size_t i = 0; i < g.mesh().size(); ++i) h5_read(gr, g.mesh().domain().names()[i], g._data[i]);
+    // h5_read(gr,"symmetry",g._symmetry);
+   }
+  };
+
+  /// ---------------------------  data access  ---------------------------------
+
+  template <typename Variable, typename Opt> struct data_proxy<Variable, local::tail, Opt> : data_proxy_vector<local::tail> {};
+
+  // -------------------------------   Factories  --------------------------------------------------
+
+  template <typename Variable, typename Opt> struct factories<Variable, local::tail, nothing, Opt> {
+   using mesh_t=gf_mesh<Variable, Opt> ;
+   using gf_t=gf<Variable, local::tail> ;
+   //using gf_view_t=gf_view<block_index, Target> ;
+   struct target_shape_t {};
+
+   static typename gf_t::data_t make_data(mesh_t const &m, target_shape_t) { return std::vector<local::tail>(m.size()); }
+   static nothing make_singularity(mesh_t const &m, target_shape_t) {
+    return {};
+   }
+  };
+
+  // partial_eval
+  template <typename Variable, typename Opt, bool IsConst> struct partial_eval_impl<Variable, local::tail, nothing, Opt, IsConst> {
+
+   using gv_t = gf_view<Variable, local::tail, nothing, Opt, IsConst>;
+   template <int... pos, typename... T> static auto invoke(gv_t g, T const &... x) {
+    return invoke_impl(g, std14::index_sequence<pos...>(), x...);
+   }
+
+   template <typename T> static local::tail_view invoke_impl(gv_t g, std14::index_sequence<0>, T const &x) {
+    return g[g.mesh().index_to_linear(x)];
+   }
+
+   template <typename T> static nothing invoke_impl(gv_t g, std14::index_sequence<1>, T const &x) {
+    return {};
+   }
+  };
+
+ } // gfs_implementation
+
+
+}}
diff --git a/triqs/gfs/meshes/discrete.hpp b/triqs/gfs/meshes/discrete.hpp
index 6279512f..e661781f 100644
--- a/triqs/gfs/meshes/discrete.hpp
+++ b/triqs/gfs/meshes/discrete.hpp
@@ -29,12 +29,18 @@ namespace gfs {
 
   using domain_t = Domain;
   using index_t = long;
+  using linear_index_t = long;
 
   discrete_mesh(domain_t dom) : _dom(std::move(dom)) {}
   discrete_mesh() = default;
 
   domain_t const &domain() const { return _dom; }
-  long size() const { return _dom.size(); }
+  size_t size() const { return _dom.size(); }
+ 
+  ///
+  utility::mini_vector<size_t, 1> size_of_components() const {
+   return {size()};
+  }
 
   /// Conversions point <-> index <-> discrete_index
   long index_to_point(index_t ind) const { return ind; }
diff --git a/triqs/gfs/meshes/linear.hpp b/triqs/gfs/meshes/linear.hpp
index 0fe1834e..2d8646b0 100644
--- a/triqs/gfs/meshes/linear.hpp
+++ b/triqs/gfs/meshes/linear.hpp
@@ -37,6 +37,7 @@ namespace gfs {
 
   using domain_t = Domain;
   using index_t = long;
+  using linear_index_t = long;
   using domain_pt_t = typename domain_t::point_t;
 
   static_assert(!std::is_base_of<std::complex<double>, domain_pt_t>::value,
@@ -64,7 +65,12 @@ namespace gfs {
   }
 
   domain_t const &domain() const { return _dom; }
-  long size() const { return L; }
+  size_t size() const { return L; }
+  
+  utility::mini_vector<size_t, 1> size_of_components() const {
+   return {size()};
+  }
+
   double delta() const { return del; }
   double x_max() const { return xmax; }
   double x_min() const { return xmin; }
diff --git a/triqs/gfs/meshes/matsubara_freq.hpp b/triqs/gfs/meshes/matsubara_freq.hpp
index 90245eac..04ea9cbd 100644
--- a/triqs/gfs/meshes/matsubara_freq.hpp
+++ b/triqs/gfs/meshes/matsubara_freq.hpp
@@ -21,6 +21,7 @@
 #pragma once
 #include "./mesh_tools.hpp"
 #include "../domains/matsubara.hpp"
+#include <triqs/mpi.hpp>
 
 namespace triqs {
 namespace gfs {
@@ -31,6 +32,7 @@ namespace gfs {
   using domain_t = matsubara_domain<true>;
   ///
   using index_t = long;
+  using linear_index_t = long;
   ///
   using domain_pt_t = typename domain_t::point_t;
 
@@ -98,6 +100,11 @@ namespace gfs {
   /// Size (linear) of the mesh of the window
   long size() const { return _last_index_window - _first_index_window + 1; }
 
+  ///
+  utility::mini_vector<size_t, 1> size_of_components() const {
+   return {size_t(size())};
+  }
+
   /// From an index of a point in the mesh, returns the corresponding point in the domain
   domain_pt_t index_to_point(index_t ind) const { return 1_j * M_PI * (2 * ind + (_dom.statistic == Fermion)) / _dom.beta; }
 
diff --git a/triqs/gfs/meshes/product.c14.hpp b/triqs/gfs/meshes/product.c14.hpp
index 3c90620f..dd0f3e28 100644
--- a/triqs/gfs/meshes/product.c14.hpp
+++ b/triqs/gfs/meshes/product.c14.hpp
@@ -2,7 +2,7 @@
  *
  * TRIQS: a Toolbox for Research in Interacting Quantum Systems
  *
- * Copyright (C) 2012-2013 by O. Parcollet
+ * Copyright (C) 2012-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
@@ -34,6 +34,7 @@ namespace gfs {
   using m_tuple_t = std::tuple<Meshes...>;
   using m_pt_tuple_t = std::tuple<typename Meshes::mesh_point_t...>;
   using domain_pt_t = typename domain_t::point_t;
+  using linear_index_t = std::tuple<typename Meshes::linear_index_t...>;
 
   static constexpr int dim = sizeof...(Meshes);
 
@@ -49,7 +50,7 @@ namespace gfs {
   size_t size() const {
    return triqs::tuple::fold([](auto const &m, size_t R) { return R * m.size(); }, m_tuple, 1);
   }
-
+/*
   /// Scatter the first mesh over the communicator c
   friend mesh_product mpi_scatter(mesh_product const &m, mpi::communicator c, int root) {
    auto r = m; // same domain, but mesh with a window. Ok ?
@@ -58,10 +59,18 @@ namespace gfs {
   }
 
   /// Opposite of scatter : rebuild the original mesh, without a window
-  friend matsubara_freq_mesh mpi_gather(matsubara_freq_mesh m, mpi::communicator c, int root) {
+  friend mesh_product mpi_gather(mesh_product m, mpi::communicator c, int root) {
    auto r = m; // same domain, but mesh with a window. Ok ?
    std::get<0>(r.m_tuple) = mpi_gather(std::get<0>(r.m_tuple), c, root);
    return r;
+  */
+
+   /// The sizes of all mesh components
+  utility::mini_vector<size_t, dim> size_of_components() const {
+   utility::mini_vector<size_t, dim> res;
+   auto l = [&res](int i, auto const &m) mutable { res[i] = m.size(); };
+   triqs::tuple::for_each_enumerate(m_tuple, l);
+   return res;
   }
 
   /// Conversions point <-> index <-> linear_index
@@ -72,32 +81,23 @@ namespace gfs {
    return res;
   }
 
-  /// Flattening index to linear :  index[0] + component[0].size * (index[1] + component[1].size* (index[2] + ....))
-  size_t index_to_linear(index_t const &ii) const {
-   auto l = [](auto const &m, auto const &i, size_t R) { return m.index_to_linear(i) + R * m.size(); };
-   return triqs::tuple::fold(l, reverse(m_tuple), reverse(ii), size_t(0));
+  /// The linear_index is the tuple of the linear_index of the components
+  linear_index_t index_to_linear(index_t const &ind) const {
+   auto l = [](auto const &m, auto const &i) { return m.index_to_linear(i); };
+   return triqs::tuple::map_on_zip(l, m_tuple, ind);
   }
   
-  /// Flattening index to linear :  index[0] + component[0].size * (index[1] + component[1].size* (index[2] + ....))
-  size_t mp_to_linear(m_pt_tuple_t const &mp) const {
-   auto l = [](auto const &m, auto const &p, size_t R) { return p.linear_index() + R * m.size(); };
-   return triqs::tuple::fold(l, reverse(m_tuple), reverse(mp), size_t(0));
-  }
-
-  utility::mini_vector<size_t, dim> shape() const {
-   utility::mini_vector<size_t, dim> res;
-   auto l = [&res](int i, auto const &m) mutable { res[i] = m.size(); };
-   triqs::tuple::for_each_enumerate(m_tuple, l);
-   return res;
+  /// 
+  linear_index_t mp_to_linear(m_pt_tuple_t const &mp) const {
+   auto l = [](auto const &p) { return p.linear_index(); };
+   return triqs::tuple::map(l, mp);
   }
 
   // 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 ");
-   // static_assert(std::is_same< std::tuple<typename std::remove_cv<typename std::remove_reference<MP>::type>::type...>,
-   // m_pt_tuple_t>::value, "Call incorrect ");
    return mp_to_linear(std::forward_as_tuple(mp...));
-  } // speed test ? or make a variadic fold...
+  } 
 
   /// The wrapper for the mesh point
   class mesh_point_t : tag::mesh_point {
@@ -116,7 +116,7 @@ namespace gfs {
       , _atend(false) {}
    mesh_point_t(mesh_product const &m_) : m(&m_), _c(triqs::tuple::map(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); }
+   linear_index_t linear_index() const { return m->mp_to_linear(_c); }
    const mesh_product *mesh() const { return m; }
 
    using cast_t = domain_pt_t;
@@ -193,14 +193,5 @@ namespace gfs {
 
  template <int pos, typename P> decltype(auto) get_component(P const &p) { return 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, ull_t To, int R, bool B, bool C>
- arrays::array_view<T, sizeof...(Meshes) + R - 1, OptionsFlags, To, true, C>
- reinterpret_linear_array(mesh_product<Meshes...> const &m, arrays::array_view<T, R, OptionsFlags, To, B, C> A) {
-  return {{join(m.shape(), get_shape(A).front_pop())}, A.storage()};
- }
 }
 }
diff --git a/triqs/gfs/product.hpp b/triqs/gfs/product.hpp
index 51a8f559..69830510 100644
--- a/triqs/gfs/product.hpp
+++ b/triqs/gfs/product.hpp
@@ -35,6 +35,7 @@ namespace gfs {
  // use alias
  template <typename... Ms> struct cartesian_product<std::tuple<Ms...>> : cartesian_product<Ms...> {};
 
+ /// TODO : Put inheriting constructor, simpler...
  // the mesh is simply a cartesian product
  template <typename Opt, typename... Ms> struct gf_mesh<cartesian_product<Ms...>, Opt> : mesh_product<gf_mesh<Ms, Opt>...> {
   // using mesh_product< gf_mesh<Ms,Opt> ... >::mesh_product< gf_mesh<Ms,Opt> ... > ;
@@ -45,57 +46,37 @@ namespace gfs {
 
  namespace gfs_implementation {
 
-  /// ---------------------------  hdf5 ---------------------------------
-
-  // h5 name : name1_x_name2_.....
-  template <typename Opt, typename... Ms> struct h5_name<cartesian_product<Ms...>, matrix_valued, Opt> {
-   static std::string invoke() {
-    return triqs::tuple::fold([](std::string a, std::string b) { return a + std::string(b.empty() ? "" : "_x_") + b; },
-                              std::make_tuple(h5_name<Ms, matrix_valued, Opt>::invoke()...), std::string());
-   }
-  };
-  template <typename Opt, int R, typename... Ms>
-  struct h5_name<cartesian_product<Ms...>, tensor_valued<R>, Opt> : h5_name<cartesian_product<Ms...>, matrix_valued, Opt> {};
-
-  // a slight difference with the generic case : reinterpret the data array to avoid flattening the variables
-  template <typename Opt, int R, typename... Ms> struct h5_rw<cartesian_product<Ms...>, tensor_valued<R>, Opt> {
-   using g_t = gf<cartesian_product<Ms...>, tensor_valued<R>, Opt>;
-
-   static void write(h5::group gr, typename g_t::const_view_type g) {
-    h5_write(gr, "data", reinterpret_linear_array(g.mesh(), g().data()));
-    h5_write(gr, "singularity", g._singularity);
-    h5_write(gr, "mesh", g._mesh);
-    h5_write(gr, "symmetry", g._symmetry);
-   }
-
-   template <bool IsView>
-   static void read(h5::group gr, gf_impl<cartesian_product<Ms...>, tensor_valued<R>, Opt, IsView, false> &g) {
-    using G_t = gf_impl<cartesian_product<Ms...>, tensor_valued<R>, Opt, IsView, false>;
-    h5_read(gr, "mesh", g._mesh);
-    auto arr = arrays::array<typename G_t::data_t::value_type, sizeof...(Ms) + R>{};
-    h5_read(gr, "data", arr);
-    auto sh = arr.shape();
-    arrays::mini_vector<size_t, R + 1> sh2;
-    sh2[0] = g._mesh.size();
-    for (int u = 1; u < R + 1; ++u) sh2[u] = sh[sizeof...(Ms) - 1 + u];
-    g._data = arrays::array<typename G_t::data_t::value_type, R + 1>{sh2, std::move(arr.storage())};
-    h5_read(gr, "singularity", g._singularity);
-    h5_read(gr, "symmetry", g._symmetry);
-   }
-  };
-
   /// ---------------------------  data access  ---------------------------------
 
   template <typename Opt, typename... Ms>
-  struct data_proxy<cartesian_product<Ms...>, scalar_valued, Opt> : data_proxy_array<std::complex<double>, 1> {};
+  struct data_proxy<cartesian_product<Ms...>, scalar_valued, Opt> : data_proxy_array_multivar<std::complex<double>,
+                                                                                              sizeof...(Ms)> {};
+
   template <typename Opt, typename... Ms>
-  struct data_proxy<cartesian_product<Ms...>, matrix_valued, Opt> : data_proxy_array<std::complex<double>, 3> {};
+  struct data_proxy<cartesian_product<Ms...>, matrix_valued, Opt> : data_proxy_array_multivar_matrix_valued<std::complex<double>,
+                                                                                                            2 + sizeof...(Ms)> {};
+
   template <int R, typename Opt, typename... Ms>
-  struct data_proxy<cartesian_product<Ms...>, tensor_valued<R>, Opt> : data_proxy_array<std::complex<double>, R + 1> {};
- 
+  struct data_proxy<cartesian_product<Ms...>, tensor_valued<R>, Opt> : data_proxy_array_multivar<std::complex<double>,
+                                                                                                 R + sizeof...(Ms)> {};
+
+  // special case ? Or make a specific container....
   template <typename Opt, typename M0>
-  struct data_proxy<cartesian_product<M0,imtime>, matrix_valued, Opt> : data_proxy_array<double, 3> {};
-  
+  struct data_proxy<cartesian_product<M0, imtime>, matrix_valued, Opt> : data_proxy_array_multivar_matrix_valued<double, 2 + 2> {
+  };
+
+  /// ---------------------------  hdf5 ---------------------------------
+
+  // h5 name : name1_x_name2_.....
+  template <typename Opt, typename... Ms> struct h5_name<cartesian_product<Ms...>, matrix_valued, nothing, Opt> {
+   static std::string invoke() {
+    return triqs::tuple::fold([](std::string a, std::string b) { return a + std::string(b.empty() ? "" : "_x_") + b; },
+                              std::make_tuple(h5_name<Ms, matrix_valued, nothing, Opt>::invoke()...), std::string());
+   }
+  };
+  template <typename Opt, int R, typename... Ms>
+  struct h5_name<cartesian_product<Ms...>, tensor_valued<R>, nothing, Opt> : h5_name<cartesian_product<Ms...>, matrix_valued, nothing, Opt> {};
+
   /// ---------------------------  evaluator ---------------------------------
 
   /**
@@ -146,7 +127,7 @@ namespace gfs {
   };
 
   // now the multi d evaluator itself.
-  template <typename Target, typename Opt, typename... Ms> struct evaluator<cartesian_product<Ms...>, Target, Opt> {
+  template <typename Target, typename Opt, typename... Ms> struct evaluator<cartesian_product<Ms...>, Target, nothing, Opt> {
    static constexpr int arity = sizeof...(Ms);
    mutable std::tuple<evaluator_fnt_on_mesh<Ms>...> evals;
 
diff --git a/triqs/gfs/refreq.hpp b/triqs/gfs/refreq.hpp
index d0715940..c11a96ee 100644
--- a/triqs/gfs/refreq.hpp
+++ b/triqs/gfs/refreq.hpp
@@ -35,18 +35,18 @@ namespace gfs {
   //using segment_mesh::segment_mesh;
  };
 
- namespace gfs_implementation {
+ // singularity
+ template <> struct gf_default_singularity<refreq, matrix_valued> {
+  using type = local::tail;
+ };
+ template <> struct gf_default_singularity<refreq, scalar_valued> {
+  using type = local::tail;
+ };
 
-  // singularity
-  template <typename Opt> struct singularity<refreq, matrix_valued, Opt> {
-   using type = local::tail;
-  };
-  template <typename Opt> struct singularity<refreq, scalar_valued, Opt> {
-   using type = local::tail;
-  };
+  namespace gfs_implementation {
 
   // h5 name
-  template <typename Opt> struct h5_name<refreq, matrix_valued, Opt> {
+  template <typename Singularity, typename Opt> struct h5_name<refreq, matrix_valued, Singularity, Opt> {
    static std::string invoke() { return "ReFreq"; }
   };
 
@@ -56,7 +56,7 @@ namespace gfs {
   struct evaluator_fnt_on_mesh<refreq> TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh,
                                                                              evaluator_grid_linear_interpolation);
 
-  template <typename Opt, typename Target> struct evaluator<refreq, Target, Opt> : evaluator_one_var<refreq> {};
+  template <typename Singularity, typename Opt, typename Target> struct evaluator<refreq, Target, Singularity, Opt> : evaluator_one_var<refreq> {};
 
 
   /// ---------------------------  data access  ---------------------------------
diff --git a/triqs/gfs/retime.hpp b/triqs/gfs/retime.hpp
index ed8a7b90..d0a0d264 100644
--- a/triqs/gfs/retime.hpp
+++ b/triqs/gfs/retime.hpp
@@ -35,18 +35,18 @@ namespace gfs {
   //using segment_mesh::segment_mesh;
  };
 
+ // singularity
+ template <> struct gf_default_singularity<retime, matrix_valued> {
+  using type = local::tail;
+ };
+ template <> struct gf_default_singularity<retime, scalar_valued> {
+  using type = local::tail;
+ };
+
  namespace gfs_implementation {
 
-  // singularity
-  template <typename Opt> struct singularity<retime, matrix_valued, Opt> {
-   using type = local::tail;
-  };
-  template <typename Opt> struct singularity<retime, scalar_valued, Opt> {
-   using type = local::tail;
-  };
-
   // h5 name
-  template <typename Opt> struct h5_name<retime, matrix_valued, Opt> {
+  template <typename Singularity, typename Opt> struct h5_name<retime, matrix_valued, Singularity, Opt> {
    static std::string invoke() { return "ReTime"; }
   };
 
@@ -55,7 +55,7 @@ namespace gfs {
   struct evaluator_fnt_on_mesh<retime> TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh,
                                                                              evaluator_grid_linear_interpolation);
 
-  template <typename Opt, typename Target> struct evaluator<retime, Target, Opt> : evaluator_one_var<retime> {};
+  template <typename Singularity, typename Opt, typename Target> struct evaluator<retime, Target, Singularity, Opt> : evaluator_one_var<retime> {};
 
   /// ---------------------------  data access  ---------------------------------
   template <typename Opt> struct data_proxy<retime, matrix_valued, Opt> : data_proxy_array<std::complex<double>, 3> {};
diff --git a/triqs/gfs/tools.hpp b/triqs/gfs/tools.hpp
index b29f9121..f9a58e63 100644
--- a/triqs/gfs/tools.hpp
+++ b/triqs/gfs/tools.hpp
@@ -186,6 +186,7 @@ namespace gfs {
   bool is_empty() const { return false;}
  };
 
+ template<int ... pos, typename ...T> nothing partial_eval(nothing, T&&...) { return {};}
  inline nothing transpose(nothing) { return {};}
  inline nothing inverse(nothing) { return {};}
  inline nothing conj(nothing) { return {};}
diff --git a/triqs/lattice/bz_mesh.hpp b/triqs/lattice/bz_mesh.hpp
index 585262ac..bca6c980 100644
--- a/triqs/lattice/bz_mesh.hpp
+++ b/triqs/lattice/bz_mesh.hpp
@@ -30,6 +30,7 @@ namespace lattice {
 
   using domain_t = brillouin_zone;
   using index_t = long;
+  using linear_index_t = long;
   using domain_pt_t = typename domain_t::point_t;
 
   bz_mesh() = default;
@@ -37,7 +38,12 @@ namespace lattice {
   bz_mesh(brillouin_zone const &bz, std::vector<k_t> k_pt_stack) : bz(bz), k_pt_stack(std::move(k_pt_stack)) {}
 
   domain_t const &domain() const { return bz; }
-  long size() const { return k_pt_stack.size(); }
+  size_t size() const { return k_pt_stack.size(); }
+
+  ///
+  utility::mini_vector<size_t, 1> size_of_components() const {
+   return {size()};
+  }
 
   /// Conversions point <-> index <-> linear_index
   domain_pt_t const &index_to_point(index_t i) const { return k_pt_stack[i]; }
diff --git a/triqs/utility/mini_vector.hpp b/triqs/utility/mini_vector.hpp
index ebe72e8d..816925e0 100644
--- a/triqs/utility/mini_vector.hpp
+++ b/triqs/utility/mini_vector.hpp
@@ -38,7 +38,7 @@ namespace triqs { namespace utility {
 
  template <typename T, int Rank> 
   class mini_vector { 
-   T _data[Rank];
+   T _data[(Rank!=0 ? Rank : 1)];
    friend class boost::serialization::access;
    template<class Archive> void serialize(Archive & ar, const unsigned int version) { ar & TRIQS_MAKE_NVP("_data",_data); }
    void init() { for (int i=0;i<Rank; ++i) _data[i] = 0;}
@@ -126,6 +126,16 @@ namespace triqs { namespace utility {
    friend std::stringstream & operator << ( std::stringstream & out, mini_vector const & v ) { out<<v.to_string(); return out;}
   }; // class mini_vector
 
+ // ------------ specialization for size 0 ------------- 
+  template <typename T> class mini_vector<T, 0> {
+   T _data[1];
+   public:
+   T & operator[](size_t i) { return _data[i];}
+   const T & operator[](size_t i) const { return _data[i];}
+   };
+
+ // ------------- Comparison -------------------------------------- 
+ 
  template <typename T, int R> 
   bool operator ==(mini_vector<T,R> const & v1, mini_vector<T,R> const & v2) {
    for (int i=0;i<R; ++i) { if (v1[i]!=v2[i]) return false;}
@@ -134,14 +144,16 @@ namespace triqs { namespace utility {
 
  template <typename T, int R> bool operator !=(mini_vector<T,R> const & v1, mini_vector<T,R> const & v2) { return (!(v1==v2));}
 
- template <typename T, int R1, int R2> 
-  mini_vector<T, R1+R2> join (mini_vector<T,R1> const & v1, mini_vector<T,R2> const & v2) {
-   mini_vector<T, R1+R2> res;
+ // ------------- join  -------------------------------------- 
+ template <typename T1, typename T2,  int R1, int R2> 
+  mini_vector<T1, R1+R2> join (mini_vector<T1,R1> const & v1, mini_vector<T2,R2> const & v2) {
+   mini_vector<T1, R1+R2> res;
    for (int i=0;i<R1; ++i)  res[i]=v1[i];
    for (int i=0;i<R2; ++i)  res[R1+i]=v2[i];
    return res;
   }
 
+ // ------------- dot -------------------------------------- 
  template <typename T1, typename T2, int Rank> 
   T1 dot_product(mini_vector<T1,Rank> const & v1, mini_vector<T2,Rank> const & v2) {
    T1 res=0;
@@ -149,25 +161,21 @@ namespace triqs { namespace utility {
    return res;
   }
 
- struct tuple_to_mini_vector_aux { template<typename M, typename V> V * operator()(M const & m,  V * v) { *v = m; return ++v;}};
+#ifndef TRIQS_C11
+  // ------------- transform a tuple into a minivector --------------------------------------
 
- // change : the first version crash clang 3.3, but not svn version.
- // must be a bug, corrected since then
-/* 
- template<typename T, typename ... U> 
-  mini_vector<T,sizeof...(U)> tuple_to_mini_vector(std::tuple<U...> const & t) { 
-   mini_vector<T,sizeof...(U)> res;
-   triqs::tuple::fold(tuple_to_mini_vector_aux(),t,&res[0]); 
-   return res;
+  template <typename T, typename TU> mini_vector<T, std::tuple_size<TU>::value> tuple_to_mini_vector(TU const &t) {
+   return tuple::apply_construct_parenthesis<mini_vector<T, std::tuple_size<TU>::value>>(t);
   }
-*/
- template<typename T, typename TU> 
-  mini_vector<T,std::tuple_size<TU>::value> tuple_to_mini_vector(TU const & t) { 
-   mini_vector<T,std::tuple_size<TU>::value> res;
-   triqs::tuple::fold(tuple_to_mini_vector_aux(),t,&res[0]); 
-   return res;
-  }
-
-}}//namespace triqs::arrays 
+#endif
+
+}}//namespace triqs::arrays 
+
+namespace std { // overload std::get to work with it
+
+template <int i, typename T, int R> AUTO_DECL get(triqs::utility::mini_vector<T, R> &v) RETURN(v[i]);
+template <int i, typename T, int R> AUTO_DECL get(triqs::utility::mini_vector<T, R> const &v) RETURN(v[i]);
+}
+
 #endif
 
diff --git a/triqs/utility/tuple_tools.hpp b/triqs/utility/tuple_tools.hpp
index 4ce8061f..9d20401f 100644
--- a/triqs/utility/tuple_tools.hpp
+++ b/triqs/utility/tuple_tools.hpp
@@ -68,6 +68,25 @@ namespace std {
 
 namespace triqs { namespace tuple {
 
+ /// Repeat an element
+ template<typename T, int R> struct make_tuple_repeat_impl;
+ 
+ template <typename T> struct make_tuple_repeat_impl<T,1> {
+  static std::tuple<T> invoke(T&&x) { return {x};}
+ };
+ template <typename T> struct make_tuple_repeat_impl<T,2> {
+  static std::tuple<T,T> invoke(T&&x) { return {x,x};}
+ };
+ template <typename T> struct make_tuple_repeat_impl<T,3> {
+  static std::tuple<T,T,T> invoke(T&&x) { return {x,x,x};}
+ };
+
+ template <typename T> struct make_tuple_repeat_impl<T, 4> {
+  static std::tuple<T,T,T,T> invoke(T &&x) { return {x, x, x, x}; }
+ };
+
+ template <int R, typename T> auto make_tuple_repeat(T &&x) { return make_tuple_repeat_impl<T, R>::invoke(std::forward<T>(x)); }
+
  /// _get_seq<T>() : from a tuple T, return the index sequence of the tuple length
  template <typename T> std14::make_index_sequence<std::tuple_size<std14::decay_t<T>>::value> _get_seq() {
   return {};
@@ -120,6 +139,18 @@ namespace triqs { namespace tuple {
  template <typename C, typename T>
  AUTO_DECL apply_construct(T &&t) RETURN(apply_construct_impl<C>(std::forward<T>(t), _get_seq<T>()));
 
+ /**
+  * apply_construct_parenthesis<C>(t)
+  * C : a class
+  * t a tuple
+  * Returns : C(t0, t1, ....)
+  */
+ template <typename C, typename T, size_t... Is>
+ AUTO_DECL apply_construct_parenthesis_impl(T &&t, std14::index_sequence<Is...>) RETURN(C(std::get<Is>(std::forward<T>(t))...));
+
+ template <typename C, typename T>
+ AUTO_DECL apply_construct_parenthesis(T &&t) RETURN(apply_construct_parenthesis_impl<C>(std::forward<T>(t), _get_seq<T>()));
+
  /**
   * called_on_tuple(f)
   * f : a callable object
@@ -455,7 +486,7 @@ namespace triqs { namespace tuple {
  /*
  * print a tuple
  */
- inline void _triqs_print_tuple_impl(std::ostream &os) {}
+ /*inline void _triqs_print_tuple_impl(std::ostream &os) {}
  template <typename T0, typename... T> void _triqs_print_tuple_impl(std::ostream &os, T0 const &x0, T const &... x) {
   os << x0;
   if (sizeof...(T) > 0) os << ',';
@@ -466,6 +497,8 @@ namespace triqs { namespace tuple {
  void _triqs_print_tuple(std::ostream &os, std::tuple<T...> const &t, std14::index_sequence<Is...>) {
   _triqs_print_tuple_impl(os, std::get<Is>(t)...);
  }
+*/
+
 }
 }
 
@@ -473,7 +506,8 @@ namespace std {
 
 template <typename... T> std::ostream &operator<<(std::ostream &os, std::tuple<T...> const &t) {
  os << "(";
- triqs::tuple::_triqs_print_tuple(os, t, std14::make_index_sequence<sizeof...(T)>());
+ triqs::tuple::for_each(t, [&os, c=0](auto & x) mutable {if (c++) os << ','; os << x;});
+ //triqs::tuple::_triqs_print_tuple(os, t, std14::make_index_sequence<sizeof...(T)>());
  return os << ")";
 }