Four floating point exceptions in CPP launch of pp_ttW (IEEE_DIVIDE_BY_ZERO FPE in vxxxxx function in SIMD mode)

[valassi]

Floating point exceptions in CPP launch of pp_ttW

This is in my branch nobm

[avalassi@itscrd80 gcc11.2/cvmfs] /data/avalassi/GPU2023/madgraph4gpuX/epochX/cudacpp> HRDCOD=1 tlau/lauX.sh -CPP nobm_pp_ttW.mad/
Execute lauX.sh for process nobm_pp_ttW.mad/ in directory /data/avalassi/GPU2023/madgraph4gpuX/epochX/cudacpp
...
Generating 10000 events with run name run_01
survey  run_01 
INFO: compile directory 
compile Source Directory
Using random number seed offset = 21
INFO: Running Survey 
Creating Jobs
Working on SubProcesses
INFO:     P1_gu_ttxwpd  
INFO: Building madevent in madevent_interface.py with 'CPP' matrix elements 
INFO:     P1_gd_ttxwmu  
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
...

[hageboeck]

I investigated this a little for
epochX/cudacpp/nobm_pp_ttW.mad/SubProcesses/P1_gd_ttxwmu

The first floating-point exception is a devision by zero here:

      const fptype_v pzpt = pvec3 / ( pp * pt ) * sqh * hel;

Caused by pt = {0, 0, 0, 0}.

It's a different process, but here's the relevant code:

madgraph4gpu/epochX/cudacpp/gg_ttgg.mad/src/HelAmps_sm.h

Lines 496 to 529 in 603ff90

	const fptype_sv& pp = pvec0; // NB: rewrite the following as in Fortran, using pp instead of pvec0
	const fptype_sv pt = fpsqrt( ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ) );
	vc[2] = cxzero_sv();
	vc[5] = cxmake( hel * pt / pp * sqh, 0. );
	#ifndef MGONGPU_CPPSIMD
	if( pt != 0. )
	{
	const fptype pzpt = pvec3 / ( pp * pt ) * sqh * hel;
	vc[3] = cxmake( -pvec1 * pzpt, -nsv * pvec2 / pt * sqh );
	vc[4] = cxmake( -pvec2 * pzpt, nsv * pvec1 / pt * sqh );
	}
	else
	{
	vc[3] = cxmake( -hel * sqh, 0. );
	// NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs!
	//vc[4] = cxmake( 0, nsv * ( pvec3 < 0. ? -std::abs( sqh ) : std::abs( sqh ) ) ); // AV why abs here?
	vc[4] = cxmake( 0., nsv * ( pvec3 < 0. ? -sqh : sqh ) ); // AV: removed an abs here
	}
	#else
	// Branch A: pt != 0.
	const fptype_v pzpt = pvec3 / ( pp * pt ) * sqh * hel;
	const cxtype_v vcA_3 = cxmake( -pvec1 * pzpt, -(fptype)nsv * pvec2 / pt * sqh );
	const cxtype_v vcA_4 = cxmake( -pvec2 * pzpt, (fptype)nsv * pvec1 / pt * sqh );
	// Branch B: pt == 0.
	const cxtype vcB_3 = cxmake( -(fptype)hel * sqh, 0 );
	const cxtype_v vcB_4 = cxmake( 0, (fptype)nsv * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here
	// Choose between the results from branch A and branch B
	const bool_v mask = ( pt != 0. );
	vc[3] = cxternary( mask, vcA_3, vcB_3 );
	vc[4] = cxternary( mask, vcA_4, vcB_4 );
	#endif
	}
	mgDebug( 1, __FUNCTION__ );
	return;

The fp exception is in line 516, caused by line 497, where pt = 0:

const fptype_sv pt = fpsqrt( ( pvec1 * pvec1 ) + ( pvec2 * pvec2 ) );
pvec1 = @0x7fffffbd98e0: {0, 0, 0, 0}
pvec2 = @0x7fffffbd9900: {0, 0, 0, 0}

[valassi]

Hi @hageboeck thanks a lot!

Ok this looks like "my" bug, something that I introduced. I am a bit surprised because I thought I had tested this extensively including for values where pt=0 (testxxx.cc), but it seems there is an issue.

Definitely this code is ugly, as it makes BOTH the calculation with pT=0 and that with pT>0 and then chooses a posteriori... with the tiny problem that if pt+0 then the calculation assuming pt>0 gives a FP exception

madgraph4gpu/epochX/cudacpp/gg_ttgg.mad/src/HelAmps_sm.h

Line 515 in 603ff90

// Branch A: pt != 0.

#ifndef MGONGPU_CPPSIMD
      if( pt != 0. )
      {
        const fptype pzpt = pvec3 / ( pp * pt ) * sqh * hel;
        vc[3] = cxmake( -pvec1 * pzpt, -nsv * pvec2 / pt * sqh );
        vc[4] = cxmake( -pvec2 * pzpt, nsv * pvec1 / pt * sqh );
      }
      else
      {
        vc[3] = cxmake( -hel * sqh, 0. );
        // NB: Do not use "abs" for floats! It returns an integer with no build warning! Use std::abs!
        //vc[4] = cxmake( 0, nsv * ( pvec3 < 0. ? -std::abs( sqh ) : std::abs( sqh ) ) ); // AV why abs here?
        vc[4] = cxmake( 0., nsv * ( pvec3 < 0. ? -sqh : sqh ) ); // AV: removed an abs here
      }
#else
      // Branch A: pt != 0.
      const fptype_v pzpt = pvec3 / ( pp * pt ) * sqh * hel;
      const cxtype_v vcA_3 = cxmake( -pvec1 * pzpt, -(fptype)nsv * pvec2 / pt * sqh );
      const cxtype_v vcA_4 = cxmake( -pvec2 * pzpt, (fptype)nsv * pvec1 / pt * sqh );
      // Branch B: pt == 0.
      const cxtype vcB_3 = cxmake( -(fptype)hel * sqh, 0 );
      const cxtype_v vcB_4 = cxmake( 0, (fptype)nsv * fpternary( ( pvec3 < 0 ), -sqh, sqh ) ); // AV: removed an abs here
      // Choose between the results from branch A and branch B
      const bool_v mask = ( pt != 0. );
      vc[3] = cxternary( mask, vcA_3, vcB_3 );
      vc[4] = cxternary( mask, vcA_4, vcB_4 );
#endif

There are TWO (or three) things TODO

• first extens the testxxx test and understand why it did not pick this up
• second, fix the issue
• (third, check if there are other places in ixx/oxx/vxx functions that use this trick of masks and that may give issues)

To fix it, one way may be to simply use a loop over ieppV=0...neppV-1 and forget about SIMD.

Otherwise, a truly SIMD solution would require having two branches both of which do not send exceptions. In this case it is probably impossible. Anyway, the performance penalty is probably small: the vxxx functions are caled ONCE even if you have 1000 Feynamn diagrams...

[valassi]

PS About the checking other masks... I guess that my idea whan I implemented this was that the two "branches" were doing different calculations completely, not that if statement was there as a protection against FP exceptions. To be checked...

[valassi]

I remembered a discussion with @hageboeck where he mentioned that fortran enables FPE traps.

I am doing some tests with FP traps enabled in C++. I do see errors in testxxx.

Debugging is complex, things are not straightforward. For instance I get some FPEs in ixxxxx in debug mode, but ixxxxx goes through if not in debug mode...