matsubara.hh

Go to the documentation of this file.

#pragma once
 
// DiFfRG
#include <DiFfRG/common/kokkos.hh>
#include <DiFfRG/common/math.hh>
 
// C++ standard library
#include <vector>
 
namespace DiFfRG
{
  template <typename NT> class MatsubaraQuadrature
  {
  public:
    static int predict_size(const NT T, const NT typical_E = 1., const double precision_factor = 1.,
                            const int step = 2);
 
    MatsubaraQuadrature(const NT T, const NT typical_E = 1., const int step = 2, const int min_size = 8,
                        const int max_size = 128, const int vacuum_quad_size = 64, const double precision_factor = 1);
 
    MatsubaraQuadrature();
 
    void reinit(const NT T, const NT typical_E = 1., const int step = 2, const int min_size = 8,
                const int max_size = 128, const int vacuum_quad_size = 64, const double precision_factor = 1);
 
    void reinit_with_size(const int size, const NT T, const NT typical_E = 1.);
 
    static int modes_below(const NT T, const NT freq_cutoff);
 
    void reinit_exact_sum(const NT T, const NT freq_cutoff);
 
    void reinit_finite_interval(const NT cutoff, const Kokkos::View<const NT *, CPU_memory> gl_nodes,
                                const Kokkos::View<const NT *, CPU_memory> gl_weights);
 
    size_t size() const;
 
    size_t sum_size() const;
 
    NT get_T() const;
 
    NT get_typical_E() const;
 
    template <typename F> auto sum(const F &f) const
    {
      auto sum = T * f(static_cast<NT>(0));
      for (int i = 0; i < m_size; ++i)
        sum += host_weights[i] * (f(host_nodes[i]) + f(-host_nodes[i]));
      return sum;
    }
 
    template <typename MemorySpace> Kokkos::View<const NT *, MemorySpace> nodes() const
    {
      return Kokkos::subview(raw_nodes<MemorySpace>(), Kokkos::make_pair(size_t(0), size_t(m_size)));
    }
 
    template <typename MemorySpace> Kokkos::View<const NT *, MemorySpace> weights() const
    {
      return Kokkos::subview(raw_weights<MemorySpace>(), Kokkos::make_pair(size_t(0), size_t(m_size)));
    }
 
    template <typename MemorySpace> Kokkos::View<const NT *, MemorySpace> sum_nodes() const
    {
      return Kokkos::subview(raw_nodes<MemorySpace>(), Kokkos::make_pair(size_t(0), sum_size()));
    }
 
    template <typename MemorySpace> Kokkos::View<const NT *, MemorySpace> sum_weights() const
    {
      return Kokkos::subview(raw_weights<MemorySpace>(), Kokkos::make_pair(size_t(0), sum_size()));
    }
 
  private:
    template <typename MemorySpace> Kokkos::View<const NT *, MemorySpace> raw_nodes() const
    {
      if constexpr (std::is_same_v<MemorySpace, Kokkos::DefaultExecutionSpace::memory_space>) {
        return device_nodes;
      } else if constexpr (std::is_same_v<MemorySpace, Kokkos::DefaultHostExecutionSpace::memory_space>) {
        return host_nodes;
      } else {
        throw std::runtime_error("Invalid memory space");
      }
    }
 
    template <typename MemorySpace> Kokkos::View<const NT *, MemorySpace> raw_weights() const
    {
      if constexpr (std::is_same_v<MemorySpace, Kokkos::DefaultExecutionSpace::memory_space>) {
        return device_weights;
      } else if constexpr (std::is_same_v<MemorySpace, Kokkos::DefaultHostExecutionSpace::memory_space>) {
        return host_weights;
      } else {
        throw std::runtime_error("Invalid memory space");
      }
    }
 
    NT T, typical_E;
 
    Kokkos::View<NT *, GPU_memory> device_nodes;
    Kokkos::View<NT *, GPU_memory> device_weights;
 
    Kokkos::View<NT *, CPU_memory> host_nodes;
    Kokkos::View<NT *, CPU_memory> host_weights;
 
    int m_size;
 
    void build_monien();
 
    void reinit_0();
 
    void write_data(const std::vector<NT> &x, const std::vector<NT> &w);
 
    int vacuum_quad_size;
    double precision_factor;
  };
} // namespace DiFfRG