diff --git a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cuh b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cuh index 8da2658ac..1e3984d59 100644 --- a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cuh +++ b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cuh @@ -37,7 +37,9 @@ public: */ static void set(int nb, int nk, int nw, const std::vector& delta_k, const std::vector& 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& q_minus_k_phase_real, + const std::vector& q_minus_k_phase_imag); __device__ auto get_bands() const { return nb_; @@ -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) @@ -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_; @@ -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, diff --git a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_cpu.hpp b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_cpu.hpp index 77e34f9ab..cb712c2cc 100644 --- a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_cpu.hpp +++ b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_cpu.hpp @@ -173,6 +173,16 @@ class TpAccumulator : public TpAccumulatorBase 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 double computeM(const std::array, 2>& M_pair, const std::array& configs); @@ -206,13 +216,20 @@ class TpAccumulator : public TpAccumulatorBase> + M_r_r_w_w_; }; template TpAccumulator::TpAccumulator( const func::function>& 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"; @@ -257,20 +274,20 @@ double TpAccumulator::computeM( const std::array, 2>& M_pair, const std::array& configs) { double flops = 0.; - - func::function> - 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::execute(M_r_r_w_w, G_); + auto M_r_r_w_w_transformed = M_r_r_w_w_; + math::transform::SpaceTransform2D::execute( + M_r_r_w_w_transformed, G_); return flops; } @@ -385,6 +402,45 @@ auto TpAccumulator::getGSingleband(const int s, con return G_(0, 0, s, k1, k2, w1_ext, w2_ext); } +template +auto TpAccumulator::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 +void TpAccumulator::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 +void TpAccumulator::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 void TpAccumulator::getGMultiband(int s, int k1, int k2, int w1, int w2, Matrix& G, @@ -395,7 +451,7 @@ void TpAccumulator::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 << "," @@ -419,6 +475,7 @@ double TpAccumulator::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; }; @@ -636,15 +693,11 @@ double TpAccumulator::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) @@ -678,8 +731,10 @@ double TpAccumulator::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) @@ -689,8 +744,12 @@ double TpAccumulator::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) diff --git a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_base.hpp b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_base.hpp index d71728c64..5c25322fb 100644 --- a/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_base.hpp +++ b/include/dca/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_base.hpp @@ -165,10 +165,38 @@ void TpAccumulatorGpuBase::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 q_minus_k_phase_real(q_indices.size() * KDmn::dmn_size() * n_bands_); + std::vector 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(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); } diff --git a/src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cu b/src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cu index ce25f9368..48113c6c8 100644 --- a/src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cu +++ b/src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/g4_helper.cu @@ -14,11 +14,9 @@ #include #include +#include #include #include -#include - -#include "dca/platform/dca_gpu.h" namespace dca { namespace phys { @@ -33,7 +31,9 @@ __CONSTANT__ int* k_ex_indices_actual; void G4Helper::set(int nb, int nk, int nw, const std::vector& delta_k, const std::vector& 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& q_minus_k_phase_real, + const std::vector& q_minus_k_phase_imag) { // Initialize the reciprocal cluster if not done already. solver::details::ClusterHelper::setMomentum(nk, add_k, lda, sub_k, lds); @@ -72,6 +72,22 @@ void G4Helper::set(int nb, int nk, int nw, const std::vector& delta_k, checkRC(cudaMemcpy(const_cast(host_helper.k_ex_indices_), delta_k.data(), sizeof(int) * delta_k.size(), cudaMemcpyHostToDevice)); + const std::size_t phase_table_size = + static_cast(nb) * static_cast(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(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(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 diff --git a/src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_kernels.cu b/src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_kernels.cu index 03b8d44a8..e9f47b929 100644 --- a/src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_kernels.cu +++ b/src/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_kernels.cu @@ -45,6 +45,14 @@ using dca::util::RealAlias; using phys::FourPointType; using dca::util::SignType; +template +__device__ inline GPUComplex> makeQMinusKPhase(const int k, const int k_ex, + const int band) { + return GPUComplex>{ + static_cast>(g4_helper.qMinusKPhaseReal(k, k_ex, band)), + static_cast>(g4_helper.qMinusKPhaseImag(k, k_ex, band))}; +} + std::string toString(const std::array& dims) { std::ostringstream oss; oss << "{{" << static_cast((dims[0]).x) << "," << (dims[0]).y << "},{" << dims[1].x << "," @@ -167,24 +175,26 @@ __global__ void computeGMultibandKernel(GPUComplex* __restrict__ G, int ld __syncthreads(); GPUComplex G_val_store = G[id_i + ldg * id_j]; - const GPUComplex* const G0_w1 = G0 + nb * k2 + no * w2; - const GPUComplex* const G0_w2 = G0 + nb * k1 + no * w1; + const GPUComplex* const G0_w1 = G0 + nb * k1 + no * w1; + const GPUComplex* const G0_w2 = G0 + nb * k2 + no * w2; G_val_store.x = 0; G_val_store.y = 0; for (int j = 0; j < nb; ++j) { for (int i = 0; i < nb; ++i) { - const GPUComplex G_band = -G0_w1[i + ldg0 * b1] * M[j + ldm * i] * G0_w2[b2 + ldg0 * j]; + const GPUComplex G_band = + -G0_w1[b1 + ldg0 * i] * M[i + ldm * j] * G0_w2[j + ldg0 * b2]; G_val_store += G_band; } } if (k1 == k2 && w1 == w2) // G0_w1 == G0_w2) - G_val_store += G0_w1[b2 + ldg0 * b1] * beta; + G_val_store += G0_w1[b1 + ldg0 * b2] * beta; #ifdef DEBUG_G4_GPU printf("%lf %lf %lf %lf %lf %lf -- %d %d %d %d %d %d %f,%f\n", M[b1 + ldm * b2].x, - M[b1 + ldm * b2].y, G0_w1[b2 + ldg0 * b1].x, G0_w1[b2 + ldg0 * b1].y, - G0_w2[b1 + ldg0 * b2].x, G0_w2[b1 + ldg0 * b2].y, b1, b2, k1, k2, w1, w2, G_val.x, G_val.y); + M[b1 + ldm * b2].y, G0_w1[b1 + ldg0 * b2].x, G0_w1[b1 + ldg0 * b2].y, + G0_w2[b2 + ldg0 * b1].x, G0_w2[b2 + ldg0 * b1].y, b1, b2, k1, k2, w1, w2, G_val.x, + G_val.y); #endif G_val = G_val_store; } @@ -650,14 +660,13 @@ __global__ void updateG4Kernel(GPUComplex>* __restrict__ G4, g4_helper.extendGIndicesMultiBand(k1_a, k2_a, w1_a, w2_a); int i_a = nb * k1_a + no * w1_a + b1; int j_a = nb * k2_a + no * w2_a + b3; - const GPUComplex> Ga_1 = G_up[i_a + ldgu * j_a]; const GPUComplex> Ga_2 = G_down[i_a + ldgd * j_a]; int w1_b(g4_helper.wexMinus(w1, w_ex)); int w2_b(g4_helper.wexMinus(w2, w_ex)); - int k1_b = g4_helper.kexMinus(k1, k_ex); - int k2_b = g4_helper.kexMinus(k2, k_ex); + int k1_b(g4_helper.kexMinus(k1, k_ex)); + int k2_b(g4_helper.kexMinus(k2, k_ex)); if (g4_helper.get_bands() == 1) g4_helper.extendGIndices(k1_b, k2_b, w1_b, w2_b); @@ -670,7 +679,9 @@ __global__ void updateG4Kernel(GPUComplex>* __restrict__ G4, const GPUComplex> Gb_1 = G_down[i_b + ldgd * j_b]; const GPUComplex> Gb_2 = G_up[i_b + ldgu * j_b]; - contribution = sign_over_2 * (Ga_1 * Gb_1 + Ga_2 * Gb_2); + const auto phase = + makeQMinusKPhase(k1, k_ex, b2) * conj(makeQMinusKPhase(k2, k_ex, b4)); + contribution = sign_over_2 * (Ga_1 * Gb_1 + Ga_2 * Gb_2) * phase; } decltype(G4) const result_ptr = G4 + local_g4_index; @@ -775,6 +786,8 @@ __global__ void updateG4KernelNoSpin(GPUComplex>* __restrict__ // the exchange momentum, implies the same operation is performed with the exchange frequency. // See tp_accumulator.hpp for more details. if constexpr (type == FourPointType::PARTICLE_PARTICLE_UP_DOWN) { + const auto phase = + makeQMinusKPhase(k1, k_ex, b2) * conj(makeQMinusKPhase(k2, k_ex, b4)); { int w1_a(w1); int w2_a(w2); @@ -797,7 +810,7 @@ __global__ void updateG4KernelNoSpin(GPUComplex>* __restrict__ const GPUComplex> Ga_1 = G_dn[i_a + ldgd * j_a]; const GPUComplex> Gb_1 = G_dn[i_b + ldgd * j_b]; - contribution = complex_factor * (Ga_1 * Gb_1); + contribution = complex_factor * (Ga_1 * Gb_1) * phase; } { int w1_a(w1); @@ -821,7 +834,7 @@ __global__ void updateG4KernelNoSpin(GPUComplex>* __restrict__ const GPUComplex> Ga_1 = G_dn[i_a + ldgd * j_a]; const GPUComplex> Gb_1 = G_dn[i_b + ldgd * j_b]; - contribution -= complex_factor * (Ga_1 * Gb_1); + contribution -= complex_factor * (Ga_1 * Gb_1) * phase; } } decltype(G4) const result_ptr = G4 + local_g4_index; diff --git a/test/integration/cluster_solver/stdthread_qmci/gpu/stdthread_ctaux_gpu_tp_test.cpp b/test/integration/cluster_solver/stdthread_qmci/gpu/stdthread_ctaux_gpu_tp_test.cpp index 9d78666f3..30975fcf2 100644 --- a/test/integration/cluster_solver/stdthread_qmci/gpu/stdthread_ctaux_gpu_tp_test.cpp +++ b/test/integration/cluster_solver/stdthread_qmci/gpu/stdthread_ctaux_gpu_tp_test.cpp @@ -11,7 +11,10 @@ // nearest-neighbour hopping and on site interaction. The results are expected to be the // same up to numerical error. +#include +#include #include +#include #include #include "dca/config/cmake_options.hpp" @@ -69,6 +72,43 @@ using QmcSolverGpu = dca::phys::solver::StdThreadQmciClusterSolver; using QmcSolverCpu = dca::phys::solver::StdThreadQmciClusterSolver; +template +dca::func::util::Difference differenceForTransfer(const Function& cpu, const Function& gpu, + const int q_index, const int w_index) { + double l1 = 0.; + double l2 = 0.; + double linf = 0.; + + double l1_error = 0.; + double l2_error = 0.; + double linf_error = 0.; + + for (int i = 0; i < cpu.size(); ++i) { + const auto subind = cpu.linind_2_subind(i); + if (subind[8] != std::size_t(q_index) || subind[9] != std::size_t(w_index)) + continue; + + const double ref = std::abs(cpu(i)); + l1 += ref; + l2 += ref * ref; + linf = std::max(linf, ref); + + const double err = std::abs(cpu(i) - gpu(i)); + l1_error += err; + l2_error += err * err; + linf_error = std::max(linf_error, err); + } + + const auto relative = [](const double numerator, const double denominator) { + if (denominator > 0) + return numerator / denominator; + return numerator == 0 ? 0. : std::numeric_limits::infinity(); + }; + + return {relative(l1_error, l1), relative(std::sqrt(l2_error), std::sqrt(l2)), + relative(linf_error, linf)}; +} + TEST(PosixCtauxClusterSolverTest, G_k_w) { dca::linalg::util::initializeMagma(); TestConcurrency concurrency(0, nullptr); @@ -78,7 +118,11 @@ TEST(PosixCtauxClusterSolverTest, G_k_w) { } Parameters parameters(dca::util::GitVersion::string(), concurrency); - parameters.read_input_and_broadcast(input_dir + "threaded_input.json"); + const char* input_override = std::getenv("DCA_CPU_GPU_TEST_INPUT"); + const std::string input_file = + input_override ? input_override : input_dir + "threaded_input.json"; + std::cout << "Reading input " << input_file << "\n"; + parameters.read_input_and_broadcast(input_file); parameters.update_model(); parameters.update_domains(); @@ -105,10 +149,45 @@ TEST(PosixCtauxClusterSolverTest, G_k_w) { QmcSolverGpu qmc_solver_gpu(parameters, data_gpu, nullptr); perform_integration(qmc_solver_gpu); + // dca::func::util::difference reports relative CPU/GPU norm errors: + // err_i = abs(cpu_i - gpu_i), ref_i = abs(cpu_i) + // l1 = sum_i err_i / sum_i ref_i + // l2 = sqrt(sum_i err_i^2 / sum_i ref_i^2) + // l_inf = max_i err_i / max_i ref_i const auto err_g = dca::func::util::difference(data_cpu.G_k_w, data_gpu.G_k_w); - const auto err_g4 = dca::func::util::difference(data_cpu.get_G4()[0], data_gpu.get_G4()[0]); + std::cout << "CPU/GPU G_k_w relative differences:" + << " l1=" << err_g.l1 << " l2=" << err_g.l2 << " l_inf=" << err_g.l_inf + << "\n"; EXPECT_GE(5e-7, err_g.l_inf); - // Is this too large? - EXPECT_GE(5e-5, err_g4.l_inf); + + const auto& g4_cpu = data_cpu.get_G4(); + const auto& g4_gpu = data_gpu.get_G4(); + constexpr double g4_tolerance = 5e-5; + ASSERT_EQ(g4_cpu.size(), g4_gpu.size()); + for (std::size_t channel = 0; channel < g4_cpu.size(); ++channel) { + const auto err_g4 = dca::func::util::difference(g4_cpu[channel], g4_gpu[channel]); + std::cout << "CPU/GPU G4 channel " << channel << " all-transfer relative differences:" + << " l1=" << err_g4.l1 << " l2=" << err_g4.l2 << " l_inf=" << err_g4.l_inf + << "\n"; + EXPECT_GE(g4_tolerance, err_g4.l_inf) << "G4 channel: " << channel; + + const auto& sizes = g4_cpu[channel].getDomainSizes(); + ASSERT_EQ(10, sizes.size()); + ASSERT_EQ(sizes, g4_gpu[channel].getDomainSizes()); + + for (int w_index = 0; w_index < sizes[9]; ++w_index) { + for (int q_index = 0; q_index < sizes[8]; ++q_index) { + const auto err_g4_transfer = + differenceForTransfer(g4_cpu[channel], g4_gpu[channel], q_index, w_index); + if (err_g4_transfer.l_inf > g4_tolerance) { + ADD_FAILURE() << "G4 channel: " << channel << ", q_index: " << q_index + << ", w_index: " << w_index + << " relative differences: l1=" << err_g4_transfer.l1 + << " l2=" << err_g4_transfer.l2 + << " l_inf=" << err_g4_transfer.l_inf; + } + } + } + } } diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test_complex_baseline.hdf5 b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test_complex_baseline.hdf5 index 99e7f9a7d..e08bfc931 100644 Binary files a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test_complex_baseline.hdf5 and b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test_complex_baseline.hdf5 differ