Skip to content

external partons with branch-free helicity #1095

New issue

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

Sign up for GitHub

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

Already on GitHub? Sign in to your account

Merged
zeniheisser merged 2 commits into from
Jul 17, 2026
Merged

external partons with branch-free helicity #1095

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7fff20d
made external fermions branch-free in helicity, no discernible perfor…
zeniheisser Jul 16, 2026
b5f5e34
squashed implementation of branch-free helicities in external partons
zeniheisser Jul 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
164 changes: 51 additions & 113 deletions epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.
!==========================================================================
//--------------------------------------------------------------------------

Expand Down Expand Up @@ -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
Expand All @@ -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. )
{
Expand Down Expand Up @@ -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;
Expand All @@ -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];

@Qubitol Qubitol Jul 17, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
fi[3] = ip * fi[1];
fi[3] = ip * fi[1] + im * sqp0p3;

fi[4] = im * fi[1];

@Qubitol Qubitol Jul 17, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
fi[4] = im * fi[1];
fi[4] = im * fi[1] + ip * sqp0p3;

fi[5] = fi[1];
Comment on lines 299 to 308

@stloufra stloufra Jul 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This does not seem to be right (even though we never exercise this function).

mgDebug( 1, __FUNCTION__ );
return;
Expand All @@ -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;
}
Expand Down Expand Up @@ -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;

@stloufra stloufra Jul 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not to be pedantic but bit shift would be perfect here instead of /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;
}
Expand Down Expand Up @@ -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
Expand All @@ -607,21 +579,19 @@
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++

@stloufra stloufra Jul 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we reuse ip and im ? :

const int ipp = -im * nhel;

or better yet

const int ipp = ip * nhel;

(not sure if the - adds some clarity from physics point of view)

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
{
const fptype sf[2] = { fptype( 1 + nsf + ( 1 - nsf ) * nh ) * (fptype)0.5,
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. ),
Expand Down Expand Up @@ -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;
Expand All @@ -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;
}
Expand All @@ -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!

@stloufra stloufra Jul 16, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do we need to keep bugs documented in the code?

@Qubitol Qubitol Jul 17, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreee, we can remove them.

fo[5] = cxzero_sv(); // AV: BUG FIX
mgDebug( 1, __FUNCTION__ );
return;
}
Expand All @@ -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;
}
Expand Down
Loading
Loading