Commit d3564261 authored by Doak, Peter W.'s avatar Doak, Peter W.
Browse files

Merge branch 'ct_int-configuration' of github.com:CompFUSE/DCA into ct_int-configuration

parents f991ac0e 6ffbf390

include/dca/phys/dca_step/cluster_solver/ctint/structs/ct_int_matrix_configuration.hpp

+1 −1
Original line number Diff line number Diff line
@@ -110,7 +110,7 @@ std::array<int, 2> MatrixConfiguration::addVertex(const Vertex& v) { 
  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.partner_id != -1)

    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.

include/dca/phys/dca_step/cluster_solver/ctint/structs/device_configuration_manager.hpp

+1 −1
Original line number Diff line number Diff line
@@ -17,7 +17,7 @@ 
#endif



#include "dca/linalg/device_type.hpp"

#include "dca/phys/dca_step/cluster_solver/ctint/structs/ct_int_configuration.hpp"

#include "dca/phys/dca_step/cluster_solver/ctint/structs/solver_configuration.hpp"

#include "dca/phys/dca_step/cluster_solver/ctint/structs/device_configuration.hpp"

#include "dca/phys/dca_step/cluster_solver/ctint/structs/sector_entry.hpp"

#include "dca/util/cuda_definitions.hpp"

include/dca/phys/dca_step/cluster_solver/ctint/structs/interaction_vertices.hpp

+71 −9
Original line number Diff line number Diff line
@@ -28,10 +28,11 @@ namespace ctint { 
// Represent a matrix element W(i,j,k,l) of the interaction hamiltonian.

// Each index represent a cluster position and a spin-band index.

struct InteractionElement {

  ushort r[4];

  ushort nu[4];

  std::array<ushort, 4> r;

  std::array<ushort, 4> nu;

  double w;

  short partner_id = -1;

  // TODO: write proper constructor.

  std::vector<ushort> partners_id = std::vector<ushort>();

  // Returns true if r and nu members are equal.

  bool operator==(const InteractionElement& other) const;

};
@@ -51,11 +52,17 @@ public: 


  void insertElement(const std::vector<double>& vec);



  template <class Nu, class Rdmn, class TDmn>

  void checkForInterbandPropagators(

      const func::function<double, func::dmn_variadic<Nu, Nu, Rdmn, TDmn>>& G_r_t);



  void reset();

  // In: random number generator object.

  // Returns: first: random vertex sampled with probability proportional to |vertex.w|.

  //          second: first vertex partner's id if it exists, -1 otherwise.

  std::pair<short, short> getInsertionIndices(double rand) const;

  template <class Rng>

  std::pair<short, short> getInsertionIndices(Rng& rng, double double_update_prob) const;



  // Returns: the sum of the absolute values of the interaction strengths.

  double integratedInteraction() const {

    return total_weigth_;
@@ -69,12 +76,51 @@ public: 
    return elements_[idx];

  }



private:

  // Returns the number of possible partners for each non density-density interaction.

  int possiblePartners() const {

    const int partners = elements_.back().partners_id.size();

    // TODO: generalize if number of possible pairings is not constant or at the back.

    return partners;

  }



  std::vector<InteractionElement> elements_;



private:

  std::vector<double> cumulative_weigths_;

  double total_weigth_ = 0;

  bool interband_propagator_ = false;

};



template <class Rng>

std::pair<short, short> InteractionVertices::getInsertionIndices(Rng& rng,

                                                                 double double_update_prob) const {

  const double random = rng() * total_weigth_;

  // search in reverse order.

  const auto it_to_vertex =

      std::upper_bound(cumulative_weigths_.rbegin(), cumulative_weigths_.rend(), random);

  const int index = cumulative_weigths_.rend() - it_to_vertex - 1;

  assert(index >= 0 && index < size());



  assert(double_update_prob >= 0 && double_update_prob <= 1);

  auto do_double = [&]() -> bool {

    if (double_update_prob == 0)

      return 0;

    else if (double_update_prob == 1)

      return 1;

    else

      return rng() < double_update_prob;

  };



  if (elements_[index].partners_id.size() && do_double()) {  // double insertion

    const auto& partners = elements_[index].partners_id;

    auto partner_id = partners[rng() * partners.size()];

    return std::make_pair(index, partner_id);

  }

  else {

    return std::make_pair(index, -1);

  }

}



template <class Nu, class Rdmn>

void InteractionVertices::initializeFromHamiltonian(

    const func::function<double, func::dmn_variadic<Nu, Nu, Rdmn>>& H_int) {
@@ -116,7 +162,7 @@ void InteractionVertices::initializeFromHamiltonian( 
        already_inserted(nu1, nu2, delta_r) = true;

        for (ushort r1 = 0; r1 < Rdmn::dmn_size(); r1++) {

          const ushort r2 = Rdmn::parameter_type::subtract(delta_r, r1);  // delta_r = r1 - r2

          insertElement(InteractionElement{{r1, r1, r2, r2}, {nu1, nu1, nu2, nu2}, value, short(-1)});

          insertElement(InteractionElement{{r1, r1, r2, r2}, {nu1, nu1, nu2, nu2}, value});

        }

      }

  }
@@ -127,7 +173,7 @@ void InteractionVertices::initializeFromNonDensityHamiltonian( 
    const func::function<double, func::dmn_variadic<Nu, Nu, Nu, Nu, Rdmn>>& H_int) {

  auto spin = [](ushort nu) { return nu >= Nu::dmn_size() / 2; };

  auto check_spins = [&](ushort nu1, ushort nu2, ushort nu3, ushort nu4) -> bool {

    return spin(nu1) == spin(nu2) and spin(nu3) == spin(nu4);

    return spin(nu1) == spin(nu2) && spin(nu3) == spin(nu4);

  };



  for (ushort nu1 = 0; nu1 < Nu::dmn_size(); nu1++)
@@ -144,8 +190,24 @@ void InteractionVertices::initializeFromNonDensityHamiltonian( 
                  "operators have same spin."));

            for (ushort r1 = 0; r1 < Rdmn::dmn_size(); r1++) {

              const ushort r2 = Rdmn::parameter_type::subtract(delta_r, r1);

              insertElement(

                  InteractionElement{{r1, r1, r2, r2}, {nu1, nu2, nu3, nu4}, value, short(-1)});

              insertElement(InteractionElement{{r1, r1, r2, r2}, {nu1, nu2, nu3, nu4}, value});

            }

          }

}



template <class Nu, class RDmn, class TDmn>

void InteractionVertices::checkForInterbandPropagators(

    const func::function<double, func::dmn_variadic<Nu, Nu, RDmn, TDmn>>& G_r_t) {

  interband_propagator_ = false;

  const int t0 = TDmn::dmn_size() / 2;

  const int nb = Nu::dmn_size() / 2;

  const int r0 = RDmn::parameter_type::origin_index();

  for (int r = 0; r < RDmn::dmn_size(); ++r)

    for (int b1 = 0; b1 < nb; ++b1)

      for (int b2 = 0; b2 < nb; ++b2) {

        if (r != r0 && b1 != b2 && std::abs(G_r_t(b1, 0, b2, 0, r, t0)) > 1e-8) {

          interband_propagator_ = true;

          return;

        }

      }

}

include/dca/phys/dca_step/cluster_solver/ctint/structs/ct_int_configuration.hpp→include/dca/phys/dca_step/cluster_solver/ctint/structs/solver_configuration.hpp

+110 −98
Original line number Diff line number Diff line
@@ -9,15 +9,17 @@ 
//

// This class organizes the vertex configuration for ct-int.



#ifndef DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_CONFIGURATION_HPP

#define DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_CONFIGURATION_HPP

#ifndef DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_SOLVER_CONFIGURATION_HPP

#define DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_SOLVER_CONFIGURATION_HPP



#include <array>

#include <numeric>

#include <vector>



#include "dca/io/buffer.hpp"

#include "dca/phys/dca_step/cluster_solver/ctint/structs/ct_int_matrix_configuration.hpp"

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

#include "dca/phys/dca_step/cluster_solver/ctint/structs/utils.hpp"

#include "dca/linalg/device_type.hpp"



namespace dca {
@@ -30,21 +32,25 @@ class SolverConfiguration : public MatrixConfiguration { 
public:

  using BaseClass = MatrixConfiguration;



  inline SolverConfiguration() {}

  inline SolverConfiguration(const SolverConfiguration& other) = default;

  inline SolverConfiguration(double beta, int n_bands, const InteractionVertices& H_int,

  SolverConfiguration() {}

  SolverConfiguration(const SolverConfiguration& other) = default;

  SolverConfiguration(SolverConfiguration&& other) = default;

  SolverConfiguration(double beta, int n_bands, const InteractionVertices& H_int,

                      double double_move_prob = 0);



  inline SolverConfiguration& operator=(const SolverConfiguration& other);

  inline SolverConfiguration& operator=(SolverConfiguration&& other);

  SolverConfiguration& operator=(const SolverConfiguration& other) = default;

  SolverConfiguration& operator=(SolverConfiguration&& other) = default;



  bool operator==(const SolverConfiguration& rhs) const;



  template <class RngType>

  void insertRandom(RngType& rng);



  // TODO: pass output as argument.

  template <class RngType>

  std::pair<short, short> randomRemovalCandidate(RngType& rng) const;

  std::vector<int> randomRemovalCandidate(RngType& rng, double removal_rand);

  template <class RngType>

  std::vector<int> randomRemovalCandidate(RngType& rng);



  // Remove elements with id remove by copying elements from the end.

  // Postcondition: from and sector_from contains the indices that where moved into the place
@@ -63,8 +69,31 @@ public: 
  inline void pop(int n = 1);

  inline void push_back(const Vertex& v);



  inline int occupationNumber(int vertex_index) const;

  inline int occupationNumber(const Vertex& vertex) const;

  inline int nPartners(int vertex_index) const;



  void prepareForSubmatrixUpdate() {

    removable_.resize(vertices_.size());

    std::iota(removable_.begin(), removable_.end(), 0);

  }



  void commitInsertion(int idx) {

    assert(idx < size() && idx >= removable_.size());

    removable_.push_back(idx);

    if (double_insertion_prob_) {

      const auto tag = vertices_[idx].tag;

      auto& list = existing_[vertices_[idx].interaction_id];

      list.push_back(tag);

    }

  }



  void markForRemoval(int idx) {

    removable_.erase(std::find(removable_.begin(), removable_.end(), idx));

    if (double_insertion_prob_) {

      const auto tag = vertices_[idx].tag;

      auto& list = existing_[vertices_[idx].interaction_id];

      list.erase(std::find(list.begin(), list.end(), tag));

    }

  }



  std::size_t size() const {

    return vertices_.size();
@@ -79,25 +108,10 @@ public: 
    return vertices_.back();

  }



  const InteractionElement& coordinates(int v_idx) const {

    assert(v_idx >= 0 and v_idx < (int)size());

    return (*H_int_)[vertices_[v_idx].interaction_id];

  }



  int get_bands() const {

    return n_bands_;

  }



  const auto& getExisting(const short id) const {

    return existing_[id];

  }



  inline ushort getLeftNu(const int matrix_index) const;

  inline ushort getRightNu(const int matrix_index) const;

  inline ushort getLeftR(const int matrix_index) const;

  inline ushort getRightR(const int matrix_index) const;

  inline double getTau(const int matrix_index) const;

  inline bool getAuxSpin(int matrix_index) const;

  auto getTag(int i) const {

    return vertices_[i].tag;

  }
@@ -115,13 +129,30 @@ public: 
  // Return index corresponding to tag.

  int findTag(std::uint64_t tag) const;



  auto possiblePartners() const {

      return H_int_->possiblePartners();

  }



  friend io::Buffer& operator<<(io::Buffer& buff, const SolverConfiguration& config);

  friend io::Buffer& operator>>(io::Buffer& buff, SolverConfiguration& config);



private:

  inline void addSectorSizes(int idx, std::array<int, 2>& sizes) const;

  const InteractionElement& coordinates(int v_idx) const {

    assert(v_idx >= 0 and v_idx < (int)size());

    return (*H_int_)[vertices_[v_idx].interaction_id];

  }

  void addSectorSizes(int idx, std::array<int, 2>& sizes) const;



  template <class Rng>

  bool doDoubleUpdate(Rng& rng) const {

    if (double_insertion_prob_ == 0)

      return false;

    else if (double_insertion_prob_ == 1)

      return true;

    else

      return rng() < double_insertion_prob_;

  }



protected:

  // List of points entering into the first or second member of g0

  const double double_insertion_prob_ = 0;
@@ -129,6 +160,8 @@ protected: 


  const InteractionVertices* H_int_ = nullptr;

  std::vector<std::vector<std::uint64_t>> existing_;

  std::vector<std::vector<std::size_t>*> partners_lists_;

  std::vector<int> removable_;

  ushort last_insertion_size_ = 1;

  const double max_tau_ = 0;

  const int n_bands_ = 0;
@@ -136,22 +169,7 @@ protected: 
  std::uint64_t current_tag_ = 0;

};



SolverConfiguration& SolverConfiguration::operator=(const SolverConfiguration& other) {

  BaseClass::operator=(other);

  H_int_ = other.H_int_;

  vertices_ = other.vertices_;

  existing_ = other.existing_;

  return *this;

}



SolverConfiguration& SolverConfiguration::operator=(SolverConfiguration&& other) {

  H_int_ = other.H_int_;

  BaseClass::operator=(other);

  vertices_ = std::move(other.vertices_);

  return *this;

}



SolverConfiguration::SolverConfiguration(const double beta, const int n_bands,

inline SolverConfiguration::SolverConfiguration(const double beta, const int n_bands,

                                                const InteractionVertices& H_int,

                                                const double double_insertion)

    : MatrixConfiguration(&H_int, n_bands),
@@ -164,19 +182,19 @@ SolverConfiguration::SolverConfiguration(const double beta, const int n_bands, 


template <class RngType>

void SolverConfiguration::insertRandom(RngType& rng) {

  const std::pair<short, short> indices = H_int_->getInsertionIndices(rng());

  const std::pair<short, short> indices = H_int_->getInsertionIndices(rng, double_insertion_prob_);

  const double tau = rng() * max_tau_;

  const bool aux_spin = rng() > 0.5;



  push_back(Vertex{aux_spin, ushort(indices.first), ++current_tag_, tau});

  push_back(Vertex{aux_spin, ushort(indices.first), current_tag_++, tau});



  if (double_insertion_prob_) {

    // TODO: generalize to multiband n_bands > 2

    if (indices.second != -1 and rng() < double_insertion_prob_) {

    if (indices.second != -1 && double_insertion_prob_) {

      const double tau2 = rng() * max_tau_;

      const bool aux_spin2 = rng() > 0.5;

      // TODO: decide if tag is same or not.

      push_back(Vertex{aux_spin2, ushort(indices.second), ++current_tag_, tau2});

      push_back(Vertex{aux_spin2, ushort(indices.second), current_tag_++, tau2});

      last_insertion_size_ = 2;

    }

    else
@@ -186,51 +204,47 @@ void SolverConfiguration::insertRandom(RngType& rng) { 
}



template <class RngType>

std::pair<short, short> SolverConfiguration::randomRemovalCandidate(RngType& rng) const {

  std::pair<short, short> candidates(rng() * size(), -1);

  if (double_insertion_prob_) {

    const int partner_id = (*H_int_)[vertices_[candidates.first].interaction_id].partner_id;

    if (partner_id != -1 and rng() < double_insertion_prob_) {

      const std::size_t n_partners = existing_[partner_id].size();

      if (n_partners) {

        auto tag = existing_[partner_id][int(rng() * n_partners)];

        candidates.second = findTag(tag);

      }

std::vector<int> SolverConfiguration::randomRemovalCandidate(RngType& rng)  {

  return randomRemovalCandidate(rng, rng());

}



template <class RngType>

std::vector<int> SolverConfiguration::randomRemovalCandidate(RngType& rng, double removal_rand) {

  std::vector<int> candidates;

  if (removable_.size() == 0)

    return candidates;



  candidates.push_back(removable_[removal_rand * removable_.size()]);



  if ((*H_int_)[vertices_[candidates[0]].interaction_id].partners_id.size() &&

      doDoubleUpdate(rng)) {  // Double removal.

    partners_lists_.clear();

    for (const auto& partner_id : (*H_int_)[vertices_[candidates[0]].interaction_id].partners_id)

      partners_lists_.push_back(&existing_[partner_id]);



    const auto tag = details::getRandomElement(partners_lists_, rng());

    candidates.push_back(findTag(tag));

    assert(candidates[1] < int(size()) && candidates[1] >= 0);

  }

  assert(candidates.first < size() and candidates.second < int(size()));



  assert(candidates[0] < int(size()));

  return candidates;

}



void SolverConfiguration::push_back(const Vertex& v) {

  vertices_.push_back(v);

  BaseClass::addVertex(v);



  if (double_insertion_prob_)

    existing_[v.interaction_id].push_back(v.tag);

}



void SolverConfiguration::pop(const int n) {

  assert(n <= (int)size());

  assert(n == 1 or (*H_int_)[vertices_[size() - 2].interaction_id].partner_id ==

                       vertices_[size() - 1].interaction_id);

  if (double_insertion_prob_) {

    for (int i = 1; i <= n; ++i) {

      const int type = vertices_[size() - i].interaction_id;

      auto& list = existing_[type];

      const auto iter = std::find(list.begin(), list.end(), vertices_[size() - i].tag);

      list.erase(iter);

    }

  }



  const int first_idx = size() - n;

  std::array<int, 2> removal_size{0, 0};

  for (std::size_t i = first_idx; i < size(); ++i)

    addSectorSizes(i, removal_size);

  vertices_.erase(vertices_.begin() + first_idx, vertices_.end());

  BaseClass::pop(removal_size[0], removal_size[1]);



  // ctint base walker leaves configuration temporary inconsistent by (poor) design.

  // assert(checkConsistency());

}



std::array<int, 2> SolverConfiguration::sizeIncrease() const {
@@ -242,27 +256,29 @@ std::array<int, 2> SolverConfiguration::sizeIncrease() const { 
  return result;

}



void SolverConfiguration::addSectorSizes(int idx, std::array<int, 2>& sizes) const {

inline void SolverConfiguration::addSectorSizes(int idx, std::array<int, 2>& sizes) const {

  auto spin = [=](ushort nu) { return nu >= n_bands_; };

  const auto& nu = coordinates(idx).nu;

  ++sizes[spin(nu[0])];

  ++sizes[spin(nu[2])];

}



short SolverConfiguration::getSign(const int vertex_index) const {

  return (*H_int_)[vertices_[vertex_index].interaction_id].w >= 0 ? 1 : -1;

inline short SolverConfiguration::getSign(const int vertex_index) const {

  return getStrength(vertex_index) > 0 ? 1 : -1;

}



double SolverConfiguration::getStrength(int vertex_index) const {

inline double SolverConfiguration::getStrength(int vertex_index) const {

  assert(vertex_index >= 0 && vertex_index < (int)size());

  return (*H_int_)[vertices_[vertex_index].interaction_id].w;

}



int SolverConfiguration::occupationNumber(int vertex_index) const {

  return existing_[vertices_[vertex_index].interaction_id].size();

}

int SolverConfiguration::occupationNumber(const Vertex& vertex) const {

  return existing_[vertex.interaction_id].size();

int SolverConfiguration::nPartners(int vertex_index) const {

  assert(vertex_index < vertices_.size());

  const auto& partners_id = (*H_int_)[vertices_[vertex_index].interaction_id].partners_id;

  int n_partners = 0;

  for (const auto partner_id : partners_id)

    n_partners += existing_[partner_id].size();

  return n_partners;

}



template <class Alloc>
@@ -303,16 +319,6 @@ inline void SolverConfiguration::moveAndShrink(std::array<std::vector<int, Alloc 
      from.push_back(source_idx);

  }



  // Remove from list of partners.

  if (double_insertion_prob_) {

    for (int i = 0; i < to.size(); ++i) {

      const int type = vertices_[to[i]].interaction_id;

      auto& list = existing_[type];

      const auto iter = std::find(list.begin(), list.end(), vertices_[to[i]].tag);

      list.erase(iter);

    }

  }



  // Move and shrink configuration.

  for (int i = 0; i < from.size(); ++i)

    vertices_[to[i]] = vertices_[from[i]];
@@ -359,6 +365,12 @@ inline bool SolverConfiguration::checkConsistency() const { 
      const auto& list = existing_[v.interaction_id];

      if (std::find(list.begin(), list.end(), v.tag) == list.end())

        return false;

      // Check total size.

      int size_sum = 0;

      for (const auto& list : existing_)

        size_sum += list.size();

      if (size_sum != vertices_.size())

        return false;

    }

  }

  return true;
@@ -376,4 +388,4 @@ inline int SolverConfiguration::findTag(std::uint64_t tag) const { 
}  // namespace phys

}  // namespace dca



#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_CT_INT_CONFIGURATION_HPP
 
#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_SOLVER_CONFIGURATION_HPP

include/dca/phys/dca_step/cluster_solver/ctint/structs/utils.hpp

0 → 100644
+51 −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 class organizes the vertex configuration for ct-int.



#ifndef DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_UTIL_HPP

#define DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_UTIL_HPP



#include <stdexcept>

#include <vector>



namespace dca {

namespace phys {

namespace solver {

namespace ctint {

namespace details {

// dca::phys::solver::ctint::details::



template<class T>

T getRandomElement(const std::vector<std::vector<T>*>& v_ptrs, const double rand){

    assert(rand >= 0 && rand <= 1);



    unsigned size = 0;

    for(int i = v_ptrs.size() - 1; i >= 0; --i)

        size += v_ptrs[i]->size();



    // TODO: use binary search or other efficient scheme.

    unsigned idx = size * rand;

    for(auto v_ptr : v_ptrs){

        if (idx < v_ptr->size())

            return (*v_ptr)[idx];

        idx -= v_ptr->size();

    }



    throw(std::logic_error("Something went wrong."));

}





}  // namespace details

}  // namespace ctint

}  // namespace solver

}  // namespace phys

}  // namespace dca



#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTINT_STRUCTS_UTIL_HPP