Unverified Commit b874cf7c authored by Peter Doak's avatar Peter Doak Committed by GitHub
Browse files

Merge pull request #147 from CompFUSE/ct_int-configuration

Ct int configuration
parents e0eed7eb ba9b17bf
Loading
Loading
Loading
Loading
+1 −0
Original line number Diff line number Diff line
@@ -120,6 +120,7 @@ set(DCA_LIBS
  timer
  cmake_options
  ctaux
  ctint
  ss_ct_hyb
  function_transform
  gaussian_quadrature
+19 −2
Original line number Diff line number Diff line
@@ -26,9 +26,16 @@ namespace util {
// dca::linalg::util::

template <typename T>
class PinnedAllocator : public std::allocator<T> {
class PinnedAllocator {
public:
  T* allocate(std::size_t n) {
  PinnedAllocator() = default;

  using size_type = std::size_t;
  using pointer = T*;
  using const_pointer = const T*;
  using value_type = T;

  T* allocate(std::size_t n, const void* /*hint*/ = nullptr) {
    if (n == 0)
      return nullptr;
    T* ptr;
@@ -49,8 +56,18 @@ public:
    }
    ptr = nullptr;
  }

};

// These are part of the requirements for a C++ Allocator however these are insufficient
// and they do not appear to be relevant to our codebase yet.
// They are part of what's needed to do a std::move on a PinnedAllocator backed object
// and avoid deallocation and reallocation.
template <class T, class U>
bool operator==(const PinnedAllocator<T>&, const PinnedAllocator<U>&) { return true; }
template <class T, class U>
bool operator!=(const PinnedAllocator<T>&, const PinnedAllocator<U>&) { return false; }

}  // util
}  // linalg
}  // dca
+3 −0
Original line number Diff line number Diff line
@@ -221,6 +221,9 @@ public: // Optional members getters.
    return G4_err_;
  }

  // The non density-density Hamiltonian is given by:
  // H = \sum(nu1, nu2, nu3, nu4, r1, r2) c^+(nu1, r1) c(nu2, r1) c^+(nu3, r2) c(nu4, r2) *
  //     non_density_interactions_(nu1, nu2, nu3, nu4, r1 - r2)
  // Note: this contribution to the Hamiltonian is not double counted.
  auto& get_non_density_interactions() {
    if (not non_density_interactions_)
+176 −0
Original line number Diff line number Diff line
// Copyright (C) 2018 ETH Zurich
// Copyright (C) 2018 UT-Battelle, LLC
// All rights reserved.
//
// See LICENSE.txt for terms of usage.
// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications.
//
// Authors: Giovanni Balduzzi (gbalduzz@phys.ethz.ch)
//
// This class process the information contained in the SolverConfiguration class into a
// representation of the M matrix.

#ifndef DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_MATRIX_CONFIGURATION_HPP
#define DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_MATRIX_CONFIGURATION_HPP

#include <array>
#include <vector>

#include "dca/linalg/device_type.hpp"
#include "dca/phys/dca_step/cluster_solver/ctint/structs/ct_int_sector.hpp"
#include "dca/phys/dca_step/cluster_solver/ctint/structs/interaction_vertices.hpp"

namespace dca {
namespace phys {
namespace solver {
namespace ctint {
// dca::phys::solver::ctint::

// Expansion term.
struct Vertex {
  bool aux_spin;
  ushort interaction_id;
  std::uint64_t tag;
  double tau;

  bool operator==(const Vertex& b) const {
    return aux_spin == b.aux_spin && interaction_id == b.interaction_id && tau == b.tau;
  }
};

class MatrixConfiguration {
public:
  MatrixConfiguration() = default;
  MatrixConfiguration(const MatrixConfiguration& rhs) = default;
  inline MatrixConfiguration& operator=(const MatrixConfiguration& rhs);
  inline MatrixConfiguration& operator=(MatrixConfiguration&& rhs);

  inline void swapSectorLabels(int a, int b, int s);

  const Sector& getSector(int s) const {
    return sectors_[s];
  }

  std::size_t size() const {
    return (size(0) + size(1)) / 2;
  }

  std::size_t size(int s) const {
    return sectors_[s].size();
  }

  const std::array<Sector, 2>& get_sectors() const {
    return sectors_;
  }

  bool operator==(const MatrixConfiguration& rhs) const {
    return sectors_ == rhs.sectors_;
  }

  // TODO try to return on the stack or preallocated mem.
  inline std::vector<int> findIndices(const uint64_t tag, const int s) const;

protected:
  inline MatrixConfiguration(const InteractionVertices* H_int, int bands);

  inline std::array<int, 2> addVertex(const Vertex& v);

  inline void pop(ushort idx_up, ushort idx_down);

  const auto& getEntries(const int s) const {
    return sectors_[s].entries_;
  }
  auto& getEntries(const int s) {
    return sectors_[s].entries_;
  }

protected:
  const InteractionVertices* H_int_ = nullptr;
  const int n_bands_ = -1;
  std::array<Sector, 2> sectors_;
};

MatrixConfiguration::MatrixConfiguration(const InteractionVertices* H_int, const int bands)
    : H_int_(H_int), n_bands_(bands), sectors_{Sector(), Sector()} {}

MatrixConfiguration& MatrixConfiguration::operator=(const MatrixConfiguration& rhs) {
  sectors_ = rhs.sectors_;
  return *this;
}

MatrixConfiguration& MatrixConfiguration::operator=(MatrixConfiguration&& rhs) {
  sectors_ = std::move(rhs.sectors_);
  return *this;
}

std::array<int, 2> MatrixConfiguration::addVertex(const Vertex& v) {
  auto spin = [=](const int nu) { return nu >= n_bands_; };
  auto band = [=](const int nu) -> ushort { return nu - n_bands_ * spin(nu); };

  auto field_type = [&](const Vertex& v, const int leg) -> short {
    const short sign = v.aux_spin ? 1 : -1;
    const InteractionElement& elem = (*H_int_)[v.interaction_id];
    if (elem.partners_id.size())
      return leg == 1 ? -3 * sign : 3 * sign;  // non density-density.
    else if (elem.w > 0)
      return leg == 1 ? -1 * sign : 1 * sign;  // positive dd interaction.
    else
      return 2 * sign;  // negative dd interaction.
  };

  std::array<int, 2> sizes{0, 0};
  const auto& nu = (*H_int_)[v.interaction_id].nu;
  const auto& r = (*H_int_)[v.interaction_id].r;

  for (ushort leg = 0; leg < 2; ++leg) {
    assert(spin(nu[0 + 2 * leg]) == spin(nu[1 + 2 * leg]));
    const short s = spin(nu[0 + 2 * leg]);
    Sector& sector = sectors_[s];
    ++sizes[s];
    sector.entries_.emplace_back(details::SectorEntry{band(nu[0 + 2 * leg]), r[0 + 2 * leg],
                                                      band(nu[1 + 2 * leg]), r[1 + 2 * leg], v.tau,
                                                      field_type(v, leg)});
    sector.tags_.push_back(v.tag);
  }
  assert(sectors_[0].entries_.size() == sectors_[0].tags_.size());
  assert(sectors_[1].entries_.size() == sectors_[1].tags_.size());
  return sizes;
}

void MatrixConfiguration::pop(ushort n_up, ushort n_down) {
  assert(n_up <= sectors_[0].size() and n_down <= sectors_[1].size());
  sectors_[0].entries_.erase(sectors_[0].entries_.end() - n_up, sectors_[0].entries_.end());
  sectors_[0].tags_.erase(sectors_[0].tags_.end() - n_up, sectors_[0].tags_.end());
  sectors_[1].entries_.erase(sectors_[1].entries_.end() - n_down, sectors_[1].entries_.end());
  sectors_[1].tags_.erase(sectors_[1].tags_.end() - n_down, sectors_[1].tags_.end());
}

std::vector<int> MatrixConfiguration::findIndices(const uint64_t tag, const int s) const {
  std::vector<int> indices;

  auto start = sectors_[s].tags_.begin();
  const auto end = sectors_[s].tags_.end();
  while (true) {
    start = std::find(start, end, tag);
    if (start != end) {
      indices.push_back(start - sectors_[s].tags_.begin());
      ++start;
    }
    else
      break;
  }

  return indices;
}

void MatrixConfiguration::swapSectorLabels(const int a, const int b, const int s) {
  std::swap(sectors_[s].entries_[a], sectors_[s].entries_[b]);
  std::swap(sectors_[s].tags_[a], sectors_[s].tags_[b]);
}

}  // namespace ctint
}  // namespace solver
}  // namespace phys
}  // namespace dca

#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_MATRIX_CONFIGURATION_HPP
+89 −0
Original line number Diff line number Diff line
// Copyright (C) 2018 ETH Zurich
// Copyright (C) 2018 UT-Battelle, LLC
// All rights reserved.
//
// See LICENSE.txt for terms of usage.
// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications.
//
// Author: Giovanni Balduzzi (gbalduzz@phys.ethz.ch)
//
// This structure organize the data of MatrixConfiguration for each sector in {up, down}

#ifndef DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_SECTOR_HPP
#define DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_SECTOR_HPP

#include <array>
#include <cassert>
#include <vector>

#include "dca/linalg/util/allocators/vectors_typedefs.hpp"
#include "dca/phys/dca_step/cluster_solver/ctint/structs/interaction_vertices.hpp"
#include "dca/phys/dca_step/cluster_solver/ctint/structs/sector_entry.hpp"

namespace dca {
namespace phys {
namespace solver {
namespace ctint {
// dca::phys::solver::ctint::

// TODO: replace by vector
class Sector {
public:
  Sector() : entries_(){};

  friend class SolverConfiguration;
  friend class MatrixConfiguration;
  friend class DeviceConfigurationManager;

  ushort size() const {
    return entries_.size();
  }

  const details::SectorEntry& operator[](const std::size_t idx) const {
    assert(idx < size());
    return entries_[idx];
  }

  inline ushort getLeftB(const int i) const {
    return entries_[i].b_left_;
  }

  inline ushort getRightB(const int i) const {
    return entries_[i].b_right_;
  }

  inline ushort getLeftR(const int i) const {
    return entries_[i].r_left_;
  }

  inline ushort getRightR(const int i) const {
    return entries_[i].r_right_;
  }

  inline double getTau(const int i) const {
    return entries_[i].tau_;
  }

  inline short getAuxFieldType(int i) const {
    return entries_[i].aux_field_type_;
  }

  auto getTag(int i) const {
    return tags_[i];
  }

  bool operator==(const Sector& rhs) const {
    return entries_ == rhs.entries_;
  }

private:
  linalg::util::HostVector<details::SectorEntry> entries_;
  std::vector<std::uint64_t> tags_;
};

}  // namespace ctint
}  // namespace solver
}  // namespace phys
}  // namespace dca

#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_SECTOR_HPP
Loading