gf.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2023 Simons Foundation
//
// This program 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.
//
// This program 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 may obtain a copy of the License at
//     https://www.gnu.org/licenses/gpl-3.0.txt
//
// Authors: Michel Ferrero, Olivier Parcollet, Nils Wentzell

 
#pragma once
 
#include "./defs.hpp"
#include "./targets.hpp"
#include "../../utility/concept_tools.hpp"
#include "../../utility/exceptions.hpp" // IWYU pragma: keep (used by ./_gf_view_common.hpp)
#include "../../utility/macros.hpp"     // IWYU pragma: keep (used by ./_gf_view_common.hpp)
 
#include <itertools/itertools.hpp>
#include <mpi/mpi.hpp>
#include <nda/nda.hpp>
 
#include <array>
#include <concepts>
#include <numeric> // IWYU pragma: keep (used by ./_gf_view_common.hpp)
#include <ostream> // IWYU pragma: keep (used by ./_gf_view_common.hpp)
#include <string>  // IWYU pragma: keep (used by ./_gf_view_common.hpp)
#include <tuple>   // IWYU pragma: keep (used by ./_gf_view_common.hpp)
#include <type_traits>
#include <utility>
 
namespace triqs::gfs {
 
  // Elevate certain types into the triqs::gfs namespace.
  using nda::C_layout;
  using nda::C_stride_layout;
 
  // Forward declarations.
  template <Mesh M, typename Target = matrix_valued, typename Layout = C_layout> class gf;
  template <Mesh M, typename Target = matrix_valued, typename Layout = C_stride_layout> class gf_view;
  template <Mesh M, typename Target = matrix_valued, typename Layout = C_stride_layout> class gf_const_view;

 
  // Trait and concept tag marking a type as a Green's function.
  TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT(GreenFunction);


  template <template <typename...> class TMPLT, typename T> struct is_instantiation_of : std::false_type {};
 
  // Specialization of is_instantiation_of for a matching instantiation.
  template <template <typename...> class TMPLT, typename... U> struct is_instantiation_of<TMPLT, TMPLT<U...>> : std::true_type {};

  template <template <typename...> class gf, typename T>
  inline constexpr bool is_instantiation_of_v = is_instantiation_of<gf, std::decay_t<T>>::value;
 
  //template <typename G> using mesh_t_of = std::decay_t<decltype(std::declval<G>().mesh())>;

   * 
   * When a second template argument `M` is supplied, it additionally requires the mesh type of `G` to be `M`.
   *
   * @deprecated In favour of the triqs::gfs::MemoryGf concept.
   *
   * @tparam G Type to test.
   * @tparam M Optional mesh type to require.
   */
  template <typename G, typename M = void> inline constexpr bool is_gf_v = false;
 
  // Specialization of is_gf_v for cvref-qualified types, which decays G first.
  template <typename G, typename M>
    requires(!std::is_same_v<G, std::remove_cvref_t<G>>)
  inline constexpr bool is_gf_v<G, M> = is_gf_v<std::remove_cvref_t<G>, M>;
 
  // Specialization of is_gf_v without a mesh constraint.
  template <typename G>
  inline constexpr bool is_gf_v<G, void> =
     is_instantiation_of_v<gf, G> or is_instantiation_of_v<gf_view, G> or is_instantiation_of_v<gf_const_view, G>;
 
  // Specialization of is_gf_v requiring the mesh type of G.
  template <typename G> inline constexpr bool is_gf_v<G, typename std::remove_cvref_t<G>::mesh_t> = is_gf_v<G, void>;


  template <typename G, typename M = typename std::remove_cvref_t<G>::mesh_t>
  concept MemoryGf = mesh::Mesh<M> and requires(G g) {
    { g.data() } -> nda::MemoryArray;
    requires std::same_as<std::decay_t<decltype(g.mesh())>, M>;
  };
 
  // Forward declarations with default arguments. The corresponding friend declarations inside `gf` and
  // the definitions in `./mpi.hpp` may then re-declare these templates without re-introducing defaults.
  template <MemoryGf G> void mpi_broadcast(G &&g, mpi::communicator c = {}, int root = 0);
  template <MemoryGf G1, MemoryGf G2>
  void mpi_reduce_into(G1 const &g_in, G2 &&g_out, mpi::communicator c = {}, int root = 0, bool all = false, MPI_Op op = MPI_SUM);


  template <typename G> using get_target_t = typename std::decay_t<G>::target_t;

  template <int N = 0, MemoryGf G> auto const &get_mesh(G const &g) {
    if constexpr (mesh::is_product<typename G::mesh_t>)
      return std::get<N>(g.mesh());
    else
      return g.mesh();
  }

  template <MemoryGf G> using target_value_t = decltype(std::declval<G>()[std::declval<typename G::mesh_point_t>()]);

  /** @} */
 
  // Tags used internally to select implementation-detail constructors.
  struct impl_tag {};
  struct impl_tag2 {};

   * @brief The owning Green's function container.
   *
   * @details triqs::gfs::gf is the central, value-semantic container of the Green's function module. It owns a mesh
   * and a data array, and represents a function defined on the mesh whose value at each mesh point is fixed by the
   * `Target` type (see triqs::gfs::scalar_valued, triqs::gfs::matrix_valued, triqs::gfs::tensor_valued).
   *
   * The data array has rank `arity + Target::rank`, where `arity` is the number of mesh dimensions. The function can
   * be accessed at a mesh point with `operator[]` and evaluated at an arbitrary point with `operator()` (using the
   * interpolation scheme of the mesh). 
   * 
   * Non-owning views are provided by triqs::gfs::gf_view and triqs::gfs::gf_const_view.
   *
   * @tparam M Mesh type, modeling triqs::mesh::Mesh; fixes the domain of definition.
   * @tparam Target Target type; fixes the value stored at each mesh point.
   * @tparam Layout Memory layout policy of the data array.
   */
  template <Mesh M, typename Target, typename Layout> class gf : TRIQS_CONCEPT_TAG_NAME(GreenFunction) {
 
    static_assert(not std::is_same_v<M, triqs::lattice::brillouin_zone>,
                  "Since TRIQS 2.3, brillouin_zone is replaced by mesh::brzone as a mesh name. Cf Doc, changelog");
 
    using this_t = gf<M, Target, Layout>; // used in common code
 
    public:
    static constexpr bool is_view = false;

    static constexpr bool is_const = false;

    using mutable_view_type = gf_view<M, Target, typename Layout::with_lowest_guarantee_t>;

    using const_view_type = gf_const_view<M, Target, typename Layout::with_lowest_guarantee_t>;

    using view_type = gf_view<M, Target, typename Layout::with_lowest_guarantee_t>;

    using regular_type = gf<M, Target, Layout>;

    using real_t = gf<M, typename Target::real_t, Layout>;

    using complex_t = gf<M, typename Target::complex_t, Layout>;

    using target_t = Target;

    using mesh_t = M;

    using mesh_point_t = typename mesh_t::mesh_point_t;

    using mesh_index_t = typename mesh_t::index_t;

    using scalar_t = typename Target::scalar_t;

    static constexpr int arity = n_variables<M>;

    static constexpr int target_rank = Target::rank;

    static constexpr int data_rank = arity + Target::rank;

    using data_t = nda::basic_array<scalar_t, data_rank, Layout, 'A', nda::heap<>>;

    using target_shape_t = std::array<long, Target::rank>;

    struct target_and_shape_t {
      target_shape_t _shape;

      using target_t = Target;

      target_shape_t const &shape() const { return _shape; }
    };
 
    // ------------- Accessors -----------------------------

    mesh_t const &mesh() const { return _mesh; }

    data_t &data() & { return _data; }

    data_t const &data() const & { return _data; }

    data_t data() && { return std::move(_data); }

    auto const &data_shape() const { return _data.shape(); }

    target_and_shape_t target() const { return target_and_shape_t{stdutil::front_mpop<arity>(_data.shape())}; } // drop arity dims

    std::array<long, Target::rank> target_shape() const { return target().shape(); } // drop arity dims

    auto target_indices() const { return itertools::product_range(target().shape()); }
 
    private:
    mesh_t _mesh;
    data_t _data;
 
    // -------------------------------- impl. details common to all classes -----------------------------------------------
 
    public:
    gf() = default; // {}

    gf(gf const &x) = default;

    gf(gf &&) = default;
 
    private:
    void swap_impl(gf &b) noexcept {
      using std::swap;
      swap(this->_mesh, b._mesh);
      swap(this->_data, b._data);
    }
 
    private:
    // Compute the shape of the data array from the mesh and the target shape.
    static auto make_data_shape(mesh_t const &m, target_shape_t const &shape) {
      if constexpr (mesh::is_product<mesh_t>)
        return stdutil::join(m.size_of_components(), shape);
      else
        return stdutil::front_append(shape, m.size());
    }
 
    public:
    gf(mesh_t m, data_t dat) : _mesh(std::move(m)), _data(std::move(dat)) {}

    gf(mesh_t m, target_shape_t shape = {}) : _mesh(std::move(m)), _data(make_data_shape(_mesh, shape)) {}

    explicit gf(gf_view<M, Target> const &g) : _mesh(g.mesh()), _data(g.data()) {}

    explicit gf(gf_const_view<M, Target> const &g) : _mesh(g.mesh()), _data(g.data()) {}

    template <typename G>
    explicit gf(G const &g)
      requires(GreenFunction<G>::value and std::is_same_v<mesh_t, typename G::mesh_t>)
       : gf() {
      *this = g;
    } // explicit is very important here.
    // TODO: We would like to refine this, G should have the same mesh, target, at least ...
 
    // ---------------  Operator = --------------------

    gf &operator=(gf const &rhs) = default;

    gf &operator=(gf &&rhs) noexcept {
      this->swap_impl(rhs);
      return *this;
    }

    template <typename RHS>
    gf &operator=(RHS &&rhs) // NOLINT
      requires(GreenFunction<RHS>::value and not std::is_same_v<std::decay_t<RHS>, gf>)
    {
      _mesh = rhs.mesh();
      _data.resize(rhs.data_shape());
      for (auto w : _mesh) (*this)[w] = rhs[w];
      return *this;
    }
 
    // other = late, cf MPI
 
    public:
    // ------------- apply_on_data -----------------------------

    template <typename Fdata> auto apply_on_data(Fdata &&fd) {
      auto d2    = std::forward<Fdata>(fd)(_data);
      using t2   = target_from_array<decltype(d2), arity>;
      using gv_t = gf_view<M, t2>;
      return gv_t{mesh(), d2};
    }

    template <typename Fdata> auto apply_on_data(Fdata &&fd) const {
      auto d2    = std::forward<Fdata>(fd)(_data);
      using t2   = target_from_array<decltype(d2), arity>;
      using gv_t = gf_const_view<M, t2>;
      return gv_t{mesh(), d2};
    }
 
    // Friend declarations (hidden from doxygen; documented as free functions in gf/mpi.hpp).
    template <MemoryGf G> friend void mpi_broadcast(G &&, mpi::communicator, int root);
    template <MemoryGf G1, MemoryGf G2> friend void mpi_reduce_into(G1 const &, G2 &&, mpi::communicator, int, bool, MPI_Op);
 
    // Common code for gf, gf_view, gf_const_view
#include "./_gf_view_common.hpp"
  };
 
  /*------------------------------------------------------------------------
   *   Deduction guides
   *-----------------------------------------------------------------------*/

  template <Mesh M, nda::MemoryArray DataArray> gf(M, DataArray) -> gf<M, target_from_array<DataArray, n_variables<M>>>;

  template <Mesh M, std::integral I, size_t R> gf(M, std::array<I, R>) -> gf<M, typename _target_from_type_rank<dcomplex, R>::type>;

  template <Mesh M> gf(M) -> gf<M, scalar_valued>;
 
  // Forward declare gf_expr
  template <typename Tag, typename L, typename R> struct gf_expr;

  template <typename Tag, typename L, typename R>
  gf(gf_expr<Tag, L, R> &&) -> gf<typename gf_expr<Tag, L, R>::mesh_t, typename gf_expr<Tag, L, R>::target_t>;

  template <typename Tag, typename L, typename R>
  gf(gf_expr<Tag, L, R> const &) -> gf<typename gf_expr<Tag, L, R>::mesh_t, typename gf_expr<Tag, L, R>::target_t>;
 
} // namespace triqs::gfs