Remove neppV and fptype_sv from API and decouple neppV from neppM (epochX)

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.cu

@@ -1,7 +1,7 @@
   //--------------------------------------------------------------------------
 
-#ifdef __CUDACC__
-  // Return by reference
+  // Decode momentum AOSOA: compute address of fptype for the given particle, 4-momentum component and event
+  // Return the fptype by reference (equivalent to returning its memory address)
   __device__ inline
   const fptype& pIparIp4Ievt( const fptype* momenta, // input: momenta as AOSOA[npagM][npar][4][neppM]
                               const int ipar,
@@ -10,46 +10,157 @@
   {
     using mgOnGpu::np4;
     using mgOnGpu::npar;
-    const int neppM = mgOnGpu::neppM; // AOSOA layout: constant at compile-time
-    const int ipagM = ievt/neppM; // #eventpage in this iteration
-    const int ieppM = ievt%neppM; // #event in the current eventpage in this iteration
-    //printf( "%f\n", momenta[ipagM*npar*np4*neppM + ipar*np4*neppM + ip4*neppM + ieppM] );
+    constexpr int neppM = mgOnGpu::neppM; // AOSOA layout: constant at compile-time
+    const int ipagM = ievt/neppM; // #event "M-page"
+    const int ieppM = ievt%neppM; // #event in the current event M-page
+    //printf( "%2d %2d %8d %8.3f\n", ipar, ip4, ievt, momenta[ipagM*npar*np4*neppM + ipar*np4*neppM + ip4*neppM + ieppM] );
     return momenta[ipagM*npar*np4*neppM + ipar*np4*neppM + ip4*neppM + ieppM]; // AOSOA[ipagM][ipar][ip4][ieppM]
     //fptype (*momenta)[npar][np4][neppM] = (fptype (*)[npar][np4][neppM]) momenta; // cast to multiD array pointer (AOSOA)
     //return momenta[ipagM][ipar][ip4][ieppM]; // this seems ~1-2% faster in eemumu C++?
   }
-#else
-  // Return by value: it seems a tiny bit faster than returning a reference (both for scalar and vector), not clear why
+
+#ifndef __CUDACC__
+  // Return a SIMD vector of fptype's for neppV events (for the given particle, 4-momentum component and event "V-page")
+  // Return the vector by value: strictly speaking, this is only unavoidable for neppM<neppV
+  // For neppM>=neppV (both being powers of 2), the momentum neppM-AOSOA is reinterpreted in terms of neppV-vectors:
+  // it could also be returned by reference, but no performance degradation is observed when returning by value
   inline
-  fptype_sv pIparIp4Ipag( const fptype_sv* momenta, // input: momenta as AOSOA[npagM][npar][4][neppM]
+  fptype_sv pIparIp4Ipag( const fptype* momenta, // input: momenta as AOSOA[npagM][npar][4][neppM]
                           const int ipar,
                           const int ip4,
-                          const int ipagM )
+                          const int ipagV )
   {
-    /*
-    //#ifndef MGONGPU_CPPSIMD
-    // TEMPORARY during epochX step3! Eventually remove this section (start)
-    // TEMPORARY during epochX step3! THERE IS NO SIMD YET: HENCE ipagM=ievt
-    // NB: this is needed for neppM>1 while neppV==1 (no SIMD) in eemumu.auto
-    // NB: this remains valid for neppM==1 with neppV==1 in eemumu (scalar)
-    return pIparIp4Ievt( momenta, ipar, ip4, ipagM );
-    // TEMPORARY during epochX step3! Eventually remove this section (end)
-    //#else
-    */
-    // NB: this assumes that neppV == neppM!
-    // NB: this is the same as "pIparIp4Ievt( momenta, ipar, ip4, ipagM )" for neppM==1 with neppV==1
-    using mgOnGpu::np4;
-    using mgOnGpu::npar;
-    //printf( "%f\n", momenta[ipagM*npar*np4 + ipar*np4 + ip4] );
-    return momenta[ipagM*npar*np4 + ipar*np4 + ip4]; // AOSOA[ipagM][ipar][ip4][ieppM]
+    const int ievt0 = ipagV*neppV; // virtual event V-page ipagV contains neppV events [ievt0...ievt0+neppV-1]
+#ifdef MGONGPU_CPPSIMD
+    constexpr int neppM = mgOnGpu::neppM; // AOSOA layout: constant at compile-time
+    // Use c++17 "if constexpr": compile-time branching
+    if constexpr ( ( neppM >= neppV ) && ( neppM%neppV == 0 ) )
+    {
+      constexpr bool useReinterpretCastIfPossible = true; // DEFAULT
+      //constexpr bool useReinterpretCastIfPossible = false; // FOR PERFORMANCE TESTS
+      constexpr bool skipAlignmentCheck = true; // DEFAULT (MAY SEGFAULT, NEEDS A SANITY CHECK ELSEWHERE!)
+      //constexpr bool skipAlignmentCheck = false; // SLOWER BUT SAFER
+      if constexpr ( useReinterpretCastIfPossible && skipAlignmentCheck )
+      {
+        //static bool first=true; if( first ){ std::cout << "WARNING! skip alignment check" << std::endl; first=false; } // SLOWS DOWN...
+        // Fastest (4.92E6 in eemumu 512y)
+        // This requires alignment for momenta - segmentation fault otherwise!
+        return *reinterpret_cast<const fptype_sv*>( &( pIparIp4Ievt( momenta, ipar, ip4, ievt0 ) ) );
+      }
+      else if ( useReinterpretCastIfPossible && ( (size_t)(momenta) % mgOnGpu::cppAlign == 0 ) )
+      {
+        //static bool first=true; if( first ){ std::cout << "WARNING! alignment ok, use reinterpret cast" << std::endl; first=false; } // SLOWS DOWN...
+        // A bit (6%) slower (4.62E6 in eemumu 512y) because of the alignment check
+        // This explicitly checks alignment for momenta to avoid segmentation faults
+        return *reinterpret_cast<const fptype_sv*>( &( pIparIp4Ievt( momenta, ipar, ip4, ievt0 ) ) );
+      }
+      else
+      {
+        //static bool first=true; if( first ){ std::cout << "WARNING! AOSOA but no reinterpret cast" << std::endl; first=false; } // SLOWS DOWN...
+        // A bit (2%) slower (4.86E6 in eemumu 512y)
+        // This does not require alignment for momenta, but it requires AOSOA with neppM>=neppV and neppM%neppV==0
+        const fptype* out0 = &( pIparIp4Ievt( momenta, ipar, ip4, ievt0) );
+#if MGONGPU_CPPSIMD == 2
+        return fptype_v{ *( out0 ),
+                         *( out0+1 ) };
+#elif MGONGPU_CPPSIMD == 4
+        return fptype_v{ *( out0 ),
+                         *( out0+1 ),
+                         *( out0+2 ),
+                         *( out0+3 ) };
+#elif MGONGPU_CPPSIMD == 8
+        return fptype_v{ *( out0 ),
+                         *( out0+1 ),
+                         *( out0+2 ),
+                         *( out0+3 ),
+                         *( out0+4 ),
+                         *( out0+5 ),
+                         *( out0+6 ),
+                         *( out0+7 ) };
+#elif MGONGPU_CPPSIMD == 16
+        return fptype_v{ *( out0 ),
+                         *( out0+1 ),
+                         *( out0+2 ),
+                         *( out0+3 ),
+                         *( out0+4 ),
+                         *( out0+5 ),
+                         *( out0+6 ),
+                         *( out0+7 ),
+                         *( out0+8 ),
+                         *( out0+9 ),
+                         *( out0+10 ),
+                         *( out0+11 ),
+                         *( out0+12 ),
+                         *( out0+13 ),
+                         *( out0+14 ),
+                         *( out0+15 ) };
+#else
+#warning Internal error? This code should never be reached!
+        // A bit (2%) slower (4.86E6 in eemumu 512y)
+        // This does not require alignment for momenta, but it requires AOSOA with neppM>=neppV and neppM%neppV==0
+        fptype_v out;
+        for ( int ieppV=0; ieppV<neppV; ieppV++ ) out[ieppV] = *( out0 + ieppV );
+        return out;
+#endif
+      }
+    }
+    else
+    {
+      // Much (20%) slower (4.07E6 in eemumu 512y)
+      // This does not even require AOSOA with neppM>=neppV and neppM%neppV==0 (e.g. can be used with AOS neppM==1)
+#if MGONGPU_CPPSIMD == 2
+      return fptype_v{ pIparIp4Ievt( momenta, ipar, ip4, ievt0 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+1 ) };
+#elif MGONGPU_CPPSIMD == 4
+      return fptype_v{ pIparIp4Ievt( momenta, ipar, ip4, ievt0 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+1 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+2 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+3 ) };
+#elif MGONGPU_CPPSIMD == 8
+      return fptype_v{ pIparIp4Ievt( momenta, ipar, ip4, ievt0 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+1 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+2 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+3 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+4 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+5 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+6 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+7 ) };
+#elif MGONGPU_CPPSIMD == 16
+      return fptype_v{ pIparIp4Ievt( momenta, ipar, ip4, ievt0 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+1 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+2 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+3 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+4 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+5 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+6 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+7 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+8 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+9 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+10 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+11 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+12 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+13 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+14 ),
+                       pIparIp4Ievt( momenta, ipar, ip4, ievt0+15 ) };
+#else
+#warning Internal error? This code should never be reached!
+      // Much much (40%) slower (3.02E6 in eemumu 512y)
+      fptype_v out;
+      for ( int ieppV=0; ieppV<neppV; ieppV++ ) out[ieppV] = pIparIp4Ievt( momenta, ipar, ip4, ievt0+ieppV );
+      return out;
+#endif
+    }
+#else
+    return pIparIp4Ievt( momenta, ipar, ip4, ievt0 );
+#endif
   }
 #endif
 
   //--------------------------------------------------------------------------
 
   // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt]
   __device__
-  void ixxxxx( const fptype_sv* momenta,
+  void ixxxxx( const fptype* momenta,
                const fptype fmass,             // input: fermion mass
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -185,7 +296,7 @@
   // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == +PZ > 0)
   __device__
-  void ipzxxx( const fptype_sv* momenta,
+  void ipzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -234,7 +345,7 @@
   // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == -PZ > 0)
   __device__
-  void imzxxx( const fptype_sv* momenta,
+  void imzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -283,7 +394,7 @@
   // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PT > 0)
   __device__
-  void ixzxxx( const fptype_sv* momenta,
+  void ixzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -344,7 +455,7 @@
 
   // Compute the output wavefunction vc[6] from the input momenta[npar*4*nevt]
   __device__
-  void vxxxxx( const fptype_sv* momenta,
+  void vxxxxx( const fptype* momenta,
                const fptype vmass,             // input: vector boson mass
                const int nhel,                 // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson)
                const int nsv,                  // input: +1 (final) or -1 (initial)
@@ -480,7 +591,7 @@
 
   // Compute the output wavefunction sc[3] from the input momenta[npar*4*nevt]
   __device__
-  void sxxxxx( const fptype_sv* momenta,
+  void sxxxxx( const fptype* momenta,
                const fptype,                   // WARNING: input "smass" unused (missing in Fortran) - scalar boson mass
                const int,                      // WARNING: input "nhel" unused (missing in Fortran) - scalar has no helicity!
                const int nss,                  // input: +1 (final) or -1 (initial)
@@ -521,7 +632,7 @@
 
   // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt]
   __device__
-  void oxxxxx( const fptype_sv* momenta,
+  void oxxxxx( const fptype* momenta,
                const fptype fmass,             // input: fermion mass
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -658,7 +769,7 @@
   // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == +PZ > 0)
   __device__
-  void opzxxx( const fptype_sv* momenta,
+  void opzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -707,7 +818,7 @@
   // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == -PZ > 0)
   __device__
-  void omzxxx( const fptype_sv* momenta,
+  void omzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -759,7 +870,7 @@
   // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PT > 0)
   __device__
-  void oxzxxx( const fptype_sv* momenta,
+  void oxzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)

epochX/cudacpp/CODEGEN/PLUGIN/CUDACPP_SA_OUTPUT/aloha/template_files/gpu/helas.h

@@ -11,7 +11,7 @@
 
   // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt]
   __device__ INLINE
-  void ixxxxx( const fptype_sv* momenta,
+  void ixxxxx( const fptype* momenta,
                const fptype fmass,             // input: fermion mass
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -26,7 +26,7 @@
   // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == +PZ > 0)
   __device__ INLINE
-  void ipzxxx( const fptype_sv* momenta,
+  void ipzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -41,7 +41,7 @@
   // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == -PZ > 0)
   __device__ INLINE
-  void imzxxx( const fptype_sv* momenta,
+  void imzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -56,7 +56,7 @@
   // Compute the output wavefunction fi[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PT > 0)
   __device__ INLINE
-  void ixzxxx( const fptype_sv* momenta,
+  void ixzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -70,7 +70,7 @@
 
   // Compute the output wavefunction vc[6] from the input momenta[npar*4*nevt]
   __device__ INLINE
-  void vxxxxx( const fptype_sv* momenta,
+  void vxxxxx( const fptype* momenta,
                const fptype vmass,             // input: vector boson mass
                const int nhel,                 // input: -1, 0 (only if vmass!=0) or +1 (helicity of vector boson)
                const int nsv,                  // input: +1 (final) or -1 (initial)
@@ -84,7 +84,7 @@
 
   // Compute the output wavefunction sc[3] from the input momenta[npar*4*nevt]
   __device__ INLINE
-  void sxxxxx( const fptype_sv* momenta,
+  void sxxxxx( const fptype* momenta,
                const fptype,                   // WARNING: input "smass" unused (missing in Fortran) - scalar boson mass
                const int,                      // WARNING: input "nhel" unused (missing in Fortran) - scalar has no helicity!
                const int nss,                  // input: +1 (final) or -1 (initial)
@@ -98,7 +98,7 @@
 
   // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt]
   __device__ INLINE
-  void oxxxxx( const fptype_sv* momenta,
+  void oxxxxx( const fptype* momenta,
                const fptype fmass,             // input: fermion mass
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -113,7 +113,7 @@
   // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == +PZ > 0)
   __device__ INLINE
-  void opzxxx( const fptype_sv* momenta,
+  void opzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -128,7 +128,7 @@
   // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt]
   // ASSUMPTIONS: (FMASS == 0) and (PX == PY == 0 and E == -PZ > 0)
   __device__ INLINE
-  void omzxxx( const fptype_sv* momenta,
+  void omzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)
@@ -142,7 +142,7 @@
 
   // Compute the output wavefunction fo[6] from the input momenta[npar*4*nevt]
   __device__ INLINE
-  void oxzxxx( const fptype_sv* momenta,
+  void oxzxxx( const fptype* momenta,
                //const fptype fmass,           // ASSUME fermion mass==0
                const int nhel,                 // input: -1 or +1 (helicity of fermion)
                const int nsf,                  // input: +1 (particle) or -1 (antiparticle)