Skip to content

Correct finite-Q multiband G4 accumulation #373

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
PeizhiMai wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Correct finite-Q multiband G4 accumulation #373

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 18 additions & 1 deletion include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cuh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,9 @@ public:
*/
static void set(int nb, int nk, int nw,
const std::vector<int>& delta_k, const std::vector<int>& delta_w,
const int extension_offset, const int* add_k, int lda, const int* sub_k, int lds);
const int extension_offset, const int* add_k, int lda, const int* sub_k, int lds,
const std::vector<double>& q_minus_k_phase_real,
const std::vector<double>& q_minus_k_phase_imag);

__device__ auto get_bands() const {
return nb_;
Expand All @@ -57,6 +59,10 @@ public:
// Returns the index of w_ex - w.
__device__ inline int wexMinus(int w_idx, int w_ex_idx) const;

// Returns the phase factor exp(i * (K_ex - K - fold(K_ex - K)) . a_band).
__device__ inline double qMinusKPhaseReal(int k_idx, int k_ex_idx, int band) const;
__device__ inline double qMinusKPhaseImag(int k_idx, int k_ex_idx, int band) const;

// Maps the indices w1 w2 from the compact frequency domain of G4,
// to the extended (positive for w1) domain used by G.
// k1, k2 mapped to minusK(k1), minus(k2)
Expand All @@ -81,6 +87,8 @@ protected:

const int* w_ex_indices_;
const int* k_ex_indices_;
const double* q_minus_k_phase_real_;
const double* q_minus_k_phase_imag_;
int ext_size_;

int nw_;
Expand Down Expand Up @@ -122,6 +130,15 @@ inline __device__ int G4Helper::kMinus(const int k_idx) const {
return solver::details::cluster_momentum_helper.minus(k_idx);
}

inline __device__ double G4Helper::qMinusKPhaseReal(const int k_idx, const int k_ex_idx,
const int band) const {
return q_minus_k_phase_real_[band + nb_ * (k_idx + nc_ * k_ex_idx)];
}

inline __device__ double G4Helper::qMinusKPhaseImag(const int k_idx, const int k_ex_idx,
const int band) const {
return q_minus_k_phase_imag_[band + nb_ * (k_idx + nc_ * k_ex_idx)];
}

inline __device__ void G4Helper::unrollIndex(std::size_t index, int& b1, int& b2,
int& b3, int& b4, int& k1, int& w1,
Expand Down
95 changes: 77 additions & 18 deletions include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_cpu.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -173,6 +173,16 @@ class TpAccumulator<Parameters, DT, linalg::CPU> : public TpAccumulatorBase<Para

auto getGSingleband(int s, int k1, int k2, int w1, int w2) -> TpComplex const;

using Momentum = typename KDmn::parameter_type::element_type;

static Momentum qMinusKVector(int k, int q);

void getGMultibandAtMomenta(int s, const Momentum& k1, int k1_folded, const Momentum& k2,
int k2_folded, int w1, int w2, Matrix& G);

void applyFoldedMomentumPhase(Matrix& G, const Momentum& k1, int k1_folded,
const Momentum& k2, int k2_folded) const;

template <class Configuration, typename SpScalar>
double computeM(const std::array<linalg::Matrix<SpScalar, linalg::CPU>, 2>& M_pair,
const std::array<Configuration, 2>& configs);
Expand Down Expand Up @@ -206,13 +216,20 @@ class TpAccumulator<Parameters, DT, linalg::CPU> : public TpAccumulatorBase<Para
private:
// work spaces for computeGMultiband.
Matrix G0_M_, G_a_, G_b_;

func::function<TpComplex, func::dmn_variadic<RDmn, RDmn, BDmn, BDmn, SDmn, WTpExtPosDmn, WTpExtDmn>>
M_r_r_w_w_;
};

template <class Parameters, DistType DT>
TpAccumulator<Parameters, DT, linalg::CPU>::TpAccumulator(
const func::function<TpComplex, func::dmn_variadic<NuDmn, NuDmn, KDmn, WDmn>>& G0,
const Parameters& pars, const int thread_id)
: Base(G0, pars, thread_id), G0_M_(n_bands_), G_a_(n_bands_), G_b_(n_bands_) {
: Base(G0, pars, thread_id),
G0_M_(n_bands_),
G_a_(n_bands_),
G_b_(n_bands_),
M_r_r_w_w_("M_r_r_w_w") {
if constexpr (DT == DistType::BLOCKED) {
std::cerr << "Blocked distribution is not supported in the CPU accumulator. "
<< "Reverting to no distribution.\n";
Expand Down Expand Up @@ -257,20 +274,20 @@ double TpAccumulator<Parameters, DT, linalg::CPU>::computeM(
const std::array<linalg::Matrix<SpScalar, linalg::CPU>, 2>& M_pair,
const std::array<Configuration, 2>& configs) {
double flops = 0.;

func::function<TpComplex, func::dmn_variadic<RDmn, RDmn, BDmn, BDmn, SDmn, WTpExtPosDmn, WTpExtDmn>>
M_r_r_w_w;
M_r_r_w_w_ = TpComplex(0., 0.);

for (int spin = 0; spin < SDmn::dmn_size(); ++spin) {
Profiler prf_a("Frequency FT", "tp-accumulation", __LINE__, thread_id_);
if (not configs[spin].size())
continue;
flops += ndft_obj_.execute(configs[spin], M_pair[spin], M_r_r_w_w, spin);
flops += ndft_obj_.execute(configs[spin], M_pair[spin], M_r_r_w_w_, spin);
}

Profiler prf_b("Space FT", "tp-accumulation", __LINE__, thread_id_);
// TODO: add the gflops here.
math::transform::SpaceTransform2D<RDmn, KDmn, BDmn, SDmn, TpPrecision>::execute(M_r_r_w_w, G_);
auto M_r_r_w_w_transformed = M_r_r_w_w_;
math::transform::SpaceTransform2D<RDmn, KDmn, BDmn, SDmn, TpPrecision>::execute(
M_r_r_w_w_transformed, G_);

return flops;
}
Expand Down Expand Up @@ -385,6 +402,45 @@ auto TpAccumulator<Parameters, DT, linalg::CPU>::getGSingleband(const int s, con
return G_(0, 0, s, k1, k2, w1_ext, w2_ext);
}

template <class Parameters, DistType DT>
auto TpAccumulator<Parameters, DT, linalg::CPU>::qMinusKVector(const int k, const int q)
-> Momentum {
const auto& k_elements = KDmn::parameter_type::get_elements();
const auto& q_vec = k_elements[q];
const auto& k_vec = k_elements[k];

Momentum q_minus_k_vec(q_vec.size(), 0);
for (std::size_t d = 0; d < q_vec.size(); ++d)
q_minus_k_vec[d] = q_vec[d] - k_vec[d];

return q_minus_k_vec;
}

template <class Parameters, DistType DT>
void TpAccumulator<Parameters, DT, linalg::CPU>::applyFoldedMomentumPhase(
Matrix& G, const Momentum& k1, const int k1_folded, const Momentum& k2, const int k2_folded) const {
const auto& k_elements = KDmn::parameter_type::get_elements();
const auto& bands = BDmn::get_elements();
const auto& K1 = k_elements[k1_folded];
const auto& K2 = k_elements[k2_folded];

for (int b2 = 0; b2 < n_bands_; ++b2)
for (int b1 = 0; b1 < n_bands_; ++b1) {
Real phase = 0;
for (std::size_t d = 0; d < k1.size(); ++d)
phase += (k1[d] - K1[d]) * bands[b1].a_vec[d] - (k2[d] - K2[d]) * bands[b2].a_vec[d];
G(b1, b2) *= std::exp(TpComplex(0, phase));
}
}

template <class Parameters, DistType DT>
void TpAccumulator<Parameters, DT, linalg::CPU>::getGMultibandAtMomenta(
const int s, const Momentum& k1, const int k1_folded, const Momentum& k2, const int k2_folded,
const int w1, const int w2, Matrix& G) {
getGMultiband(s, k1_folded, k2_folded, w1, w2, G);
applyFoldedMomentumPhase(G, k1, k1_folded, k2, k2_folded);
}

template <class Parameters, DistType DT>
void TpAccumulator<Parameters, DT, linalg::CPU>::getGMultiband(int s, int k1, int k2, int w1,
int w2, Matrix& G,
Expand All @@ -395,7 +451,7 @@ void TpAccumulator<Parameters, DT, linalg::CPU>::getGMultiband(int s, int k1, in
const auto* const G_ptr = &G_(0, 0, s, k1, k2, w1_ext, w2_ext);
for (int b2 = 0; b2 < n_bands_; ++b2)
for (int b1 = 0; b1 < n_bands_; ++b1) {
G(b1, b2) = sign * G(b1, b2) + G_ptr[b2 + b1 * n_bands_];
G(b1, b2) = sign * G(b1, b2) + G_ptr[b1 + b2 * n_bands_];
#ifndef NDEBUG
// if (std::abs(G(b1, b2).imag()) > 10) // std::isnan(imag(G_(b1, b2, s, k1, k2, w1_ext, w2_ext))))
// std::cout << w1 << "," << w2 << "," << k1 << "," << k2 << "," << b1 << "," << b2 << ","
Expand All @@ -419,6 +475,7 @@ double TpAccumulator<Parameters, DT, linalg::CPU>::updateG4(const int channel_id

auto momentum_sum = [](const int k, const int q) { return KDmn::parameter_type::add(k, q); };
auto q_minus_k = [](const int k, const int q) { return KDmn::parameter_type::subtract(k, q); };
auto q_minus_k_vector = [](const int k, const int q) { return qMinusKVector(k, q); };
auto q_plus_k = [](const int k, const int q) { return KDmn::parameter_type::add(k, q); };
// Returns the index of the exchange frequency w_ex plus the Matsubara frequency with index w.
auto w_plus_w_ex = [](const int w, const int w_ex) { return w + w_ex; };
Expand Down Expand Up @@ -636,15 +693,11 @@ double TpAccumulator<Parameters, DT, linalg::CPU>::updateG4(const int channel_id
for (int k2 = 0; k2 < KDmn::dmn_size(); ++k2)
for (int w1 = 0; w1 < WTpDmn::dmn_size(); ++w1)
for (int k1 = 0; k1 < KDmn::dmn_size(); ++k1) {
// TpComplex* const G4_ptr = &G4(0, 0, 0, 0, k1, w1, k2, w2, k_ex_idx, w_ex_idx);
for (int s = 0; s < 2; ++s) {
// updateG4Atomic(G4_ptr, s, k1, k2, w1, w2, !s, q_minus_k(k1, k_ex),
// q_minus_k(k2, k_ex), w_ex_minus_w(w1, w_ex),
// w_ex_minus_w(w2, w_ex), sign_over_2, false);
// contraction: G_{b1,b3}(k1, k2) * G_{b2,b4}(q-k1, q-k2).
getGMultiband(s, k1, k2, w1, w2, G_a_);
getGMultiband(!s, q_minus_k(k1, k_ex), q_minus_k(k2, k_ex),
w_ex_minus_w(w1, w_ex), w_ex_minus_w(w2, w_ex), G_b_);
getGMultibandAtMomenta(!s, q_minus_k_vector(k1, k_ex), q_minus_k(k1, k_ex),
q_minus_k_vector(k2, k_ex), q_minus_k(k2, k_ex),
w_ex_minus_w(w1, w_ex), w_ex_minus_w(w2, w_ex), G_b_);
for (int b4 = 0; b4 < BDmn::dmn_size(); ++b4)
for (int b3 = 0; b3 < BDmn::dmn_size(); ++b3)
for (int b2 = 0; b2 < BDmn::dmn_size(); ++b2)
Expand Down Expand Up @@ -678,8 +731,10 @@ double TpAccumulator<Parameters, DT, linalg::CPU>::updateG4(const int channel_id
for (int k1 = 0; k1 < KDmn::dmn_size(); ++k1) {
// contraction: G(k2, k1, b1, b3) * G(k_ex - k2, k_ex - k1, b2, b4).
getGMultiband(0, k1, k2, w1, w2, G_a_);
getGMultiband(0, q_minus_k(k1, k_ex), q_minus_k(k2, k_ex),
w_ex_minus_w(w1, w_ex), w_ex_minus_w(w2, w_ex), G_b_);
getGMultibandAtMomenta(0, q_minus_k_vector(k1, k_ex),
q_minus_k(k1, k_ex), q_minus_k_vector(k2, k_ex),
q_minus_k(k2, k_ex), w_ex_minus_w(w1, w_ex),
w_ex_minus_w(w2, w_ex), G_b_);
for (int b4 = 0; b4 < BDmn::dmn_size(); ++b4)
for (int b3 = 0; b3 < BDmn::dmn_size(); ++b3)
for (int b2 = 0; b2 < BDmn::dmn_size(); ++b2)
Expand All @@ -689,8 +744,12 @@ double TpAccumulator<Parameters, DT, linalg::CPU>::updateG4(const int channel_id
}

// contraction: -G(k1, k_ex - k2, b1, b4) * G(k_ex - k1, k2, b2, b3).
getGMultiband(0, k1, q_minus_k(k2, k_ex), w1, w_ex_minus_w(w2, w_ex), G_a_);
getGMultiband(0, q_minus_k(k1, k_ex), k2, w_ex_minus_w(w1, w_ex), w2, G_b_);
getGMultibandAtMomenta(0, KDmn::parameter_type::get_elements()[k1], k1,
q_minus_k_vector(k2, k_ex), q_minus_k(k2, k_ex), w1,
w_ex_minus_w(w2, w_ex), G_a_);
getGMultibandAtMomenta(0, q_minus_k_vector(k1, k_ex), q_minus_k(k1, k_ex),
KDmn::parameter_type::get_elements()[k2], k2,
w_ex_minus_w(w1, w_ex), w2, G_b_);
for (int b4 = 0; b4 < BDmn::dmn_size(); ++b4)
for (int b3 = 0; b3 < BDmn::dmn_size(); ++b3)
for (int b2 = 0; b2 < BDmn::dmn_size(); ++b2)
Expand Down
30 changes: 29 additions & 1 deletion ...dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_base.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -165,10 +165,38 @@ void TpAccumulatorGpuBase<Parameters, DT>::initializeG4Helpers() const {
const auto& q_indices = domains::MomentumExchangeDomain::get_elements();
const auto extension_offset = (WTpExtDmn::dmn_size() - WTpDmn::dmn_size()) / 2;

const auto& k_elements = KDmn::parameter_type::get_elements();
const auto& bands = BDmn::get_elements();
std::vector<double> q_minus_k_phase_real(q_indices.size() * KDmn::dmn_size() * n_bands_);
std::vector<double> q_minus_k_phase_imag(q_minus_k_phase_real.size());

auto phaseIndex = [this](const int q_pos, const int k, const int band) {
return band + n_bands_ * (k + KDmn::dmn_size() * q_pos);
};

for (int q_pos = 0; q_pos < static_cast<int>(q_indices.size()); ++q_pos) {
const int q = q_indices[q_pos];
const auto& q_vec = k_elements[q];
for (int k = 0; k < KDmn::dmn_size(); ++k) {
const auto& k_vec = k_elements[k];
const int folded_q_minus_k = KDmn::parameter_type::subtract(k, q);
const auto& folded_vec = k_elements[folded_q_minus_k];
for (int band = 0; band < n_bands_; ++band) {
double phase = 0.;
for (std::size_t d = 0; d < q_vec.size(); ++d)
phase += (q_vec[d] - k_vec[d] - folded_vec[d]) * bands[band].a_vec[d];

const auto idx = phaseIndex(q_pos, k, band);
q_minus_k_phase_real[idx] = std::cos(phase);
q_minus_k_phase_imag[idx] = std::sin(phase);
}
}
}

// CurrentlyA WTpPosDmn should always be == WTpDmn
details::G4Helper::set(n_bands_, KDmn::dmn_size(), WTpDmn::dmn_size(), q_indices, w_indices,
extension_offset, add_mat.ptr(), add_mat.leadingDimension(), sub_mat.ptr(),
sub_mat.leadingDimension());
sub_mat.leadingDimension(), q_minus_k_phase_real, q_minus_k_phase_imag);
assert(cudaPeekAtLastError() == cudaSuccess);
}

Expand Down
24 changes: 20 additions & 4 deletions src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,9 @@

#include <algorithm>
#include <array>
#include <iostream>
#include <mutex>
#include <stdexcept>
#include <iostream>

#include "dca/platform/dca_gpu.h"

namespace dca {
namespace phys {
Expand All @@ -33,7 +31,9 @@ __CONSTANT__ int* k_ex_indices_actual;

void G4Helper::set(int nb, int nk, int nw, const std::vector<int>& delta_k,
const std::vector<int>& delta_w, const int extension_offset, const int* add_k,
int lda, const int* sub_k, int lds) {
int lda, const int* sub_k, int lds,
const std::vector<double>& q_minus_k_phase_real,
const std::vector<double>& q_minus_k_phase_imag) {
// Initialize the reciprocal cluster if not done already.
solver::details::ClusterHelper::setMomentum(nk, add_k, lda, sub_k, lds);

Expand Down Expand Up @@ -72,6 +72,22 @@ void G4Helper::set(int nb, int nk, int nw, const std::vector<int>& delta_k,
checkRC(cudaMemcpy(const_cast<int*>(host_helper.k_ex_indices_), delta_k.data(),
sizeof(int) * delta_k.size(), cudaMemcpyHostToDevice));

const std::size_t phase_table_size =
static_cast<std::size_t>(nb) * static_cast<std::size_t>(nk) * delta_k.size();
if (q_minus_k_phase_real.size() != phase_table_size ||
q_minus_k_phase_imag.size() != phase_table_size)
throw(std::logic_error("Unexpected G4 q-k phase table size."));

cudaMalloc(&host_helper.q_minus_k_phase_real_, sizeof(double) * q_minus_k_phase_real.size());
checkRC(cudaMemcpy(const_cast<double*>(host_helper.q_minus_k_phase_real_),
q_minus_k_phase_real.data(), sizeof(double) * q_minus_k_phase_real.size(),
cudaMemcpyHostToDevice));

cudaMalloc(&host_helper.q_minus_k_phase_imag_, sizeof(double) * q_minus_k_phase_imag.size());
checkRC(cudaMemcpy(const_cast<double*>(host_helper.q_minus_k_phase_imag_),
q_minus_k_phase_imag.data(), sizeof(double) * q_minus_k_phase_imag.size(),
cudaMemcpyHostToDevice));

#ifndef NDEBUG
checkRC(cudaMalloc(&host_helper.bad_indicies_, sizeof(int) * 1024));
#endif
Expand Down
Loading
Loading