external partons with branch-free helicity by zeniheisser · Pull Request #1095 · madgraph5/madgraph4gpu

zeniheisser · 2026-07-16T12:39:18Z

Removed if-else blocks for helicity states in external fermions with branch-free arithmetic; always define floating point ip and im denoting whether it has plus or minus polarisation (1.0 if yes, 0.0 if no) and multiply outputs by this number to automatically set correct spinors to zero. (This trick is already done for massive fermions as well as external vector bosons)

Minimal checks to see if this had any measureable performance impact: Running standalone_cudacpp with no SIMD and commenting out everything but the IXXXXX and OXXXXX calls in sigmaKin, there is no difference in evaluation runtime. Using the internal timers for check_cpp.exe with -p and running ~a billion events (./check_cpp.exe 1024 1024 1024):
(Branch-free) 3a SigmaKin : 1.169649 sec
(Original ver) 3a SigmaKin : 1.172651 sec

I ran a couple of times to see, and even running just external fermions there was no measureable runtime difference.

…mance drop

stloufra

Do we care about the formatting CI not passing?

stloufra · 2026-07-16T12:59:08Z

    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];


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

stloufra · 2026-07-16T13:01:40Z

+    fo[2] = cxzero_sv();
+    fo[3] = ip * chi1;
+    fo[4] = im * chi1;
+    //fo[5] = chi1; // AV: BUG!


Do we need to keep bugs documented in the code?

Agreee, we can remove them.

stloufra · 2026-07-16T13:04:47Z

+    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];


stloufra · 2026-07-16T13:10:05Z

-        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++


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)

stloufra · 2026-07-16T13:10:48Z

    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;


Not to be pedantic but bit shift would be perfect here instead of /2

Qubitol

Thank you for modifying this, it looks much more clean.
I fixed few things which also @stloufra pointed out, please check them.
Also, I guess we can remove the stray comments.

Qubitol · 2026-07-17T08:32:36Z

-      fi[3] = sqp0p3;
-      fi[4] = fi[1];
-    }
+    fi[3] = ip * fi[1];


Suggested change

fi[3] = ip * fi[1];

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

Qubitol · 2026-07-17T08:33:04Z

-      fi[4] = fi[1];
-    }
+    fi[3] = ip * fi[1];
+    fi[4] = im * fi[1];


Suggested change

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

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

Qubitol · 2026-07-17T08:39:26Z

+    fo[2] = cxzero_sv();
+    fo[3] = ip * chi1;
+    fo[4] = im * chi1;
+    //fo[5] = chi1; // AV: BUG!


Agreee, we can remove them.

Qubitol · 2026-07-17T08:40:46Z

-      fi[3] = sqp0p3;
-      fi[4] = fi[1];
-    }
+    fi[3] = ip * fi[1];


Suggested change

fi[3] = ip * fi[1];

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

Qubitol · 2026-07-17T08:41:05Z

-      fi[4] = fi[1];
-    }
+    fi[3] = ip * fi[1];
+    fi[4] = im * fi[1];


Suggested change

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

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

stloufra · 2026-07-17T09:13:00Z

Thanks for the simplification - LGTM

Qubitol · 2026-07-17T09:29:26Z

Do we care about the formatting CI not passing?

Yes, that needs to be fixed.

Qubitol · 2026-07-17T09:32:20Z

+    const int ip = ( 1 + nh ) >> 1; // 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 ) >> 1; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2


I reckon bit-shift makes it less readable, is there a real performance gain into using that?

Probably none whatsoever, most any compiler will cast integer division by two to a bit shift --- but it does look cooler!

FWIW I would go for more readable rather than cooler ;-)

zeniheisser · 2026-07-17T09:39:44Z

@stloufra @Qubitol I've added the bug fixes + another one which didn't appear in the tests due to formatting (cxtype * int implicitly cast to float rather than double) and fixed formatting (also applied to helas_fd.h to make it consistent, although that's not checked by the CI). Also changed integer division to a bit shift --- doesn't matter, but it does look a bit cooler ;-) I'm fine either way, though.

@stloufra , I took your suggestion and recycled the evaluated ip and im for the fermion at rest. I don't know why the plus and minus polarisations are inverted for those, but that's also how the Fortran version are written so I've kept the ipp and imp variables. I presume it's just a way to consistently keep track of rest-state helicity (which isn't well-defined) between initial and final state particles.

roiser · 2026-07-17T10:00:32Z

Do we care about the formatting CI not passing?

Good point, let's discuss at the next occasion

zeniheisser · 2026-07-17T10:14:32Z

For reference, even with the cast to double there is no noticeable runtime difference at O3. Running standalone branchless and original code with everything but IXXXXX and OXXXXXX commented out consecutively for 1024^3 PS points:
(Branch-free) 3a SigmaKin : 1.504317 sec
(Original ver) 3a SigmaKin : 1.526404 sec
(apparently my machine is running a bit hotter today than yesterday hehe)

made external fermions branch-free in helicity, no discernible perfor…

7fff20d

…mance drop

zeniheisser requested review from Qubitol and stloufra July 16, 2026 12:39

stloufra requested changes Jul 16, 2026

View reviewed changes

Qubitol requested changes Jul 17, 2026

View reviewed changes

stloufra approved these changes Jul 17, 2026

View reviewed changes

Qubitol reviewed Jul 17, 2026

View reviewed changes

squashed implementation of branch-free helicities in external partons

b5f5e34

zeniheisser force-pushed the helicity_branching branch from 7eb6f0d to b5f5e34 Compare July 17, 2026 10:22

zeniheisser merged commit 578406f into madgraph5:master Jul 17, 2026
205 checks passed

		const int ip = ( 1 + nh ) >> 1; // 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 ) >> 1; // NB: same as in Fortran pp==0, differs from Fortran pp>0, which is (3-nh)/2 because omega(2) has indexes 1,2

Uh oh!

Conversation

zeniheisser commented Jul 16, 2026

Uh oh!

stloufra left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Qubitol left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

stloufra commented Jul 17, 2026

Uh oh!

Qubitol commented Jul 17, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

zeniheisser commented Jul 17, 2026

Uh oh!

roiser commented Jul 17, 2026

Uh oh!

zeniheisser commented Jul 17, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

5 participants

Qubitol left a comment •

edited

Loading