From 7fff20dec22d2a33010b88f6d017bc4d3a873227 Mon Sep 17 00:00:00 2001 From: Zenny Date: Thu, 16 Jul 2026 14:31:04 +0200 Subject: [PATCH 1/2] made external fermions branch-free in helicity, no discernible performance drop --- .../aloha/template_files/gpu/helas.h | 164 ++++++------------ .../aloha/template_files/gpu/helas_fd.h | 164 ++++++------------ 2 files changed, 102 insertions(+), 226 deletions(-) diff --git a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h index fcfc4b3153..5b7db11740 100644 --- a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h +++ b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h @@ -2,10 +2,10 @@ ! Copyright (C) 2010 The MadGraph5_aMC@NLO development team and contributors. ! Created by: J. Alwall (Sep 2010) for the MG5aMC CPP backend. !========================================================================== -! Copyright (C) 2020-2024 CERN and UCLouvain. +! Copyright (C) 2020-2026 CERN and UCLouvain. ! Licensed under the GNU Lesser General Public License (version 3 or later). ! Modified by: O. Mattelaer (Mar 2020) for the MG5aMC CUDACPP plugin. -! Further modified by: O. Mattelaer, A. Valassi (2020-2024) for the MG5aMC CUDACPP plugin. +! Further modified by: O. Mattelaer, A. Valassi, Z. Wettersten (2020-2026) for the MG5aMC CUDACPP plugin. !========================================================================== //-------------------------------------------------------------------------- @@ -177,6 +177,8 @@ fi[0] = cxmake( -pvec0 * (fptype)nsf, -pvec3 * (fptype)nsf ); fi[1] = cxmake( -pvec1 * (fptype)nsf, -pvec2 * (fptype)nsf ); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3+nh)/2 because omega(2) has indexes 1,2 + const int im = ( 1 - nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2 if( fmass != 0. ) { #ifndef MGONGPU_CPPSIMD @@ -188,8 +190,6 @@ // In C++ ixxxxx, use a single ip/im numbering that is valid both for pp==0 and pp>0, which have two numbering schemes in Fortran ixxxxx: // for pp==0, Fortran sqm(0:1) has indexes 0,1 as in C++; but for Fortran pp>0, omega(2) has indexes 1,2 and not 0,1 // NB: this is only possible in ixxxx, but in oxxxxx two different numbering schemes must be used - const int ip = ( 1 + nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3+nh)/2 because omega(2) has indexes 1,2 - const int im = ( 1 - nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2 #ifndef MGONGPU_CPPSIMD if( pp == 0. ) { @@ -271,20 +271,10 @@ const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), ( sqp0p3 == 0. ? cxmake( -(fptype)nhel * fpsqrt( 2. * pvec0 ), 0. ) : cxmake( (fptype)nh * pvec1, pvec2 ) / sqp0p3 ) }; #endif - if( nh == 1 ) - { - fi[2] = cxzero_sv(); - fi[3] = cxzero_sv(); - fi[4] = chi[0]; - fi[5] = chi[1]; - } - else - { - fi[2] = chi[1]; - fi[3] = chi[0]; - fi[4] = cxzero_sv(); - fi[5] = cxzero_sv(); - } + fi[2] = im * chi[1]; + fi[3] = im * chi[0]; + fi[4] = ip * chi[0]; + fi[5] = ip * chi[1]; } mgDebug( 1, __FUNCTION__ ); return; @@ -309,18 +299,12 @@ fi[0] = cxmake( -pvec3 * (fptype)nsf, -pvec3 * (fptype)nsf ); fi[1] = cxzero_sv(); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; const cxtype_sv sqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); fi[2] = fi[1]; - if( nh == 1 ) - { - fi[3] = fi[1]; - fi[4] = sqp0p3; - } - else - { - fi[3] = sqp0p3; - fi[4] = fi[1]; - } + fi[3] = ip * fi[1]; + fi[4] = im * fi[1]; fi[5] = fi[1]; mgDebug( 1, __FUNCTION__ ); return; @@ -345,19 +329,13 @@ fi[0] = cxmake( pvec3 * (fptype)nsf, -pvec3 * (fptype)nsf ); fi[1] = cxzero_sv(); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; const cxtype_sv chi = cxmake( -(fptype)nhel * fpsqrt( -2. * pvec3 ), 0. ); fi[3] = cxzero_sv(); fi[4] = cxzero_sv(); - if( nh == 1 ) - { - fi[2] = cxzero_sv(); - fi[5] = chi; - } - else - { - fi[2] = chi; - fi[5] = cxzero_sv(); - } + fi[2] = im * chi; + fi[5] = ip * chi; mgDebug( 1, __FUNCTION__ ); return; } @@ -386,24 +364,16 @@ fi[0] = cxmake( -pvec0 * (fptype)nsf, -pvec3 * (fptype)nsf ); // AV: BUG FIX fi[1] = cxmake( -pvec1 * (fptype)nsf, -pvec2 * (fptype)nsf ); // AV: BUG FIX const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); - if( nh == 1 ) - { - fi[2] = cxzero_sv(); - fi[3] = cxzero_sv(); - fi[4] = chi0; - fi[5] = chi1; - } - else - { - fi[2] = chi1; - fi[3] = chi0; - fi[4] = cxzero_sv(); - fi[5] = cxzero_sv(); - } + fi[2] = im * chi1; + fi[3] = im * chi0; + fi[4] = ip * chi0; + fi[5] = ip * chi1; mgDebug( 1, __FUNCTION__ ); return; } @@ -597,6 +567,8 @@ fo[0] = cxmake( pvec0 * (fptype)nsf, pvec3 * (fptype)nsf ); fo[1] = cxmake( pvec1 * (fptype)nsf, pvec2 * (fptype)nsf ); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; // NB: Fortran is (3+nh)/2 because omega(2) has indexes 1,2 and not 0,1 + const int im = ( 1 - nh ) / 2; // NB: Fortran is (3-nh)/2 because omega(2) has indexes 1,2 and not 0,1 if( fmass != 0. ) { #ifndef MGONGPU_CPPSIMD @@ -607,12 +579,12 @@ fptype sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here - const int ip = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - const int im = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - fo[2] = cxmake( im * sqm[std::abs( ip )], 0 ); - fo[3] = cxmake( ip * nsf * sqm[std::abs( ip )], 0 ); - fo[4] = cxmake( im * nsf * sqm[std::abs( im )], 0 ); - fo[5] = cxmake( ip * sqm[std::abs( im )], 0 ); + const int ipp = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + const int imp = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + fo[2] = cxmake( imp * sqm[std::abs( ipp )], 0 ); + fo[3] = cxmake( ipp * nsf * sqm[std::abs( ipp )], 0 ); + fo[4] = cxmake( imp * nsf * sqm[std::abs( imp )], 0 ); + fo[5] = cxmake( ipp * sqm[std::abs( imp )], 0 ); } else { @@ -620,8 +592,6 @@ fptype( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype)0.5 }; fptype omega[2] = { fpsqrt( pvec0 + pp ), 0. }; omega[1] = fmass / omega[0]; - const int ip = ( 1 + nh ) / 2; // NB: Fortran is (3+nh)/2 because omega(2) has indexes 1,2 and not 0,1 - const int im = ( 1 - nh ) / 2; // NB: Fortran is (3-nh)/2 because omega(2) has indexes 1,2 and not 0,1 const fptype sfomeg[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; const fptype pp3 = fpmax( pp + pvec3, 0. ); const cxtype chi[2] = { cxmake( fpsqrt( pp3 * (fptype)0.5 / pp ), 0. ), @@ -692,20 +662,10 @@ const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( 2. * pvec0 ) : cxmake( (fptype)nh * pvec1, -pvec2 ) / sqp0p3 ) }; #endif - if( nh == 1 ) - { - fo[2] = chi[0]; - fo[3] = chi[1]; - fo[4] = cxzero_sv(); - fo[5] = cxzero_sv(); - } - else - { - fo[2] = cxzero_sv(); - fo[3] = cxzero_sv(); - fo[4] = chi[1]; - fo[5] = chi[0]; - } + fo[2] = ip * chi[0]; + fo[3] = ip * chi[1]; + fo[4] = im * chi[1]; + fo[5] = im * chi[0]; } mgDebug( 1, __FUNCTION__ ); return; @@ -730,19 +690,13 @@ fo[0] = cxmake( pvec3 * (fptype)nsf, pvec3 * (fptype)nsf ); fo[1] = cxzero_sv(); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; const cxtype_sv csqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); fo[3] = cxzero_sv(); fo[4] = cxzero_sv(); - if( nh == 1 ) - { - fo[2] = csqp0p3; - fo[5] = cxzero_sv(); - } - else - { - fo[2] = cxzero_sv(); - fo[5] = csqp0p3; - } + fo[2] = ip * csqp0p3; + fo[5] = im * csqp0p3; mgDebug( 1, __FUNCTION__ ); return; } @@ -766,22 +720,14 @@ fo[0] = cxmake( -pvec3 * (fptype)nsf, pvec3 * (fptype)nsf ); // remember pvec0 == -pvec3 fo[1] = cxzero_sv(); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; const cxtype_sv chi1 = cxmake( -nhel, 0. ) * fpsqrt( -2. * pvec3 ); - if( nh == 1 ) - { - fo[2] = cxzero_sv(); - fo[3] = chi1; - fo[4] = cxzero_sv(); - fo[5] = cxzero_sv(); - } - else - { - fo[2] = cxzero_sv(); - fo[3] = cxzero_sv(); - fo[4] = chi1; - //fo[5] = chi1; // AV: BUG! - fo[5] = cxzero_sv(); // AV: BUG FIX - } + fo[2] = cxzero_sv(); + fo[3] = ip * chi1; + fo[4] = im * chi1; + //fo[5] = chi1; // AV: BUG! + fo[5] = cxzero_sv(); // AV: BUG FIX mgDebug( 1, __FUNCTION__ ); return; } @@ -808,24 +754,16 @@ fo[0] = cxmake( pvec0 * (fptype)nsf, pvec3 * (fptype)nsf ); fo[1] = cxmake( pvec1 * (fptype)nsf, pvec2 * (fptype)nsf ); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); - if( nh == 1 ) - { - fo[2] = chi0; - fo[3] = chi1; - fo[4] = cxzero_sv(); - fo[5] = cxzero_sv(); - } - else - { - fo[2] = cxzero_sv(); - fo[3] = cxzero_sv(); - fo[4] = chi1; - fo[5] = chi0; - } + fo[2] = ip * chi0; + fo[3] = ip * chi1; + fo[4] = im * chi1; + fo[5] = im * chi0; mgDebug( 1, __FUNCTION__ ); return; } diff --git a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas_fd.h b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas_fd.h index 205708a5ba..4e7da2197b 100644 --- a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas_fd.h +++ b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas_fd.h @@ -2,10 +2,10 @@ ! Copyright (C) 2010 The MadGraph5_aMC@NLO development team and contributors. ! Created by: J. Alwall (Sep 2010) for the MG5aMC CPP backend. !========================================================================== -! Copyright (C) 2020-2024 CERN and UCLouvain. +! Copyright (C) 2020-2026 CERN and UCLouvain. ! Licensed under the GNU Lesser General Public License (version 3 or later). ! Modified by: O. Mattelaer (Mar 2020) for the MG5aMC CUDACPP plugin. -! Further modified by: O. Mattelaer, A. Valassi (2020-2024) for the MG5aMC CUDACPP plugin. +! Further modified by: O. Mattelaer, A. Valassi, F. Stloukal, Z. Wettersten (2020-2026) for the MG5aMC CUDACPP plugin. !========================================================================== //-------------------------------------------------------------------------- @@ -204,6 +204,8 @@ fi[0] = cxmake( -pvec0 * (fptype)nsf, -pvec3 * (fptype)nsf ); fi[1] = cxmake( -pvec1 * (fptype)nsf, -pvec2 * (fptype)nsf ); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3+nh)/2 because omega(2) has indexes 1,2 + const int im = ( 1 - nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2 if( fmass != 0. ) { #ifndef MGONGPU_CPPSIMD @@ -215,8 +217,6 @@ // In C++ ixxxxx, use a single ip/im numbering that is valid both for pp==0 and pp>0, which have two numbering schemes in Fortran ixxxxx: // for pp==0, Fortran sqm(0:1) has indexes 0,1 as in C++; but for Fortran pp>0, omega(2) has indexes 1,2 and not 0,1 // NB: this is only possible in ixxxx, but in oxxxxx two different numbering schemes must be used - const int ip = ( 1 + nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3+nh)/2 because omega(2) has indexes 1,2 - const int im = ( 1 - nh ) / 2; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2 #ifndef MGONGPU_CPPSIMD if( pp == 0. ) { @@ -298,20 +298,10 @@ const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), ( sqp0p3 == 0. ? cxmake( -(fptype)nhel * fpsqrt( 2. * pvec0 ), 0. ) : cxmake( (fptype)nh * pvec1, pvec2 ) / sqp0p3 ) }; #endif - if( nh == 1 ) - { - fi[2] = cxzero_sv(); - fi[3] = cxzero_sv(); - fi[4] = chi[0]; - fi[5] = chi[1]; - } - else - { - fi[2] = chi[1]; - fi[3] = chi[0]; - fi[4] = cxzero_sv(); - fi[5] = cxzero_sv(); - } + fi[2] = im * chi[1]; + fi[3] = im * chi[0]; + fi[4] = ip * chi[0]; + fi[5] = ip * chi[1]; } mgDebug( 1, __FUNCTION__ ); return; @@ -336,18 +326,12 @@ fi[0] = cxmake( -pvec3 * (fptype)nsf, -pvec3 * (fptype)nsf ); fi[1] = cxzero_sv(); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; const cxtype_sv sqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); fi[2] = fi[1]; - if( nh == 1 ) - { - fi[3] = fi[1]; - fi[4] = sqp0p3; - } - else - { - fi[3] = sqp0p3; - fi[4] = fi[1]; - } + fi[3] = ip * fi[1]; + fi[4] = im * fi[1]; fi[5] = fi[1]; mgDebug( 1, __FUNCTION__ ); return; @@ -372,19 +356,13 @@ fi[0] = cxmake( pvec3 * (fptype)nsf, -pvec3 * (fptype)nsf ); fi[1] = cxzero_sv(); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; const cxtype_sv chi = cxmake( -(fptype)nhel * fpsqrt( -2. * pvec3 ), 0. ); fi[3] = cxzero_sv(); fi[4] = cxzero_sv(); - if( nh == 1 ) - { - fi[2] = cxzero_sv(); - fi[5] = chi; - } - else - { - fi[2] = chi; - fi[5] = cxzero_sv(); - } + fi[2] = im * chi; + fi[5] = ip * chi; mgDebug( 1, __FUNCTION__ ); return; } @@ -413,24 +391,16 @@ fi[0] = cxmake( -pvec0 * (fptype)nsf, -pvec3 * (fptype)nsf ); // AV: BUG FIX fi[1] = cxmake( -pvec1 * (fptype)nsf, -pvec2 * (fptype)nsf ); // AV: BUG FIX const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); - if( nh == 1 ) - { - fi[2] = cxzero_sv(); - fi[3] = cxzero_sv(); - fi[4] = chi0; - fi[5] = chi1; - } - else - { - fi[2] = chi1; - fi[3] = chi0; - fi[4] = cxzero_sv(); - fi[5] = cxzero_sv(); - } + fi[2] = im * chi1; + fi[3] = im * chi0; + fi[4] = ip * chi0; + fi[5] = ip * chi1; mgDebug( 1, __FUNCTION__ ); return; } @@ -712,6 +682,8 @@ fo[0] = cxmake( pvec0 * (fptype)nsf, pvec3 * (fptype)nsf ); fo[1] = cxmake( pvec1 * (fptype)nsf, pvec2 * (fptype)nsf ); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; // NB: Fortran is (3+nh)/2 because omega(2) has indexes 1,2 and not 0,1 + const int im = ( 1 - nh ) / 2; // NB: Fortran is (3-nh)/2 because omega(2) has indexes 1,2 and not 0,1 if( fmass != 0. ) { #ifndef MGONGPU_CPPSIMD @@ -722,12 +694,12 @@ fptype sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here - const int ip = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - const int im = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - fo[2] = cxmake( im * sqm[std::abs( ip )], 0 ); - fo[3] = cxmake( ip * nsf * sqm[std::abs( ip )], 0 ); - fo[4] = cxmake( im * nsf * sqm[std::abs( im )], 0 ); - fo[5] = cxmake( ip * sqm[std::abs( im )], 0 ); + const int ipp = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + const int imp = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + fo[2] = cxmake( imp * sqm[std::abs( ipp )], 0 ); + fo[3] = cxmake( ipp * nsf * sqm[std::abs( ipp )], 0 ); + fo[4] = cxmake( imp * nsf * sqm[std::abs( imp )], 0 ); + fo[5] = cxmake( ipp * sqm[std::abs( imp )], 0 ); } else { @@ -735,8 +707,6 @@ fptype( 1 + nsf - ( 1 - nsf ) * nh ) * (fptype)0.5 }; fptype omega[2] = { fpsqrt( pvec0 + pp ), 0. }; omega[1] = fmass / omega[0]; - const int ip = ( 1 + nh ) / 2; // NB: Fortran is (3+nh)/2 because omega(2) has indexes 1,2 and not 0,1 - const int im = ( 1 - nh ) / 2; // NB: Fortran is (3-nh)/2 because omega(2) has indexes 1,2 and not 0,1 const fptype sfomeg[2] = { sf[0] * omega[ip], sf[1] * omega[im] }; const fptype pp3 = fpmax( pp + pvec3, 0. ); const cxtype chi[2] = { cxmake( fpsqrt( pp3 * (fptype)0.5 / pp ), 0. ), @@ -807,20 +777,10 @@ const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( 2. * pvec0 ) : cxmake( (fptype)nh * pvec1, -pvec2 ) / sqp0p3 ) }; #endif - if( nh == 1 ) - { - fo[2] = chi[0]; - fo[3] = chi[1]; - fo[4] = cxzero_sv(); - fo[5] = cxzero_sv(); - } - else - { - fo[2] = cxzero_sv(); - fo[3] = cxzero_sv(); - fo[4] = chi[1]; - fo[5] = chi[0]; - } + fo[2] = ip * chi[0]; + fo[3] = ip * chi[1]; + fo[4] = im * chi[1]; + fo[5] = im * chi[0]; } mgDebug( 1, __FUNCTION__ ); return; @@ -845,19 +805,13 @@ fo[0] = cxmake( pvec3 * (fptype)nsf, pvec3 * (fptype)nsf ); fo[1] = cxzero_sv(); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; const cxtype_sv csqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); + fo[2] = ip * csqp0p3; fo[3] = cxzero_sv(); fo[4] = cxzero_sv(); - if( nh == 1 ) - { - fo[2] = csqp0p3; - fo[5] = cxzero_sv(); - } - else - { - fo[2] = cxzero_sv(); - fo[5] = csqp0p3; - } + fo[5] = im * csqp0p3; mgDebug( 1, __FUNCTION__ ); return; } @@ -881,22 +835,14 @@ fo[0] = cxmake( -pvec3 * (fptype)nsf, pvec3 * (fptype)nsf ); // remember pvec0 == -pvec3 fo[1] = cxzero_sv(); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; const cxtype_sv chi1 = cxmake( -nhel, 0. ) * fpsqrt( -2. * pvec3 ); - if( nh == 1 ) - { - fo[2] = cxzero_sv(); - fo[3] = chi1; - fo[4] = cxzero_sv(); - fo[5] = cxzero_sv(); - } - else - { - fo[2] = cxzero_sv(); - fo[3] = cxzero_sv(); - fo[4] = chi1; - //fo[5] = chi1; // AV: BUG! - fo[5] = cxzero_sv(); // AV: BUG FIX - } + fo[2] = cxzero_sv(); + fo[3] = ip * chi1; + fo[4] = im * chi1; + //fo[5] = chi1; // AV: BUG! + fo[5] = cxzero_sv(); // AV: BUG FIX mgDebug( 1, __FUNCTION__ ); return; } @@ -923,24 +869,16 @@ fo[0] = cxmake( pvec0 * (fptype)nsf, pvec3 * (fptype)nsf ); fo[1] = cxmake( pvec1 * (fptype)nsf, pvec2 * (fptype)nsf ); const int nh = nhel * nsf; + const int ip = ( 1 + nh ) / 2; + const int im = ( 1 - nh ) / 2; //const float sqp0p3 = sqrtf( pvec0 + pvec3 ) * nsf; // AV: why force a float here? const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); - if( nh == 1 ) - { - fo[2] = chi0; - fo[3] = chi1; - fo[4] = cxzero_sv(); - fo[5] = cxzero_sv(); - } - else - { - fo[2] = cxzero_sv(); - fo[3] = cxzero_sv(); - fo[4] = chi1; - fo[5] = chi0; - } + fo[2] = ip * chi0; + fo[3] = ip * chi1; + fo[4] = im * chi1; + fo[5] = im * chi0; mgDebug( 1, __FUNCTION__ ); return; } From b5f5e341c9f08adacc4a69558668d306dc871324 Mon Sep 17 00:00:00 2001 From: Zenny Date: Fri, 17 Jul 2026 12:22:06 +0200 Subject: [PATCH 2/2] squashed implementation of branch-free helicities in external partons --- .../aloha/template_files/gpu/helas.h | 55 ++-- .../aloha/template_files/gpu/helas_fd.h | 236 ++++++++---------- 2 files changed, 137 insertions(+), 154 deletions(-) diff --git a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h index 5b7db11740..120d6f0fb0 100644 --- a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h +++ b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h @@ -271,10 +271,10 @@ const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), ( sqp0p3 == 0. ? cxmake( -(fptype)nhel * fpsqrt( 2. * pvec0 ), 0. ) : cxmake( (fptype)nh * pvec1, pvec2 ) / sqp0p3 ) }; #endif - fi[2] = im * chi[1]; - fi[3] = im * chi[0]; - fi[4] = ip * chi[0]; - fi[5] = ip * chi[1]; + fi[2] = (fptype)im * chi[1]; + fi[3] = (fptype)im * chi[0]; + fi[4] = (fptype)ip * chi[0]; + fi[5] = (fptype)ip * chi[1]; } mgDebug( 1, __FUNCTION__ ); return; @@ -303,8 +303,8 @@ const int im = ( 1 - nh ) / 2; const cxtype_sv sqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); fi[2] = fi[1]; - fi[3] = ip * fi[1]; - fi[4] = im * fi[1]; + fi[3] = (fptype)ip * fi[1] + (fptype)im * sqp0p3; + fi[4] = (fptype)im * fi[1] + (fptype)ip * sqp0p3; fi[5] = fi[1]; mgDebug( 1, __FUNCTION__ ); return; @@ -334,8 +334,8 @@ const cxtype_sv chi = cxmake( -(fptype)nhel * fpsqrt( -2. * pvec3 ), 0. ); fi[3] = cxzero_sv(); fi[4] = cxzero_sv(); - fi[2] = im * chi; - fi[5] = ip * chi; + fi[2] = (fptype)im * chi; + fi[5] = (fptype)ip * chi; mgDebug( 1, __FUNCTION__ ); return; } @@ -370,10 +370,10 @@ const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); - fi[2] = im * chi1; - fi[3] = im * chi0; - fi[4] = ip * chi0; - fi[5] = ip * chi1; + fi[2] = (fptype)im * chi1; + fi[3] = (fptype)im * chi0; + fi[4] = (fptype)ip * chi0; + fi[5] = (fptype)ip * chi1; mgDebug( 1, __FUNCTION__ ); return; } @@ -579,8 +579,8 @@ fptype sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here - const int ipp = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - const int imp = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + const int ipp = -im * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + const int imp = ip * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ fo[2] = cxmake( imp * sqm[std::abs( ipp )], 0 ); fo[3] = cxmake( ipp * nsf * sqm[std::abs( ipp )], 0 ); fo[4] = cxmake( imp * nsf * sqm[std::abs( imp )], 0 ); @@ -662,10 +662,10 @@ const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( 2. * pvec0 ) : cxmake( (fptype)nh * pvec1, -pvec2 ) / sqp0p3 ) }; #endif - fo[2] = ip * chi[0]; - fo[3] = ip * chi[1]; - fo[4] = im * chi[1]; - fo[5] = im * chi[0]; + fo[2] = (fptype)ip * chi[0]; + fo[3] = (fptype)ip * chi[1]; + fo[4] = (fptype)im * chi[1]; + fo[5] = (fptype)im * chi[0]; } mgDebug( 1, __FUNCTION__ ); return; @@ -695,8 +695,8 @@ const cxtype_sv csqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); fo[3] = cxzero_sv(); fo[4] = cxzero_sv(); - fo[2] = ip * csqp0p3; - fo[5] = im * csqp0p3; + fo[2] = (fptype)ip * csqp0p3; + fo[5] = (fptype)im * csqp0p3; mgDebug( 1, __FUNCTION__ ); return; } @@ -724,10 +724,9 @@ const int im = ( 1 - nh ) / 2; const cxtype_sv chi1 = cxmake( -nhel, 0. ) * fpsqrt( -2. * pvec3 ); fo[2] = cxzero_sv(); - fo[3] = ip * chi1; - fo[4] = im * chi1; - //fo[5] = chi1; // AV: BUG! - fo[5] = cxzero_sv(); // AV: BUG FIX + fo[3] = (fptype)ip * chi1; + fo[4] = (fptype)im * chi1; + fo[5] = cxzero_sv(); mgDebug( 1, __FUNCTION__ ); return; } @@ -760,10 +759,10 @@ const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); - fo[2] = ip * chi0; - fo[3] = ip * chi1; - fo[4] = im * chi1; - fo[5] = im * chi0; + fo[2] = (fptype)ip * chi0; + fo[3] = (fptype)ip * chi1; + fo[4] = (fptype)im * chi1; + fo[5] = (fptype)im * chi0; mgDebug( 1, __FUNCTION__ ); return; } diff --git a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas_fd.h b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas_fd.h index 4e7da2197b..aac2742c06 100644 --- a/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas_fd.h +++ b/epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas_fd.h @@ -152,8 +152,7 @@ const int ipar // input: particle# out of npar ) ALWAYS_INLINE; - -//-------------------------------------------------------------------------- + //-------------------------------------------------------------------------- // Compute the direction n[5] of the gauge q[5] __host__ __device__ INLINE void @@ -161,25 +160,24 @@ fptype n[] // output: direction ) ALWAYS_INLINE; - //-------------------------------------------------------------------------- // Compute a propagator factor d out of gauge q[5] and a mass __host__ __device__ INLINE void calculate_propagator_factor( const fptype_sv q[5], // input: gauge const fptype_sv mass, // input: mass - fptype_sv *d // output: propagator factor + fptype_sv* d // output: propagator factor ) ALWAYS_INLINE; //-------------------------------------------------------------------------- // multiply by propagation factor from m and wawefunctionsin[] and output them // as wavefunctionout[] - template< class W_ACCESS> + template __host__ __device__ INLINE void multiply_propagator_factor( const fptype wavefunctionsin[], // input: wavefunctions const fptype m, // input: mass fptype wavefunctionsout[] // output: wavefunctions ) ALWAYS_INLINE; -//========================================================================== + //========================================================================== // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt] template @@ -298,10 +296,10 @@ const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), ( sqp0p3 == 0. ? cxmake( -(fptype)nhel * fpsqrt( 2. * pvec0 ), 0. ) : cxmake( (fptype)nh * pvec1, pvec2 ) / sqp0p3 ) }; #endif - fi[2] = im * chi[1]; - fi[3] = im * chi[0]; - fi[4] = ip * chi[0]; - fi[5] = ip * chi[1]; + fi[2] = (fptype)im * chi[1]; + fi[3] = (fptype)im * chi[0]; + fi[4] = (fptype)ip * chi[0]; + fi[5] = (fptype)ip * chi[1]; } mgDebug( 1, __FUNCTION__ ); return; @@ -330,8 +328,8 @@ const int im = ( 1 - nh ) / 2; const cxtype_sv sqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); fi[2] = fi[1]; - fi[3] = ip * fi[1]; - fi[4] = im * fi[1]; + fi[3] = (fptype)ip * fi[1] + (fptype)im * sqp0p3; + fi[4] = (fptype)im * fi[1] + (fptype)ip * sqp0p3; fi[5] = fi[1]; mgDebug( 1, __FUNCTION__ ); return; @@ -361,8 +359,8 @@ const cxtype_sv chi = cxmake( -(fptype)nhel * fpsqrt( -2. * pvec3 ), 0. ); fi[3] = cxzero_sv(); fi[4] = cxzero_sv(); - fi[2] = im * chi; - fi[5] = ip * chi; + fi[2] = (fptype)im * chi; + fi[5] = (fptype)ip * chi; mgDebug( 1, __FUNCTION__ ); return; } @@ -397,10 +395,10 @@ const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, pvec2 / sqp0p3 ); - fi[2] = im * chi1; - fi[3] = im * chi0; - fi[4] = ip * chi0; - fi[5] = ip * chi1; + fi[2] = (fptype)im * chi1; + fi[3] = (fptype)im * chi0; + fi[4] = (fptype)ip * chi0; + fi[5] = (fptype)ip * chi1; mgDebug( 1, __FUNCTION__ ); return; } @@ -427,21 +425,21 @@ const fptype_sv& pvec2 = M_ACCESS::kernelAccessIp4IparConst( momenta, 2, ipar ); const fptype_sv& pvec3 = M_ACCESS::kernelAccessIp4IparConst( momenta, 3, ipar ); cxtype_sv* vc = W_ACCESS::kernelAccess( wavefunctions ); // needs to be size 7 instead of 6 - const fptype sqh = fpsqrt( 0.5 ); // AV this is > 0! + const fptype sqh = fpsqrt( 0.5 ); // AV this is > 0! const fptype hel = nhel; vc[0] = cxmake( pvec0 * (fptype)nsv, pvec3 * (fptype)nsv ); vc[1] = cxmake( pvec1 * (fptype)nsv, pvec2 * (fptype)nsv ); // FD gauge - const cxtype_sv cI = cxmake( 0 + fptype_sv{ 0 }, 1 + fptype_sv{ 0 } ); + const cxtype_sv cI = cxmake( 0 + fptype_sv{ 0 }, 1 + fptype_sv{ 0 } ); #ifdef MGONGPU_CPPSIMD fptype_sv nA[5]; fptype_sv nB[5]; #endif fptype_sv n[5]; fptype_sv nk; - const fptype_sv zero{0.}; - const fptype_sv one{0.}; + const fptype_sv zero{ 0. }; + const fptype_sv one{ 0. }; if( vmass != 0. ) { @@ -484,33 +482,33 @@ //FD gauge if( pp > 0. ) { - n[0] = ( pvec0 >= zero) ? one : -one; - n[1] = -pvec1/pp; - n[2] = -pvec2/pp; - n[3] = -pvec3/pp; + n[0] = ( pvec0 >= zero ) ? one : -one; + n[1] = -pvec1 / pp; + n[2] = -pvec2 / pp; + n[3] = -pvec3 / pp; n[4] = zero; } else { - n[0] = ( pvec0 >= zero) ? one : -one; + n[0] = ( pvec0 >= zero ) ? one : -one; n[1] = zero; n[2] = zero; - n[3] = ( pvec0 >= zero) ? -one : one; + n[3] = ( pvec0 >= zero ) ? -one : one; } + nk = n[0] * pvec0 - n[1] * pvec1 - n[2] * pvec2 - n[3] * pvec3; - nk = n[0]*pvec0 - n[1]*pvec1 - n[2]*pvec2 - n[3]*pvec3; - - if ( abs(nhel) == 1) + if( abs( nhel ) == 1 ) { vc[6] = cxzero_sv(); } - else{ - vc[2] = cxmake( -vmass/nk * n[0], zero ); - vc[3] = cxmake( -vmass/nk * n[1], zero ); - vc[4] = cxmake( -vmass/nk * n[2], zero ); - vc[5] = cxmake( -vmass/nk * n[3], zero ); - vc[6] = static_cast(nsv)*cI; + else + { + vc[2] = cxmake( -vmass / nk * n[0], zero ); + vc[3] = cxmake( -vmass / nk * n[1], zero ); + vc[4] = cxmake( -vmass / nk * n[2], zero ); + vc[5] = cxmake( -vmass / nk * n[3], zero ); + vc[6] = static_cast( nsv ) * cI; } #else @@ -547,10 +545,10 @@ //FD gauge //branch A - nA[0] = fpternary( pvec0 >= zero , one , -one); - nA[1] = -pvec1/pp; - nA[2] = -pvec2/pp; - nA[3] = -pvec3/pp; + nA[0] = fpternary( pvec0 >= zero, one, -one ); + nA[1] = -pvec1 / pp; + nA[2] = -pvec2 / pp; + nA[3] = -pvec3 / pp; nA[4] = zero; //branch B @@ -559,23 +557,23 @@ nB[2] = zero; nB[3] = -nA[0]; - const fptype_sv b_A = fpternary(pp > zero, one , zero); - const fptype_sv b_B = fpternary(pp <= zero , one , zero); + const fptype_sv b_A = fpternary( pp > zero, one, zero ); + const fptype_sv b_B = fpternary( pp <= zero, one, zero ); - n[0] = nA[0]*b_A + nB[0]*b_B; - n[1] = nA[1]*b_A + nB[1]*b_B; - n[2] = nA[2]*b_A + nB[2]*b_B; - n[3] = nA[3]*b_A + nB[3]*b_B; + n[0] = nA[0] * b_A + nB[0] * b_B; + n[1] = nA[1] * b_A + nB[1] * b_B; + n[2] = nA[2] * b_A + nB[2] * b_B; + n[3] = nA[3] * b_A + nB[3] * b_B; n[4] = nA[4]; - nk = n[0]*pvec0 - n[1]*pvec1 - n[2]*pvec2 - n[3]*pvec3; + nk = n[0] * pvec0 - n[1] * pvec1 - n[2] * pvec2 - n[3] * pvec3; - const bool_v mask3 = { (abs(nhel) == 1 ? 1 : 0) }; // first element replicated - vc[2] = cxternary( mask3, vc[2], cxmake( -vmass/nk * n[0], zero)); - vc[3] = cxternary( mask3, vc[3], cxmake( -vmass/nk * n[1], zero)); - vc[4] = cxternary( mask3, vc[4], cxmake( -vmass/nk * n[2], zero)); - vc[5] = cxternary( mask3, vc[5], cxmake( -vmass/nk * n[3], zero)); - vc[6] = cxternary( mask3, cxzero_sv(), -static_cast(nsv)*cI); + const bool_v mask3 = { ( abs( nhel ) == 1 ? 1 : 0 ) }; // first element replicated + vc[2] = cxternary( mask3, vc[2], cxmake( -vmass / nk * n[0], zero ) ); + vc[3] = cxternary( mask3, vc[3], cxmake( -vmass / nk * n[1], zero ) ); + vc[4] = cxternary( mask3, vc[4], cxmake( -vmass / nk * n[2], zero ) ); + vc[5] = cxternary( mask3, vc[5], cxmake( -vmass / nk * n[3], zero ) ); + vc[6] = cxternary( mask3, cxzero_sv(), -static_cast( nsv ) * cI ); #endif } else @@ -694,8 +692,8 @@ fptype sqm[2] = { fpsqrt( std::abs( fmass ) ), 0. }; // possibility of negative fermion masses //sqm[1] = ( fmass < 0. ? -abs( sqm[0] ) : abs( sqm[0] ) ); // AV: why abs here? sqm[1] = ( fmass < 0. ? -sqm[0] : sqm[0] ); // AV: removed an abs here - const int ipp = -( ( 1 - nh ) / 2 ) * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ - const int imp = ( 1 + nh ) / 2 * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + const int ipp = -im * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ + const int imp = ip * nhel; // NB: Fortran sqm(0:1) also has indexes 0,1 as in C++ fo[2] = cxmake( imp * sqm[std::abs( ipp )], 0 ); fo[3] = cxmake( ipp * nsf * sqm[std::abs( ipp )], 0 ); fo[4] = cxmake( imp * nsf * sqm[std::abs( imp )], 0 ); @@ -777,10 +775,10 @@ const cxtype_sv chi[2] = { cxmake( sqp0p3, 0. ), ( sqp0p3 == 0. ? cxmake( -nhel, 0. ) * fpsqrt( 2. * pvec0 ) : cxmake( (fptype)nh * pvec1, -pvec2 ) / sqp0p3 ) }; #endif - fo[2] = ip * chi[0]; - fo[3] = ip * chi[1]; - fo[4] = im * chi[1]; - fo[5] = im * chi[0]; + fo[2] = (fptype)ip * chi[0]; + fo[3] = (fptype)ip * chi[1]; + fo[4] = (fptype)im * chi[1]; + fo[5] = (fptype)im * chi[0]; } mgDebug( 1, __FUNCTION__ ); return; @@ -808,10 +806,10 @@ const int ip = ( 1 + nh ) / 2; const int im = ( 1 - nh ) / 2; const cxtype_sv csqp0p3 = cxmake( fpsqrt( 2. * pvec3 ) * (fptype)nsf, 0. ); - fo[2] = ip * csqp0p3; + fo[2] = (fptype)ip * csqp0p3; fo[3] = cxzero_sv(); fo[4] = cxzero_sv(); - fo[5] = im * csqp0p3; + fo[5] = (fptype)im * csqp0p3; mgDebug( 1, __FUNCTION__ ); return; } @@ -839,10 +837,9 @@ const int im = ( 1 - nh ) / 2; const cxtype_sv chi1 = cxmake( -nhel, 0. ) * fpsqrt( -2. * pvec3 ); fo[2] = cxzero_sv(); - fo[3] = ip * chi1; - fo[4] = im * chi1; - //fo[5] = chi1; // AV: BUG! - fo[5] = cxzero_sv(); // AV: BUG FIX + fo[3] = (fptype)ip * chi1; + fo[4] = (fptype)im * chi1; + fo[5] = cxzero_sv(); mgDebug( 1, __FUNCTION__ ); return; } @@ -875,10 +872,10 @@ const fptype_sv sqp0p3 = fpsqrt( pvec0 + pvec3 ) * (fptype)nsf; const cxtype_sv chi0 = cxmake( sqp0p3, 0. ); const cxtype_sv chi1 = cxmake( (fptype)nh * pvec1 / sqp0p3, -pvec2 / sqp0p3 ); - fo[2] = ip * chi0; - fo[3] = ip * chi1; - fo[4] = im * chi1; - fo[5] = im * chi0; + fo[2] = (fptype)ip * chi0; + fo[3] = (fptype)ip * chi1; + fo[4] = (fptype)im * chi1; + fo[5] = (fptype)im * chi0; mgDebug( 1, __FUNCTION__ ); return; } @@ -888,21 +885,19 @@ __host__ __device__ INLINE void define_gauge_dir( const cxtype_sv q[5], // input: gauge fptype_sv n[5] ) // output: direction - { - const fptype_sv qabs2 = q[1].real()*q[1].real() - + q[2].real()*q[2].real() - + q[3].real()*q[3].real(); + { + const fptype_sv qabs2 = q[1].real() * q[1].real() + q[2].real() * q[2].real() + q[3].real() * q[3].real(); - const fptype_sv one = 1. + fptype_sv{0}; - const fptype_sv zero = 0. + fptype_sv{0}; + const fptype_sv one = 1. + fptype_sv{ 0 }; + const fptype_sv zero = 0. + fptype_sv{ 0 }; #ifndef MGONGPU_CPPSIMD - if (qabs2 > 0.0) + if( qabs2 > 0.0 ) { - const fptype_sv qabs = fpsqrt(qabs2); + const fptype_sv qabs = fpsqrt( qabs2 ); - n[0] = fpternary( q[0].real() >= 0. , one , -one); + n[0] = fpternary( q[0].real() >= 0., one, -one ); n[1] = -q[1].real() / qabs; n[2] = -q[2].real() / qabs; n[3] = -q[3].real() / qabs; @@ -910,45 +905,44 @@ } else { - n[0] = fpternary( q[0].real() >= 0. , one , -one ); + n[0] = fpternary( q[0].real() >= 0., one, -one ); n[1] = zero; n[2] = zero; - n[3] = fpternary( q[0].real() >= 0. , one , -one); //possible error in Fortran + n[3] = fpternary( q[0].real() >= 0., one, -one ); //possible error in Fortran n[4] = zero; } #else - const fptype_sv qabs = fpsqrt(qabs2); - const bool_v qsign = (qabs2 > 0.); - n[0] = fpternary( q[0].real() >= 0. , one , -one); - n[1] = fpternary( qsign , -q[1].real() / qabs , zero ); - n[2] = fpternary( qsign , -q[2].real() / qabs , zero ); - n[3] = fpternary( qsign , -q[3].real() / qabs , fpternary( q[0].real() >= 0. , one , -one)); + const fptype_sv qabs = fpsqrt( qabs2 ); + const bool_v qsign = ( qabs2 > 0. ); + n[0] = fpternary( q[0].real() >= 0., one, -one ); + n[1] = fpternary( qsign, -q[1].real() / qabs, zero ); + n[2] = fpternary( qsign, -q[2].real() / qabs, zero ); + n[3] = fpternary( qsign, -q[3].real() / qabs, fpternary( q[0].real() >= 0., one, -one ) ); n[4] = zero; #endif - } + } -//-------------------------------------------------------------------------- -// Compute propagator factor n[5] of the gauge q[5] + //-------------------------------------------------------------------------- + // Compute propagator factor n[5] of the gauge q[5] __host__ __device__ INLINE void calculate_propagator_factor( const cxtype_sv q[5], // input: gauge const fptype mass, // input: mass - fptype_sv *d ) // output: propagator factor + fptype_sv* d ) // output: propagator factor { - const fptype_sv one = 1. + fptype_sv{0}; - const fptype_sv q2 = q[0].real()*q[0].real() - ( q[1].real()*q[1].real() + q[2].real()*q[2].real() + q[3].real()*q[3].real() ); - *d = one / (q2 - mass*mass); + const fptype_sv one = 1. + fptype_sv{ 0 }; + const fptype_sv q2 = q[0].real() * q[0].real() - ( q[1].real() * q[1].real() + q[2].real() * q[2].real() + q[3].real() * q[3].real() ); + *d = one / ( q2 - mass * mass ); } -//-------------------------------------------------------------------------- -// multiply by propagation factor from m and wawefunctionsin[] and output them -// as wavefunctionout[] - template< class W_ACCESS> + //-------------------------------------------------------------------------- + // multiply by propagation factor from m and wawefunctionsin[] and output them + // as wavefunctionout[] + template __host__ __device__ INLINE void multiply_propagator_factor( const fptype wavefunctionsin[], // input: wavefunctions const fptype m, // input: mass fptype wavefunctionsout[] ) // output: wavefunctions { - const cxtype_sv* win = W_ACCESS::kernelAccessConst( wavefunctionsin ); cxtype_sv* wout = W_ACCESS::kernelAccess( wavefunctionsout ); @@ -957,20 +951,20 @@ fptype_sv n[5]; cxtype_sv w0[5], w1[5]; - const cxtype_sv cI = cxmake( 0 + fptype_sv{ 0 }, 1. + fptype_sv{ 0 } ); + const cxtype_sv cI = cxmake( 0 + fptype_sv{ 0 }, 1. + fptype_sv{ 0 } ); // Construct q - q[0] = cxmake( -win[0].real(), 0.); - q[1] = cxmake( -win[1].real(), 0.); - q[2] = cxmake( -win[1].imag(), 0.); - q[3] = cxmake( -win[0].imag(), 0.); - q[4] = -cI*m; + q[0] = cxmake( -win[0].real(), 0. ); + q[1] = cxmake( -win[1].real(), 0. ); + q[2] = cxmake( -win[1].imag(), 0. ); + q[3] = cxmake( -win[0].imag(), 0. ); + q[4] = -cI * m; // Copy first two components wout[0] = win[0]; wout[1] = win[1]; - define_gauge_dir(q, n); + define_gauge_dir( q, n ); w0[0] = win[2]; w0[1] = win[3]; @@ -979,31 +973,21 @@ w0[4] = win[6]; fptype_sv nq = - n[0]*q[0].real() - - n[1]*q[1].real() - - n[2]*q[2].real() - - n[3]*q[3].real(); + n[0] * q[0].real() - n[1] * q[1].real() - n[2] * q[2].real() - n[3] * q[3].real(); - calculate_propagator_factor(q, m, &d); + calculate_propagator_factor( q, m, &d ); cxtype_sv js1 = - ( n[0]*w0[0] - - n[1]*w0[1] - - n[2]*w0[2] - - n[3]*w0[3] ) / nq; + ( n[0] * w0[0] - n[1] * w0[1] - n[2] * w0[2] - n[3] * w0[3] ) / nq; cxtype_sv js2 = - ( q[0]*w0[0] - - q[1]*w0[1] - - q[2]*w0[2] - - q[3]*w0[3] - - cxconj(q[4]) * w0[4] ) / nq; - - w1[0] = w0[0] - q[0]*js1 - n[0]*js2; - w1[1] = w0[1] - q[1]*js1 - n[1]*js2; - w1[2] = w0[2] - q[2]*js1 - n[2]*js2; - w1[3] = w0[3] - q[3]*js1 - n[3]*js2; - w1[4] = w0[4] - q[4]*js1 - n[4]*js2; + ( q[0] * w0[0] - q[1] * w0[1] - q[2] * w0[2] - q[3] * w0[3] - cxconj( q[4] ) * w0[4] ) / nq; + + w1[0] = w0[0] - q[0] * js1 - n[0] * js2; + w1[1] = w0[1] - q[1] * js1 - n[1] * js2; + w1[2] = w0[2] - q[2] * js1 - n[2] * js2; + w1[3] = w0[3] - q[3] * js1 - n[3] * js2; + w1[4] = w0[4] - q[4] * js1 - n[4] * js2; wout[2] = w1[0]; wout[3] = w1[1];