Skip to content

Commit

Permalink
Fix compilation of G2 measurement
Browse files Browse the repository at this point in the history
  • Loading branch information
@Wentzell
Wentzell committed Sep 2, 2023
1 parent 4fc0e54 commit 67ff398
Show file tree
Hide file tree
Showing 4 changed files with 31 additions and 36 deletions.
12 changes: 6 additions & 6 deletions c++/triqs_cthyb/measures/G2_iw_acc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,12 +86,12 @@ namespace triqs_cthyb {

timer_G2.start();
for (auto &m : G2_measures()) {
auto G2_iw_block = G2_iw(m.b1.index(), m.b2.idx);
bool diag_block = (m.b1.index() == m.b2.idx);
auto G2_iw_block = G2_iw(m.b1.idx, m.b2.idx);
bool diag_block = (m.b1.idx == m.b2.idx);
if (order == block_order::AABB || diag_block)
accumulate_impl_AABB<Channel>(G2_iw_block, s, M(m.b1.index()), M(m.b2.idx));
accumulate_impl_AABB<Channel>(G2_iw_block, s, M(m.b1.idx), M(m.b2.idx));
if (order == block_order::ABBA || diag_block)
accumulate_impl_ABBA<Channel>(G2_iw_block, s, M(m.b1.index()), M(m.b2.idx));
accumulate_impl_ABBA<Channel>(G2_iw_block, s, M(m.b1.idx), M(m.b2.idx));
}
timer_G2.stop();
}
Expand Down Expand Up @@ -170,7 +170,7 @@ namespace triqs_cthyb {
using mesh_point_t = typename std::remove_reference<decltype(iw_mesh)>::type::mesh_point_t;

for (const auto &[n1, n2, n3] : G2.mesh()) {
mesh_point_t n4{iw_mesh, n1.index() + n3.idx - n2.idx};
auto n4 = iw_mesh(n1.index() + n3.index() - n2.index());
for (const auto [i, j, k, l] : G2.target_indices())
G2[n1, n2, n3](i, j, k, l) += s * M_ij[n2, n1](j, i) * M_kl[n4, n3](l, k);
}
Expand All @@ -186,7 +186,7 @@ namespace triqs_cthyb {
using mesh_point_t = typename std::remove_reference<decltype(iw_mesh)>::type::mesh_point_t;

for (const auto &[n1, n2, n3] : G2.mesh()) {
mesh_point_t n4{iw_mesh, n1.index() + n3.idx - n2.idx};
auto n4 = iw_mesh(n1.index() + n3.index() - n2.index());
for (const auto [i, j, k, l] : G2.target_indices())
G2[n1, n2, n3](i, j, k, l) -= s * M_il[n4, n1](l, i) * M_kj[n2, n3](j, k);
}
Expand Down
37 changes: 19 additions & 18 deletions c++/triqs_cthyb/measures/G2_iwll.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,9 +30,9 @@ namespace triqs_cthyb {
const double beta = data.config.beta();

order = G2_measures.params.measure_G2_block_order;
size_t n_l = G2_measures.params.measure_G2_n_l;
long n_l = G2_measures.params.measure_G2_n_l;
int n_bosonic = G2_measures.params.measure_G2_n_bosonic;
int nfft_buf_size = G2_measures.params.measure_G2_iwll_nfft_buf_size;
int nfft_buf_size = G2_measures.params.measure_G2_iwll_nfft_buf_size;

// Allocate the two-particle Green's function
{
Expand All @@ -54,8 +54,8 @@ namespace triqs_cthyb {
for (auto const &m : G2_measures()) {
auto s = m.target_shape;
array<int, 6> buf_sizes(n_l, n_l, s[0], s[1], s[2], s[3]);
buf_sizes() = nfft_buf_size;
nfft_buf(m.b1.index(), m.b2.idx) = nfft_array_t<1, 6>{mesh_w, G2_iwll(m.b1.idx, m.b2.idx).data(), buf_sizes};
buf_sizes() = nfft_buf_size;
nfft_buf(m.b1.idx, m.b2.idx) = nfft_array_t<1, 6>{mesh_w, G2_iwll(m.b1.idx, m.b2.idx).data(), buf_sizes};
}
}
}
Expand All @@ -70,10 +70,9 @@ namespace triqs_cthyb {

for (auto const &m : G2_measures()) {

if (data.dets[m.b1.index()].size() == 0 || data.dets[m.b2.idx].size() == 0) continue;
if (data.dets[m.b1.idx].size() == 0 || data.dets[m.b2.idx].size() == 0) continue;

auto accumulate_impl = [&](op_t const &i, op_t const &j, op_t const &k, op_t const &l, mc_weight_t val) {

tilde_p_gen p_l1_gen(beta), p_l2_gen(beta);
double dtau = setup_times(p_l1_gen, p_l2_gen, i, j, k, l);

Expand All @@ -82,26 +81,26 @@ namespace triqs_cthyb {
for (int l2 : range(n_l)) {
double p_l2 = p_l2_gen.next();
std::array<int, 6> vec{l1, l2, i.second, j.second, k.second, l.second};
nfft_buf(m.b1.index(), m.b2.idx).push_back({dtau}, vec, val * p_l1 * p_l2);
nfft_buf(m.b1.idx, m.b2.idx).push_back({dtau}, vec, val * p_l1 * p_l2);
}
}
};

bool diag_block = (m.b1.index() == m.b2.idx);
bool diag_block = (m.b1.idx == m.b2.idx);

// Perform the accumulation looping over both determinants
if (order == block_order::AABB || diag_block) {
foreach (data.dets[m.b1.index()], [&](op_t const &i, op_t const &j, mc_weight_t M_ij) {
foreach (data.dets[m.b2.index()], [&](op_t const &k, op_t const &l, mc_weight_t M_kl) {
foreach (data.dets[m.b1.idx], [&](op_t const &i, op_t const &j, mc_weight_t M_ij) {
foreach (data.dets[m.b2.idx], [&](op_t const &k, op_t const &l, mc_weight_t M_kl) {
accumulate_impl(i, j, k, l, s * M_ij * M_kl); // Accumulate in legendre-nfft buffer
})
;
})
;
}
if (order == block_order::ABBA || diag_block) {
foreach (data.dets[m.b1.index()], [&](op_t const &i, op_t const &l, mc_weight_t M_il) {
foreach (data.dets[m.b2.index()], [&](op_t const &k, op_t const &j, mc_weight_t M_kj) {
foreach (data.dets[m.b1.idx], [&](op_t const &i, op_t const &l, mc_weight_t M_il) {
foreach (data.dets[m.b2.idx], [&](op_t const &k, op_t const &j, mc_weight_t M_kj) {
accumulate_impl(i, j, k, l, -s * M_il * M_kj); // Accumulate in legendre-nfft buffer
})
;
Expand All @@ -112,15 +111,17 @@ namespace triqs_cthyb {
}

template <>
double measure_G2_iwll<G2_channel::PH>::setup_times(tilde_p_gen &p_l1_gen, tilde_p_gen &p_l2_gen, op_t const &i, op_t const &j, op_t const &k, op_t const &l) {
double measure_G2_iwll<G2_channel::PH>::setup_times(tilde_p_gen &p_l1_gen, tilde_p_gen &p_l2_gen, op_t const &i, op_t const &j, op_t const &k,
op_t const &l) {
p_l1_gen.reset(i.first, j.first);
p_l2_gen.reset(k.first, l.first);
double dtau = 0.5 * double(i.first + j.first - k.first - l.first);
return dtau;
}

template <>
double measure_G2_iwll<G2_channel::PP>::setup_times(tilde_p_gen &p_l1_gen, tilde_p_gen &p_l2_gen, op_t const &i, op_t const &j, op_t const &k, op_t const &l) {
double measure_G2_iwll<G2_channel::PP>::setup_times(tilde_p_gen &p_l1_gen, tilde_p_gen &p_l2_gen, op_t const &i, op_t const &j, op_t const &k,
op_t const &l) {
p_l1_gen.reset(i.first, k.first);
p_l2_gen.reset(j.first, l.first);
double dtau = 0.5 * double(i.first + mult_by_int(j.first, 3) - mult_by_int(k.first, 3) - l.first);
Expand All @@ -129,7 +130,7 @@ namespace triqs_cthyb {

template <G2_channel Channel> void measure_G2_iwll<Channel>::collect_results(mpi::communicator const &c) {

for (auto const &m : G2_measures()) { nfft_buf(m.b1.index(), m.b2.idx).flush(); }
for (auto const &m : G2_measures()) { nfft_buf(m.b1.idx, m.b2.idx).flush(); }

G2_iwll = mpi::all_reduce(G2_iwll, c);

Expand All @@ -139,9 +140,9 @@ namespace triqs_cthyb {

for (auto l : std::get<1>(G2_iwll_block.mesh().components())) {
auto _ = all_t{};
double s = std::sqrt(2 * l + 1);
double s = std::sqrt(2 * l.index() + 1);
G2_iwll_block[_, l, _] *= s;
G2_iwll_block[_, _, l] *= s * (l % 2 ? 1 : -1);
G2_iwll_block[_, _, l] *= s * (l.index() % 2 ? 1 : -1);
}

G2_iwll_block /= (real(average_sign) * data.config.beta());
Expand All @@ -150,4 +151,4 @@ namespace triqs_cthyb {

template struct measure_G2_iwll<G2_channel::PP>;
template struct measure_G2_iwll<G2_channel::PH>;
}
} // namespace triqs_cthyb
8 changes: 4 additions & 4 deletions c++/triqs_cthyb/measures/G2_tau.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,11 +51,11 @@ namespace triqs_cthyb {
// loop only over block-combinations that should be measured
for (auto &m : G2_measures()) {

auto G2_tau_block = G2_tau(m.b1.index(), m.b2.idx);
bool diag_block = (m.b1.index() == m.b2.idx);
auto G2_tau_block = G2_tau(m.b1.idx, m.b2.idx);
bool diag_block = (m.b1.idx == m.b2.idx);

foreach (data.dets[m.b1.index()], [&](auto const &i, auto const &j, auto const M_ij) {
foreach (data.dets[m.b2.index()], [&](auto const &k, auto const &l, auto const M_kl) {
foreach (data.dets[m.b1.idx], [&](auto const &i, auto const &j, auto const M_ij) {
foreach (data.dets[m.b2.idx], [&](auto const &k, auto const &l, auto const M_kl) {

// lambda for computing a single product term of M_ij and M_kl
auto compute_M2_product = [&](auto const &i, auto const &j, auto const &k, auto const &l, mc_weight_t sign) {
Expand Down
10 changes: 2 additions & 8 deletions c++/triqs_cthyb/types.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -99,17 +99,11 @@ namespace triqs {

for (auto const &[bl1, bl1_size] : gf_struct) {
block_names.push_back(bl1);
std::vector<std::string> indices1;
for (auto idx : range(bl1_size)) indices1.push_back(std::to_string(idx));

std::vector<gf<Var_t, tensor_valued<4>>> gf_vec;
for (auto const &[bl2, bl2_size] : gf_struct) {
std::vector<std::string> indices2;
for (auto idx : range(bl2_size)) indices2.push_back(std::to_string(idx));
auto I = std::vector<std::vector<std::string>>{indices1, indices1, indices2, indices2};
switch (order) {
case triqs_cthyb::block_order::AABB: gf_vec.emplace_back(m, make_shape(bl1_size, bl1_size, bl2_size, bl2_size), I); break;
case triqs_cthyb::block_order::ABBA: gf_vec.emplace_back(m, make_shape(bl1_size, bl2_size, bl2_size, bl1_size), I); break;
case triqs_cthyb::block_order::AABB: gf_vec.emplace_back(m, make_shape(bl1_size, bl1_size, bl2_size, bl2_size)); break;
case triqs_cthyb::block_order::ABBA: gf_vec.emplace_back(m, make_shape(bl1_size, bl2_size, bl2_size, bl1_size)); break;
}
}
gf_vecvec.emplace_back(std::move(gf_vec));
Expand Down

0 comments on commit 67ff398

Please sign in to comment.