From 244c003a1bd5316e80cfdac38959d7ff7528bb0b Mon Sep 17 00:00:00 2001 From: Raoul Hodgson Date: Fri, 12 Jun 2020 11:00:25 +0100 Subject: [PATCH 1/8] Updated Baryon code --- Grid/qcd/utils/BaryonUtils.h | 40 +++++++++++++++--------------------- 1 file changed, 17 insertions(+), 23 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 23267270..3384d273 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -7,6 +7,7 @@ Copyright (C) 2019 Author: Felix Erben + Author: Raoul Hodgson This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -58,7 +59,7 @@ public: const Gamma GammaA_right, const Gamma GammaB_right, const int parity, - const int * wick_contractions, + const bool * wick_contractions, robj &result); public: static void ContractBaryons(const PropagatorField &q1_left, @@ -68,8 +69,7 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const char * quarks_left, - const char * quarks_right, + const bool* wick_contractions, const int parity, ComplexField &baryon_corr); template @@ -80,9 +80,9 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const char * quarks_left, - const char * quarks_right, + const int wick_contraction, const int parity, + const int nt, robj &result); private: template @@ -173,7 +173,7 @@ void BaryonUtils::baryon_site(const mobj &D1, const Gamma GammaA_right, const Gamma GammaB_right, const int parity, - const int * wick_contraction, + const bool * wick_contraction, robj &result) { @@ -279,8 +279,7 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const char * quarks_left, - const char * quarks_right, + const bool* wick_contractions, const int parity, ComplexField &baryon_corr) { @@ -288,7 +287,6 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); - std::cout << "Contraction <" << quarks_right[0] << quarks_right[1] << quarks_right[2] << "|" << quarks_left[0] << quarks_left[1] << quarks_left[2] << ">" << std::endl; std::cout << "GammaA (left) " << (GammaA_left.g) << std::endl; std::cout << "GammaB (left) " << (GammaB_left.g) << std::endl; std::cout << "GammaA (right) " << (GammaA_right.g) << std::endl; @@ -298,10 +296,6 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, GridBase *grid = q1_left.Grid(); - int wick_contraction[6]; - for (int ie=0; ie < 6 ; ie++) - wick_contraction[ie] = (quarks_left[0] == quarks_right[epsilon[ie][0]] && quarks_left[1] == quarks_right[epsilon[ie][1]] && quarks_left[2] == quarks_right[epsilon[ie][2]]) ? 1 : 0; - auto vbaryon_corr= baryon_corr.View(); auto v1 = q1_left.View(); auto v2 = q2_left.View(); @@ -311,10 +305,10 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, bytes += grid->oSites() * (432.*sizeof(vComplex) + 126.*sizeof(int) + 36.*sizeof(Real)); for (int ie=0; ie < 6 ; ie++){ if(ie==0 or ie==3){ - bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contraction[ie]; + bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie]; } else{ - bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contraction[ie]; + bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie]; } } Real t=0.; @@ -325,7 +319,7 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, auto D2 = v2[ss]; auto D3 = v3[ss]; vobj result=Zero(); - baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result); + baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contractions,result); vbaryon_corr[ss] = result; } );//end loop over lattice sites @@ -343,16 +337,15 @@ void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const char * quarks_left, - const char * quarks_right, + const int wick_contraction, const int parity, + const int nt, robj &result) { assert(Ns==4 && "Baryon code only implemented for N_spin = 4"); assert(Nc==3 && "Baryon code only implemented for N_colour = 3"); - std::cout << "Contraction <" << quarks_right[0] << quarks_right[1] << quarks_right[2] << "|" << quarks_left[0] << quarks_left[1] << quarks_left[2] << ">" << std::endl; std::cout << "GammaA (left) " << (GammaA_left.g) << std::endl; std::cout << "GammaB (left) " << (GammaB_left.g) << std::endl; std::cout << "GammaA (right) " << (GammaA_right.g) << std::endl; @@ -360,12 +353,13 @@ void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); - int wick_contraction[6]; + bool wick_contractions[6]; for (int ie=0; ie < 6 ; ie++) - wick_contraction[ie] = (quarks_left[0] == quarks_right[epsilon[ie][0]] && quarks_left[1] == quarks_right[epsilon[ie][1]] && quarks_left[2] == quarks_right[epsilon[ie][2]]) ? 1 : 0; + wick_contractions[ie] = (ie == wick_contraction); - result=Zero(); - baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result); + for (int t=0; t Date: Fri, 12 Jun 2020 11:35:52 +0100 Subject: [PATCH 2/8] Added Baryon3pt code --- Grid/qcd/utils/BaryonUtils.h | 408 +++++++++++++++++++++++++++++++++++ 1 file changed, 408 insertions(+) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 23267270..a6f8b78d 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -84,6 +84,55 @@ public: const char * quarks_right, const int parity, robj &result); + private: + template + static void Baryon_Gamma_3pt_Group1_Site( + const mobj &Dq1_ti, + const mobj2 &Dq2_spec, + const mobj2 &Dq3_spec, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result); + + template + static void Baryon_Gamma_3pt_Group2_Site( + const mobj2 &Dq1_spec, + const mobj &Dq2_ti, + const mobj2 &Dq3_spec, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result); + + template + static void Baryon_Gamma_3pt_Group3_Site( + const mobj2 &Dq1_spec, + const mobj2 &Dq2_spec, + const mobj &Dq3_ti, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result); + public: + template + static void Baryon_Gamma_3pt( + const PropagatorField &q_ti, + const mobj &Dq_spec1, + const mobj &Dq_spec2, + const PropagatorField &q_tf, + int group, + int wick_contraction, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + SpinMatrixField &stn_corr); private: template static void Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop, @@ -371,6 +420,365 @@ void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, /*********************************************************************** * End of Baryon 2pt-function code. * * * + * The following code is for baryonGamma3pt function * + **********************************************************************/ + +/* Dq1_ti is a quark line from t_i to t_J + * Dq2_spec is a quark line from t_i to t_f + * Dq3_spec is a quark line from t_i to t_f + * Dq4_tf is a quark line from t_f to t_J */ +template +template +void BaryonUtils::Baryon_Gamma_3pt_Group1_Site( + const mobj &Dq1_ti, + const mobj2 &Dq2_spec, + const mobj2 &Dq3_spec, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result) +{ + Gamma g5(Gamma::Algebra::Gamma5); + + auto adjD4_g_D1 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq1_ti; + auto Gf_adjD4_g_D1 = GammaBf * adjD4_g_D1; + auto D2_Gi = Dq2_spec * GammaBi; + auto Gf_D2_Gi = GammaBf * D2_Gi; + auto Gf_D3 = GammaBf * Dq3_spec; + + int a_f, b_f, c_f; + int a_i, b_i, c_i; + + Complex ee; + + auto D2_Gi_ab_aa = D2_Gi ()(0,0)(0,0); + auto Gf_D3_ab_bb = Gf_D3 ()(0,0)(0,0); + auto Gf_D2_Gi_ab_ba = Gf_D2_Gi ()(0,0)(0,0); + auto ee_adjD4_g_D1_ag_ac = ee * adjD4_g_D1 ()(0,0)(0,0); + auto ee_Gf_adjD4_g_D1_ag_bc = ee * Gf_adjD4_g_D1()(0,0)(0,0); + auto Dq3_spec_ab_ab = Dq3_spec ()(0,0)(0,0); + auto ee_adjD4_g_D1_gg_cc = ee * adjD4_g_D1 ()(0,0)(0,0); + auto Dq3_spec_gb_cb = Dq3_spec ()(0,0)(0,0); + auto D2_Gi_gb_ca = D2_Gi ()(0,0)(0,0); + + for (int ie_f=0; ie_f < 6 ; ie_f++){ + a_f = epsilon[ie_f][0]; //a + b_f = epsilon[ie_f][1]; //b + c_f = epsilon[ie_f][2]; //c + for (int ie_i=0; ie_i < 6 ; ie_i++){ + a_i = epsilon[ie_i][0]; //a' + b_i = epsilon[ie_i][1]; //b' + c_i = epsilon[ie_i][2]; //c' + + ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + + for (int alpha_f=0; alpha_f +template +void BaryonUtils::Baryon_Gamma_3pt_Group2_Site( + const mobj2 &Dq1_spec, + const mobj &Dq2_ti, + const mobj2 &Dq3_spec, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result) +{ + Gamma g5(Gamma::Algebra::Gamma5); + + auto adjD4_g_D2_Gi = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi; + auto Gf_adjD4_g_D2_Gi = GammaBf * adjD4_g_D2_Gi; + auto Gf_D1 = GammaBf * Dq1_spec; + auto Gf_D3 = GammaBf * Dq3_spec; + + int a_f, b_f, c_f; + int a_i, b_i, c_i; + + Complex ee; + + auto adjD4_g_D2_Gi_ab_aa = adjD4_g_D2_Gi ()(0,0)(0,0); + auto Gf_D3_ab_bb = Gf_D3 ()(0,0)(0,0); + auto Gf_adjD4_g_D2_Gi_ab_ba = Gf_adjD4_g_D2_Gi ()(0,0)(0,0); + auto Dq3_spec_ab_ab = Dq3_spec ()(0,0)(0,0); + auto ee_Dq1_spec_ag_ac = ee* Dq1_spec ()(0,0)(0,0); + auto ee_Gf_D1_ag_bc = ee*Gf_D1 ()(0,0)(0,0); + auto ee_Dq1_spec_gg_cc = ee*Dq1_spec ()(0,0)(0,0); + auto Dq3_spec_gb_cb = Dq3_spec ()(0,0)(0,0); + auto adjD4_g_D2_Gi_gb_ca = adjD4_g_D2_Gi ()(0,0)(0,0); + + for (int ie_f=0; ie_f < 6 ; ie_f++){ + a_f = epsilon[ie_f][0]; //a + b_f = epsilon[ie_f][1]; //b + c_f = epsilon[ie_f][2]; //c + for (int ie_i=0; ie_i < 6 ; ie_i++){ + a_i = epsilon[ie_i][0]; //a' + b_i = epsilon[ie_i][1]; //b' + c_i = epsilon[ie_i][2]; //c' + + ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + + for (int alpha_f=0; alpha_f +template +void BaryonUtils::Baryon_Gamma_3pt_Group3_Site( + const mobj2 &Dq1_spec, + const mobj2 &Dq2_spec, + const mobj &Dq3_ti, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result) +{ + Gamma g5(Gamma::Algebra::Gamma5); + + auto adjD4_g_D3 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti; + auto Gf_adjD4_g_D3 = GammaBf * adjD4_g_D3; + auto Gf_D1 = GammaBf * Dq1_spec; + auto D2_Gi = Dq2_spec * GammaBi; + auto Gf_D2_Gi = GammaBf * D2_Gi; + + int a_f, b_f, c_f; + int a_i, b_i, c_i; + + Complex ee; + + auto D2_Gi_ab_aa = D2_Gi ()(0,0)(0,0); + auto Gf_adjD4_g_D3_ab_bb = Gf_adjD4_g_D3 ()(0,0)(0,0); + auto Gf_D2_Gi_ab_ba = Gf_D2_Gi ()(0,0)(0,0); + auto adjD4_g_D3_ab_ab = adjD4_g_D3 ()(0,0)(0,0); + auto ee_Dq1_spec_ag_ac = ee * Dq1_spec ()(0,0)(0,0); + auto ee_Gf_D1_ag_bc = ee * Gf_D1 ()(0,0)(0,0); + auto ee_Dq1_spec_gg_cc = ee * Dq1_spec ()(0,0)(0,0); + auto adjD4_g_D3_gb_cb = adjD4_g_D3 ()(0,0)(0,0); + auto D2_Gi_gb_ca = D2_Gi ()(0,0)(0,0); + + for (int ie_f=0; ie_f < 6 ; ie_f++){ + a_f = epsilon[ie_f][0]; //a + b_f = epsilon[ie_f][1]; //b + c_f = epsilon[ie_f][2]; //c + for (int ie_i=0; ie_i < 6 ; ie_i++){ + a_i = epsilon[ie_i][0]; //a' + b_i = epsilon[ie_i][1]; //b' + c_i = epsilon[ie_i][2]; //c' + + ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + + for (int alpha_f=0; alpha_f +template +void BaryonUtils::Baryon_Gamma_3pt( + const PropagatorField &q_ti, + const mobj &Dq_spec1, + const mobj &Dq_spec2, + const PropagatorField &q_tf, + int group, + int wick_contraction, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + SpinMatrixField &stn_corr) +{ + GridBase *grid = q_tf.Grid(); + + auto vcorr= stn_corr.View(); + auto vq_ti = q_ti.View(); + auto vq_tf = q_tf.View(); + + if (group ==1) { + thread_for(ss,grid->oSites(),{ + auto Dq_ti = vq_ti[ss]; + auto Dq_tf = vq_tf[ss]; + sobj result=Zero(); + Baryon_Gamma_3pt_Group1_Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + vcorr[ss] += result; + });//end loop over lattice sites + } else if (group == 2) { + thread_for(ss,grid->oSites(),{ + auto Dq_ti = vq_ti[ss]; + auto Dq_tf = vq_tf[ss]; + sobj result=Zero(); + Baryon_Gamma_3pt_Group2_Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + vcorr[ss] += result; + });//end loop over lattice sites + } else if (group == 3) { + thread_for(ss,grid->oSites(),{ + auto Dq_ti = vq_ti[ss]; + auto Dq_tf = vq_tf[ss]; + sobj result=Zero(); + Baryon_Gamma_3pt_Group3_Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + + vcorr[ss] += result; + });//end loop over lattice sites + } + +} + + +/*********************************************************************** + * End of BaryonGamma3pt-function code. * + * * * The following code is for Sigma -> N rare hypeon decays * **********************************************************************/ From 599f28f6efa5deae9a3fa639ea9b00e4ceba44e1 Mon Sep 17 00:00:00 2001 From: Raoul Hodgson Date: Tue, 23 Jun 2020 11:10:26 +0100 Subject: [PATCH 3/8] Baryon bug fixes --- Grid/qcd/utils/BaryonUtils.h | 168 +++++++++++++++++------------------ 1 file changed, 84 insertions(+), 84 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 3384d273..297f6a2e 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -168,107 +168,107 @@ template void BaryonUtils::baryon_site(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA_left, - const Gamma GammaB_left, - const Gamma GammaA_right, - const Gamma GammaB_right, + const Gamma GammaA_i, + const Gamma GammaB_i, + const Gamma GammaA_f, + const Gamma GammaB_f, const int parity, const bool * wick_contraction, robj &result) { Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) - auto gD1a = GammaA_left * GammaA_right * D1; - auto gD1b = GammaA_left * g4 * GammaA_right * D1; - auto pD1 = 0.5* (gD1a + (Real)parity * gD1b); - auto gD3 = GammaB_right * D3; - auto D2g = D2 * GammaB_left; - auto pD1g = pD1 * GammaB_left; - auto gD3g = gD3 * GammaB_left; + + auto D1_GAi = D1 * GammaA_i; + auto D1_GAi_g4 = D1_GAi * g4; + auto D1_GAi_P = 0.5*(D1_GAi + (Real)parity * D1_GAi_g4); + auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P; + auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P; - for (int ie_left=0; ie_left < 6 ; ie_left++){ - int a_left = epsilon[ie_left][0]; //a - int b_left = epsilon[ie_left][1]; //b - int c_left = epsilon[ie_left][2]; //c - for (int ie_right=0; ie_right < 6 ; ie_right++){ - int a_right = epsilon[ie_right][0]; //a' - int b_right = epsilon[ie_right][1]; //b' - int c_right = epsilon[ie_right][2]; //c' - Real ee = epsilon_sgn[ie_left] * epsilon_sgn[ie_right]; + auto D2_GBi = D2 * GammaB_i; + auto GBf_D2_GBi = GammaB_f * D2_GBi; + auto GAf_D2_GBi = GammaA_f * D2_GBi; + + auto GBf_D3 = GammaB_f * D3; + auto GAf_D3 = GammaA_f * D3; + + for (int ie_f=0; ie_f < 6 ; ie_f++){ + int a_f = epsilon[ie_f][0]; //a + int b_f = epsilon[ie_f][1]; //b + int c_f = epsilon[ie_f][2]; //c + for (int ie_i=0; ie_i < 6 ; ie_i++){ + int a_i = epsilon[ie_i][0]; //a' + int b_i = epsilon[ie_i][1]; //b' + int c_i = epsilon[ie_i][2]; //c' + + Real ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; //This is the \delta_{456}^{123} part - if (wick_contraction[0]){ - for (int gamma_left=0; gamma_left From 3e97a26f90b66876fa44fd72fd444efb8ea61cb9 Mon Sep 17 00:00:00 2001 From: Raoul Hodgson Date: Tue, 23 Jun 2020 11:35:32 +0100 Subject: [PATCH 4/8] BaryonGamm3pt threads -> accelerator --- Grid/qcd/utils/BaryonUtils.h | 96 +++++++++++++----------------------- 1 file changed, 33 insertions(+), 63 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index a6f8b78d..d36c3cb8 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -451,17 +451,7 @@ void BaryonUtils::Baryon_Gamma_3pt_Group1_Site( int a_f, b_f, c_f; int a_i, b_i, c_i; - Complex ee; - - auto D2_Gi_ab_aa = D2_Gi ()(0,0)(0,0); - auto Gf_D3_ab_bb = Gf_D3 ()(0,0)(0,0); - auto Gf_D2_Gi_ab_ba = Gf_D2_Gi ()(0,0)(0,0); - auto ee_adjD4_g_D1_ag_ac = ee * adjD4_g_D1 ()(0,0)(0,0); - auto ee_Gf_adjD4_g_D1_ag_bc = ee * Gf_adjD4_g_D1()(0,0)(0,0); - auto Dq3_spec_ab_ab = Dq3_spec ()(0,0)(0,0); - auto ee_adjD4_g_D1_gg_cc = ee * adjD4_g_D1 ()(0,0)(0,0); - auto Dq3_spec_gb_cb = Dq3_spec ()(0,0)(0,0); - auto D2_Gi_gb_ca = D2_Gi ()(0,0)(0,0); + Real ee; for (int ie_f=0; ie_f < 6 ; ie_f++){ a_f = epsilon[ie_f][0]; //a @@ -476,18 +466,18 @@ void BaryonUtils::Baryon_Gamma_3pt_Group1_Site( for (int alpha_f=0; alpha_f::Baryon_Gamma_3pt_Group2_Site( int a_f, b_f, c_f; int a_i, b_i, c_i; - Complex ee; - - auto adjD4_g_D2_Gi_ab_aa = adjD4_g_D2_Gi ()(0,0)(0,0); - auto Gf_D3_ab_bb = Gf_D3 ()(0,0)(0,0); - auto Gf_adjD4_g_D2_Gi_ab_ba = Gf_adjD4_g_D2_Gi ()(0,0)(0,0); - auto Dq3_spec_ab_ab = Dq3_spec ()(0,0)(0,0); - auto ee_Dq1_spec_ag_ac = ee* Dq1_spec ()(0,0)(0,0); - auto ee_Gf_D1_ag_bc = ee*Gf_D1 ()(0,0)(0,0); - auto ee_Dq1_spec_gg_cc = ee*Dq1_spec ()(0,0)(0,0); - auto Dq3_spec_gb_cb = Dq3_spec ()(0,0)(0,0); - auto adjD4_g_D2_Gi_gb_ca = adjD4_g_D2_Gi ()(0,0)(0,0); + Real ee; for (int ie_f=0; ie_f < 6 ; ie_f++){ a_f = epsilon[ie_f][0]; //a @@ -577,18 +557,18 @@ void BaryonUtils::Baryon_Gamma_3pt_Group2_Site( for (int alpha_f=0; alpha_f::Baryon_Gamma_3pt_Group3_Site( int a_f, b_f, c_f; int a_i, b_i, c_i; - Complex ee; - - auto D2_Gi_ab_aa = D2_Gi ()(0,0)(0,0); - auto Gf_adjD4_g_D3_ab_bb = Gf_adjD4_g_D3 ()(0,0)(0,0); - auto Gf_D2_Gi_ab_ba = Gf_D2_Gi ()(0,0)(0,0); - auto adjD4_g_D3_ab_ab = adjD4_g_D3 ()(0,0)(0,0); - auto ee_Dq1_spec_ag_ac = ee * Dq1_spec ()(0,0)(0,0); - auto ee_Gf_D1_ag_bc = ee * Gf_D1 ()(0,0)(0,0); - auto ee_Dq1_spec_gg_cc = ee * Dq1_spec ()(0,0)(0,0); - auto adjD4_g_D3_gb_cb = adjD4_g_D3 ()(0,0)(0,0); - auto D2_Gi_gb_ca = D2_Gi ()(0,0)(0,0); + Real ee; for (int ie_f=0; ie_f < 6 ; ie_f++){ a_f = epsilon[ie_f][0]; //a @@ -678,18 +648,18 @@ void BaryonUtils::Baryon_Gamma_3pt_Group3_Site( for (int alpha_f=0; alpha_f::Baryon_Gamma_3pt( auto vq_tf = q_tf.View(); if (group ==1) { - thread_for(ss,grid->oSites(),{ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti[ss]; auto Dq_tf = vq_tf[ss]; sobj result=Zero(); @@ -755,7 +725,7 @@ void BaryonUtils::Baryon_Gamma_3pt( vcorr[ss] += result; });//end loop over lattice sites } else if (group == 2) { - thread_for(ss,grid->oSites(),{ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti[ss]; auto Dq_tf = vq_tf[ss]; sobj result=Zero(); @@ -763,7 +733,7 @@ void BaryonUtils::Baryon_Gamma_3pt( vcorr[ss] += result; });//end loop over lattice sites } else if (group == 3) { - thread_for(ss,grid->oSites(),{ + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { auto Dq_ti = vq_ti[ss]; auto Dq_tf = vq_tf[ss]; sobj result=Zero(); From 4ef50ba31f2b39f39492c35b57ee9e4c10457f3b Mon Sep 17 00:00:00 2001 From: Raoul Hodgson Date: Tue, 23 Jun 2020 11:44:20 +0100 Subject: [PATCH 5/8] Baryon speedup --- Grid/qcd/utils/BaryonUtils.h | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 297f6a2e..9e9a6957 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -205,11 +205,12 @@ void BaryonUtils::baryon_site(const mobj &D1, //This is the \delta_{456}^{123} part if (wick_contraction[0]){ for (int rho=0; rho::baryon_site(const mobj &D1, if (wick_contraction[1]){ for (int rho=0; rho::baryon_site(const mobj &D1, if (wick_contraction[2]){ for (int rho=0; rho::baryon_site(const mobj &D1, //This is the \delta_{456}^{132} part if (wick_contraction[3]){ for (int rho=0; rho::baryon_site(const mobj &D1, if (wick_contraction[4]){ for (int rho=0; rho::baryon_site(const mobj &D1, if (wick_contraction[5]){ for (int rho=0; rho Date: Thu, 25 Jun 2020 16:41:58 +0100 Subject: [PATCH 6/8] BaryonUtils: update to autoView --- Grid/qcd/utils/BaryonUtils.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index c7a72812..11744d16 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -712,9 +712,9 @@ void BaryonUtils::Baryon_Gamma_3pt( { GridBase *grid = q_tf.Grid(); - auto vcorr= stn_corr.View(); - auto vq_ti = q_ti.View(); - auto vq_tf = q_tf.View(); + autoView( vcorr, stn_corr, CpuWrite); + autoView( vq_ti , q_ti, CpuRead); + autoView( vq_tf , q_tf, CpuRead); if (group ==1) { accelerator_for(ss, grid->oSites(), grid->Nsimd(), { From 77af9a3ddca85a34097a5282c290601b95c1a345 Mon Sep 17 00:00:00 2001 From: Raoul Hodgson Date: Fri, 26 Jun 2020 10:08:42 +0100 Subject: [PATCH 7/8] Baryon revert sign --- Grid/qcd/utils/BaryonUtils.h | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 11744d16..df16db45 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -257,7 +257,7 @@ void BaryonUtils::baryon_site(const mobj &D1, auto GAf_D1_GAi_P_rr_cc = GAf_D1_GAi_P()(rho,rho)(c_f,c_i); for (int alpha_f=0; alpha_f::baryon_site(const mobj &D1, for (int alpha_f=0; alpha_f::baryon_site(const mobj &D1, for (int alpha_f=0; alpha_f::baryon_site(const mobj &D1, auto GAf_D1_GAi_P_rr_cc = GAf_D1_GAi_P()(rho,rho)(c_f,c_i); for (int alpha_f=0; alpha_f::baryon_site(const mobj &D1, for (int alpha_f=0; alpha_f::baryon_site(const mobj &D1, for (int alpha_f=0; alpha_f Date: Mon, 29 Jun 2020 09:43:01 +0100 Subject: [PATCH 8/8] Update to baryon and added comments/fix whitespace --- Grid/qcd/utils/BaryonUtils.h | 827 ++++++++++++++++++----------------- 1 file changed, 428 insertions(+), 399 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index df16db45..b268b684 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -62,6 +62,9 @@ public: const bool * wick_contractions, robj &result); public: + static void Wick_Contractions(std::string qi, + std::string qf, + bool* wick_contractions); static void ContractBaryons(const PropagatorField &q1_left, const PropagatorField &q2_left, const PropagatorField &q3_left, @@ -80,59 +83,59 @@ public: const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const int wick_contraction, + const bool* wick_contractions, const int parity, const int nt, robj &result); - private: - template - static void Baryon_Gamma_3pt_Group1_Site( - const mobj &Dq1_ti, - const mobj2 &Dq2_spec, - const mobj2 &Dq3_spec, - const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, - int wick_contraction, - robj &result); + private: + template + static void Baryon_Gamma_3pt_Group1_Site( + const mobj &Dq1_ti, + const mobj2 &Dq2_spec, + const mobj2 &Dq3_spec, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result); - template - static void Baryon_Gamma_3pt_Group2_Site( - const mobj2 &Dq1_spec, - const mobj &Dq2_ti, - const mobj2 &Dq3_spec, - const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, - int wick_contraction, - robj &result); + template + static void Baryon_Gamma_3pt_Group2_Site( + const mobj2 &Dq1_spec, + const mobj &Dq2_ti, + const mobj2 &Dq3_spec, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result); - template - static void Baryon_Gamma_3pt_Group3_Site( - const mobj2 &Dq1_spec, - const mobj2 &Dq2_spec, - const mobj &Dq3_ti, - const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, - int wick_contraction, - robj &result); - public: - template - static void Baryon_Gamma_3pt( - const PropagatorField &q_ti, - const mobj &Dq_spec1, - const mobj &Dq_spec2, - const PropagatorField &q_tf, - int group, - int wick_contraction, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, - SpinMatrixField &stn_corr); + template + static void Baryon_Gamma_3pt_Group3_Site( + const mobj2 &Dq1_spec, + const mobj2 &Dq2_spec, + const mobj &Dq3_ti, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result); + public: + template + static void Baryon_Gamma_3pt( + const PropagatorField &q_ti, + const mobj &Dq_spec1, + const mobj &Dq_spec2, + const PropagatorField &q_tf, + int group, + int wick_contraction, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + SpinMatrixField &stn_corr); private: template static void Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop, @@ -215,15 +218,15 @@ const Real BaryonUtils::epsilon_sgn[6] = {1.,1.,1.,-1.,-1.,-1.}; template template void BaryonUtils::baryon_site(const mobj &D1, - const mobj &D2, - const mobj &D3, - const Gamma GammaA_i, - const Gamma GammaB_i, - const Gamma GammaA_f, - const Gamma GammaB_f, - const int parity, - const bool * wick_contraction, - robj &result) + const mobj &D2, + const mobj &D3, + const Gamma GammaA_i, + const Gamma GammaB_i, + const Gamma GammaA_f, + const Gamma GammaB_f, + const int parity, + const bool * wick_contraction, + robj &result) { Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) @@ -231,101 +234,121 @@ void BaryonUtils::baryon_site(const mobj &D1, auto D1_GAi = D1 * GammaA_i; auto D1_GAi_g4 = D1_GAi * g4; auto D1_GAi_P = 0.5*(D1_GAi + (Real)parity * D1_GAi_g4); - auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P; - auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P; + auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P; + auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P; - auto D2_GBi = D2 * GammaB_i; - auto GBf_D2_GBi = GammaB_f * D2_GBi; - auto GAf_D2_GBi = GammaA_f * D2_GBi; + auto D2_GBi = D2 * GammaB_i; + auto GBf_D2_GBi = GammaB_f * D2_GBi; + auto GAf_D2_GBi = GammaA_f * D2_GBi; - auto GBf_D3 = GammaB_f * D3; - auto GAf_D3 = GammaA_f * D3; + auto GBf_D3 = GammaB_f * D3; + auto GAf_D3 = GammaA_f * D3; for (int ie_f=0; ie_f < 6 ; ie_f++){ - int a_f = epsilon[ie_f][0]; //a - int b_f = epsilon[ie_f][1]; //b - int c_f = epsilon[ie_f][2]; //c + int a_f = epsilon[ie_f][0]; //a + int b_f = epsilon[ie_f][1]; //b + int c_f = epsilon[ie_f][2]; //c for (int ie_i=0; ie_i < 6 ; ie_i++){ - int a_i = epsilon[ie_i][0]; //a' - int b_i = epsilon[ie_i][1]; //b' - int c_i = epsilon[ie_i][2]; //c' + int a_i = epsilon[ie_i][0]; //a' + int b_i = epsilon[ie_i][1]; //b' + int c_i = epsilon[ie_i][2]; //c' - Real ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + Real ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; //This is the \delta_{456}^{123} part - if (wick_contraction[0]){ - for (int rho=0; rho +void BaryonUtils::Wick_Contractions(std::string qi, std::string qf, bool* wick_contractions) { + const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + for (int ie=0; ie < 6 ; ie++) { + wick_contractions[ie] = (qi.size() == 3 && qf.size() == 3 + && qi[0] == qf[epsilon[ie][0]] + && qi[1] == qf[epsilon[ie][1]] + && qi[2] == qf[epsilon[ie][2]]); + } +} + +/* The array wick_contractions must be of length 6. The order * + * corresponds to the to that shown in the Hadrons documentation * + * at https://aportelli.github.io/Hadrons-doc/#/mcontraction * + * This can be computed from the quark flavours using the * + * Wick_Contractions function above */ template void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, const PropagatorField &q2_left, @@ -383,6 +406,12 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, std::cout << std::setw(10) << bytes/t*1.0e6/1024/1024/1024 << " GB/s " << std::endl; } + +/* The array wick_contractions must be of length 6. The order * + * corresponds to the to that shown in the Hadrons documentation * + * at https://aportelli.github.io/Hadrons-doc/#/mcontraction * + * This can also be computed from the quark flavours using the * + * Wick_Contractions function above */ template template void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, @@ -392,7 +421,7 @@ void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, const Gamma GammaB_left, const Gamma GammaA_right, const Gamma GammaB_right, - const int wick_contraction, + const bool* wick_contractions, const int parity, const int nt, robj &result) @@ -408,10 +437,6 @@ void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); - bool wick_contractions[6]; - for (int ie=0; ie < 6 ; ie++) - wick_contractions[ie] = (ie == wick_contraction); - for (int t=0; t::ContractBaryons_Sliced(const mobj &D1, template template void BaryonUtils::Baryon_Gamma_3pt_Group1_Site( - const mobj &Dq1_ti, - const mobj2 &Dq2_spec, - const mobj2 &Dq3_spec, - const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, - int wick_contraction, - robj &result) + const mobj &Dq1_ti, + const mobj2 &Dq2_spec, + const mobj2 &Dq3_spec, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); + Gamma g5(Gamma::Algebra::Gamma5); - auto adjD4_g_D1 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq1_ti; - auto Gf_adjD4_g_D1 = GammaBf * adjD4_g_D1; - auto D2_Gi = Dq2_spec * GammaBi; - auto Gf_D2_Gi = GammaBf * D2_Gi; - auto Gf_D3 = GammaBf * Dq3_spec; + auto adjD4_g_D1 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq1_ti; + auto Gf_adjD4_g_D1 = GammaBf * adjD4_g_D1; + auto D2_Gi = Dq2_spec * GammaBi; + auto Gf_D2_Gi = GammaBf * D2_Gi; + auto Gf_D3 = GammaBf * Dq3_spec; - int a_f, b_f, c_f; - int a_i, b_i, c_i; + int a_f, b_f, c_f; + int a_i, b_i, c_i; - Real ee; + Real ee; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c - for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + for (int ie_f=0; ie_f < 6 ; ie_f++){ + a_f = epsilon[ie_f][0]; //a + b_f = epsilon[ie_f][1]; //b + c_f = epsilon[ie_f][2]; //c + for (int ie_i=0; ie_i < 6 ; ie_i++){ + a_i = epsilon[ie_i][0]; //a' + b_i = epsilon[ie_i][1]; //b' + c_i = epsilon[ie_i][2]; //c' - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; - for (int alpha_f=0; alpha_f::Baryon_Gamma_3pt_Group1_Site( template template void BaryonUtils::Baryon_Gamma_3pt_Group2_Site( - const mobj2 &Dq1_spec, - const mobj &Dq2_ti, - const mobj2 &Dq3_spec, - const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, - int wick_contraction, - robj &result) + const mobj2 &Dq1_spec, + const mobj &Dq2_ti, + const mobj2 &Dq3_spec, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); + Gamma g5(Gamma::Algebra::Gamma5); - auto adjD4_g_D2_Gi = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi; - auto Gf_adjD4_g_D2_Gi = GammaBf * adjD4_g_D2_Gi; - auto Gf_D1 = GammaBf * Dq1_spec; - auto Gf_D3 = GammaBf * Dq3_spec; + auto adjD4_g_D2_Gi = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi; + auto Gf_adjD4_g_D2_Gi = GammaBf * adjD4_g_D2_Gi; + auto Gf_D1 = GammaBf * Dq1_spec; + auto Gf_D3 = GammaBf * Dq3_spec; - int a_f, b_f, c_f; - int a_i, b_i, c_i; + int a_f, b_f, c_f; + int a_i, b_i, c_i; - Real ee; + Real ee; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c - for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + for (int ie_f=0; ie_f < 6 ; ie_f++){ + a_f = epsilon[ie_f][0]; //a + b_f = epsilon[ie_f][1]; //b + c_f = epsilon[ie_f][2]; //c + for (int ie_i=0; ie_i < 6 ; ie_i++){ + a_i = epsilon[ie_i][0]; //a' + b_i = epsilon[ie_i][1]; //b' + c_i = epsilon[ie_i][2]; //c' - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; - for (int alpha_f=0; alpha_f::Baryon_Gamma_3pt_Group2_Site( template template void BaryonUtils::Baryon_Gamma_3pt_Group3_Site( - const mobj2 &Dq1_spec, - const mobj2 &Dq2_spec, - const mobj &Dq3_ti, - const mobj &Dq4_tf, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, - int wick_contraction, - robj &result) + const mobj2 &Dq1_spec, + const mobj2 &Dq2_spec, + const mobj &Dq3_ti, + const mobj &Dq4_tf, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + int wick_contraction, + robj &result) { - Gamma g5(Gamma::Algebra::Gamma5); + Gamma g5(Gamma::Algebra::Gamma5); - auto adjD4_g_D3 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti; - auto Gf_adjD4_g_D3 = GammaBf * adjD4_g_D3; - auto Gf_D1 = GammaBf * Dq1_spec; - auto D2_Gi = Dq2_spec * GammaBi; - auto Gf_D2_Gi = GammaBf * D2_Gi; + auto adjD4_g_D3 = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti; + auto Gf_adjD4_g_D3 = GammaBf * adjD4_g_D3; + auto Gf_D1 = GammaBf * Dq1_spec; + auto D2_Gi = Dq2_spec * GammaBi; + auto Gf_D2_Gi = GammaBf * D2_Gi; - int a_f, b_f, c_f; - int a_i, b_i, c_i; + int a_f, b_f, c_f; + int a_i, b_i, c_i; - Real ee; + Real ee; - for (int ie_f=0; ie_f < 6 ; ie_f++){ - a_f = epsilon[ie_f][0]; //a - b_f = epsilon[ie_f][1]; //b - c_f = epsilon[ie_f][2]; //c - for (int ie_i=0; ie_i < 6 ; ie_i++){ - a_i = epsilon[ie_i][0]; //a' - b_i = epsilon[ie_i][1]; //b' - c_i = epsilon[ie_i][2]; //c' + for (int ie_f=0; ie_f < 6 ; ie_f++){ + a_f = epsilon[ie_f][0]; //a + b_f = epsilon[ie_f][1]; //b + c_f = epsilon[ie_f][2]; //c + for (int ie_i=0; ie_i < 6 ; ie_i++){ + a_i = epsilon[ie_i][0]; //a' + b_i = epsilon[ie_i][1]; //b' + c_i = epsilon[ie_i][2]; //c' - ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; + ee = epsilon_sgn[ie_f] * epsilon_sgn[ie_i]; - for (int alpha_f=0; alpha_f template void BaryonUtils::Baryon_Gamma_3pt( - const PropagatorField &q_ti, - const mobj &Dq_spec1, - const mobj &Dq_spec2, - const PropagatorField &q_tf, - int group, - int wick_contraction, - const Gamma GammaJ, - const Gamma GammaBi, - const Gamma GammaBf, - SpinMatrixField &stn_corr) + const PropagatorField &q_ti, + const mobj &Dq_spec1, + const mobj &Dq_spec2, + const PropagatorField &q_tf, + int group, + int wick_contraction, + const Gamma GammaJ, + const Gamma GammaBi, + const Gamma GammaBf, + SpinMatrixField &stn_corr) { - GridBase *grid = q_tf.Grid(); + GridBase *grid = q_tf.Grid(); - autoView( vcorr, stn_corr, CpuWrite); - autoView( vq_ti , q_ti, CpuRead); - autoView( vq_tf , q_tf, CpuRead); + autoView( vcorr, stn_corr, CpuWrite); + autoView( vq_ti , q_ti, CpuRead); + autoView( vq_tf , q_tf, CpuRead); - if (group ==1) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - Baryon_Gamma_3pt_Group1_Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - vcorr[ss] += result; - });//end loop over lattice sites - } else if (group == 2) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - Baryon_Gamma_3pt_Group2_Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - vcorr[ss] += result; - });//end loop over lattice sites - } else if (group == 3) { - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - auto Dq_ti = vq_ti[ss]; - auto Dq_tf = vq_tf[ss]; - sobj result=Zero(); - Baryon_Gamma_3pt_Group3_Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); - - vcorr[ss] += result; - });//end loop over lattice sites - } + if (group == 1) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti[ss]; + auto Dq_tf = vq_tf[ss]; + sobj result=Zero(); + Baryon_Gamma_3pt_Group1_Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + vcorr[ss] += result; + });//end loop over lattice sites + } else if (group == 2) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti[ss]; + auto Dq_tf = vq_tf[ss]; + sobj result=Zero(); + Baryon_Gamma_3pt_Group2_Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + vcorr[ss] += result; + });//end loop over lattice sites + } else if (group == 3) { + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto Dq_ti = vq_ti[ss]; + auto Dq_tf = vq_tf[ss]; + sobj result=Zero(); + Baryon_Gamma_3pt_Group3_Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); + vcorr[ss] += result; + });//end loop over lattice sites + } }