Merge pull request #162 from gbalduzz/optimize_ctint_configuration (ce7cf572) · Commits · NDIP / Tool Sources / Direct-Geometry Spectroscopy / DCA / DCA Main

cmake/dca_compiler_options.cmake

+1 −1

Original line number	Diff line number	Diff line
		@@ -21,7 +21,7 @@
		set(DCA_WARNINGS -Wall -Wextra -Wpedantic -Wno-sign-compare -Wno-dangling-else)

		# Languange standard
		set(DCA_STD_FLAG -std=c++14)
		set(DCA_STD_FLAG -std=c++17)

		# Set C and CXX flags.
		add_compile_options(${DCA_WARNINGS})

include/dca/math/random/std_random_wrapper.hpp

+1 −1

Original line number	Diff line number	Diff line
		@@ -63,7 +63,7 @@ public:
		counter_ = 0;
		}

		// Returns a uniformly distributied pseudo-random number in the interval [0, 1).
		// Returns a uniformly distributed pseudo-random number in the interval [0, 1).
		inline double operator()() {
		return distro_(engine_);
		}

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

+30 −87

Original line number	Diff line number	Diff line
		@@ -28,22 +28,44 @@ namespace ctint {

		// Expansion term.
		struct Vertex {
		Vertex() = default;
		Vertex(bool _aux_spin, ushort _interaction_id, std::uint64_t _tag, double _tau)
		: aux_spin(_aux_spin), interaction_id(_interaction_id), tag(_tag), tau(_tau) {
		spins.fill(-1);
		matrix_config_indices.fill(-1);
		}

		bool aux_spin;
		ushort interaction_id;
		std::uint64_t tag;
		double tau;

		// Reference to the matrix config
		std::array<std::uint8_t, 2> spins;
		std::array<unsigned, 2> matrix_config_indices;

		// Marks if this vertex is part of the accepted configuration.
		bool annihilatable = false;

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

		struct ConfigRef {
		ConfigRef() = default;
		ConfigRef(unsigned _config_id, std::uint8_t _leg_id) : config_id(_config_id), leg_id(_leg_id) {}

		unsigned config_id; // Index of the interaction in the SolverConfiguration.
		std::uint8_t leg_id; // In {0, 1}. Stores if this is the first or second leg of an interaction vertex.
		};

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

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

		@@ -67,15 +89,15 @@ public:
		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);
		MatrixConfiguration(const InteractionVertices* H_int, int bands);

		inline std::array<int, 2> addVertex(const Vertex& v);
		// In/Out: v. The vertex to be added, the matrix configuration spins and indices are updated.
		// In. config_id. Position of the vertex in the solver configuration.
		// Out: config_refs. References from matrix configuration to solver configuration to be updated.
		void addVertex(Vertex& v, unsigned config_id, std::array<std::vector<ConfigRef>, 2>& config_refs);

		inline void pop(ushort idx_up, ushort idx_down);
		// inline void pop(ushort idx_up, ushort idx_down);

		const auto& getEntries(const int s) const {
		return sectors_[s].entries_;
		@@ -84,90 +106,11 @@ protected:
		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

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

+6 −6

Original line number	Diff line number	Diff line
		@@ -26,7 +26,8 @@ namespace solver {
		namespace ctint {
		// dca::phys::solver::ctint::

		// TODO: replace by vector


		class Sector {
		public:
		Sector() : entries_(){};
		@@ -39,6 +40,10 @@ public:
		return entries_.size();
		}

		void pop_back(){
		entries_.pop_back();
		}

		const details::SectorEntry& operator[](const std::size_t idx) const {
		assert(idx < size());
		return entries_[idx];
		@@ -68,17 +73,12 @@ public:
		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

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

+91 −224

Original line number	Diff line number	Diff line
		@@ -15,6 +15,7 @@
		#include <array>
		#include <numeric>
		#include <vector>
		#include <vector>

		#include "dca/io/buffer.hpp"
		#include "dca/phys/dca_step/cluster_solver/ctint/structs/ct_int_matrix_configuration.hpp"
		@@ -31,6 +32,8 @@ namespace ctint {
		class SolverConfiguration : public MatrixConfiguration {
		public:
		using BaseClass = MatrixConfiguration;
		template <class T>
		using HostVector = dca::linalg::util::HostVector<T>;

		SolverConfiguration() {}
		SolverConfiguration(const SolverConfiguration& other) = default;
		@@ -46,54 +49,29 @@ public:
		template <class RngType>
		void insertRandom(RngType& rng);

		// TODO: pass output as argument.
		// Returns the indices of the removal candidates. -1 stands for a missing candidate.
		template <class RngType>
		std::vector<int> randomRemovalCandidate(RngType& rng, double removal_rand);
		std::array<int, 2> 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
		// of the old elements.
		// Postcondition: all vectors are ordered, resized, and ready to be used for moving the rows of
		// the M matrix
		// In/Out: sector_remove, remove.
		// Out: sector_from, from.
		template <class Alloc>
		void moveAndShrink(std::array<std::vector<int, Alloc>, 2>& sector_from,
		std::array<std::vector<int, Alloc>, 2>& sector_remove, std::vector<int>& from,
		std::vector<int>& remove);
		std::array<int, 2> randomRemovalCandidate(RngType& rng);

		inline void swapVertices(short i, short j);
		// Out: indices. Appends the result of the search to indices.
		template <class Alloc>
		void findIndices(std::vector<int, Alloc>& indices, unsigned config_index, int s) const;

		inline void pop(int n = 1);
		inline void push_back(const Vertex& v);
		// Remove elements with index in 'remove' by copying elements from the end.
		// In: remove. List of vertices to be removed. Sorted on exit.
		// Out: sector remove, sector_from. Lists of matrix indices to be moved to maintain consistency.
		void moveAndShrink(std::array<HostVector<int>, 2>& sector_from,
		std::array<HostVector<int>, 2>& sector_remove, std::vector<int>& remove);

		inline int nPartners(int vertex_index) const;
		// In/Out: v. Matrix indices are updated.
		void push_back(Vertex& v);

		void prepareForSubmatrixUpdate() {
		removable_.resize(vertices_.size());
		std::iota(removable_.begin(), removable_.end(), 0);
		}
		int nPartners(int vertex_index) const;

		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));
		}
		}
		void commitInsertion(int idx);
		void markForRemoval(int idx);

		std::size_t size() const {
		return vertices_.size();
		@@ -108,10 +86,6 @@ public:
		return vertices_.back();
		}

		int get_bands() const {
		return n_bands_;
		}

		auto getTag(int i) const {
		return vertices_[i].tag;
		}
		@@ -123,7 +97,7 @@ public:
		return last_insertion_size_;
		}

		inline std::array<int, 2> sizeIncrease() const;
		std::array<int, 2> sizeIncrease() const;

		bool checkConsistency() const;
		// Return index corresponding to tag.
		@@ -137,64 +111,48 @@ public:
		friend io::Buffer& operator>>(io::Buffer& buff, SolverConfiguration& config);

		private:
		const InteractionElement& coordinates(int v_idx) const {
		assert(v_idx >= 0 and v_idx < (int)size());
		return (*H_int_)[vertices_[v_idx].interaction_id];
		}
		const InteractionElement& coordinates(int v_idx) const;
		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_;
		}
		bool doDoubleUpdate(Rng& rng) const;

		// List of points entering into the first or second member of g0
		const double double_insertion_prob_ = 0;

		std::vector<Vertex> vertices_;
		// Connection from Vertices to MatrixConfiguration elements.
		std::array<std::vector<ConfigRef>, 2> matrix_config_indices_;

		const InteractionVertices* H_int_ = nullptr;
		using BaseClass::H_int_;
		std::vector<std::vector<std::uint64_t>> existing_;

		// TODO: use a structure with fast (log N?) removal/insertion and random access.
		// Or sample randomly from std::unordered_set using its hash function, if it's good enough.
		std::vector<const 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;

		unsigned n_annihilatable_ = 0;
		std::uint64_t current_tag_ = 0;
		};

		inline SolverConfiguration::SolverConfiguration(const double beta, const int n_bands,
		const InteractionVertices& H_int,
		const double double_insertion)
		: MatrixConfiguration(&H_int, n_bands),
		double_insertion_prob_(double_insertion),

		H_int_(&H_int),
		existing_(double_insertion ? H_int_->size() : 0),
		max_tau_(beta),
		n_bands_(n_bands) {}

		template <class RngType>
		void SolverConfiguration::insertRandom(RngType& rng) {
		template <class Rng>
		void SolverConfiguration::insertRandom(Rng& 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});
		Vertex v(aux_spin, ushort(indices.first), current_tag_++, tau);
		push_back(v);

		if (double_insertion_prob_) {
		// TODO: generalize to multiband n_bands > 2
		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});
		Vertex v2(aux_spin2, ushort(indices.second), current_tag_++, tau2);
		push_back(v2);
		last_insertion_size_ = 2;
		}
		else
		@@ -204,63 +162,83 @@ void SolverConfiguration::insertRandom(RngType& rng) {
		}

		template <class RngType>
		std::vector<int> SolverConfiguration::randomRemovalCandidate(RngType& rng) {
		std::array<int, 2> SolverConfiguration::randomRemovalCandidate(RngType& rng) {
		return randomRemovalCandidate(rng, rng());
		}

		// TODO: possibly use only the above signature.
		template <class RngType>
		std::vector<int> SolverConfiguration::randomRemovalCandidate(RngType& rng, double removal_rand) {
		std::vector<int> candidates;
		if (removable_.size() == 0)
		std::array<int, 2> SolverConfiguration::randomRemovalCandidate(RngType& rng, double removal_rand) {
		// TODO: generalize to n > 2.
		std::array<int, 2> candidates{-1, -1};
		if (n_annihilatable_ == 0)
		return candidates;

		candidates.push_back(removable_[removal_rand * removable_.size()]);
		// Note:
		// When sampling by retrying in case of failure, the probability of success is p_s n /
		// n_annihlatable, with n = size(). This translates to an expected cost of \sum_l l p_s (1 -
		// p_s)^(l - 1) = 1 / p_s. Therefore this algorithm is faster than a read on all the
		// vertices when n_annihlatable > cost(random _number) / cost(vertex_read).

		// Lets assume cost(random _number) / cost(vertex_read) =~ 10
		// This also helps when testing on a small configuration as we need only one random number.
		constexpr unsigned threshold = 10;

		if (n_annihilatable_ >= threshold) {
		candidates[0] = removal_rand * size();
		while (!vertices_[candidates[0]].annihilatable) {
		candidates[0] = rng() * size();
		}
		}
		else {
		unsigned annihilatable_idx = removal_rand * n_annihilatable_;
		unsigned annihilatable_found = 0;
		for (int i = 0; i < vertices_.size(); ++i) {
		if (vertices_[i].annihilatable) {
		if (annihilatable_found == annihilatable_idx) {
		candidates[0] = i;
		break;
		}
		++annihilatable_found;
		}
		}
		}

		if ((*H_int_)[vertices_[candidates[0]].interaction_id].partners_id.size() &&
		doDoubleUpdate(rng)) { // Double removal.
		if (doDoubleUpdate(rng) &&
		(*H_int_)[vertices_[candidates[0]].interaction_id].partners_id.size()) { // 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));
		candidates[1] = findTag(tag);
		assert(candidates[1] < int(size()) && candidates[1] >= 0);
		assert(vertices_[candidates[1]].annihilatable);
		}

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

		void SolverConfiguration::push_back(const Vertex& v) {
		vertices_.push_back(v);
		BaseClass::addVertex(v);
		template <class Rng>
		bool SolverConfiguration::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_;
		}

		void SolverConfiguration::pop(const int n) {
		assert(n <= (int)size());

		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]);
		template <class Alloc>
		void SolverConfiguration::findIndices(std::vector<int, Alloc>& indices, unsigned config_index,
		int s) const {
		const auto& v = vertices_[config_index];
		for (int leg = 0; leg < 2; ++leg) {
		if (v.spins[leg] == s) {
		indices.push_back(v.matrix_config_indices[leg]);
		}

		std::array<int, 2> SolverConfiguration::sizeIncrease() const {
		assert(last_insertion_size_ > 0);
		std::array<int, 2> result{0, 0};
		for (std::size_t i = size() - last_insertion_size_; i < size(); ++i)
		addSectorSizes(i, result);

		return result;
		}

		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])];
		}

		inline short SolverConfiguration::getSign(const int vertex_index) const {
		@@ -272,117 +250,6 @@ inline double SolverConfiguration::getStrength(int vertex_index) const {
		return (*H_int_)[vertices_[vertex_index].interaction_id].w;
		}

		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>
		inline void SolverConfiguration::moveAndShrink(std::array<std::vector<int, Alloc>, 2>& sector_from,
		std::array<std::vector<int, Alloc>, 2>& sector_remove,
		std::vector<int>& from, std::vector<int>& remove) {
		for (int s = 0; s < 2; ++s) {
		// Sort and prepare source array.
		auto& sector = BaseClass::sectors_[s].entries_;
		std::sort(sector_remove[s].begin(), sector_remove[s].end());
		auto& tags = BaseClass::sectors_[s].tags_;
		sector_from[s].clear();
		int source_idx = sector.size() - sector_remove[s].size();
		for (; source_idx < sector.size(); ++source_idx) {
		while (std::binary_search(sector_remove[s].begin(), sector_remove[s].end(), source_idx))
		++source_idx;
		if (source_idx < sector.size())
		sector_from[s].push_back(source_idx);
		}

		// Move configuration elements.
		for (int i = 0; i < sector_from[s].size(); ++i) {
		sector[sector_remove[s][i]] = sector[sector_from[s][i]];
		tags[sector_remove[s][i]] = tags[sector_from[s][i]];
		}
		// Shrink sector configuration.
		sector.erase(sector.end() - sector_remove[s].size(), sector.end());
		tags.erase(tags.end() - sector_remove[s].size(), tags.end());
		}

		std::sort(remove.begin(), remove.end());
		from.clear();
		int source_idx = size() - remove.size();
		for (; source_idx < size(); ++source_idx) {
		while (std::binary_search(remove.begin(), remove.end(), source_idx))
		++source_idx;
		if (source_idx < size())
		from.push_back(source_idx);
		}

		// Move and shrink configuration.
		for (int i = 0; i < from.size(); ++i)
		vertices_[remove[i]] = vertices_[from[i]];
		vertices_.erase(vertices_.end() - remove.size(), vertices_.end());

		assert(checkConsistency());
		}

		inline bool SolverConfiguration::operator==(const SolverConfiguration& rhs) const {
		bool result = true;
		result &= vertices_ == rhs.vertices_;
		result &= existing_ == rhs.existing_;
		result &= max_tau_ == rhs.max_tau_;
		result &= n_bands_ == rhs.n_bands_;
		result &= double_insertion_prob_ == rhs.double_insertion_prob_;

		result &= static_cast<BaseClass>(*this) == static_cast<BaseClass>(rhs);

		return result;
		}

		void SolverConfiguration::swapVertices(const short i, const short j) {
		std::swap(vertices_[i], vertices_[j]);
		}

		inline bool SolverConfiguration::checkConsistency() const {
		for (auto v : vertices_) {
		int count = 0;
		for (int s = 0; s < 2; ++s) {
		auto indices = findIndices(v.tag, s);
		count += indices.size();
		for (auto idx : indices) {
		const auto& matrix_el = sectors_[s].entries_[idx];
		if (matrix_el.tau_ != v.tau)
		return false;
		}
		}
		if ((double_insertion_prob_ == 0 && count != 2) \|\| !count)
		return false;
		}

		if (double_insertion_prob_) {
		for (const auto& v : vertices_) {
		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;
		}

		inline int SolverConfiguration::findTag(std::uint64_t tag) const {
		for (int i = 0; i < vertices_.size(); ++i)
		if (vertices_[i].tag == tag)
		return i;
		return -1;
		}

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