Fix the summation of the equal time accumulation results. (de330d64) · Commits · NDIP / Tool Sources / Direct-Geometry Spectroscopy / DCA / DCA Main

include/dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp

+272 −175

File changed.

Preview size limit exceeded, changes collapsed.

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

+141 −102

Original line number	Diff line number	Diff line
		@@ -3,7 +3,8 @@
		// All rights reserved.
		//
		// See LICENSE for terms of usage.
		// See CITATION.md for citation guidelines, if DCA++ is used for scientific publications.
		// See CITATION.md for citation guidelines, if DCA++ is used for scientific
		// publications.
		//
		// Author: Peter Staar (taa@zurich.ibm.com)
		// Raffaele Solca' (rasolca@itp.phys.ethz.ch)
		@@ -23,14 +24,14 @@
		#include "dca/function/function.hpp"
		#include "dca/linalg/matrix.hpp"
		#include "dca/linalg/util/cuda_event.hpp"
		#include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator.hpp"
		#include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp"
		#include "dca/phys/dca_step/cluster_solver/ctaux/accumulator/tp/tp_equal_time_accumulator.hpp"
		#include "dca/phys/dca_step/cluster_solver/ctaux/domains/feynman_expansion_order_domain.hpp"
		#include "dca/phys/dca_step/cluster_solver/ctaux/structs/ct_aux_hs_configuration.hpp"
		#include "dca/phys/dca_step/cluster_solver/ctaux/structs/vertex_pair.hpp"
		#include "dca/phys/dca_step/cluster_solver/ctaux/structs/vertex_singleton.hpp"
		#include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/mc_accumulator_data.hpp"
		#include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/sp/sp_accumulator.hpp"
		#include "dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator.hpp"
		#include "dca/phys/domains/cluster/cluster_domain.hpp"
		#include "dca/phys/domains/quantum/electron_band_domain.hpp"
		#include "dca/phys/domains/quantum/electron_spin_domain.hpp"
		@@ -63,7 +64,8 @@ public:

		using t = func::dmn_0<domains::time_domain>;
		using w = func::dmn_0<domains::frequency_domain>;
		using w_VERTEX = func::dmn_0<domains::vertex_frequency_domain<domains::COMPACT>>;
		using w_VERTEX =
		func::dmn_0<domains::vertex_frequency_domain<domains::COMPACT>>;

		using b = func::dmn_0<domains::electron_band_domain>;
		using s = func::dmn_0<domains::electron_spin_domain>;
		@@ -83,29 +85,31 @@ public:

		CtauxAccumulator(Parameters &parameters_ref, Data &data_ref, int id);

		template <typename Writer>
		void write(Writer& writer);
		template <typename Writer> void write(Writer &writer);

		void initialize(int dca_iteration);

		template <typename walker_type>
		void updateFrom(walker_type& walker);
		template <typename walker_type> void updateFrom(walker_type &walker);

		void measure();

		// Sums all accumulated objects of this accumulator to the equivalent objects of the 'other'
		// Sums all accumulated objects of this accumulator to the equivalent objects
		// of the 'other'
		// accumulator.
		void sumTo(this_type &other);

		void finalize();

		std::vector<vertex_singleton_type>& get_configuration(e_spin_states_type e_spin = e_UP);
		std::vector<vertex_singleton_type> &
		get_configuration(e_spin_states_type e_spin = e_UP);

		func::function<double, func::dmn_0<domains::numerical_error_domain>>& get_error_distribution() {
		func::function<double, func::dmn_0<domains::numerical_error_domain>> &
		get_error_distribution() {
		return error;
		}

		func::function<double, func::dmn_0<Feynman_expansion_order_domain>>& get_visited_expansion_order_k() {
		func::function<double, func::dmn_0<Feynman_expansion_order_domain>> &
		get_visited_expansion_order_k() {
		return visited_expansion_order_k;
		}

		@@ -113,17 +117,21 @@ public:
		func::function<double, func::dmn_variadic<nu, nu, r_dmn_t, t>> &get_G_r_t() {
		return G_r_t;
		}
		func::function<double, func::dmn_variadic<nu, nu, r_dmn_t, t>>& get_G_r_t_stddev() {
		func::function<double, func::dmn_variadic<nu, nu, r_dmn_t, t>> &
		get_G_r_t_stddev() {
		return G_r_t_stddev;
		}

		func::function<double, func::dmn_variadic<b, r_dmn_t>>& get_charge_cluster_moment() {
		func::function<double, func::dmn_variadic<b, r_dmn_t>> &
		get_charge_cluster_moment() {
		return charge_cluster_moment;
		}
		func::function<double, func::dmn_variadic<b, r_dmn_t>>& get_magnetic_cluster_moment() {
		func::function<double, func::dmn_variadic<b, r_dmn_t>> &
		get_magnetic_cluster_moment() {
		return magnetic_cluster_moment;
		}
		func::function<double, func::dmn_variadic<b, r_dmn_t>>& get_dwave_pp_correlator() {
		func::function<double, func::dmn_variadic<b, r_dmn_t>> &
		get_dwave_pp_correlator() {
		return dwave_pp_correlator;
		}

		@@ -141,9 +149,7 @@ public:
		return two_particle_accumulator_.get_sign_times_G4();
		}

		bool compute_std_deviation() const {
		return compute_std_deviation_;
		}
		bool compute_std_deviation() const { return compute_std_deviation_; }

		#ifdef MEASURE_ERROR_BARS
		void store_standard_deviation(int nr_measurements, std::ofstream &points_file,
		@@ -158,7 +164,8 @@ public:
		}

		static std::size_t staticDeviceFingerprint() {
		return accumulator::TpAccumulator<Parameters, device_t>::staticDeviceFingerprint();
		return accumulator::TpAccumulator<Parameters,
		device_t>::staticDeviceFingerprint();
		}

		private:
		@@ -166,8 +173,10 @@ private:

		void accumulate_equal_time_quantities();
		using AccumType = typename Parameters::MC_measurement_scalar_type;
		void accumulate_equal_time_quantities(const std::array<linalg::Matrix<AccumType, linalg::GPU>, 2>& M);
		void accumulate_equal_time_quantities(const std::array<linalg::Matrix<AccumType, linalg::CPU>, 2>& M);
		void accumulate_equal_time_quantities(
		const std::array<linalg::Matrix<AccumType, linalg::GPU>, 2> &M);
		void accumulate_equal_time_quantities(
		const std::array<linalg::Matrix<AccumType, linalg::CPU>, 2> &M);

		void accumulate_two_particle_quantities();

		@@ -194,20 +203,25 @@ protected:
		std::array<dca::linalg::Matrix<AccumType, linalg::CPU>, 2> M_host_;

		func::function<double, func::dmn_0<domains::numerical_error_domain>> error;
		func::function<double, func::dmn_0<Feynman_expansion_order_domain>> visited_expansion_order_k;
		func::function<double, func::dmn_0<Feynman_expansion_order_domain>>
		visited_expansion_order_k;

		func::function<double, func::dmn_variadic<nu, nu, r_dmn_t, t>> G_r_t;
		func::function<double, func::dmn_variadic<nu, nu, r_dmn_t, t>> G_r_t_stddev;

		func::function<double, func::dmn_variadic<b, r_dmn_t>> charge_cluster_moment;
		func::function<double, func::dmn_variadic<b, r_dmn_t>> magnetic_cluster_moment;
		func::function<double, func::dmn_variadic<b, r_dmn_t>>
		magnetic_cluster_moment;
		func::function<double, func::dmn_variadic<b, r_dmn_t>> dwave_pp_correlator;

		func::function<std::complex<double>, func::dmn_variadic<nu, nu, r_dmn_t, w>> M_r_w_stddev;
		func::function<std::complex<double>, func::dmn_variadic<nu, nu, r_dmn_t, w>>
		M_r_w_stddev;

		accumulator::SpAccumulator<Parameters, device_t> single_particle_accumulator_obj;
		accumulator::SpAccumulator<Parameters, device_t>
		single_particle_accumulator_obj;

		ctaux::TpEqualTimeAccumulator<Parameters, Data> MC_two_particle_equal_time_accumulator_obj;
		ctaux::TpEqualTimeAccumulator<Parameters, Data>
		MC_two_particle_equal_time_accumulator_obj;

		accumulator::TpAccumulator<Parameters, device_t> two_particle_accumulator_;

		@@ -215,12 +229,11 @@ protected:
		};

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		CtauxAccumulator<device_t, Parameters, Data>::CtauxAccumulator(Parameters& parameters_ref,
		Data& data_ref, int id)
		CtauxAccumulator<device_t, Parameters, Data>::CtauxAccumulator(
		Parameters &parameters_ref, Data &data_ref, int id)
		: MC_accumulator_data(),

		parameters_(parameters_ref),
		data_(data_ref),
		parameters_(parameters_ref), data_(data_ref),
		concurrency(parameters_.get_concurrency()),

		thread_id(id),
		@@ -231,8 +244,7 @@ CtauxAccumulator<device_t, Parameters, Data>::CtauxAccumulator(Parameters& param
		error("numerical-error-distribution-of-N-matrices"),
		visited_expansion_order_k("<k>"),

		G_r_t("G_r_t_measured"),
		G_r_t_stddev("G_r_t_stddev"),
		G_r_t("G_r_t_measured"), G_r_t_stddev("G_r_t_stddev"),

		charge_cluster_moment("charge-cluster-moment"),
		magnetic_cluster_moment("magnetic-cluster-moment"),
		@@ -247,11 +259,14 @@ CtauxAccumulator<device_t, Parameters, Data>::CtauxAccumulator(Parameters& param
		two_particle_accumulator_(data_.G0_k_w_cluster_excluded, parameters_) {}

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		void CtauxAccumulator<device_t, Parameters, Data>::initialize(int dca_iteration) {
		// Note: profiling this function breaks the PAPI profiler as both the master thread and the first
		void CtauxAccumulator<device_t, Parameters, Data>::initialize(
		int dca_iteration) {
		// Note: profiling this function breaks the PAPI profiler as both the master
		// thread and the first
		// worker call this with the same thread_id.
		// TODO: fix thread id assignment.
		// profiler_type profiler(__FUNCTION__, "CT-AUX accumulator", __LINE__, thread_id);
		// profiler_type profiler(__FUNCTION__, "CT-AUX accumulator", __LINE__,
		// thread_id);

		MC_accumulator_data::initialize(dca_iteration);

		@@ -288,9 +303,11 @@ void CtauxAccumulator<device_t, Parameters, Data>::finalize() {

		if (compute_std_deviation_) {
		const auto &M_r_w = single_particle_accumulator_obj.get_sign_times_M_r_w();
		const auto& M_r_w_squared = single_particle_accumulator_obj.get_sign_times_M_r_w_sqr();
		const auto &M_r_w_squared =
		single_particle_accumulator_obj.get_sign_times_M_r_w_sqr();
		for (int l = 0; l < M_r_w_stddev.size(); l++)
		M_r_w_stddev(l) = std::sqrt(abs(M_r_w_squared(l)) - std::pow(abs(M_r_w(l)), 2));
		M_r_w_stddev(l) =
		std::sqrt(abs(M_r_w_squared(l)) - std::pow(abs(M_r_w(l)), 2));

		double factor = 1. / std::sqrt(parameters_.get_measurements() - 1);

		@@ -298,15 +315,18 @@ void CtauxAccumulator<device_t, Parameters, Data>::finalize() {
		}

		if (parameters_.additional_time_measurements()) {
		MC_two_particle_equal_time_accumulator_obj.finalize(); // G_r_t, G_r_t_stddev);
		MC_two_particle_equal_time_accumulator_obj.finalize();

		G_r_t = MC_two_particle_equal_time_accumulator_obj.get_G_r_t();
		G_r_t_stddev = MC_two_particle_equal_time_accumulator_obj.get_G_r_t_stddev();

		charge_cluster_moment = MC_two_particle_equal_time_accumulator_obj.get_charge_cluster_moment();
		magnetic_cluster_moment =
		MC_two_particle_equal_time_accumulator_obj.get_magnetic_cluster_moment();
		dwave_pp_correlator = MC_two_particle_equal_time_accumulator_obj.get_dwave_pp_correlator();
		G_r_t_stddev =
		MC_two_particle_equal_time_accumulator_obj.get_G_r_t_stddev();

		charge_cluster_moment =
		MC_two_particle_equal_time_accumulator_obj.get_charge_cluster_moment();
		magnetic_cluster_moment = MC_two_particle_equal_time_accumulator_obj
		.get_magnetic_cluster_moment();
		dwave_pp_correlator =
		MC_two_particle_equal_time_accumulator_obj.get_dwave_pp_correlator();
		}

		if (perform_tp_accumulation_)
		@@ -314,7 +334,8 @@ void CtauxAccumulator<device_t, Parameters, Data>::finalize() {
		}

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		std::vector<vertex_singleton>& CtauxAccumulator<device_t, Parameters, Data>::get_configuration(
		std::vector<vertex_singleton> &
		CtauxAccumulator<device_t, Parameters, Data>::get_configuration(
		e_spin_states_type e_spin) {
		if (e_spin == e_UP)
		return hs_configuration_[0];
		@@ -349,7 +370,8 @@ void CtauxAccumulator<device_t, Parameters, Data>::write(Writer& writer) {
		}

		/*!
		* \brief Get all the information from the walker in order to start a measurement.
		* \brief Get all the information from the walker in order to start a
		* measurement.
		*
		* \f{eqnarray}{
		* M_{i,j} &=& (e^{V_i}-1) N_{i,j}
		@@ -357,8 +379,10 @@ void CtauxAccumulator<device_t, Parameters, Data>::write(Writer& writer) {
		*/
		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		template <typename walker_type>
		void CtauxAccumulator<device_t, Parameters, Data>::updateFrom(walker_type& walker) {
		profiler_type profiler("update from", "CT-AUX accumulator", __LINE__, thread_id);
		void CtauxAccumulator<device_t, Parameters, Data>::updateFrom(
		walker_type &walker) {
		profiler_type profiler("update from", "CT-AUX accumulator", __LINE__,
		thread_id);

		GFLOP += walker.get_Gflop();

		@@ -406,21 +430,26 @@ void CtauxAccumulator<device_t, Parameters, Data>::measure() {
		/*!
		* \brief Output and store standard deviation and error.
		*
		* It computes and write to the given files the standard deviation of the measurements of the one
		* It computes and write to the given files the standard deviation of the
		* measurements of the one
		* particle accumulator.
		* It outputs the L1-Norm, i.e. \f$\sum_{i=1}^N \left\|x_i\right\|/N\f$, the L2-Norm, i.e.
		* It outputs the L1-Norm, i.e. \f$\sum_{i=1}^N \left\|x_i\right\|/N\f$, the
		* L2-Norm, i.e.
		* \f$\sqrt{\sum_{i=1}^N \left\|x_i\right\|^2/N}\f$,
		* and the Linf-Norm, i.e. \f$\max_{i=1}^N \left\|x_i\right\|\f$ of the standard deviation and of the
		* and the Linf-Norm, i.e. \f$\max_{i=1}^N \left\|x_i\right\|\f$ of the standard
		* deviation and of the
		* error.
		*/
		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		void CtauxAccumulator<device_t, Parameters, Data>::store_standard_deviation(
		int nr_measurements, std::ofstream &points_file, std::ofstream &norm_file) {
		single_particle_accumulator_obj.store_standard_deviation(nr_measurements, points_file, norm_file);
		single_particle_accumulator_obj.store_standard_deviation(
		nr_measurements, points_file, norm_file);
		}

		/*!
		* \brief Update the sum of the squares of the measurements of the single particle accumulator.
		* \brief Update the sum of the squares of the measurements of the single
		* particle accumulator.
		* It has to be called after each measurement.
		*/
		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		@@ -436,10 +465,13 @@ void CtauxAccumulator<device_t, Parameters, Data>::update_sum_squares() {
		*************************************************************/

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		void CtauxAccumulator<device_t, Parameters, Data>::accumulate_single_particle_quantities() {
		profiler_type profiler("sp-accumulation", "CT-AUX accumulator", __LINE__, thread_id);
		void CtauxAccumulator<device_t, Parameters,
		Data>::accumulate_single_particle_quantities() {
		profiler_type profiler("sp-accumulation", "CT-AUX accumulator", __LINE__,
		thread_id);

		single_particle_accumulator_obj.accumulate(M_, hs_configuration_, current_sign);
		single_particle_accumulator_obj.accumulate(M_, hs_configuration_,
		current_sign);

		GFLOP += 2. * 8. * M_[1].size().first * M_[1].size().first * (1.e-9);
		GFLOP += 2. * 8. * M_[0].size().first * M_[0].size().first * (1.e-9);
		@@ -452,14 +484,17 @@ void CtauxAccumulator<device_t, Parameters, Data>::accumulate_single_particle_qu
		*************************************************************/

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		void CtauxAccumulator<device_t, Parameters, Data>::accumulate_equal_time_quantities() {
		profiler_type profiler("equal-time-measurements", "CT-AUX accumulator", __LINE__, thread_id);
		void CtauxAccumulator<device_t, Parameters,
		Data>::accumulate_equal_time_quantities() {
		profiler_type profiler("equal-time-measurements", "CT-AUX accumulator",
		__LINE__, thread_id);

		return accumulate_equal_time_quantities(M_);
		}

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		void CtauxAccumulator<device_t, Parameters, Data>::accumulate_equal_time_quantities(
		void CtauxAccumulator<device_t, Parameters, Data>::
		accumulate_equal_time_quantities(
		const std::array<linalg::Matrix<AccumType, linalg::GPU>, 2> &M) {
		for (int s = 0; s < 2; ++s)
		M_host_[s].setAsync(M[s], thread_id, s);
		@@ -470,16 +505,18 @@ void CtauxAccumulator<device_t, Parameters, Data>::accumulate_equal_time_quantit
		}

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		void CtauxAccumulator<device_t, Parameters, Data>::accumulate_equal_time_quantities(
		void CtauxAccumulator<device_t, Parameters, Data>::
		accumulate_equal_time_quantities(
		const std::array<linalg::Matrix<AccumType, linalg::CPU>, 2> &M) {
		MC_two_particle_equal_time_accumulator_obj.compute_G_r_t(hs_configuration_[0], M[0],
		hs_configuration_[1], M[1]);
		MC_two_particle_equal_time_accumulator_obj.compute_G_r_t(
		hs_configuration_[0], M[0], hs_configuration_[1], M[1]);

		MC_two_particle_equal_time_accumulator_obj.accumulate_G_r_t(current_sign);

		MC_two_particle_equal_time_accumulator_obj.accumulate_moments(current_sign);

		MC_two_particle_equal_time_accumulator_obj.accumulate_dwave_pp_correlator(current_sign);
		MC_two_particle_equal_time_accumulator_obj.accumulate_dwave_pp_correlator(
		current_sign);

		GFLOP += MC_two_particle_equal_time_accumulator_obj.get_GFLOP();
		}
		@@ -491,9 +528,12 @@ void CtauxAccumulator<device_t, Parameters, Data>::accumulate_equal_time_quantit
		*************************************************************/

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		void CtauxAccumulator<device_t, Parameters, Data>::accumulate_two_particle_quantities() {
		profiler_type profiler("tp-accumulation", "CT-AUX accumulator", __LINE__, thread_id);
		/GFLOP +=/two_particle_accumulator_.accumulate(M_, hs_configuration_, current_sign);
		void CtauxAccumulator<device_t, Parameters,
		Data>::accumulate_two_particle_quantities() {
		profiler_type profiler("tp-accumulation", "CT-AUX accumulator", __LINE__,
		thread_id);
		/GFLOP +=/two_particle_accumulator_.accumulate(M_, hs_configuration_,
		current_sign);
		}

		template <dca::linalg::DeviceType device_t, class Parameters, class Data>
		@@ -506,16 +546,15 @@ void CtauxAccumulator<device_t, Parameters, Data>::sumTo(this_type& other) {
		other.get_visited_expansion_order_k() += visited_expansion_order_k;
		other.get_error_distribution() += error;

		// equal time measurements
		other.get_G_r_t() += G_r_t;
		other.get_G_r_t_stddev() += G_r_t_stddev;
		other.get_charge_cluster_moment() += charge_cluster_moment;
		other.get_magnetic_cluster_moment() += magnetic_cluster_moment;
		other.get_dwave_pp_correlator() += dwave_pp_correlator;

		// sp-measurements
		single_particle_accumulator_obj.sumTo(other.single_particle_accumulator_obj);

		// equal time measurements
		if (DCA_iteration == parameters_.get_dca_iterations() - 1 &&
		parameters_.additional_time_measurements())
		MC_two_particle_equal_time_accumulator_obj.sumTo(
		other.MC_two_particle_equal_time_accumulator_obj);

		// tp-measurements
		if (perform_tp_accumulation_)
		two_particle_accumulator_.sumTo(other.two_particle_accumulator_);