Commit 25997e8d authored by gbalduzz's avatar gbalduzz
Browse files

O(1) e_spin_configuration shrinkage.

parent 90f7ab89

include/dca/phys/dca_step/cluster_solver/ctaux/ctaux_walker.hpp

+2 −2
Original line number Diff line number Diff line
@@ -1030,8 +1030,8 @@ void CtauxWalker<device_t, Parameters, Data, Real>::add_delayed_spins_to_the_con 
        configuration.add_delayed_HS_spin(configuration_index, delayed_spins[i].new_HS_spin_value);

      }

      else {

        configuration[configuration_index].is_creatable() = false;

        configuration[configuration_index].is_annihilatable() = false;

        configuration[configuration_index].set_creatable(false);

        configuration[configuration_index].set_annihilatable(false);

      }

    }

  }

include/dca/phys/dca_step/cluster_solver/ctaux/structs/ct_aux_hs_configuration.hpp

+71 −66
Original line number Diff line number Diff line
@@ -49,6 +49,9 @@ public: 
  vertex_pair_type& operator[](int index);



  std::vector<vertex_singleton_type>& get(e_spin_states_type e_spin_type);

  auto& get() {

    return configuration;

  }



  void reset();
@@ -85,6 +88,10 @@ public: 
  // a later annihilation proposal.

  int get_random_noninteracting_vertex(bool mark_annihilatable);



  // Remove element at index 'idx' in O(1). Vertices order is not preserved. The electron spin

  // configuration needs to be updated independently.

  void erase(unsigned idx);



  // debug tools

  void print();

  void print(e_spin_states_type e_spin_type);
@@ -99,7 +106,7 @@ public: 
  //           temporarily been marked as annihilatable by get_random_noninteracting_vertex.

  void unmarkAsAnnihilatable(const int vertex_index) {

    assert(configuration[vertex_index].is_annihilatable() == true);

    configuration[vertex_index].is_annihilatable() = false;

    configuration[vertex_index].set_annihilatable(false);

    --current_Nb_of_annihilatable_spins;

  }
@@ -219,8 +226,7 @@ void CT_AUX_HS_configuration<parameters_type>::initialize() { 
  reset();



  for (int i = 0; i < parameters.get_initial_configuration_size(); ++i) {

    vertex_pair_type vertex(parameters, rng, configuration.size(), configuration_e_DN.size(),

                            configuration_e_UP.size(), next_vertex_id_++);

    vertex_pair_type vertex(parameters, rng, configuration.size(), next_vertex_id_++);

    vertex.set_random_interacting();



    ++current_Nb_of_annihilatable_spins;
@@ -236,9 +242,9 @@ void CT_AUX_HS_configuration<parameters_type>::shuffle_noninteracting_vertices() 
  assert(changed_spin_indices.size() == 0);



  for (size_t i = 0; i < configuration.size(); i++) {

    configuration[i].is_shuffled() = false;

    configuration[i].is_successfully_flipped() = false;

    configuration[i].is_Bennett() = false;

    configuration[i].set_shuffled(false);

    configuration[i].set_successfully_flipped(false);

    configuration[i].set_Bennett(false);



    assert(configuration[i].is_annihilatable() || configuration[i].is_creatable());

    assert(configuration[i].is_annihilatable() != configuration[i].is_creatable());
@@ -246,8 +252,7 @@ void CT_AUX_HS_configuration<parameters_type>::shuffle_noninteracting_vertices() 


  // Add non-interacting spins.

  while (current_Nb_of_creatable_spins < max_num_noninteracting_spins_) {

    vertex_pair_type vertex(parameters, rng, configuration.size(), configuration_e_DN.size(),

                            configuration_e_UP.size(), next_vertex_id_++);

    vertex_pair_type vertex(parameters, rng, configuration.size(), next_vertex_id_++);

    vertex.set_random_noninteracting();



    ++current_Nb_of_creatable_spins;
@@ -283,64 +288,34 @@ void CT_AUX_HS_configuration<parameters_type>::update_configuration_e_spin( 
template <class parameters_type>

int CT_AUX_HS_configuration<parameters_type>::get_first_non_interacting_spin_index(

    e_spin_states_type e_spin) {

  std::vector<vertex_singleton_type>& configuration_e_spin = get(e_spin);

  int configuration_size = configuration_e_spin.size();

  std::vector<vertex_singleton_type>& spin_configuration = get(e_spin);



  // FIXME: What and when to return if configuration_size = 0?

  // Note: should be handled by the caller



  int first_non_interacting_index, configuration_index;

  for (first_non_interacting_index = 0; first_non_interacting_index < configuration_size;

       first_non_interacting_index++) {

    configuration_index = configuration_e_spin[first_non_interacting_index].get_configuration_index();



    if (!configuration[configuration_index].is_annihilatable())

      break;

  }



  assert(

      (first_non_interacting_index == 0 &&

       !configuration[configuration_e_spin[first_non_interacting_index].get_configuration_index()]

            .is_annihilatable()) ||

      (first_non_interacting_index > 0 && first_non_interacting_index < configuration_size &&

       !configuration[configuration_e_spin[first_non_interacting_index].get_configuration_index()]

            .is_annihilatable() &&

       configuration[configuration_e_spin[first_non_interacting_index - 1].get_configuration_index()]

           .is_annihilatable()) ||

      (first_non_interacting_index == configuration_size &&

       configuration[configuration_e_spin[first_non_interacting_index - 1].get_configuration_index()]

           .is_annihilatable()));

  const int first_non_interacting_index = std::distance(

      spin_configuration.begin(),

      std::find_if(spin_configuration.begin(), spin_configuration.end(), [&](const auto& elem) {

        const int configuration_index = elem.get_configuration_index();

        return !configuration[configuration_index].is_annihilatable();

      }));



  return first_non_interacting_index;

}



template <class parameters_type>

int CT_AUX_HS_configuration<parameters_type>::get_first_shuffled_spin_index(e_spin_states_type e_spin) {

  std::vector<vertex_singleton_type>& configuration_e_spin = get(e_spin);

  int configuration_size = configuration_e_spin.size();



  std::vector<vertex_singleton_type>& spin_configuration = get(e_spin);

  // FIXME: What and when to return if configuration_size = 0?



  assert(assert_block_form(e_spin));



  int first_shuffled_index, configuration_index;

  for (first_shuffled_index = 0; first_shuffled_index < configuration_size; first_shuffled_index++) {

    configuration_index = configuration_e_spin[first_shuffled_index].get_configuration_index();



    if (configuration[configuration_index].is_shuffled())

      break;

  }



  assert((first_shuffled_index == 0 &&

          configuration[configuration_e_spin[first_shuffled_index].get_configuration_index()]

              .is_shuffled()) ||

         (first_shuffled_index > 0 && first_shuffled_index < configuration_size &&

          configuration[configuration_e_spin[first_shuffled_index].get_configuration_index()]

              .is_shuffled() &&

          !configuration[configuration_e_spin[first_shuffled_index - 1].get_configuration_index()]

               .is_shuffled()) ||

         (first_shuffled_index == configuration_size &&

          !configuration[configuration_e_spin[first_shuffled_index - 1].get_configuration_index()]

               .is_shuffled()));

  const int first_shuffled_index = std::distance(

      spin_configuration.begin(),

      std::find_if(spin_configuration.begin(), spin_configuration.end(), [&](const auto& elem) {

        const int index = elem.get_configuration_index();

        return configuration[index].is_shuffled();

      }));



  return first_shuffled_index;

}
@@ -382,8 +357,8 @@ void CT_AUX_HS_configuration<parameters_type>::add_delayed_HS_spin(int configura 


    current_Nb_of_annihilatable_spins -= 1;



    configuration[configuration_index].is_annihilatable() = false;

    configuration[configuration_index].is_Bennett() = true;

    configuration[configuration_index].set_annihilatable(false);

    configuration[configuration_index].set_Bennett(true);



    return;

  }
@@ -396,17 +371,17 @@ void CT_AUX_HS_configuration<parameters_type>::add_delayed_HS_spin(int configura 


    current_Nb_of_annihilatable_spins -= 1;



    configuration[configuration_index].is_annihilatable() = false;

    configuration[configuration_index].set_annihilatable(false);

  }

  else {

    // cout << "\t--> create spin : " << configuration_index << std::endl;



    current_Nb_of_annihilatable_spins += 1;



    configuration[configuration_index].is_annihilatable() = true;

    configuration[configuration_index].set_annihilatable(true);

  }



  configuration[configuration_index].is_successfully_flipped() = true;

  configuration[configuration_index].set_successfully_flipped(true);

}



template <class parameters_type>
@@ -502,20 +477,51 @@ int CT_AUX_HS_configuration<parameters_type>::get_random_noninteracting_vertex(b 
  assert(!configuration[vertex_index].is_Bennett());



  // Make sure that this spin will not be proposed for insertion again.

  configuration[vertex_index].is_creatable() = false;

  configuration[vertex_index].set_creatable(false);

  --current_Nb_of_creatable_spins;



  if (mark_annihilatable) {

    // However, this "virtual" interacting spin is eligible for removal.

    // INTERNAL: CtauxWalker::generateDelayedSpinsAbortAtBennett unmarks the "virtual" interacting

    //           spins as annihilatable when all delayed spins have been generated.

    configuration[vertex_index].is_annihilatable() = true;

    configuration[vertex_index].set_annihilatable(true);

    ++current_Nb_of_annihilatable_spins;

  }



  return vertex_index;

}



template <class parameters_type>

void CT_AUX_HS_configuration<parameters_type>::erase(unsigned idx) {

  assert(configuration[idx].get_HS_spin() == HS_ZERO);

  assert(configuration[idx].get_configuration_e_spin_indices().first == -1);

  assert(configuration[idx].get_configuration_e_spin_indices().second == -1);



  if (configuration[idx].is_annihilatable())

    --current_Nb_of_annihilatable_spins;

  if (configuration[idx].is_creatable())

    --current_Nb_of_creatable_spins;



  configuration[idx] = configuration.back();



  const auto& e_spin_indices = configuration[idx].get_configuration_e_spin_indices();

  const auto& e_spins = configuration[idx].get_e_spins();



  const auto update_index = [&](const auto& e_spin, const int spin_index) {

    if (spin_index >= 0) {

      assert(spin_index < get(e_spin).size());

      get(e_spin)[spin_index].get_configuration_index() = idx;

    }

  };



  update_index(e_spins.first, e_spin_indices.first);

  update_index(e_spins.second, e_spin_indices.second);



  configuration[idx].get_configuration_index() = idx;



  configuration.pop_back();

}



template <class parameters_type>

int CT_AUX_HS_configuration<parameters_type>::get_number_of_creatable_HS_spins() {

  assert(assert_counters());
@@ -576,9 +582,8 @@ bool CT_AUX_HS_configuration<parameters_type>::assert_counters() { 


template <class parameters_type>

bool CT_AUX_HS_configuration<parameters_type>::assert_consistency() {

  // std::cout << __FUNCTION__ <<std::endl;



  assert(2 * configuration.size() == (configuration_e_UP.size() + configuration_e_DN.size()));



  assert_counters();



  // assert configuration-index
@@ -861,9 +866,9 @@ bool CT_AUX_HS_configuration<parameters_type>::operator==( 
         current_Nb_of_creatable_spins == rhs.current_Nb_of_creatable_spins;

}



}  // ctaux

}  // solver

}  // phys

}  // dca

}  // namespace ctaux

}  // namespace solver

}  // namespace phys

}  // namespace dca



#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTAUX_STRUCTS_CT_AUX_HS_CONFIGURATION_HPP

include/dca/phys/dca_step/cluster_solver/ctaux/structs/read_write_config.hpp

+1 −2
Original line number Diff line number Diff line
@@ -66,8 +66,7 @@ io::Buffer& operator>>(io::Buffer& buff, CT_AUX_HS_configuration<Parameters>& co 


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

    vertex_pair<Parameters> vertex(config.parameters, config.rng, config.configuration.size(),

                                   config.configuration_e_DN.size(),

                                   config.configuration_e_UP.size(), config.next_vertex_id_++);

                                   config.next_vertex_id_++);



    buff >> vertex;

include/dca/phys/dca_step/cluster_solver/ctaux/structs/vertex_pair.hpp

+48 −51
Original line number Diff line number Diff line
@@ -46,8 +46,8 @@ public: 
  typedef vertex_pair<parameters_type> this_type;



public:

  vertex_pair(parameters_type& parameters_ref, rng_type& rng_ref, int configuration_index_in,

              int configuration_e_DN_index_in, int configuration_e_UP_index_in, uint64_t id);

  vertex_pair(parameters_type& parameters_ref, rng_type& rng_ref, int configuration_index,

              uint64_t id);



  this_type& operator=(const this_type& other_vertex_pair);
@@ -74,14 +74,46 @@ public: 


  // algorithm-information



  int& get_configuration_index();

  int& get_configuration_index() {

    return configuration_index;

  }

  const int& get_configuration_index() const {

    return configuration_index;

  }

  std::pair<int, int>& get_configuration_e_spin_indices();



  bool& is_creatable();

  bool& is_annihilatable();

  bool& is_shuffled();

  bool& is_Bennett();

  bool& is_successfully_flipped();

  void set_creatable(bool val) {

    creatable = val;

  }

  void set_annihilatable(bool val) {

    annihilatable = val;

  }

  void set_shuffled(bool val) {

    shuffled = val;

  }

  void set_successfully_flipped(bool val) {

    successfully_flipped = val;

  }

  void set_Bennett(bool val) {

    Bennett = val;

  }



  bool is_creatable() const {

    return creatable;

  }

  bool is_annihilatable() const {

    return annihilatable;

  }

  bool is_shuffled() const {

    return shuffled;

  }



  bool is_Bennett() const {

    return Bennett;

  }

  bool is_successfully_flipped() const {

    return successfully_flipped;

  }



  bool operator==(const vertex_pair<parameters_type>& rhs) const;
@@ -104,7 +136,6 @@ private: 
  std::pair<int, int> bands;

  std::pair<e_spin_states_type, e_spin_states_type> e_spins;

  std::pair<int, int> spin_orbitals;

  std::pair<HS_field_sign_type, HS_field_sign_type> HS_fields;

  std::pair<int, int> r_sites;



  HS_spin_states_type HS_spin;
@@ -112,7 +143,7 @@ private: 
  double tau;



  // algorithm-information

  int configuration_index;

  int configuration_index; // TODO: remove redundant info.

  std::pair<int, int> configuration_e_spin_indices;



  bool creatable;
@@ -124,12 +155,9 @@ private: 


template <class parameters_type>

vertex_pair<parameters_type>::vertex_pair(parameters_type& parameters_ref, rng_type& rng_ref,

                                          int configuration_index_in,

                                          int /*configuration_e_DN_index_in*/,

                                          int /*configuration_e_UP_index_in*/, uint64_t id)

                                          int configuration_index_in, uint64_t id)

    : parameters(parameters_ref),

      rng(rng_ref),

      //     concurrency(parameters.get_concurrency()),



      id_(id),
@@ -139,7 +167,6 @@ vertex_pair<parameters_type>::vertex_pair(parameters_type& parameters_ref, rng_t 
      bands(0, 0),

      e_spins(e_DN, e_UP),

      spin_orbitals(0, 1),

      HS_fields(HS_FIELD_DN, HS_FIELD_UP),

      r_sites(0, 0),



      HS_spin(HS_ZERO),
@@ -285,41 +312,11 @@ double& vertex_pair<parameters_type>::get_tau() { 
  return tau;

}



template <class parameters_type>

int& vertex_pair<parameters_type>::get_configuration_index() {

  return configuration_index;

}



template <class parameters_type>

std::pair<int, int>& vertex_pair<parameters_type>::get_configuration_e_spin_indices() {

  return configuration_e_spin_indices;

}



template <class parameters_type>

bool& vertex_pair<parameters_type>::is_creatable() {

  return creatable;

}



template <class parameters_type>

bool& vertex_pair<parameters_type>::is_annihilatable() {

  return annihilatable;

}



template <class parameters_type>

bool& vertex_pair<parameters_type>::is_successfully_flipped() {

  return successfully_flipped;

}



template <class parameters_type>

bool& vertex_pair<parameters_type>::is_Bennett() {

  return Bennett;

}



template <class parameters_type>

bool& vertex_pair<parameters_type>::is_shuffled() {

  return shuffled;

}



template <class parameters_type>

bool vertex_pair<parameters_type>::operator==(

    const dca::phys::solver::ctaux::vertex_pair<parameters_type>& rhs) const {
@@ -327,17 +324,17 @@ bool vertex_pair<parameters_type>::operator==( 


  return id_ == rhs.id_ && interacting_bands == rhs.interacting_bands && BANDS == rhs.BANDS &&

         bands == rhs.bands && e_spins == rhs.e_spins && spin_orbitals == rhs.spin_orbitals &&

         HS_fields == rhs.HS_fields && r_sites == rhs.r_sites && HS_spin == rhs.HS_spin &&

         delta_r == rhs.delta_r && tau == rhs.tau && configuration_index == rhs.configuration_index &&

         r_sites == rhs.r_sites && HS_spin == rhs.HS_spin && delta_r == rhs.delta_r &&

         tau == rhs.tau && configuration_index == rhs.configuration_index &&

         configuration_e_spin_indices == rhs.configuration_e_spin_indices &&

         creatable == rhs.creatable && annihilatable == rhs.annihilatable &&

         successfully_flipped == rhs.successfully_flipped && Bennett == rhs.Bennett &&

         shuffled == rhs.shuffled;

}



}  // ctaux

}  // solver

}  // phys

}  // dca

}  // namespace ctaux

}  // namespace solver

}  // namespace phys

}  // namespace dca



#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTAUX_STRUCTS_VERTEX_PAIR_HPP

include/dca/phys/dca_step/cluster_solver/ctaux/structs/vertex_singleton.hpp

+6 −6
Original line number Diff line number Diff line
@@ -106,8 +106,8 @@ private: 
  int delta_r;

  double tau;



  HS_spin_states_type HS_spin;

  HS_field_sign_type HS_field;

  HS_spin_states_type HS_spin;  // Auxiliary spin.

  HS_field_sign_type HS_field;  // Marker to denote if it is the first or second singleton of a pair.



  int configuration_index;

};
@@ -131,9 +131,9 @@ vertex_singleton& vertex_singleton::get_partner(configuration_type& configuratio 
  return configuration.get(e_spin)[configuration_e_spin];

}



}  // ctaux

}  // solver

}  // phys

}  // dca

}  // namespace ctaux

}  // namespace solver

}  // namespace phys

}  // namespace dca



#endif  // DCA_PHYS_DCA_STEP_CLUSTER_SOLVER_CTAUX_STRUCTS_VERTEX_SINGLETON_HPP