From b8378226e40eb3a226d8e1890b421189b62e4224 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 11 Jun 2026 15:25:32 -0400 Subject: [PATCH 01/32] Separate the allocation and deallocation of data and terms for amplitude calculation. --- AmpTools/GPUManager/GPUManager.cc | 186 +++++++++++++++++++----------- AmpTools/GPUManager/GPUManager.h | 19 +-- 2 files changed, 128 insertions(+), 77 deletions(-) diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc index 2660dddd..02dfa448 100644 --- a/AmpTools/GPUManager/GPUManager.cc +++ b/AmpTools/GPUManager/GPUManager.cc @@ -70,7 +70,7 @@ GPUManager::GPUManager() m_iNAmps=0; m_iNUserVars=0; - m_iGDoubleEventArrSize=0; + m_iGDoubleDataArrSize=0; m_iDoubleEventArrSize=0; m_iAmpArrSize=0; m_iVArrSize=0; @@ -173,27 +173,71 @@ GPUManager::~GPUManager() // Initialization routines: -void -GPUManager::init( const AmpVecs& a, bool use4Vectors ) +void +GPUManager::initData( const AmpVecs& a, bool use4Vectors ) { - clearAll(); - + clearData(); + // copy over some info from the AmpVecs object for array dimensions m_iNTrueEvents = a.m_iNTrueEvents; m_iNEvents = a.m_iNEvents; m_iNParticles = a.m_iNParticles; + + // the rest of the data are derived: + m_iGDoubleDataArrSize = sizeof(GDouble) * m_iNEvents; + + double totalMemory = 0; + + totalMemory += m_iGDoubleDataArrSize; + if( use4Vectors ) totalMemory += 4 * m_iNParticles * m_iGDoubleDataArrSize; + totalMemory += m_iNUserVars * m_iGDoubleDataArrSize; + totalMemory += m_iNParticles * sizeof( int ); + + totalMemory /= (1024*1024); + + report( INFO, kModule ) << "Attempting to allocate " << (int)totalMemory << " MB of global GPU memory." << endl; + + gpuErrChk( cudaMalloc( (void**)&m_pfDevWeights , m_iGDoubleDataArrSize ) ) ; + + // allocate device memory needed for amplitude calculations + if( use4Vectors ) gpuErrChk( cudaMalloc( (void**)&m_pfDevData , + 4 * m_iNParticles * m_iGDoubleDataArrSize ) ) ; + gpuErrChk( cudaMalloc( (void**)&m_pfDevUserVars, m_iNUserVars * m_iGDoubleDataArrSize ) ) ; + + report( INFO, kModule ) << "GPU memory allocated for storage of " + << m_iNEvents << " events (" << m_iNTrueEvents << " actual events)." + << endl; + + //CUDA Dims + calcCUDADims(); +} + + +void +GPUManager::initTerms( const AmpVecs& a, unsigned long chunkSize ) +{ + clearTerms(); + + // chunkSize must be power of 2 when provided + assert( chunkSize == 0 || ( (chunkSize & (chunkSize - 1)) == 0 ) ); + + // use all events or the specified size of a chunk of events + long int iNEvents = ( chunkSize == 0 ? m_iNevents : chunkSize ); + m_iNAmps = a.m_iNTerms; m_iNUserVars = a.m_userVarsPerEvent; - // the rest of the data are derived: - m_iGDoubleEventArrSize = sizeof(GDouble) * m_iNEvents; - - // double precision - m_iDoubleTrueEventArrSize = sizeof(double) * m_iNTrueEvents; - m_iDoubleEventArrSize = sizeof(double) * m_iNEvents; + // double precision intensity calculation + m_iDoubleTrueIntenArrSize = sizeof(double) * m_iNTrueEvents; + m_iDoubleIntenArrSize = sizeof(double) * m_iNEvents; + + // GDouble precision for factors and amplitudes + // ... and use the chunkSize for these (iNEvents) when it is provided + m_iGDoubleFactArrSize = sizeof(GDouble) * a.m_maxFactPerEvent * iNEvents; + // size needed to store amplitudes for each event - m_iAmpArrSize = 2 * sizeof(GDouble) * m_iNEvents * m_iNAmps; + m_iAmpArrSize = 2 * sizeof(GDouble) * iNEvents * m_iNAmps; // size of upper half of ViVj* matrix m_iVArrSize = 2 * sizeof(GDouble) * m_iNAmps * ( m_iNAmps + 1 ) / 2; @@ -212,40 +256,29 @@ GPUManager::init( const AmpVecs& a, bool use4Vectors ) totalMemory += m_iVArrSize; totalMemory += m_iNICalcSize; - totalMemory += m_iGDoubleEventArrSize; - totalMemory += 2*m_iDoubleEventArrSize; - if( use4Vectors ) totalMemory += 4 * m_iNParticles * m_iGDoubleEventArrSize; - totalMemory += m_iNUserVars * m_iGDoubleEventArrSize; - totalMemory += m_iNParticles * sizeof( int ); - totalMemory += a.m_maxFactPerEvent * m_iGDoubleEventArrSize; + totalMemory += 2*m_iDoubleIntenArrSize; + totalMemory += m_iGDoubleFactArrSize; totalMemory += m_iAmpArrSize; + totalMemory += m_iNParticles * sizeof( int ); totalMemory /= (1024*1024); - report( INFO, kModule ) << "Attempting to allocate " << (int)totalMemory << " MB of global GPU memory." << endl; + report( INFO, kModule ) << "Attempting to allocate " << (int)totalMemory + << " MB of global GPU memory." << endl; // device memory needed for intensity or integral calculation and sum - gpuErrChk( cudaMalloc( (void**)&m_pfDevVVStar , m_iVArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pdDevNICalc , m_iNICalcSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pfDevWeights , m_iGDoubleEventArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pdDevRes , m_iDoubleEventArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pdDevREDUCE , m_iDoubleEventArrSize ) ) ; - - // allocate device memory needed for amplitude calculations - if( use4Vectors ) gpuErrChk( cudaMalloc( (void**)&m_pfDevData , - 4 * m_iNParticles * m_iGDoubleEventArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pfDevUserVars, m_iNUserVars * m_iGDoubleEventArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_piDevPerm , m_iNParticles * sizeof( int ) ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pcDevCalcAmp , - a.m_maxFactPerEvent * m_iGDoubleEventArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pfDevAmps , m_iAmpArrSize ) ) ; - - report( INFO, kModule ) << "GPU memory allocated for " << m_iNAmps << " amplitudes and " - << m_iNEvents << " events (" << m_iNTrueEvents << " actual events)." + gpuErrChk( cudaMalloc( (void**)&m_pfDevVVStar , m_iVArrSize ) ) ; + gpuErrChk( cudaMalloc( (void**)&m_pdDevNICalc , m_iNICalcSize ) ) ; + gpuErrChk( cudaMalloc( (void**)&m_pdDevRes , m_iDoubleIntenArrSize ) ) ; + gpuErrChk( cudaMalloc( (void**)&m_pdDevREDUCE , m_iDoubleIntenArrSize ) ) ; + gpuErrChk( cudaMalloc( (void**)&m_pcDevCalcAmp , m_iGDoubleFactArrSize ) ) ; + gpuErrChk( cudaMalloc( (void**)&m_pfDevAmps , m_iAmpArrSize ) ) ; + + gpuErrChk( cudaMalloc( (void**)&m_piDevPerm , m_iNParticles * sizeof( int ) ) ) ; + + report( INFO, kModule ) << "GPU memory allocated for " << m_iNAmps << " amplitudes for " + << iNEvents << " events (" << m_iNTrueEvents << " total events)." << endl; - - //CUDA Dims - calcCUDADims(); } @@ -259,7 +292,7 @@ GPUManager::copyDataToGPU( const AmpVecs& a, bool use4Vectors ) // copy the weights to the GPU gpuErrChk( cudaMemcpy( m_pfDevWeights, a.m_pdWeights, - m_iGDoubleEventArrSize, cudaMemcpyHostToDevice ) ); + m_iGDoubleDataArrSize, cudaMemcpyHostToDevice ) ); // we only need to copy the four vectors to the GPU if the // amplitude evaulation kernels need them -- this is done by @@ -304,7 +337,7 @@ GPUManager::copyDataToGPU( const AmpVecs& a, bool use4Vectors ) // copy the data into the device gpuErrChk( cudaMemcpy( m_pfDevData, tmpStorage, - 4 * m_iNParticles * m_iGDoubleEventArrSize, + 4 * m_iNParticles * m_iGDoubleDataArrSize, cudaMemcpyHostToDevice ) ); delete[] tmpStorage; @@ -327,8 +360,8 @@ GPUManager::copyUserVarsToGPU( const AmpVecs& a ) // copy the data into the device gpuErrChk( cudaMemcpy( m_pfDevUserVars, a.m_pdUserVars, - m_iNUserVars * m_iGDoubleEventArrSize, - cudaMemcpyHostToDevice ) ); + m_iNUserVars * m_iGDoubleDataArrSize, + cudaMemcpyHostToDevice ) ); #ifdef SCOREP SCOREP_USER_REGION_END( copyUserVarsToGPU ) #endif @@ -352,7 +385,7 @@ GPUManager::copyAmpsFromGPU( AmpVecs& a ) cudaMemcpyDeviceToHost ) ); gpuErrChk( cudaMemcpy( a.m_pdAmpFactors, m_pcDevCalcAmp, - a.m_maxFactPerEvent * m_iGDoubleEventArrSize, + m_iGDoubleFactArrSize, cudaMemcpyDeviceToHost ) ); #ifdef SCOREP SCOREP_USER_REGION_END( copyAmpsFromGPU ) @@ -360,9 +393,9 @@ GPUManager::copyAmpsFromGPU( AmpVecs& a ) } void -GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long long offset, +GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long offset, const vector< vector< int > >* pvPermutations, - unsigned long long iUserVarsOffset ) + unsigned long iUserVarsOffset ) { #ifdef SCOREP SCOREP_USER_REGION_DEFINE( calcAmplitudeAll_gpuMgr ) @@ -379,8 +412,8 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long long offset, // if this is not true, AmplitudeManager hasn't been setup properly assert( permItr->size() == m_iNParticles ); - unsigned long long udLocalOffset = 0; - unsigned long long permOffset = 0; + unsigned long udLocalOffset = 0; + unsigned long permOffset = 0; for( ; permItr != pvPermutations->end(); ++permItr ){ // copy the permutation to global memory @@ -541,7 +574,7 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef, void GPUManager::calcIntegrals( double* result, int nElements, - const vector& iIndex, + const vector& iIndex, const vector& jIndex ){ unsigned int resultSize = 2*sizeof(double)*nElements; @@ -594,15 +627,47 @@ GPUManager::calcIntegrals( double* result, int nElements, // Methods to clear memory: -void GPUManager::clearAll() -{ +void +GPUManager::clearAll() { + + clearData(); + clearTerms(); +} + +void +GPUManager::clearData(){ + m_iNParticles=0; m_iNEvents=0; m_iNTrueEvents=0; - m_iGDoubleEventArrSize=0; + + m_iGDoubleDataArrSize=0; m_iDoubleEventArrSize=0; m_iDoubleTrueEventArrSize=0; + if(m_pfDevData) + cudaFree(m_pfDevData); + m_pfDevData=0; + + if(m_pfDevUserVars) + cudaFree(m_pfDevUserVars); + m_pfDevUserVars=0; + + if(m_pfDevWeights) + cudaFree(m_pfDevWeights); + m_pfDevWeights=0; + + + //CUDA Thread and Grid sizes + m_iDimGridX=0; + m_iDimGridY=0; + m_iDimThreadX=0; + m_iDimThreadY=0; +} + +void +GPUManager::clearTerms(){ + m_iNAmps=0; m_iAmpArrSize=0; m_iVArrSize=0; @@ -621,14 +686,6 @@ void GPUManager::clearAll() //Device Memory - if(m_pfDevData) - cudaFree(m_pfDevData); - m_pfDevData=0; - - if(m_pfDevUserVars) - cudaFree(m_pfDevUserVars); - m_pfDevUserVars=0; - if(m_pcDevCalcAmp) cudaFree(m_pcDevCalcAmp); m_pcDevCalcAmp=0; @@ -649,9 +706,6 @@ void GPUManager::clearAll() cudaFree(m_pdDevNICalc); m_pdDevNICalc=0; - if(m_pfDevWeights) - cudaFree(m_pfDevWeights); - m_pfDevWeights=0; if(m_pdDevRes) cudaFree(m_pdDevRes); @@ -660,12 +714,6 @@ void GPUManager::clearAll() if(m_pdDevREDUCE) cudaFree(m_pdDevREDUCE); m_pdDevREDUCE=0; - - //CUDA Thread and Grid sizes - m_iDimGridX=0; - m_iDimGridY=0; - m_iDimThreadX=0; - m_iDimThreadY=0; } // Internal utilities: diff --git a/AmpTools/GPUManager/GPUManager.h b/AmpTools/GPUManager/GPUManager.h index 82a8f2a1..ea5088bd 100644 --- a/AmpTools/GPUManager/GPUManager.h +++ b/AmpTools/GPUManager/GPUManager.h @@ -78,19 +78,22 @@ class GPUManager ~GPUManager(); void clearAll(); - void clearAmpCalc(); - void clearLikeCalc(); + void clearData(); + void clearTerms(); void init( const AmpVecs& a, bool use4Vectors = true ); + void initData( const AmpVecs& a, bool use4Vectors = true ); + void initTerms( const AmpVecs& a, unsigned long chunkSize = 0 ); + // Interface Utils // First Amplitude calculation interface void copyDataToGPU( const AmpVecs& a, bool use4Vectors = true ); void copyUserVarsToGPU( const AmpVecs& a ); - void calcAmplitudeAll( const Amplitude* amp, unsigned long long offset, + void calcAmplitudeAll( const Amplitude* amp, unsigned long offset, const vector< vector< int > >* pvPermutations, - unsigned long long userVarsOffset ); + unsigned long userVarsOffset ); void assembleTerms( int iAmpInd, int nFact, int nPerm ); @@ -124,10 +127,10 @@ class GPUManager unsigned int m_iNUserVars; // array sizes - unsigned long long m_iGDoubleEventArrSize; - unsigned long long m_iDoubleEventArrSize; - unsigned long long m_iDoubleTrueEventArrSize; - unsigned long long m_iAmpArrSize; + unsigned long m_iGDoubleDataArrSize; + unsigned long m_iDoubleEventArrSize; + unsigned long m_iDoubleTrueEventArrSize; + unsigned long m_iAmpArrSize; unsigned int m_iVArrSize; unsigned int m_iNICalcSize; From 4d5a86e4035e4c2aee32681373fab8e82432a5af Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 11 Jun 2026 15:26:45 -0400 Subject: [PATCH 02/32] revisions to allow allocating a chunk of memory for amplitude calclations; in many places reduce size of event indices from long or long long to int, which should be sufficient --- AmpTools/IUAmpTools/AmpVecs.cc | 34 ++++++++++++++++++++++------------ AmpTools/IUAmpTools/AmpVecs.h | 17 +++++++++++------ 2 files changed, 33 insertions(+), 18 deletions(-) diff --git a/AmpTools/IUAmpTools/AmpVecs.cc b/AmpTools/IUAmpTools/AmpVecs.cc index 28423bf4..de4fa15d 100644 --- a/AmpTools/IUAmpTools/AmpVecs.cc +++ b/AmpTools/IUAmpTools/AmpVecs.cc @@ -198,8 +198,8 @@ AmpVecs::clearFourVecs(){ } void -AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned long long iEvent, - unsigned long long iNTrueEvents ){ +AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent, + unsigned int iNTrueEvents ){ // allocate memory and set variables // if this is the first call to this method @@ -272,7 +272,7 @@ AmpVecs::loadData( DataReader* pDataReader ){ // Loop over events and load each one individually Kinematics* pKinematics; - for(unsigned long long iEvent = 0; iEvent < m_iNTrueEvents; iEvent++){ + for(unsigned int iEvent = 0; iEvent < m_iNTrueEvents; iEvent++){ pKinematics = pDataReader->getEvent(); loadEvent(pKinematics, iEvent, m_iNTrueEvents ); @@ -292,13 +292,20 @@ AmpVecs::loadData( DataReader* pDataReader ){ // Fill any remaining space in the data array with the last event's kinematics - for (unsigned long long iEvent = m_iNTrueEvents; iEvent < m_iNEvents; iEvent++){ + for (unsigned int iEvent = m_iNTrueEvents; iEvent < m_iNEvents; iEvent++){ loadEvent(pKinematics, iEvent, m_iNTrueEvents ); } if( m_iNTrueEvents ) delete pKinematics; +#ifdef GPU_ACCELERATION + + m_gpuMan.initData( *this, !intenMan.needsUserVarsOnly() ); + m_gpuMan.copyDataToGPU( *this, !intenMan.needsUserVarsOnly() ); + +#endif + m_termsValid = false; m_integralValid = false; m_dataLoaded = true; @@ -307,8 +314,10 @@ AmpVecs::loadData( DataReader* pDataReader ){ void -AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity ){ +AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity, unsigned int chunkSize ){ + unsigned int iNEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize ); + m_iNTerms = intenMan.getTermNames().size(); m_maxFactPerEvent = intenMan.maxFactorStoragePerEvent(); m_userVarsPerEvent = intenMan.userVarsPerEvent(); @@ -346,13 +355,12 @@ AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity ) #ifndef GPU_ACCELERATION - m_pdAmps = new GDouble[m_iNEvents * intenMan.termStoragePerEvent()]; - m_pdAmpFactors = new GDouble[m_iNEvents * m_maxFactPerEvent]; + m_pdAmps = new GDouble[iNEvents * intenMan.termStoragePerEvent()]; + m_pdAmpFactors = new GDouble[iNEvents * m_maxFactPerEvent]; #else - m_gpuMan.init( *this, !intenMan.needsUserVarsOnly() ); - m_gpuMan.copyDataToGPU( *this, !intenMan.needsUserVarsOnly() ); + m_gpuMan.initTerms( *this, iNEvents ); #endif // GPU_ACCELERATION @@ -362,14 +370,16 @@ AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity ) #ifdef GPU_ACCELERATION void -AmpVecs::allocateCPUAmpStorage( const IntensityManager& intenMan ){ +AmpVecs::allocateCPUAmpStorage( const IntensityManager& intenMan, unsigned int chunkSize ){ + unsigned int iNEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize ); + // we should start with unallocated memory assert( m_pdAmps == NULL && m_pdAmpFactors == NULL ); // allocate as "pinned memory" for fast CPU<->GPU memcopies - cudaMallocHost( (void**)&m_pdAmps, m_iNEvents * intenMan.termStoragePerEvent() * sizeof(GDouble) ); - cudaMallocHost( (void**)&m_pdAmpFactors, m_iNEvents * m_maxFactPerEvent * sizeof(GDouble)); + cudaMallocHost( (void**)&m_pdAmps, iNEvents * intenMan.termStoragePerEvent() * sizeof(GDouble) ); + cudaMallocHost( (void**)&m_pdAmpFactors, iNEvents * m_maxFactPerEvent * sizeof(GDouble)); cudaError_t cudaErr = cudaGetLastError(); if( cudaErr != cudaSuccess ){ diff --git a/AmpTools/IUAmpTools/AmpVecs.h b/AmpTools/IUAmpTools/AmpVecs.h index c2f61a87..646ad6a7 100644 --- a/AmpTools/IUAmpTools/AmpVecs.h +++ b/AmpTools/IUAmpTools/AmpVecs.h @@ -71,14 +71,14 @@ struct AmpVecs /** * An integer that stores the number of events. */ - unsigned long m_iNEvents; + unsigned int m_iNEvents; /** * An integer that stores the true number of events. For GPU calculations * it is necessary to pad iNEvents up to the next power of 2. This integer * stores the actual number of unique events. */ - unsigned long m_iNTrueEvents; + unsigned int m_iNTrueEvents; /** * A double that stores the absolute value of the sum of the weights. (For @@ -182,7 +182,7 @@ struct AmpVecs * utilized by the AmplitudeManager, but the values are tied * to the data set so it resides in the AmpVecs struct. */ - map< string, unsigned long long > m_userVarsOffset; + map< string, unsigned int > m_userVarsOffset; /** * These booleans track features of the set of weights are are adjusted @@ -226,9 +226,14 @@ struct AmpVecs * * \param[in] bAllocIntensity if set to true this will allocate space for * the intensity calculation + * + * \param[in] chunkSize is useful if one wants to compute amplitude + * integrals in chunks rather the entire data set -- this functionality allows + * for reduced memory consumption in some applications */ void allocateTerms( const IntensityManager& intenMan, - bool bAllocIntensity = false ); + bool bAllocIntensity = false, + unsigned int chunkSize = 0 ); #ifdef GPU_ACCELERATION /** @@ -271,8 +276,8 @@ struct AmpVecs * * \see loadData */ - void loadEvent( const Kinematics* pKinematics, unsigned long long iEvent = 0, - unsigned long long iNTrueEvents = 1 ); + void loadEvent( const Kinematics* pKinematics, unsigned int iEvent = 0, + unsigned int iNTrueEvents = 1 ); /** * A helper routine to get an event i from the array of data and weights. From 2210e0bcd4ef2fa55d094fc1c05706aacd7efe56 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 11 Jun 2026 15:27:29 -0400 Subject: [PATCH 03/32] revise so that calAmplitudeAll can calculate amplitudes for a chunk of events rather than all of them --- AmpTools/IUAmpTools/Amplitude.cc | 23 ++++++++++++++--------- AmpTools/IUAmpTools/Amplitude.h | 5 +++-- 2 files changed, 17 insertions(+), 11 deletions(-) diff --git a/AmpTools/IUAmpTools/Amplitude.cc b/AmpTools/IUAmpTools/Amplitude.cc index 053583e8..ca52d6c0 100644 --- a/AmpTools/IUAmpTools/Amplitude.cc +++ b/AmpTools/IUAmpTools/Amplitude.cc @@ -132,9 +132,10 @@ SCOREP_USER_REGION_END( calcUserVarsAll ) } void -Amplitude::calcAmplitudeAll( GDouble* pdData, GDouble* pdAmps, int iNEvents, - const vector< vector< int > >* pvPermutations, - GDouble* pdUserVars ) const +Amplitude::calcAmplitudeAll( GDouble* pdData, GDouble* pdAmps, unsigned int iNEvents, + const vector< vector< int > >* pvPermutations, + GDouble* pdUserVars, unsigned int startEvent, + unsigned int chunkSize ) const { #ifdef SCOREP @@ -154,19 +155,23 @@ SCOREP_USER_REGION_BEGIN( calcAmplitudeAll, "calcAmplitudeAll", SCOREP_USER_REGI GDouble** pKin = new GDouble*[iNParticles]; - int i, j, iEvent; + // this is the size of the chunk in the pdAmps array that we + // are going to compute amplitdues for + unsigned int nEvents = ( chunkSize == 0 ? iNEvents : chunkSize ); + + int i, j, iEvent; for( iPermutation = 0; iPermutation < iNPermutations; iPermutation++ ){ m_currentPermutation = (*pvPermutations)[iPermutation]; - for( iEvent=0; iEvent >* pvPermutations, - GDouble* pdUserVars = 0 ) const; + GDouble* pdUserVars = 0, unsigned int startEvent = 0, + unsigned int chunkSize = 0 ) const; /** From fd9118b3d18c3c43345aabca7404c39dd797ce85 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 11 Jun 2026 15:28:56 -0400 Subject: [PATCH 04/32] first attempt at a revision to allow calcIntegrals to be done in chunks instead of for the entire data set at once; this is intended to reduce memory consumption --- AmpTools/IUAmpTools/AmplitudeManager.cc | 276 ++++++++++++++---------- AmpTools/IUAmpTools/AmplitudeManager.h | 39 +++- AmpTools/IUAmpTools/IntensityManager.h | 3 +- 3 files changed, 200 insertions(+), 118 deletions(-) diff --git a/AmpTools/IUAmpTools/AmplitudeManager.cc b/AmpTools/IUAmpTools/AmplitudeManager.cc index 69e796d7..72229d97 100644 --- a/AmpTools/IUAmpTools/AmplitudeManager.cc +++ b/AmpTools/IUAmpTools/AmplitudeManager.cc @@ -273,7 +273,7 @@ SCOREP_USER_REGION_BEGIN( calcUserVars, "calcUserVars", SCOREP_USER_REGION_TYPE_ int iNAmps = ampNames.size(); int iAmpIndex; - unsigned long long iUserVarsOffset = 0; + unsigned int iUserVarsOffset = 0; for( iAmpIndex = 0; iAmpIndex < iNAmps; iAmpIndex++ ) { @@ -300,7 +300,7 @@ SCOREP_USER_REGION_BEGIN( calcUserVars, "calcUserVars", SCOREP_USER_REGION_TYPE_ int iNData = iNVars * a.m_iNEvents * iNPerms; // we will set this based on the algorithm below - unsigned long long thisOffset = 0; + unsigned int thisOffset = 0; if( pCurrAmp->areUserVarsStatic() ){ @@ -309,7 +309,7 @@ SCOREP_USER_REGION_BEGIN( calcUserVars, "calcUserVars", SCOREP_USER_REGION_TYPE_ // object and see if there is one associated with this // amplitude name - map< string, unsigned long long >::const_iterator offsetItr = + map< string, unsigned int >::const_iterator offsetItr = a.m_userVarsOffset.find( pCurrAmp->name() ); if( offsetItr == a.m_userVarsOffset.end() ){ @@ -337,7 +337,7 @@ SCOREP_USER_REGION_BEGIN( calcUserVars, "calcUserVars", SCOREP_USER_REGION_TYPE_ // the variables are not static, repeat the algorithm // above but search based on identifier of the amplitude - map< string, unsigned long long >::const_iterator offsetItr = + map< string, unsigned int >::const_iterator offsetItr = a.m_userVarsOffset.find( pCurrAmp->identifier() ); if( offsetItr == a.m_userVarsOffset.end() ){ @@ -384,9 +384,9 @@ SCOREP_USER_REGION_BEGIN( calcUserVars, "calcUserVars", SCOREP_USER_REGION_TYPE_ for( int iEvt = 0; iEvt < a.m_iNEvents; ++iEvt ){ for( int iVar = 0; iVar < iNVars; ++iVar ){ - unsigned long long cpuIndex = + unsigned int cpuIndex = thisOffset + iPerm*a.m_iNEvents*iNVars + iEvt*iNVars + iVar; - unsigned long long gpuIndex = + unsigned int gpuIndex = iPerm*a.m_iNEvents*iNVars + iVar*a.m_iNEvents + iEvt; tmpVarStorage[gpuIndex] = a.m_pdUserVars[cpuIndex]; @@ -417,21 +417,37 @@ SCOREP_USER_REGION_END( calcUserVars ) return; } - vector AmplitudeManager::calcTerms( AmpVecs& a ) const { + // this function satisfies the requirement of the base class -- + // in cases where the user is not calculating a chunk of terms use the + // specialized function below with default arguments + + vector result = calcTerms( a, 0, 0 ); + + return result; +} + +vector +AmplitudeManager::calcTerms( AmpVecs& a, unsigned int startEvent, unsigned int chunkSize ) const +{ + #ifdef SCOREP SCOREP_USER_REGION_DEFINE( calcTerms ) SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON ) #endif + + unsigned int nEvents = ( chunkSize == 0 ? a.m_iNEvents : chunkSize ); report( DEBUG, kModule ) << "Calculating terms... termsValid = " << a.m_termsValid << endl; // on the first pass through this data set be sure to calculate // the user data first, if needed, before doing term calculations - if( !a.m_termsValid && a.m_userVarsPerEvent > 0 ){ + // the last criteria will make sure that userData is only computed + // once for the first chunk in sequential calls of calcTerms + if( !a.m_termsValid && a.m_userVarsPerEvent > 0 && startEvent == 0 ){ calcUserVars( a ); if( m_needsUserVarsOnly && !m_forceUserVarRecalculation @@ -446,7 +462,7 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON #ifndef GPU_ACCELERATION assert( a.m_pdAmps && a.m_pdAmpFactors); #endif - + int iAmpIndex; for( iAmpIndex = 0; iAmpIndex < iNAmps; iAmpIndex++ ) { @@ -508,15 +524,15 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON // calculate all the factors that make up an amplitude for // for all events serially on CPU or in parallel on GPU - unsigned long long iLocalOffset = 0; + unsigned int iLocalOffset = 0; for( iFactor=0; iFactor < iNFactors; - iFactor++, iLocalOffset += 2 * a.m_iNEvents * iNPermutations ){ + iFactor++, iLocalOffset += 2 * nEvents * iNPermutations ){ pCurrAmp = vAmps.at( iFactor ); // if we have static user data, look up the location in the data array // if not, then look up by identifier - unsigned long long uOffset = + unsigned int uOffset = ( pCurrAmp->areUserVarsStatic() ? a.m_userVarsOffset[pCurrAmp->name()] : a.m_userVarsOffset[pCurrAmp->identifier()] ); @@ -526,7 +542,8 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON calcAmplitudeAll( a.m_pdData, a.m_pdAmpFactors + iLocalOffset, a.m_iNEvents, &vvPermuations, - a.m_pdUserVars + uOffset ); + a.m_pdUserVars + uOffset, + startEvent, nEvents ); #else a.m_gpuMan.calcAmplitudeAll( pCurrAmp, iLocalOffset, &vvPermuations, @@ -543,29 +560,28 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON GDouble dSymmFactor = 1.0f/sqrt( iNPermutations ); GDouble dAmpFacRe, dAmpFacIm, dTRe, dTIm; int iEvent, iPerm; - unsigned long long iOffsetA, iOffsetP, iOffsetF; + unsigned int iOffsetA, iOffsetP, iOffsetF; // re-ordering of data will be useful to not fall out of (CPU) memory cache!!! // zeroing out the entire range - memset( (void*)( a.m_pdAmps + 2 * a.m_iNEvents * iAmpIndex ), 0, - 2 * a.m_iNEvents * sizeof(GDouble) ); + memset( (void*)( a.m_pdAmps + 2 * nEvents * iAmpIndex ), 0, + 2 * nEvents * sizeof(GDouble) ); - // only sum over the true events from data and skip paddings - for( iEvent=0; iEvent < a.m_iNTrueEvents; iEvent++ ) + for( iEvent=0; iEvent < nEvents; iEvent++ ) { - iOffsetA = 2 * a.m_iNEvents * iAmpIndex + 2 * iEvent; + iOffsetA = 2 * nEvents * iAmpIndex + 2 * iEvent; for( iPerm = 0; iPerm < iNPermutations; iPerm++ ) { - iOffsetP = 2 * a.m_iNEvents * iPerm + 2 * iEvent; + iOffsetP = 2 * nEvents * iPerm + 2 * iEvent; dAmpFacRe = a.m_pdAmpFactors[iOffsetP]; dAmpFacIm = a.m_pdAmpFactors[iOffsetP+1]; for( iFactor = 1; iFactor < iNFactors; iFactor++ ) { - iOffsetF = iOffsetP + 2 * a.m_iNEvents * iNPermutations * iFactor; + iOffsetF = iOffsetP + 2 * nEvents * iNPermutations * iFactor; dTRe = dAmpFacRe; dTIm = dAmpFacIm; @@ -592,7 +608,7 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON #endif } - + a.m_termsValid = true; // finally make a loop and record the iteration of parameters @@ -796,7 +812,8 @@ SCOREP_USER_REGION_END( calcSumLogIntensity ) void -AmplitudeManager::calcIntegrals( AmpVecs& a, int iNGenEvents ) const +AmplitudeManager::calcIntegrals( AmpVecs& a, unsigned int iNGenEvents, + unsigned int chunkSize ) const { #ifdef SCOREP SCOREP_USER_REGION_DEFINE( calcIntegralsA ) @@ -808,32 +825,14 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T double* integralMatrix = a.m_pdIntegralMatrix; - // amp -> amp* -> value assert( iNGenEvents ); - // this returns a vector indicating which terms have changed - const vector& termChanged = calcTerms( a ); - -#ifdef SCOREP - SCOREP_USER_REGION_END( calcIntegralsA ) -#endif - bool anyTermChanged = - ( find( termChanged.begin(), termChanged.end(), true ) != termChanged.end() ); - - // if nothing changed and it isn't the first pass, return - if( !anyTermChanged && a.m_integralValid ) return; - -#ifdef SCOREP - SCOREP_USER_REGION_BEGIN( calcIntegralsB, "calcIntegralsB", SCOREP_USER_REGION_TYPE_COMMON ) -#endif + // figure out the number of chunks needed to go through the data + unsigned int nChunk = ( chunkSize == 0 ? 1 : iNGenEvents / chunkSize ); + if( ( chunkSize !=0 ) && ( iNGenEvents % chunkSize != 0 ) ) ++nChunk; int iNAmps = a.m_iNTerms; - // this must match the dimension of the coherence matrix or - // else there is a coding mistake somewhere - assert( iNAmps == m_sumCoherently.size() ); - assert( iNAmps == termChanged.size() ); - int maxNIElements = uniqueNIElements(); // use these vectors to store the amp indices of unique terms @@ -843,85 +842,138 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T // this stores the real and imaginary parts of the term calculations vector result(maxNIElements*2); - - // now figure out which elements of the NI matrix need to be computed - int i, j; - int nCompute = 0; - for( i = 0; i < iNAmps;i++ ) { - for( j = 0; j <= i; j++ ) { - - // if the two amplitudes do not interfere then set the integral to zero - // and move on to the next pair of amplitudes - if( !m_sumCoherently[i][j] ){ + + // a variable to track how many elements need computing + int nCompute; + + for( unsigned int iChunk = 0; iChunk < nChunk; ++iChunk ){ - integralMatrix[2*i*iNAmps+2*j] = 0; - integralMatrix[2*i*iNAmps+2*j+1] = 0; - - // must also set the complex conjugate to zero because complex - // conjugation is only done for terms that are actually computed - integralMatrix[2*j*iNAmps+2*i] = 0; - integralMatrix[2*j*iNAmps+2*i+1] = 0; - - continue; - } + unsigned int startEvent = iChunk * chunkSize; + unsigned int nEvents = iNGenEvents; + // for chuncked calculation: + // the number of events to process is the chunk size unless + // it is the last chunk then it is the remainder + if( chunkSize != 0 ){ + nEvents = chunkSize; + if( iChunk - nChunk == 1 ) nEvents = iNGenEvents % chunkSize; + } + + // this returns a vector indicating which terms have changed + const vector& termChanged = calcTerms( a, startEvent, chunkSize ); + + // if we are doing a "chunked" calculation, we should invalidate the + // termsValid boolean since the terms array will only hold a fraction + // of the terms + if( chunkSize != 0 ) a.m_termsValid = false; + +#ifdef SCOREP + SCOREP_USER_REGION_END( calcIntegralsA ) +#endif + + bool anyTermChanged = + ( find( termChanged.begin(), termChanged.end(), true ) != termChanged.end() ); + + // if nothing changed and it isn't the first pass, return + if( !anyTermChanged && a.m_integralValid ) return; + +#ifdef SCOREP + SCOREP_USER_REGION_BEGIN( calcIntegralsB, "calcIntegralsB", SCOREP_USER_REGION_TYPE_COMMON ) +#endif + + // this must match the dimension of the coherence matrix or + // else there is a coding mistake somewhere + assert( iNAmps == m_sumCoherently.size() ); + + // we need a bool for every amp or else there is a mistake + assert( iNAmps == termChanged.size() ); + + // now figure out which elements of the NI matrix need to be computed + nCompute = 0; + for( int i = 0; i < iNAmps;i++ ) { + for( int j = 0; j <= i; j++ ) { - if( a.m_integralValid && !termChanged.at(i) && !termChanged.at(j) ){ - - // if the amplitude hasn't changed and it isn't the first pass - // through these data, then its integral can't change - - continue; - } - else{ + // if the two amplitudes do not interfere then set the integral to zero + // and move on to the next pair of amplitudes + if( !m_sumCoherently[i][j] ){ + + integralMatrix[2*i*iNAmps+2*j] = 0; + integralMatrix[2*i*iNAmps+2*j+1] = 0; + + // must also set the complex conjugate to zero because complex + // conjugation is only done for terms that are actually computed + integralMatrix[2*j*iNAmps+2*i] = 0; + integralMatrix[2*j*iNAmps+2*i+1] = 0; + + continue; + } - report( DEBUG, kModule ) << "Requesting recomputation of NI term ( " - << i << ", " << j << " )" << endl; + if( a.m_integralValid && !termChanged.at(i) && !termChanged.at(j) ){ + + // if the amplitude hasn't changed and it isn't the first pass + // through these data, then its integral can't change + + continue; + } + else{ + + report( DEBUG, kModule ) << "Requesting recomputation of NI term ( " + << i << ", " << j << " )" << endl; + + // if we get to this point, we need to compute a value + // for the integral matrix + + iIndex[nCompute] = i; + jIndex[nCompute] = j; + + ++nCompute; + } + } + } + +#ifndef GPU_ACCELERATION + for( int iTerm = 0; iTerm < nCompute; ++iTerm ){ + + int i = iIndex[iTerm]; + int j = jIndex[iTerm]; + + for( int iEvent = 0; iEvent < nEvents; iEvent++ ) { - // if we get to this point, we need to compute a value - // for the integral matrix + //AiAj* + result[2*iTerm] += a.m_pdWeights[iEvent+startEvent] * + ( a.m_pdAmps[2*nEvents*i+2*iEvent] * + a.m_pdAmps[2*nEvents*j+2*iEvent] + + a.m_pdAmps[2*nEvents*i+2*iEvent+1] * + a.m_pdAmps[2*nEvents*j+2*iEvent+1] ); - iIndex[nCompute] = i; - jIndex[nCompute] = j; + if( i == j ) continue; // diagonal elements are real - ++nCompute; + result[2*iTerm+1] += a.m_pdWeights[iEvent] * + ( -a.m_pdAmps[2*nEvents*i+2*iEvent] * + a.m_pdAmps[2*nEvents*j+2*iEvent+1] + + a.m_pdAmps[2*nEvents*i+2*iEvent+1] * + a.m_pdAmps[2*nEvents*j+2*iEvent] ); } } - } - -#ifndef GPU_ACCELERATION - for( int iTerm = 0; iTerm < nCompute; ++iTerm ){ - int i = iIndex[iTerm]; - int j = jIndex[iTerm]; +#else + + // this needs to be fixed: needs a temporary vector for the chunk result + // and the calcIntegral GPU method needs to be "chunked" + + // use the GPU manager to compute the result + a.m_gpuMan.calcIntegrals( &(result[0]), nCompute, iIndex, jIndex ); + + // result += tempResult +#endif - for( int iEvent = 0; iEvent < a.m_iNTrueEvents; iEvent++ ) { - - //AiAj* - result[2*iTerm] += a.m_pdWeights[iEvent] * - ( a.m_pdAmps[2*a.m_iNEvents*i+2*iEvent] * - a.m_pdAmps[2*a.m_iNEvents*j+2*iEvent] + - a.m_pdAmps[2*a.m_iNEvents*i+2*iEvent+1] * - a.m_pdAmps[2*a.m_iNEvents*j+2*iEvent+1] ); - - if( i == j ) continue; // diagonal elements are real - - result[2*iTerm+1] += a.m_pdWeights[iEvent] * - ( -a.m_pdAmps[2*a.m_iNEvents*i+2*iEvent] * - a.m_pdAmps[2*a.m_iNEvents*j+2*iEvent+1] + - a.m_pdAmps[2*a.m_iNEvents*i+2*iEvent+1] * - a.m_pdAmps[2*a.m_iNEvents*j+2*iEvent] ); - } } - -#else - // use the GPU manager to compute the result - a.m_gpuMan.calcIntegrals( &(result[0]), nCompute, iIndex, jIndex ); + // in a "chuncked" calculation nCompute, iIndex, and jIndex get reset + // in the above loop, but they are set to the same values every time so + // at loop exit they are suitable for use in the code that follows -#endif - report( DEBUG, kModule ) << "NI terms will be renormalized by 1 / " - << iNGenEvents << endl; + << iNGenEvents << endl; // now we want to renormalize and complex congugate while // filling the integral matrix data structure @@ -932,11 +984,11 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T report( DEBUG, kModule ) << "NI result for ( " << i << ", " << j << " ) = " << complex( result[2*iTerm], result[2*iTerm+1] ) << endl; - + integralMatrix[2*i*iNAmps+2*j] = result[2*iTerm] / - static_cast< double >( iNGenEvents ); + static_cast< double >( iNGenEvents ); integralMatrix[2*i*iNAmps+2*j+1] = result[2*iTerm+1] / - static_cast< double >( iNGenEvents ); + static_cast< double >( iNGenEvents ); if( i != j ) { diff --git a/AmpTools/IUAmpTools/AmplitudeManager.h b/AmpTools/IUAmpTools/AmplitudeManager.h index eca188da..c9915be6 100644 --- a/AmpTools/IUAmpTools/AmplitudeManager.h +++ b/AmpTools/IUAmpTools/AmplitudeManager.h @@ -159,7 +159,8 @@ class AmplitudeManager : public IntensityManager * This function caculates the amplitudes for each data event and stores * them in the AmpVecs structure. It returns true if it alters the * terms in the structure (helpful to see if an update actually - * happenend when recalculating). + * happenend when recalculating). This function satisfies the need + * of the base class and uses the more complex calcTerms function below. * * \param[in,out] ampVecs a reference to the AmpVecs storage structure, four vectors * will be read from this class and amplitude caculations will be written to @@ -169,8 +170,10 @@ class AmplitudeManager : public IntensityManager * \see calcSumLogIntensity * \see calcIntegrals */ + vector calcTerms( AmpVecs& ampVecs ) const; + /** * This function calculates the intensity (one number) for all events and * stores it in the AmpVecs structure. It returns the maximum intensity @@ -227,7 +230,8 @@ class AmplitudeManager : public IntensityManager * \see calcAmplitudes * \see calcIntensities */ - void calcIntegrals( AmpVecs& ampVecs, int iNGenEvents ) const; + void calcIntegrals( AmpVecs& ampVecs, unsigned int iNGenEvents, + unsigned int chunkSize = 0 ) const; /** * The function returns a list of permutations that will be performed on @@ -427,10 +431,35 @@ class AmplitudeManager : public IntensityManager void updatePar( const string& parName ) const; - - private: + + /** + * This function is private because it can likely only safely be used by + * an implementation of calcIntegrals that reduces memory. It doesn't + * seem necessary to expose this function outside the class. + * + * This function caculates the amplitudes for each data event and stores + * them in the AmpVecs structure. It returns true if it alters the + * terms in the structure (helpful to see if an update actually + * happenend when recalculating). + * + * \param[in,out] ampVecs a reference to the AmpVecs storage structure, four vectors + * will be read from this class and amplitude caculations will be written to + * this class + * \param[in] startEvent is index of the event in the data and user variables + * array to start the computaiton with if the terms are a partial set + * \param[in] chunkSize is the number of terms to compute when startEvent + * is not equal to zero + * + * \see calcIntensities + * \see calcSumLogIntensity + * \see calcIntegrals + */ + vector calcTerms( AmpVecs& ampVecs, unsigned int startEvent, + unsigned int chunkSize ) const; + + // recursive routine to symmetrize final state void generateSymmetricCombos( const vector< pair< int, int > >& prevSwaps, vector< vector< pair< int, int > > > remainingSwaps, @@ -465,7 +494,7 @@ class AmplitudeManager : public IntensityManager mutable map< const Amplitude*, int > m_ampIteration; mutable map< AmpVecs*, map< const Amplitude*, int > > m_dataAmpIteration; - mutable map< string, unsigned long long > m_staticUserVarsOffset; + mutable map< string, unsigned int > m_staticUserVarsOffset; static const char* kModule; }; diff --git a/AmpTools/IUAmpTools/IntensityManager.h b/AmpTools/IUAmpTools/IntensityManager.h index 51aeca6c..49f5fb8a 100644 --- a/AmpTools/IUAmpTools/IntensityManager.h +++ b/AmpTools/IUAmpTools/IntensityManager.h @@ -169,7 +169,8 @@ class IntensityManager { * \see calcAmplitudes * \see calcIntensities */ - virtual void calcIntegrals( AmpVecs& ampVecs, int iNGenEvents ) const = 0; + virtual void calcIntegrals( AmpVecs& ampVecs, unsigned int iNGenEvents, + unsigned int chunkSize = 0 ) const = 0; /** * This function calculates the intensity for one event using a Kinematics From 93e7f2229f7d2b7f85ce7ebe277e04ce9d72fc21 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 11 Jun 2026 15:30:12 -0400 Subject: [PATCH 05/32] revision to compute generated integrals in parts rather than for the entire data set, which reduces memory consumption --- AmpTools/IUAmpTools/NormIntInterface.cc | 32 +++++++++++++++++++++++-- AmpTools/IUAmpTools/NormIntInterface.h | 2 ++ 2 files changed, 32 insertions(+), 2 deletions(-) diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc index 7f02bf4c..7ac5729a 100644 --- a/AmpTools/IUAmpTools/NormIntInterface.cc +++ b/AmpTools/IUAmpTools/NormIntInterface.cc @@ -428,15 +428,21 @@ NormIntInterface::forceCacheUpdate( bool normIntOnly ) const // MC in order to be able to continue assert( m_genMCVecs.m_dataLoaded ); + // computing integrals of the generated MC can be done in chunks + // to avoid exhausting memory on CPU or GPU -- there is never a need + // to have the entire generated MC in memory at once + unsigned int chunkSize = genMCChunkSize(); + // do "lazy" allocation of memory here -- this is important for MPI jobs // where forceCacheUpdate is only called on follower nodes, as it // avoids big memory allocations on the lead nodes - if( m_genMCVecs.m_iNTerms == 0 ) m_genMCVecs.allocateTerms( *m_pIntenManager ); + if( m_genMCVecs.m_iNTerms == 0 ) m_genMCVecs.allocateTerms( *m_pIntenManager, chunkSize ); report( DEBUG, kModule ) << "Asking IntensityManager to calculate integrals " << "using the generated MC." << endl; - m_pIntenManager->calcIntegrals( m_genMCVecs, m_nGenEvents ); + m_pIntenManager->calcIntegrals( m_genMCVecs, m_nGenEvents, chunkSize ); + setAmpIntMatrix( m_genMCVecs.m_pdIntegralMatrix ); m_emptyAmpIntCache = false; @@ -553,6 +559,28 @@ NormIntInterface::setNormIntMatrix( const double* input ) const { } +unsigned int +NormIntInterface::genMCChunkSize() const { + + // need to know the sizes of the data sets + assert( m_accMCVecs.m_dataLoaded && m_genMCVecs.m_dataLoaded ); + + unsigned long nGen = m_genMCVecs.m_iNEvents; + unsigned long nAcc = m_accMCVecs.m_iNEvents; + + // get the chunk size by doing integer division by 2 + // until the size of the generated MC is less than the + // accepted MC -- for GPU fits this should return a + // chunk size that is a power of 2, which is required + + int iPow = 0; + while( ( nGen >> iPow ) > nAcc ) ++iPow; + + report( DEBUG, kModule ) << "Chunk size for generated MC: " << ( nGen >> iPow ) << endl; + + return( nGen >> iPow ); +} + #ifndef __ACLIC__ void NormIntInterface::invalidateTerms(){ diff --git a/AmpTools/IUAmpTools/NormIntInterface.h b/AmpTools/IUAmpTools/NormIntInterface.h index 976420d7..6f58a2e3 100644 --- a/AmpTools/IUAmpTools/NormIntInterface.h +++ b/AmpTools/IUAmpTools/NormIntInterface.h @@ -122,6 +122,8 @@ class NormIntInterface void initializeCache(); int m_cacheSize; + unsigned int genMCChunkSize() const; + vector< string > m_termNames; map< string, int > m_termIndex; From a1f0f1a1484d499b86ba64f633d8a90650add138 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Mon, 15 Jun 2026 13:49:53 -0400 Subject: [PATCH 06/32] change array name to make consistent with other classes, add debug statemnt --- AmpTools/IUAmpTools/Amplitude.cc | 15 ++++++++------- AmpTools/IUAmpTools/Amplitude.h | 4 ++-- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/AmpTools/IUAmpTools/Amplitude.cc b/AmpTools/IUAmpTools/Amplitude.cc index ca52d6c0..c4518cd3 100644 --- a/AmpTools/IUAmpTools/Amplitude.cc +++ b/AmpTools/IUAmpTools/Amplitude.cc @@ -78,8 +78,6 @@ Amplitude::calcUserVarsAll( GDouble* pdData, GDouble* pdUserVars, int iNEvents, SCOREP_USER_REGION_DEFINE( calcUserVarsAll ) #endif - - report( DEBUG, kModule ) << "Caculating user data for " << name() << endl; unsigned int numVars = numUserVars(); @@ -132,7 +130,7 @@ SCOREP_USER_REGION_END( calcUserVarsAll ) } void -Amplitude::calcAmplitudeAll( GDouble* pdData, GDouble* pdAmps, unsigned int iNEvents, +Amplitude::calcAmplitudeAll( GDouble* pdData, GDouble* pdAmpFact, unsigned int iNEvents, const vector< vector< int > >* pvPermutations, GDouble* pdUserVars, unsigned int startEvent, unsigned int chunkSize ) const @@ -155,7 +153,7 @@ SCOREP_USER_REGION_BEGIN( calcAmplitudeAll, "calcAmplitudeAll", SCOREP_USER_REGI GDouble** pKin = new GDouble*[iNParticles]; - // this is the size of the chunk in the pdAmps array that we + // this is the size of the chunk in the pdAmpFact array that we // are going to compute amplitdues for unsigned int nEvents = ( chunkSize == 0 ? iNEvents : chunkSize ); @@ -189,15 +187,18 @@ SCOREP_USER_REGION_BEGIN( calcAmplitudeAll, "calcAmplitudeAll", SCOREP_USER_REGI cRes = calcAmplitude( pKin ); } - pdAmps[2*nEvents*iPermutation+2*iEvent] = cRes.real(); - pdAmps[2*nEvents*iPermutation+2*iEvent+1] = cRes.imag(); + pdAmpFact[2*nEvents*iPermutation+2*iEvent] = cRes.real(); + pdAmpFact[2*nEvents*iPermutation+2*iEvent+1] = cRes.imag(); } } - + #ifdef SCOREP SCOREP_USER_REGION_END( calcAmplitudeAll ) #endif + report( DEBUG, kModule ) << "First event and permutation amplitude: ( " + << pdAmpFact[0] << ", " << pdAmpFact[1] << " )" << endl; + delete[] pKin; } diff --git a/AmpTools/IUAmpTools/Amplitude.h b/AmpTools/IUAmpTools/Amplitude.h index e524924e..f5d90e9e 100644 --- a/AmpTools/IUAmpTools/Amplitude.h +++ b/AmpTools/IUAmpTools/Amplitude.h @@ -274,7 +274,7 @@ class Amplitude * of GDoubles E, px, py, pz repeated sequentially for each particle in * the event and then the series repeated for each event in the data set. * - * \param[out] pdAmps a pointer to a block of memory where the calculated + * \param[out] pdAmpFact a pointer to a block of memory where the calculated * amplitudes should go. The values are stored as real and imaginary parts * repeated for each permutation and then for each event. * Total length of memory needed is 2 * sizeof( GDDouble ) * N_evts * N_perms @@ -292,7 +292,7 @@ class Amplitude * \see calcAmplitudeAll * \see AmplitudeManager::addAmpPermutation */ - virtual void calcAmplitudeAll( GDouble* pdData, GDouble* pdAmps, unsigned int iNEvents, + virtual void calcAmplitudeAll( GDouble* pdData, GDouble* pdAmpFact, unsigned int iNEvents, const vector< vector< int > >* pvPermutations, GDouble* pdUserVars = 0, unsigned int startEvent = 0, unsigned int chunkSize = 0 ) const; From 5a4f96002fae158ae3df6af22307af4a9b15adae Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Mon, 15 Jun 2026 13:50:24 -0400 Subject: [PATCH 07/32] change variable name to make use more clear --- AmpTools/IUAmpTools/AmpVecs.cc | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/AmpTools/IUAmpTools/AmpVecs.cc b/AmpTools/IUAmpTools/AmpVecs.cc index de4fa15d..a8e04ebe 100644 --- a/AmpTools/IUAmpTools/AmpVecs.cc +++ b/AmpTools/IUAmpTools/AmpVecs.cc @@ -316,7 +316,7 @@ AmpVecs::loadData( DataReader* pDataReader ){ void AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity, unsigned int chunkSize ){ - unsigned int iNEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize ); + unsigned int ampsEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize ); m_iNTerms = intenMan.getTermNames().size(); m_maxFactPerEvent = intenMan.maxFactorStoragePerEvent(); @@ -355,12 +355,14 @@ AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity, #ifndef GPU_ACCELERATION - m_pdAmps = new GDouble[iNEvents * intenMan.termStoragePerEvent()]; - m_pdAmpFactors = new GDouble[iNEvents * m_maxFactPerEvent]; + // these use reduced size arrays when a non-zero chunk size is provided: + + m_pdAmps = new GDouble[ampsEvents * intenMan.termStoragePerEvent()]; + m_pdAmpFactors = new GDouble[ampsEvents * m_maxFactPerEvent]; #else - m_gpuMan.initTerms( *this, iNEvents ); + m_gpuMan.initTerms( *this, ampsEvents ); #endif // GPU_ACCELERATION From 9abb7acfa510972cd7e5c303e89da622762ea72f Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Mon, 15 Jun 2026 13:51:22 -0400 Subject: [PATCH 08/32] fix bug where the wrong number of events was used in memory indexing --- AmpTools/IUAmpTools/AmplitudeManager.cc | 29 +++++++++++++------------ 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/AmpTools/IUAmpTools/AmplitudeManager.cc b/AmpTools/IUAmpTools/AmplitudeManager.cc index 72229d97..86be2a0a 100644 --- a/AmpTools/IUAmpTools/AmplitudeManager.cc +++ b/AmpTools/IUAmpTools/AmplitudeManager.cc @@ -424,9 +424,7 @@ AmplitudeManager::calcTerms( AmpVecs& a ) const // in cases where the user is not calculating a chunk of terms use the // specialized function below with default arguments - vector result = calcTerms( a, 0, 0 ); - - return result; + return calcTerms( a, 0, 0 ); } vector @@ -550,8 +548,7 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON uOffset ); #endif//GPU_ACCELERATION } - - + // now assemble all the factors in an amplitude into a single // symmetrized amplitude for each event @@ -561,9 +558,7 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON GDouble dAmpFacRe, dAmpFacIm, dTRe, dTIm; int iEvent, iPerm; unsigned int iOffsetA, iOffsetP, iOffsetF; - - // re-ordering of data will be useful to not fall out of (CPU) memory cache!!! - + // zeroing out the entire range memset( (void*)( a.m_pdAmps + 2 * nEvents * iAmpIndex ), 0, 2 * nEvents * sizeof(GDouble) ); @@ -600,6 +595,9 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON a.m_pdAmps[iOffsetA+1] *= dSymmFactor; } + report( DEBUG, kModule ) << "Amplitude index " << iAmpIndex << ", event 0: (" << + a.m_pdAmps[2*nEvents*iAmpIndex] << ", " << a.m_pdAmps[2*nEvents*iAmpIndex+1] << " )" << endl; + #else // on the GPU the terms are assembled and never copied out // of GPU memory @@ -828,8 +826,8 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T assert( iNGenEvents ); // figure out the number of chunks needed to go through the data - unsigned int nChunk = ( chunkSize == 0 ? 1 : iNGenEvents / chunkSize ); - if( ( chunkSize !=0 ) && ( iNGenEvents % chunkSize != 0 ) ) ++nChunk; + unsigned int nChunk = ( chunkSize == 0 ? 1 : a.m_iNEvents / chunkSize ); + if( ( chunkSize !=0 ) && ( a.m_iNEvents % chunkSize != 0 ) ) ++nChunk; int iNAmps = a.m_iNTerms; @@ -849,13 +847,13 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T for( unsigned int iChunk = 0; iChunk < nChunk; ++iChunk ){ unsigned int startEvent = iChunk * chunkSize; - unsigned int nEvents = iNGenEvents; + unsigned int nEvents = a.m_iNEvents; // for chuncked calculation: // the number of events to process is the chunk size unless // it is the last chunk then it is the remainder if( chunkSize != 0 ){ nEvents = chunkSize; - if( iChunk - nChunk == 1 ) nEvents = iNGenEvents % chunkSize; + if( iChunk - nChunk == 1 ) nEvents = a.m_iNEvents % chunkSize; } // this returns a vector indicating which terms have changed @@ -947,12 +945,15 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T if( i == j ) continue; // diagonal elements are real - result[2*iTerm+1] += a.m_pdWeights[iEvent] * + result[2*iTerm+1] += a.m_pdWeights[iEvent+startEvent] * ( -a.m_pdAmps[2*nEvents*i+2*iEvent] * a.m_pdAmps[2*nEvents*j+2*iEvent+1] + a.m_pdAmps[2*nEvents*i+2*iEvent+1] * a.m_pdAmps[2*nEvents*j+2*iEvent] ); } + + report( DEBUG, kModule ) << "NI Term Sum ( " << i << ", " << j << " ) = ( " << + result[2*iTerm] << ", " << result[2*iTerm+1] << ")" << endl; } #else @@ -968,7 +969,7 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T } - // in a "chuncked" calculation nCompute, iIndex, and jIndex get reset + // in a "chunked" calculation nCompute, iIndex, and jIndex get reset // in the above loop, but they are set to the same values every time so // at loop exit they are suitable for use in the code that follows From 76a41da9bcfbd9c93334735b78b86df0e3ec9847 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Mon, 15 Jun 2026 13:51:35 -0400 Subject: [PATCH 09/32] typo fix in comment --- AmpTools/IUAmpTools/NormIntInterface.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc index 7ac5729a..706a2463 100644 --- a/AmpTools/IUAmpTools/NormIntInterface.cc +++ b/AmpTools/IUAmpTools/NormIntInterface.cc @@ -450,7 +450,7 @@ NormIntInterface::forceCacheUpdate( bool normIntOnly ) const // first trip through or forced update... // we either copy the generated MC integrals to accepted or - // compute thm directly: + // compute them directly: if( ( m_accMCReader == m_genMCReader ) && !m_emptyAmpIntCache ) { // optimization for perfect acceptance From a36c430cf8dcd79ebfd71855e7bfc01bc5c75255 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Wed, 17 Jun 2026 18:17:47 -0400 Subject: [PATCH 10/32] modify to pass information about need for user variables to loadData --- AmpTools/IUAmpTools/NormIntInterface.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc index 706a2463..82b72c35 100644 --- a/AmpTools/IUAmpTools/NormIntInterface.cc +++ b/AmpTools/IUAmpTools/NormIntInterface.cc @@ -106,7 +106,7 @@ m_termNames( intenManager.getTermNames() ) if( genVecs == m_uniqueDataSets.end() ){ report( INFO, kModule ) << "Loading generated Monte Carlo from file..." << endl; - m_genMCVecs.loadData( m_genMCReader ); + m_genMCVecs.loadData( m_genMCReader, intenManager.needsUserVarsOnly() ); m_uniqueDataSets[m_genMCReader] = &m_genMCVecs; } @@ -123,7 +123,7 @@ m_termNames( intenManager.getTermNames() ) if( accVecs == m_uniqueDataSets.end() ){ report( INFO, kModule ) << "Loading accepted Monte Carlo from file..." << endl; - m_accMCVecs.loadData( m_accMCReader ); + m_accMCVecs.loadData( m_accMCReader, intenManager.needsUserVarsOnly() ); m_uniqueDataSets[m_accMCReader] = &m_accMCVecs; } From 40b2e38a2c05c2f2f896714ac8f63c26ec01a0d2 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Wed, 17 Jun 2026 18:19:10 -0400 Subject: [PATCH 11/32] modify to pass information about need for user variables to loadData --- AmpTools/IUAmpTools/LikelihoodCalculator.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/AmpTools/IUAmpTools/LikelihoodCalculator.cc b/AmpTools/IUAmpTools/LikelihoodCalculator.cc index c40420b1..ed05a91a 100644 --- a/AmpTools/IUAmpTools/LikelihoodCalculator.cc +++ b/AmpTools/IUAmpTools/LikelihoodCalculator.cc @@ -225,7 +225,7 @@ SCOREP_USER_REGION_BEGIN( dataTerm, "dataTerm", SCOREP_USER_REGION_TYPE_COMMON ) report( DEBUG, kModule ) << "Allocating Data and Amplitude Array in LikelihoodCalculator for " << m_intenManager.reactionName() << "..." << endl; - m_ampVecsSignal.loadData( m_dataReaderSignal ); + m_ampVecsSignal.loadData( m_dataReaderSignal, m_intenManager.needsUserVarsOnly() ); m_ampVecsSignal.allocateTerms( m_intenManager, true ); m_numDataEvents = m_ampVecsSignal.m_iNTrueEvents; @@ -247,7 +247,7 @@ SCOREP_USER_REGION_BEGIN( dataTerm, "dataTerm", SCOREP_USER_REGION_TYPE_COMMON ) if( m_hasBackground ){ - m_ampVecsBkgnd.loadData( m_dataReaderBkgnd ); + m_ampVecsBkgnd.loadData( m_dataReaderBkgnd, m_intenManager.needsUserVarsOnly() ); m_ampVecsBkgnd.allocateTerms( m_intenManager, true ); if( m_ampVecsBkgnd.m_hasMixedSignWeights ){ From 8664259498ef84f75a789e28f616fb4f451e88a8 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:19:22 -0400 Subject: [PATCH 12/32] include math header --- AmpTools/GPUManager/CUDA-Complex.cuh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/AmpTools/GPUManager/CUDA-Complex.cuh b/AmpTools/GPUManager/CUDA-Complex.cuh index ada2f7ac..9e92b915 100644 --- a/AmpTools/GPUManager/CUDA-Complex.cuh +++ b/AmpTools/GPUManager/CUDA-Complex.cuh @@ -37,6 +37,8 @@ #ifndef __CUDA_COMPLEX__H__ #define __CUDA_COMPLEX__H__ +#include + #include "GPUManager/GPUCustomTypes.h" #ifdef __CUDACC__ From b7ab5b9c4b0cd5f7b15b2b58dfa2d08d306fb9b4 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:20:09 -0400 Subject: [PATCH 13/32] change event numbers to unsigned int --- AmpTools/GPUManager/GPUAmpProductKernel.cu | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/AmpTools/GPUManager/GPUAmpProductKernel.cu b/AmpTools/GPUManager/GPUAmpProductKernel.cu index 9a8770cc..2a0716af 100644 --- a/AmpTools/GPUManager/GPUAmpProductKernel.cu +++ b/AmpTools/GPUManager/GPUAmpProductKernel.cu @@ -38,7 +38,7 @@ __global__ void amp_kernel( GDouble* pfDevAmps, GDouble* pfDevVVStar, GDouble* pfDevWeights, - int nAmps, int nEvents, double* pdDevRes ) + int nAmps, unsigned int nEvents, double* pdDevRes ) { int i = threadIdx.x + GPU_BLOCK_SIZE_X * threadIdx.y + ( blockIdx.x + blockIdx.y * gridDim.x ) * GPU_BLOCK_SIZE_SQ; @@ -100,7 +100,7 @@ amp_kernel( GDouble* pfDevAmps, GDouble* pfDevVVStar, GDouble* pfDevWeights, extern "C" void GPU_ExecAmpKernel( dim3 dimGrid, dim3 dimBlock, GDouble* pfDevAmps, GDouble* pfDevVVStar, GDouble* pfDevWeights, - int nAmps, int nEvents, double* pdDevRes ) + int nAmps, unsigned int nEvents, double* pdDevRes ) { amp_kernel<<< dimGrid, dimBlock >>>( pfDevAmps, pfDevVVStar, pfDevWeights, nAmps, nEvents, pdDevRes ); From 9b3dfae2cc113fe59d4b0fc5755c8b6404b9620a Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:21:09 -0400 Subject: [PATCH 14/32] enhance amplitude arguments for chunked calculations of amplitudes --- AmpTools/GPUManager/GPUCustomTypes.h | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/AmpTools/GPUManager/GPUCustomTypes.h b/AmpTools/GPUManager/GPUCustomTypes.h index faff02ce..d78f50ed 100755 --- a/AmpTools/GPUManager/GPUCustomTypes.h +++ b/AmpTools/GPUManager/GPUCustomTypes.h @@ -50,12 +50,15 @@ using namespace std; // a few helpers to make operations on GPU a little more user friendly -// standard arguments to kernel launches +// standard arguments to ampliutde kernel launches #define GPU_AMP_PROTO GDouble* pfDevData, GDouble* pfDevUserVars, \ - WCUComplex* pcDevAmp, int* piDevPerm, int iNParticles, int iNEvents -#define GPU_AMP_ARGS pfDevData, pfDevUserVars, pcDevAmp, piDevPerm, iNParticles, iNEvents + WCUComplex* pcDevAmp, int* piDevPerm, unsigned int iNParticles, \ + unsigned int iNEvents, unsigned int startEvent +#define GPU_AMP_ARGS pfDevData, pfDevUserVars, pcDevAmp, piDevPerm, iNParticles, \ + iNEvents, startEvent -// how to get the event we are working with: +// how to get the event we are working with -- this is within a chunk and is typically +// set to iEvent by the user #define GPU_THIS_EVENT threadIdx.x + GPU_BLOCK_SIZE_X * threadIdx.y + \ ( blockIdx.x + blockIdx.y * gridDim.x ) * GPU_BLOCK_SIZE_SQ @@ -65,11 +68,13 @@ using namespace std; #define GPU_PY 2 #define GPU_PZ 3 -#define GPU_KIN(id,val) pfDevData[ ( piDevPerm[id] * 4 + val ) * iNEvents + iEvent] +// both the data and user vars are always stored for the entire data set so the indices +// need to be adjusted by startEvent when doing a chunked amplitude calculation +#define GPU_KIN(id,val) pfDevData[ ( piDevPerm[id] * 4 + val ) * iNEvents + iEvent + startEvent ] #define GPU_P4(id) { GPU_KIN(id,GPU_E) , GPU_KIN(id,GPU_PX), \ GPU_KIN(id,GPU_PY) , GPU_KIN(id,GPU_PZ)} -#define GPU_UVARS(val) pfDevUserVars[val*iNEvents+iEvent] +#define GPU_UVARS(val) pfDevUserVars[val*iNEvents + iEvent + startEvent] #define COPY_P4(a1,a2) GDouble a2[4]; for( int zzz = 0; zzz < 4; ++zzz ) a2[zzz] = a1[zzz]; From f7f9a1db37abf909fed78dd9c7c87c46188d92f0 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:21:36 -0400 Subject: [PATCH 15/32] change variable names to be more consistent --- AmpTools/GPUManager/GPUFactPermKernel.cu | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/AmpTools/GPUManager/GPUFactPermKernel.cu b/AmpTools/GPUManager/GPUFactPermKernel.cu index 9b3dcac3..faefbc0e 100644 --- a/AmpTools/GPUManager/GPUFactPermKernel.cu +++ b/AmpTools/GPUManager/GPUFactPermKernel.cu @@ -37,8 +37,8 @@ #include "GPUCustomTypes.h" __global__ void -fact_perm_kernel( GDouble* pfDevAmps, GDouble* pcDevCalcAmp, int nFact, - int nPerm, int nEvents ) +fact_perm_kernel( GDouble* pfDevAmps, GDouble* pcDevAmpFact, int nFact, + int nPerm, unsigned int nEvents ) { int i = threadIdx.x + GPU_BLOCK_SIZE_X * threadIdx.y + @@ -46,20 +46,20 @@ fact_perm_kernel( GDouble* pfDevAmps, GDouble* pcDevCalcAmp, int nFact, for( int iPerm = 0; iPerm < nPerm; ++iPerm ){ - int offsetP = 2*nEvents*iPerm + 2*i; + unsigned int offsetP = 2*nEvents*iPerm + 2*i; - double ampRe = pcDevCalcAmp[offsetP]; - double ampIm = pcDevCalcAmp[offsetP+1]; + GDouble ampRe = pcDevAmpFact[offsetP]; + GDouble ampIm = pcDevAmpFact[offsetP+1]; for( int iFactor = 1; iFactor < nFact; ++iFactor ){ - int offsetF = 2*nEvents*nPerm*iFactor + offsetP; + unsigned int offsetF = 2*nEvents*nPerm*iFactor + offsetP; - double re = ampRe; - double im = ampIm; + GDouble re = ampRe; + GDouble im = ampIm; - ampRe = re * pcDevCalcAmp[offsetF] - im * pcDevCalcAmp[offsetF+1]; - ampIm = re * pcDevCalcAmp[offsetF+1] + im * pcDevCalcAmp[offsetF]; + ampRe = re * pcDevAmpFact[offsetF] - im * pcDevAmpFact[offsetF+1]; + ampIm = re * pcDevAmpFact[offsetF+1] + im * pcDevAmpFact[offsetF]; } pfDevAmps[2*i] += ampRe; @@ -71,8 +71,8 @@ fact_perm_kernel( GDouble* pfDevAmps, GDouble* pcDevCalcAmp, int nFact, } extern "C" void GPU_ExecFactPermKernel( dim3 dimGrid, dim3 dimBlock, - GDouble* pfDevAmps, GDouble* pcDevCalcAmp, int nFact, int nPerm, int nEvents ) + GDouble* pfDevAmps, GDouble* pcDevAmpFact, int nFact, int nPerm, unsigned int nEvents ) { - fact_perm_kernel<<< dimGrid, dimBlock >>>( pfDevAmps, pcDevCalcAmp, + fact_perm_kernel<<< dimGrid, dimBlock >>>( pfDevAmps, pcDevAmpFact, nFact, nPerm, nEvents ); } From 640db651535941dab7e8936ec21419c7ed11d2a4 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:22:34 -0400 Subject: [PATCH 16/32] variable change to unsigned int and include functionality for chunked calculation --- AmpTools/GPUManager/GPUNICalcKernel.cu | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/AmpTools/GPUManager/GPUNICalcKernel.cu b/AmpTools/GPUManager/GPUNICalcKernel.cu index a57068a1..244ab43c 100644 --- a/AmpTools/GPUManager/GPUNICalcKernel.cu +++ b/AmpTools/GPUManager/GPUNICalcKernel.cu @@ -40,7 +40,8 @@ __global__ void ni_calc_kernel( int nElements, double* pdDevNICalc, GDouble* pfDevAmps, GDouble* pfDevWeights, - int nEvents, int nTrueEvents ) + unsigned int startEvent, unsigned int nEvents, + unsigned int nTrueEvents ) { // used shared memory block for amplitude indices and results @@ -73,10 +74,13 @@ ni_calc_kernel( int nElements, double* pdDevNICalc, __syncthreads(); - // this is the overall event index - int iEvt = iThread + + // this is the event index in the chunk + unsigned int iEvtChunk = iThread + ( blockIdx.x + blockIdx.y * gridDim.x ) * GPU_BLOCK_SIZE_SQ; + // this is the event index in the overall set of events + unsigned int iEvt = iEvtChunk + startEvent; + if( iEvt < nTrueEvents ){ // do not compute for the padding events GDouble wt = pfDevWeights[iEvt]; @@ -85,8 +89,8 @@ ni_calc_kernel( int nElements, double* pdDevNICalc, // these are the indices to the // relevant amplitudes in the amplitude array - int aInd = 2*iEvt + 2*nEvents*iIndex[i]; - int bInd = 2*iEvt + 2*nEvents*jIndex[i]; + int aInd = 2*iEvtChunk + 2*nEvents*iIndex[i]; + int bInd = 2*iEvtChunk + 2*nEvents*jIndex[i]; GDouble aRe = pfDevAmps[aInd]; GDouble aIm = pfDevAmps[aInd+1]; @@ -124,15 +128,16 @@ ni_calc_kernel( int nElements, double* pdDevNICalc, extern "C" void GPU_ExecNICalcKernel( dim3 dimGrid, dim3 dimBlock, - unsigned int sharedSize, + unsigned int sharedSize, int nElements, double* pdDevNICalc, GDouble* pfDevAmps, GDouble* pfDevWeights, - int nEvents, int nTrueEvents, unsigned int maxSize ) + unsigned int startEvent, unsigned int nEvents, + unsigned int nTrueEvents, unsigned int maxSize ) { if( maxSize ) cudaFuncSetAttribute( ni_calc_kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, maxSize ); ni_calc_kernel<<< dimGrid, dimBlock, sharedSize >>> - ( nElements, pdDevNICalc, pfDevAmps, pfDevWeights, nEvents, nTrueEvents ); + ( nElements, pdDevNICalc, pfDevAmps, pfDevWeights, startEvent, nEvents, nTrueEvents ); } From a53879216ece77fcad44acd9c21c5ca93618f192 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:23:49 -0400 Subject: [PATCH 17/32] large revisions to enable chunked calculation of integrals --- AmpTools/GPUManager/GPUManager.cc | 238 ++++++++++++++++++++---------- AmpTools/GPUManager/GPUManager.h | 64 ++++---- 2 files changed, 200 insertions(+), 102 deletions(-) diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc index 02dfa448..7fe0299a 100644 --- a/AmpTools/GPUManager/GPUManager.cc +++ b/AmpTools/GPUManager/GPUManager.cc @@ -49,6 +49,7 @@ #include "GPUManager/GPUKernel.h" #include "GPUManager/GPUManager.h" +#include "GPUManager/CUDA-Complex.cuh" #include "IUAmpTools/report.h" const char* GPUManager::kModule = "GPUManager"; @@ -71,7 +72,9 @@ GPUManager::GPUManager() m_iNUserVars=0; m_iGDoubleDataArrSize=0; - m_iDoubleEventArrSize=0; + + m_iDoubleIntenArrSize=0; + m_iGDoubleFactArrSize=0; m_iAmpArrSize=0; m_iVArrSize=0; m_iNICalcSize=0; @@ -83,7 +86,7 @@ GPUManager::GPUManager() //Device Arrays m_pfDevData=0; m_pfDevUserVars=0; - m_pcDevCalcAmp=0; + m_pcDevAmpFact=0; m_piDevPerm=0; m_pfDevAmps=0; m_pfDevWeights=0; @@ -173,6 +176,13 @@ GPUManager::~GPUManager() // Initialization routines: +void +GPUManager::init( const AmpVecs& a, bool use4Vectors ) +{ + initData( a, use4Vectors ); + initTerms( a ); +} + void GPUManager::initData( const AmpVecs& a, bool use4Vectors ) { @@ -190,23 +200,26 @@ GPUManager::initData( const AmpVecs& a, bool use4Vectors ) totalMemory += m_iGDoubleDataArrSize; if( use4Vectors ) totalMemory += 4 * m_iNParticles * m_iGDoubleDataArrSize; - totalMemory += m_iNUserVars * m_iGDoubleDataArrSize; totalMemory += m_iNParticles * sizeof( int ); totalMemory /= (1024*1024); - report( INFO, kModule ) << "Attempting to allocate " << (int)totalMemory << " MB of global GPU memory." << endl; + report( INFO, kModule ) << "Attempting to allocate " << (int)totalMemory + << " MB of global GPU memory for data." << endl; + report( DEBUG, kModule ) << "AmpVecs details (initData): " << endl; + report( DEBUG, kModule ) << "\tiNEvents: " << m_iNEvents << endl; + report( DEBUG, kModule ) << "\tiNTrueEvents: " << m_iNTrueEvents << endl; + report( DEBUG, kModule ) << "\tiNParticles: " << m_iNParticles << endl; gpuErrChk( cudaMalloc( (void**)&m_pfDevWeights , m_iGDoubleDataArrSize ) ) ; // allocate device memory needed for amplitude calculations - if( use4Vectors ) gpuErrChk( cudaMalloc( (void**)&m_pfDevData , - 4 * m_iNParticles * m_iGDoubleDataArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pfDevUserVars, m_iNUserVars * m_iGDoubleDataArrSize ) ) ; + if( use4Vectors ) + gpuErrChk( cudaMalloc( (void**)&m_pfDevData, 4 * m_iNParticles * m_iGDoubleDataArrSize ) ) ; report( INFO, kModule ) << "GPU memory allocated for storage of " - << m_iNEvents << " events (" << m_iNTrueEvents << " actual events)." - << endl; + << m_iNEvents << " events (" << m_iNTrueEvents + << " actual events)." << endl; //CUDA Dims calcCUDADims(); @@ -214,30 +227,33 @@ GPUManager::initData( const AmpVecs& a, bool use4Vectors ) void -GPUManager::initTerms( const AmpVecs& a, unsigned long chunkSize ) +GPUManager::initTerms( const AmpVecs& a, unsigned int chunkSize ) { clearTerms(); + report( DEBUG, kModule ) << "Initializing GPU arrays for amplitude calculations." << endl; + // chunkSize must be power of 2 when provided assert( chunkSize == 0 || ( (chunkSize & (chunkSize - 1)) == 0 ) ); // use all events or the specified size of a chunk of events - long int iNEvents = ( chunkSize == 0 ? m_iNevents : chunkSize ); + unsigned int iNAmpEvents = ( chunkSize == 0 ? a.m_iNEvents : chunkSize ); - m_iNAmps = a.m_iNTerms; - m_iNUserVars = a.m_userVarsPerEvent; + report( DEBUG, kModule ) << "\tChunk size for amplitude calculations: " << iNAmpEvents << endl; + + m_iNUserVars = a.m_userVarsPerEvent; + m_iNAmps = a.m_iNTerms; // double precision intensity calculation - m_iDoubleTrueIntenArrSize = sizeof(double) * m_iNTrueEvents; - m_iDoubleIntenArrSize = sizeof(double) * m_iNEvents; + m_iDoubleIntenArrSize = sizeof(double) * a.m_iNEvents; // GDouble precision for factors and amplitudes - // ... and use the chunkSize for these (iNEvents) when it is provided + // ... and use the chunkSize for these (iNAmpEvents) when it is provided - m_iGDoubleFactArrSize = sizeof(GDouble) * a.m_maxFactPerEvent * iNEvents; + m_iGDoubleFactArrSize = 2 * sizeof(GDouble) * a.m_maxFactPerEvent * iNAmpEvents; // size needed to store amplitudes for each event - m_iAmpArrSize = 2 * sizeof(GDouble) * iNEvents * m_iNAmps; + m_iAmpArrSize = 2 * sizeof(GDouble) * iNAmpEvents * m_iNAmps; // size of upper half of ViVj* matrix m_iVArrSize = 2 * sizeof(GDouble) * m_iNAmps * ( m_iNAmps + 1 ) / 2; @@ -250,10 +266,11 @@ GPUManager::initTerms( const AmpVecs& a, unsigned long chunkSize ) // host memory needed for intensity or integral calculation cudaMallocHost( (void**)&m_pfVVStar , m_iVArrSize ); - cudaMallocHost( (void**)&m_pdRes , m_iDoubleEventArrSize ); + cudaMallocHost( (void**)&m_pdRes , m_iDoubleIntenArrSize ); double totalMemory = 0; + totalMemory += m_iNUserVars * m_iGDoubleDataArrSize; totalMemory += m_iVArrSize; totalMemory += m_iNICalcSize; totalMemory += 2*m_iDoubleIntenArrSize; @@ -264,21 +281,36 @@ GPUManager::initTerms( const AmpVecs& a, unsigned long chunkSize ) totalMemory /= (1024*1024); report( INFO, kModule ) << "Attempting to allocate " << (int)totalMemory - << " MB of global GPU memory." << endl; - + << " MB of global GPU memory for amplitude calculations." << endl; + + report( DEBUG, kModule ) << "GPU memory allocation details:" << endl; + report( DEBUG, kModule ) << "\tsize of user vars array: " << m_iNUserVars * m_iGDoubleDataArrSize / 1024 << " kB ( " + << m_iNUserVars * m_iGDoubleDataArrSize/sizeof(GDouble) << " elements )" << endl; + report( DEBUG, kModule ) << "\tsize of amplitude array: " << m_iAmpArrSize / 1024 << " kB ( " + << m_iAmpArrSize/sizeof(GDouble) << " elements )" << endl; + report( DEBUG, kModule ) << "\tsize of factor array: " << m_iGDoubleFactArrSize / 1024 << " kB ( " + << m_iGDoubleFactArrSize/sizeof(GDouble) << " elements )" << endl; + report( DEBUG, kModule ) << "\tsize of VV* array: " << m_iVArrSize / 1024 << " kB ( " + << m_iVArrSize/sizeof(GDouble) << " elements )" << endl; + report( DEBUG, kModule ) << "\tsize of NICalc array: " << m_iNICalcSize / 1024 << " kB ( " + << m_iNICalcSize/sizeof(double) << " elements )" << endl; + // device memory needed for intensity or integral calculation and sum + gpuErrChk( cudaMalloc( (void**)&m_pfDevUserVars, m_iNUserVars * m_iGDoubleDataArrSize ) ) ; gpuErrChk( cudaMalloc( (void**)&m_pfDevVVStar , m_iVArrSize ) ) ; gpuErrChk( cudaMalloc( (void**)&m_pdDevNICalc , m_iNICalcSize ) ) ; gpuErrChk( cudaMalloc( (void**)&m_pdDevRes , m_iDoubleIntenArrSize ) ) ; gpuErrChk( cudaMalloc( (void**)&m_pdDevREDUCE , m_iDoubleIntenArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_pcDevCalcAmp , m_iGDoubleFactArrSize ) ) ; + gpuErrChk( cudaMalloc( (void**)&m_pcDevAmpFact , m_iGDoubleFactArrSize ) ) ; gpuErrChk( cudaMalloc( (void**)&m_pfDevAmps , m_iAmpArrSize ) ) ; - gpuErrChk( cudaMalloc( (void**)&m_piDevPerm , m_iNParticles * sizeof( int ) ) ) ; + gpuErrChk( cudaMalloc( (void**)&m_piDevPerm , m_iNParticles * sizeof( int ) ) ) ; report( INFO, kModule ) << "GPU memory allocated for " << m_iNAmps << " amplitudes for " - << iNEvents << " events (" << m_iNTrueEvents << " total events)." + << iNAmpEvents << " events (" << m_iNTrueEvents << " total events)." << endl; + + calcCUDADimsAmpFact( iNAmpEvents ); } @@ -358,6 +390,11 @@ GPUManager::copyUserVarsToGPU( const AmpVecs& a ) // make sure AmpVecs has been loaded with data assert( a.m_pdUserVars ); + report( DEBUG, kModule ) << "Copying user variables to GPU for " << m_iNUserVars << " variables" << endl; + report( DEBUG, kModule ) << "\tdevice pointer: " << m_pfDevUserVars << endl; + report( DEBUG, kModule ) << "\thost pointer: " << a.m_pdUserVars << endl; + report( DEBUG, kModule ) << "\tsize: " << m_iNUserVars * m_iGDoubleDataArrSize << " bytes" << endl; + // copy the data into the device gpuErrChk( cudaMemcpy( m_pfDevUserVars, a.m_pdUserVars, m_iNUserVars * m_iGDoubleDataArrSize, @@ -380,11 +417,22 @@ GPUManager::copyAmpsFromGPU( AmpVecs& a ) // to save CPU memory -- the user must allocate it explicitly assert( a.m_pdAmpFactors != NULL && a.m_pdAmps != NULL ); + report( DEBUG, kModule ) << "Copying amplitudes and factors from GPU to CPU for " << m_iNAmps << " amplitudes and " + << m_iNEvents << " events." << endl; + report( DEBUG, kModule ) << "\tdevice pointer for amplitude: " << m_pfDevAmps << endl; + report( DEBUG, kModule ) << "\tdevice pointer for factors: " << m_pcDevAmpFact << endl; + report( DEBUG, kModule ) << "\thost pointer for amplitude: " << a.m_pdAmps << endl; + report( DEBUG, kModule ) << "\thost pointer for factors: " << a.m_pdAmpFactors << endl; + report( DEBUG, kModule ) << "\tsize of amplitude array: " << m_iAmpArrSize << " bytes" << endl; + report( DEBUG, kModule ) << "\tsize of factor array: " << m_iGDoubleFactArrSize << " bytes" << endl; + + cudaDeviceSynchronize(); + gpuErrChk( cudaMemcpy( a.m_pdAmps, m_pfDevAmps, m_iAmpArrSize, cudaMemcpyDeviceToHost ) ); - gpuErrChk( cudaMemcpy( a.m_pdAmpFactors, m_pcDevCalcAmp, + gpuErrChk( cudaMemcpy( a.m_pdAmpFactors, m_pcDevAmpFact, m_iGDoubleFactArrSize, cudaMemcpyDeviceToHost ) ); #ifdef SCOREP @@ -393,9 +441,10 @@ GPUManager::copyAmpsFromGPU( AmpVecs& a ) } void -GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long offset, +GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, const vector< vector< int > >* pvPermutations, - unsigned long iUserVarsOffset ) + unsigned int userVarsOffset, + unsigned int startEvent ) { #ifdef SCOREP SCOREP_USER_REGION_DEFINE( calcAmplitudeAll_gpuMgr ) @@ -403,7 +452,9 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long offset, #endif dim3 dimBlock( m_iDimThreadX, m_iDimThreadY ); - dim3 dimGrid( m_iDimGridX, m_iDimGridY ); + dim3 dimGrid( m_iDimGridXAmpFact, m_iDimGridYAmpFact ); + + report( DEBUG, kModule ) << "Calculating amplitude factors for amplitude " << amp->name() << endl; // do the computation for all events for each permutation in the // vector of permunations @@ -412,8 +463,8 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long offset, // if this is not true, AmplitudeManager hasn't been setup properly assert( permItr->size() == m_iNParticles ); - unsigned long udLocalOffset = 0; - unsigned long permOffset = 0; + unsigned int udLocalOffset = 0; + unsigned int permOffset = 0; for( ; permItr != pvPermutations->end(); ++permItr ){ // copy the permutation to global memory @@ -426,9 +477,9 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long offset, // operation of both real and complex parts at once amp->calcAmplitudeGPU( dimGrid, dimBlock, m_pfDevData, - &m_pfDevUserVars[iUserVarsOffset+udLocalOffset], - (WCUComplex*)&m_pcDevCalcAmp[offset+permOffset], - m_piDevPerm, m_iNParticles, m_iNEvents, + &m_pfDevUserVars[userVarsOffset+udLocalOffset], + (WCUComplex*)&m_pcDevAmpFact[uAmpFactOffset+permOffset], + m_piDevPerm, m_iNParticles, m_iNEvents, startEvent, *permItr ); // check to be sure kernel execution was OK @@ -451,24 +502,30 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long offset, } void -GPUManager::assembleTerms( int iAmpInd, int nFact, int nPerm ){ +GPUManager::assembleTerms( int iAmpInd, int nFact, int nPerm, unsigned int nEvents ){ #ifdef SCOREP SCOREP_USER_REGION_DEFINE( assembleTerms ) SCOREP_USER_REGION_BEGIN( assembleTerms, "assembleTerms", SCOREP_USER_REGION_TYPE_COMMON ) #endif dim3 dimBlock( m_iDimThreadX, m_iDimThreadY ); - dim3 dimGrid( m_iDimGridX, m_iDimGridY ); + dim3 dimGrid( m_iDimGridXAmpFact, m_iDimGridYAmpFact ); - gpuErrChk( cudaMemset( &(m_pfDevAmps[2*m_iNEvents*iAmpInd]), 0, - 2*sizeof(GDouble)*m_iNEvents ) ); - - GPU_ExecFactPermKernel( dimGrid, dimBlock, &(m_pfDevAmps[2*m_iNEvents*iAmpInd]), - m_pcDevCalcAmp, nFact, nPerm, m_iNEvents ); + gpuErrChk( cudaMemset( &(m_pfDevAmps[2*nEvents*iAmpInd]), 0, + 2*sizeof(GDouble)*nEvents ) ); + + report( DEBUG, kModule ) << "Assembling terms for amplitude index " << iAmpInd + << " with " << nFact << " factors and " << nPerm << " permutations for " + << nEvents << " events." << endl; + report( DEBUG, kModule ) << "\tdevice pointer for amplitude: " << &(m_pfDevAmps[2*nEvents*iAmpInd]) << endl; + report( DEBUG, kModule ) << "\tdevice pointer for factors: " << m_pcDevAmpFact << endl; + + GPU_ExecFactPermKernel( dimGrid, dimBlock, &(m_pfDevAmps[2*nEvents*iAmpInd]), + m_pcDevAmpFact, nFact, nPerm, nEvents ); // check to be sure kernel execution was OK - cudaError_t cerrKernel=cudaGetLastError(); - if( cerrKernel!= cudaSuccess ){ + cudaError_t cerrKernel = cudaGetLastError(); + if( cerrKernel != cudaSuccess ){ report( ERROR, kModule ) << "\nKERNEL LAUNCH ERROR [GPU_ExecFactPermKernel]: " << cudaGetErrorString( cerrKernel ) << endl; @@ -482,7 +539,7 @@ GPUManager::assembleTerms( int iAmpInd, int nFact, int nPerm ){ double GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef, - const vector< vector< bool > >& cohMtx ) + const vector< vector< bool > >& cohMtx ) { #ifdef SCOREP SCOREP_USER_REGION_DEFINE( calcSumLogIntensity_gpuMgr ) @@ -511,11 +568,16 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef, gpuErrChk( cudaMemcpy( m_pfDevVVStar, m_pfVVStar, m_iVArrSize, cudaMemcpyHostToDevice ) ); + // note here that this will never be called in a "chunked" calculation + // because the only thing that is calculated in chunks is is the integral + // or sum of the amplitudes, so no special consideration is needed here + // we are guaranteed that pfDevAmps is sized for the entire data set + // compute the logs of the intensities dim3 dimBlock( m_iDimThreadX, m_iDimThreadY ); dim3 dimGrid( m_iDimGridX, m_iDimGridY ); GPU_ExecAmpKernel( dimGrid, dimBlock, m_pfDevAmps, - m_pfDevVVStar, m_pfDevWeights, + m_pfDevVVStar, m_pfDevWeights, m_iNAmps, m_iNEvents, m_pdDevRes ); cudaError_t cerrKernel = cudaGetLastError(); @@ -535,8 +597,8 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef, if( m_iNTrueEvents <= m_iNBlocks ) { gpuErrChk( cudaMemcpy( m_pdRes, m_pdDevRes, - m_iDoubleTrueEventArrSize,cudaMemcpyDeviceToHost ) ); - for(i=0; i >& prodCoef, // zeroing out the padding as not to alter the results after the reduction gpuErrChk( cudaMemset( m_pdDevRes+m_iNTrueEvents, 0, sizeof(double)*(m_iNEvents-m_iNTrueEvents)) ); + // execute the kernel to sum partial sums from each block on CPU - reduce( m_iNEvents, m_iNThreads, m_iNBlocks, - m_pdDevRes, m_pdDevREDUCE ); + reduce( m_iNEvents, m_iNThreads, m_iNBlocks, m_pdDevRes, m_pdDevREDUCE ); cerrKernel = cudaGetLastError(); if( cerrKernel!= cudaSuccess ){ @@ -575,14 +637,15 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef, void GPUManager::calcIntegrals( double* result, int nElements, const vector& iIndex, - const vector& jIndex ){ + const vector& jIndex, + unsigned int startEvent, unsigned int nEvents ){ unsigned int resultSize = 2*sizeof(double)*nElements; unsigned int indexSize = sizeof(int)*nElements; unsigned int totalSize = resultSize + 2*indexSize; dim3 dimBlock( m_iDimThreadX, m_iDimThreadY ); - dim3 dimGrid( m_iDimGridX, m_iDimGridY ); + dim3 dimGrid( m_iDimGridXAmpFact, m_iDimGridYAmpFact ); // this is where the index arrays start on the device int* piDevIndex = (int*)&m_pdDevNICalc[2*nElements]; @@ -607,22 +670,25 @@ GPUManager::calcIntegrals( double* result, int nElements, exit( 1 ); } - GPU_ExecNICalcKernel( dimGrid, dimBlock, totalSize, nElements, - m_pdDevNICalc, m_pfDevAmps, m_pfDevWeights, - m_iNEvents, m_iNTrueEvents, - m_devProp_major >= 7 ? m_maxShared_bytes : 0 ); + // don't bother to execute the kernel if the answer is going to be zero in the end: + if( startEvent < m_iNTrueEvents ){ + + GPU_ExecNICalcKernel( dimGrid, dimBlock, totalSize, nElements, + m_pdDevNICalc, m_pfDevAmps, m_pfDevWeights, + startEvent, nEvents, m_iNTrueEvents, + m_devProp_major >= 7 ? m_maxShared_bytes : 0 ); - // check to be sure kernel execution was OK - cudaError_t cerrKernel = cudaGetLastError(); - if( cerrKernel != cudaSuccess ){ - - report( ERROR, kModule ) << "\nKERNEL LAUNCH ERROR [GPU_ExecNICalcKernel]: " - << cudaGetErrorString( cerrKernel ) << endl; - assert( false ); + // check to be sure kernel execution was OK + cudaError_t cerrKernel = cudaGetLastError(); + if( cerrKernel != cudaSuccess ){ + + report( ERROR, kModule ) << "\nKERNEL LAUNCH ERROR [GPU_ExecNICalcKernel]: " + << cudaGetErrorString( cerrKernel ) << endl; + assert( false ); + } } - gpuErrChk( cudaMemcpy( result, m_pdDevNICalc, - resultSize, cudaMemcpyDeviceToHost ) ); + gpuErrChk( cudaMemcpy( result, m_pdDevNICalc, resultSize, cudaMemcpyDeviceToHost ) ); } // Methods to clear memory: @@ -642,22 +708,15 @@ GPUManager::clearData(){ m_iNTrueEvents=0; m_iGDoubleDataArrSize=0; - m_iDoubleEventArrSize=0; - m_iDoubleTrueEventArrSize=0; if(m_pfDevData) cudaFree(m_pfDevData); m_pfDevData=0; - if(m_pfDevUserVars) - cudaFree(m_pfDevUserVars); - m_pfDevUserVars=0; - if(m_pfDevWeights) cudaFree(m_pfDevWeights); m_pfDevWeights=0; - //CUDA Thread and Grid sizes m_iDimGridX=0; m_iDimGridY=0; @@ -668,9 +727,14 @@ GPUManager::clearData(){ void GPUManager::clearTerms(){ - m_iNAmps=0; + m_iNAmps=0; + m_iNUserVars=0; + + m_iDoubleIntenArrSize=0; + m_iGDoubleFactArrSize=0; m_iAmpArrSize=0; m_iVArrSize=0; + m_iNICalcSize=0; //Host Memory @@ -686,9 +750,13 @@ GPUManager::clearTerms(){ //Device Memory - if(m_pcDevCalcAmp) - cudaFree(m_pcDevCalcAmp); - m_pcDevCalcAmp=0; + if(m_pfDevUserVars) + cudaFree(m_pfDevUserVars); + m_pfDevUserVars=0; + + if(m_pcDevAmpFact) + cudaFree(m_pcDevAmpFact); + m_pcDevAmpFact=0; if(m_piDevPerm) cudaFree(m_piDevPerm); @@ -765,3 +833,25 @@ void GPUManager::calcCUDADims() report( DEBUG, kModule ) << "\tNumber of blocks: "<(sqrt(iNBlocksAmpFact));iDivLo>=1;iDivLo--) + { + iDivHi=iNBlocksAmpFact/iDivLo; + if(iDivLo*iDivHi==iNBlocksAmpFact) + break; + } + m_iDimGridXAmpFact = iDivLo; + m_iDimGridYAmpFact = iDivHi; + } +} diff --git a/AmpTools/GPUManager/GPUManager.h b/AmpTools/GPUManager/GPUManager.h index ea5088bd..ad84566d 100644 --- a/AmpTools/GPUManager/GPUManager.h +++ b/AmpTools/GPUManager/GPUManager.h @@ -46,20 +46,12 @@ #include #include "GPUManager/GPUCustomTypes.h" +#include "IUAmpTools/report.h" #include "cuda_runtime.h" #define gpuErrChk(ans) { gpuAssert((ans), __FILE__, __LINE__); } -inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true) -{ - if (code != cudaSuccess) - { - fprintf(stderr,"GPU ERROR: %s %s %d\n", cudaGetErrorString(code), file, line); - if (abort) exit(code); - } -} - using namespace std; // using namespace __gnu_cxx; @@ -84,18 +76,17 @@ class GPUManager void init( const AmpVecs& a, bool use4Vectors = true ); void initData( const AmpVecs& a, bool use4Vectors = true ); - void initTerms( const AmpVecs& a, unsigned long chunkSize = 0 ); + void initTerms( const AmpVecs& a, unsigned int chunkSize = 0 ); - // Interface Utils - // First Amplitude calculation interface void copyDataToGPU( const AmpVecs& a, bool use4Vectors = true ); void copyUserVarsToGPU( const AmpVecs& a ); - void calcAmplitudeAll( const Amplitude* amp, unsigned long offset, + void calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, const vector< vector< int > >* pvPermutations, - unsigned long userVarsOffset ); + unsigned int userVarsOffset, + unsigned int startEvent ); - void assembleTerms( int iAmpInd, int nFact, int nPerm ); + void assembleTerms( int iAmpInd, int nFact, int nPerm, unsigned int nEvents ); void copyAmpsFromGPU( AmpVecs& a ); @@ -103,8 +94,9 @@ class GPUManager double calcSumLogIntensity( const vector< complex< double > >& prodCoef, const vector< vector< bool > >& cohMtx ); - void calcIntegrals( double* result, int nElements, const vector& iIndex, - const vector& jIndex ); + void calcIntegrals( double* result, int nElements, + const vector& iIndex, const vector& jIndex, + unsigned int startEvent, unsigned int nEvents ); // General utils: static int calcNEventsGPU( int iNEvents ){ @@ -121,19 +113,21 @@ class GPUManager // array dimensions unsigned int m_iNParticles; - unsigned long m_iNEvents; - unsigned long m_iNTrueEvents; + unsigned int m_iNEvents; + unsigned int m_iNTrueEvents; unsigned int m_iNAmps; unsigned int m_iNUserVars; - // array sizes - unsigned long m_iGDoubleDataArrSize; - unsigned long m_iDoubleEventArrSize; - unsigned long m_iDoubleTrueEventArrSize; - unsigned long m_iAmpArrSize; + // array sizes (data) + unsigned int m_iGDoubleDataArrSize; + + // array sizes (terms) + unsigned int m_iDoubleIntenArrSize; + unsigned int m_iGDoubleFactArrSize; + unsigned int m_iAmpArrSize; unsigned int m_iVArrSize; unsigned int m_iNICalcSize; - + //Host Arrays GDouble* m_pfVVStar; double* m_pdRes; @@ -142,7 +136,7 @@ class GPUManager GDouble* m_pfDevData; GDouble* m_pfDevUserVars; GDouble* m_pfDevWeights; - GDouble* m_pcDevCalcAmp; + GDouble* m_pcDevAmpFact; int* m_piDevPerm; GDouble* m_pfDevAmps; @@ -155,10 +149,12 @@ class GPUManager // CUDA Thread and Grid sizes unsigned int m_iDimGridX; unsigned int m_iDimGridY; + unsigned int m_iDimGridXAmpFact; + unsigned int m_iDimGridYAmpFact; unsigned int m_iDimThreadX; unsigned int m_iDimThreadY; - unsigned int m_iNBlocks; + unsigned int m_iNBlocks; unsigned int m_iNThreads; // Internal Utils @@ -168,8 +164,20 @@ class GPUManager unsigned int m_maxShared_bytes; void calcCUDADims(); - + void calcCUDADimsAmpFact( unsigned int chunkSize = 0 ); + static const char* kModule; }; +inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true) +{ + if (code != cudaSuccess) { + + report( ERROR, "GPUManager" ) << "CUDA ERROR: " << cudaGetErrorString( code ) << " in " + << file << " at line " << line << endl; + + if (abort) exit(code); + } +} + #endif //__GPU_MANAGER_H__ From fb8fa68f5b6d2ee6b663a75894862f92000f8859 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:24:21 -0400 Subject: [PATCH 18/32] modify arguments to calcIntegrals for generated MC --- AmpTools/IUAmpTools/NormIntInterface.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc index 82b72c35..7208314d 100644 --- a/AmpTools/IUAmpTools/NormIntInterface.cc +++ b/AmpTools/IUAmpTools/NormIntInterface.cc @@ -436,7 +436,7 @@ NormIntInterface::forceCacheUpdate( bool normIntOnly ) const // do "lazy" allocation of memory here -- this is important for MPI jobs // where forceCacheUpdate is only called on follower nodes, as it // avoids big memory allocations on the lead nodes - if( m_genMCVecs.m_iNTerms == 0 ) m_genMCVecs.allocateTerms( *m_pIntenManager, chunkSize ); + if( m_genMCVecs.m_iNTerms == 0 ) m_genMCVecs.allocateTerms( *m_pIntenManager, false, chunkSize ); report( DEBUG, kModule ) << "Asking IntensityManager to calculate integrals " << "using the generated MC." << endl; @@ -570,11 +570,11 @@ NormIntInterface::genMCChunkSize() const { // get the chunk size by doing integer division by 2 // until the size of the generated MC is less than the - // accepted MC -- for GPU fits this should return a + // 1/2 of the accepted MC -- for GPU fits this should return a // chunk size that is a power of 2, which is required int iPow = 0; - while( ( nGen >> iPow ) > nAcc ) ++iPow; + while( ( nGen >> iPow ) > nAcc/4 ) ++iPow; report( DEBUG, kModule ) << "Chunk size for generated MC: " << ( nGen >> iPow ) << endl; From 44c128de8120efa475070b742f64781aec34eb4b Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:25:08 -0400 Subject: [PATCH 19/32] revisions to manage GPU data --- AmpTools/IUAmpTools/AmpVecs.cc | 38 ++++++++++++++++++---------------- AmpTools/IUAmpTools/AmpVecs.h | 7 +++++-- 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/AmpTools/IUAmpTools/AmpVecs.cc b/AmpTools/IUAmpTools/AmpVecs.cc index a8e04ebe..ecf349f3 100644 --- a/AmpTools/IUAmpTools/AmpVecs.cc +++ b/AmpTools/IUAmpTools/AmpVecs.cc @@ -199,10 +199,10 @@ AmpVecs::clearFourVecs(){ void AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent, - unsigned int iNTrueEvents ){ + unsigned int iNTrueEvents, bool needsUserVarsOnly ){ // allocate memory and set variables - // if this is the first call to this method + // if this is the first call to this method if (m_pdData == NULL){ @@ -218,6 +218,10 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent, m_pdData = new GDouble[4*m_iNParticles*m_iNEvents]; m_pdWeights = new GDouble[m_iNEvents]; + +#ifdef GPU_ACCELERATION + m_gpuMan.initData( *this, !needsUserVarsOnly ); +#endif } // check to be sure we won't exceed the bounds of the array @@ -232,6 +236,7 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent, m_pdWeights[iEvent] = pKinematics->weight(); + m_termsValid = false; m_integralValid = false; m_dataLoaded = true; @@ -240,7 +245,7 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent, void -AmpVecs::loadData( DataReader* pDataReader ){ +AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly ){ // Make sure no data is already loaded @@ -274,7 +279,7 @@ AmpVecs::loadData( DataReader* pDataReader ){ Kinematics* pKinematics; for(unsigned int iEvent = 0; iEvent < m_iNTrueEvents; iEvent++){ pKinematics = pDataReader->getEvent(); - loadEvent(pKinematics, iEvent, m_iNTrueEvents ); + loadEvent(pKinematics, iEvent, m_iNTrueEvents, needsUserVarsOnly ); float weight = pKinematics->weight(); @@ -293,19 +298,18 @@ AmpVecs::loadData( DataReader* pDataReader ){ // Fill any remaining space in the data array with the last event's kinematics for (unsigned int iEvent = m_iNTrueEvents; iEvent < m_iNEvents; iEvent++){ - loadEvent(pKinematics, iEvent, m_iNTrueEvents ); + loadEvent(pKinematics, iEvent, m_iNTrueEvents, needsUserVarsOnly ); } if( m_iNTrueEvents ) delete pKinematics; - + #ifdef GPU_ACCELERATION - - m_gpuMan.initData( *this, !intenMan.needsUserVarsOnly() ); - m_gpuMan.copyDataToGPU( *this, !intenMan.needsUserVarsOnly() ); - -#endif - + + m_gpuMan.copyDataToGPU( *this, !needsUserVarsOnly ); + +#endif + m_termsValid = false; m_integralValid = false; m_dataLoaded = true; @@ -362,7 +366,7 @@ AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity, #else - m_gpuMan.initTerms( *this, ampsEvents ); + m_gpuMan.initTerms( *this, chunkSize ); #endif // GPU_ACCELERATION @@ -372,16 +376,14 @@ AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity, #ifdef GPU_ACCELERATION void -AmpVecs::allocateCPUAmpStorage( const IntensityManager& intenMan, unsigned int chunkSize ){ +AmpVecs::allocateCPUAmpStorage( const IntensityManager& intenMan ){ - unsigned int iNEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize ); - // we should start with unallocated memory assert( m_pdAmps == NULL && m_pdAmpFactors == NULL ); // allocate as "pinned memory" for fast CPU<->GPU memcopies - cudaMallocHost( (void**)&m_pdAmps, iNEvents * intenMan.termStoragePerEvent() * sizeof(GDouble) ); - cudaMallocHost( (void**)&m_pdAmpFactors, iNEvents * m_maxFactPerEvent * sizeof(GDouble)); + cudaMallocHost( (void**)&m_pdAmps, m_iNEvents * intenMan.termStoragePerEvent() * sizeof(GDouble) ); + cudaMallocHost( (void**)&m_pdAmpFactors, 2 * m_iNEvents * m_maxFactPerEvent * sizeof(GDouble)); cudaError_t cudaErr = cudaGetLastError(); if( cudaErr != cudaSuccess ){ diff --git a/AmpTools/IUAmpTools/AmpVecs.h b/AmpTools/IUAmpTools/AmpVecs.h index 646ad6a7..71195592 100644 --- a/AmpTools/IUAmpTools/AmpVecs.h +++ b/AmpTools/IUAmpTools/AmpVecs.h @@ -245,6 +245,7 @@ struct AmpVecs * * \param[in] intenMan a reference to the IntensityManager to use as a guide * for allocating the arrays for terms, factors, and intensities + * */ void allocateCPUAmpStorage( const IntensityManager& intenMan ); #endif @@ -256,11 +257,13 @@ struct AmpVecs * a null pointer is provided (signaling the end of the source). * * \param[in] pDataReader a pointer to a user-defined data reader + * \param[in] needsUserVarsOnly if set to true, this will allow memory + * saving by only copying user variables * * \see DataReader::resetSource * \see DataReader::getEvent */ - void loadData( DataReader* pDataReader ); + void loadData( DataReader* pDataReader, bool needsUserVarsOnly = false ); /** * This routine fills the arrays of data and weights event by event @@ -277,7 +280,7 @@ struct AmpVecs * \see loadData */ void loadEvent( const Kinematics* pKinematics, unsigned int iEvent = 0, - unsigned int iNTrueEvents = 1 ); + unsigned int iNTrueEvents = 1, bool needsUserVarsOnly = false ); /** * A helper routine to get an event i from the array of data and weights. From 08644b224e7ac50e5c9676e97ef412b429e871bd Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:26:27 -0400 Subject: [PATCH 20/32] clean up screen printing --- AmpTools/IUAmpTools/AmpToolsInterface.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/AmpTools/IUAmpTools/AmpToolsInterface.cc b/AmpTools/IUAmpTools/AmpToolsInterface.cc index b49c620f..4ae72a43 100644 --- a/AmpTools/IUAmpTools/AmpToolsInterface.cc +++ b/AmpTools/IUAmpTools/AmpToolsInterface.cc @@ -830,8 +830,8 @@ AmpToolsInterface::printAmplitudes(string reactionName, Kinematics* kin) const { vector< const Amplitude* > ampFactors = ampMan->getFactors(ampNames[iamp]); vector > permutations = ampMan->getPermutations(ampNames[iamp]); - report( INFO, kModule ) << " PRODUCT OF FACTORS" << endl - << " SUMMED OVER PERMUTATIONS = ( " + report( INFO, kModule ) << " PRODUCT OF FACTORS" << endl; + report( INFO, kModule ) << " SUMMED OVER PERMUTATIONS = ( " << aVecs.m_pdAmps[iamp*2] << ", " << aVecs.m_pdAmps[iamp*2+1] << " )" << endl << endl; @@ -845,7 +845,7 @@ AmpToolsInterface::printAmplitudes(string reactionName, Kinematics* kin) const { for (unsigned int iperm = 0; iperm < nPerm; iperm++){ - report( INFO, kModule ) << " PERMUTATION = "; + report( INFO, kModule ) << " PERMUTATION = "; for (unsigned int ipar = 0; ipar < permutations[iperm].size(); ipar++){ report( INFO, kModule ) << permutations[iperm][ipar] << " "; } From c7203cfc418d788d55c6ed23148fa681f8f63184 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:27:12 -0400 Subject: [PATCH 21/32] change some variable names and fix arguments to assembleTerms --- AmpTools/IUAmpTools/AmplitudeManager.cc | 41 +++++++++++++++---------- 1 file changed, 25 insertions(+), 16 deletions(-) diff --git a/AmpTools/IUAmpTools/AmplitudeManager.cc b/AmpTools/IUAmpTools/AmplitudeManager.cc index 86be2a0a..4d83657d 100644 --- a/AmpTools/IUAmpTools/AmplitudeManager.cc +++ b/AmpTools/IUAmpTools/AmplitudeManager.cc @@ -522,15 +522,15 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON // calculate all the factors that make up an amplitude for // for all events serially on CPU or in parallel on GPU - unsigned int iLocalOffset = 0; + unsigned int uAmpFactOffset = 0; for( iFactor=0; iFactor < iNFactors; - iFactor++, iLocalOffset += 2 * nEvents * iNPermutations ){ + iFactor++, uAmpFactOffset += 2 * nEvents * iNPermutations ){ pCurrAmp = vAmps.at( iFactor ); // if we have static user data, look up the location in the data array // if not, then look up by identifier - unsigned int uOffset = + unsigned int userVarsOffset = ( pCurrAmp->areUserVarsStatic() ? a.m_userVarsOffset[pCurrAmp->name()] : a.m_userVarsOffset[pCurrAmp->identifier()] ); @@ -538,14 +538,14 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON #ifndef GPU_ACCELERATION pCurrAmp-> calcAmplitudeAll( a.m_pdData, - a.m_pdAmpFactors + iLocalOffset, + a.m_pdAmpFactors + uAmpFactOffset, a.m_iNEvents, &vvPermuations, - a.m_pdUserVars + uOffset, + a.m_pdUserVars + userVarsOffset, startEvent, nEvents ); #else - a.m_gpuMan.calcAmplitudeAll( pCurrAmp, iLocalOffset, + a.m_gpuMan.calcAmplitudeAll( pCurrAmp, uAmpFactOffset, &vvPermuations, - uOffset ); + userVarsOffset, startEvent ); #endif//GPU_ACCELERATION } @@ -602,7 +602,7 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON // on the GPU the terms are assembled and never copied out // of GPU memory - a.m_gpuMan.assembleTerms( iAmpIndex, iNFactors, iNPermutations ); + a.m_gpuMan.assembleTerms( iAmpIndex, iNFactors, iNPermutations, nEvents ); #endif } @@ -957,14 +957,23 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T } #else - - // this needs to be fixed: needs a temporary vector for the chunk result - // and the calcIntegral GPU method needs to be "chunked" - - // use the GPU manager to compute the result - a.m_gpuMan.calcIntegrals( &(result[0]), nCompute, iIndex, jIndex ); - - // result += tempResult + + vector temp(maxNIElements*2); + + // use the GPU manager to compute the result for the chunk + a.m_gpuMan.calcIntegrals( &(temp[0]), nCompute, iIndex, jIndex, + startEvent, nEvents ); + + // add the to the result array that is accumulting the sum + for( int iTerm = 0; iTerm < nCompute; ++iTerm ){ + result[2*iTerm] += temp[2*iTerm]; + result[2*iTerm+1] += temp[2*iTerm+1]; + + report( DEBUG, kModule ) << "NI Term Sum ( " << iIndex[iTerm] << ", " + << jIndex[iTerm] << " ) = ( " + << result[2*iTerm] << ", " << result[2*iTerm+1] + << ")" << endl; + } #endif } From f2ede93f4272a0c461cc7034edb58a7ef28486e5 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Thu, 18 Jun 2026 15:27:41 -0400 Subject: [PATCH 22/32] change type to unsigned int --- AmpTools/GPUManager/GPUKernel.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/AmpTools/GPUManager/GPUKernel.h b/AmpTools/GPUManager/GPUKernel.h index b7071995..44116c18 100644 --- a/AmpTools/GPUManager/GPUKernel.h +++ b/AmpTools/GPUManager/GPUKernel.h @@ -41,16 +41,16 @@ extern "C" void GPU_ExecAmpKernel(dim3 dimGrid,dim3 dimBlock,GDouble* pfDevAmps, GDouble* pfDevVVStar, GDouble* pfDevWeights, - int nAmps, int nEvents, double* pfDevRes); + int nAmps, unsigned int nEvents, double* pfDevRes); extern "C" void GPU_ExecFactPermKernel( dim3 dimGrid, dim3 dimBlock, - GDouble* pfDevAmps, GDouble* pcDevCalcAmp, - int nFact, int nPerm, int nEvents ); - + GDouble* pfDevAmps, GDouble* pcDevAmpFact, + int nFact, int nPerm, unsigned int nEvents ); extern "C" void GPU_ExecNICalcKernel( dim3 dimGrid, dim3 dimBlock, unsigned int sharedSize, int nElements, double* pdDevNICalc, GDouble* pfDevAmps, GDouble* pfDevWeights, - int nEvents, int nTrueEvents, unsigned int maxSize ); + unsigned int startEvent, unsigned int nEvents, + unsigned int nTrueEvents, unsigned int maxSize ); #endif //__GPU_KERNEL__H__ From 0bd90bde5dad3d5b1321ddd40b5b0e21dba82ca4 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Fri, 19 Jun 2026 13:29:27 -0400 Subject: [PATCH 23/32] modify calcAmplitudeAll method in GPU manager to properly use the chunk size to increment the factor arrays; improve some debugging outputs --- AmpTools/GPUManager/GPUManager.cc | 31 +++++++++++-------------- AmpTools/GPUManager/GPUManager.h | 2 +- AmpTools/IUAmpTools/Amplitude.cc | 2 +- AmpTools/IUAmpTools/AmplitudeManager.cc | 10 ++++---- 4 files changed, 21 insertions(+), 24 deletions(-) diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc index 7fe0299a..8e01ad5e 100644 --- a/AmpTools/GPUManager/GPUManager.cc +++ b/AmpTools/GPUManager/GPUManager.cc @@ -290,10 +290,15 @@ GPUManager::initTerms( const AmpVecs& a, unsigned int chunkSize ) << m_iAmpArrSize/sizeof(GDouble) << " elements )" << endl; report( DEBUG, kModule ) << "\tsize of factor array: " << m_iGDoubleFactArrSize / 1024 << " kB ( " << m_iGDoubleFactArrSize/sizeof(GDouble) << " elements )" << endl; - report( DEBUG, kModule ) << "\tsize of VV* array: " << m_iVArrSize / 1024 << " kB ( " + report( DEBUG, kModule ) << "\tsize of VV* array: " << m_iVArrSize << " B ( " << m_iVArrSize/sizeof(GDouble) << " elements )" << endl; - report( DEBUG, kModule ) << "\tsize of NICalc array: " << m_iNICalcSize / 1024 << " kB ( " + report( DEBUG, kModule ) << "\tsize of NICalc array: " << m_iNICalcSize << " B ( " << m_iNICalcSize/sizeof(double) << " elements )" << endl; + report( DEBUG, kModule ) << "\tsize of permutation array: " << m_iNParticles * sizeof( int ) << " B ( " + << m_iNParticles << " elements )" << endl; + + // allocate device memory needed for amplitude calculations + // ... and use the chunkSize for these (iNAmpEvents) when it is provided // device memory needed for intensity or integral calculation and sum gpuErrChk( cudaMalloc( (void**)&m_pfDevUserVars, m_iNUserVars * m_iGDoubleDataArrSize ) ) ; @@ -419,14 +424,6 @@ GPUManager::copyAmpsFromGPU( AmpVecs& a ) report( DEBUG, kModule ) << "Copying amplitudes and factors from GPU to CPU for " << m_iNAmps << " amplitudes and " << m_iNEvents << " events." << endl; - report( DEBUG, kModule ) << "\tdevice pointer for amplitude: " << m_pfDevAmps << endl; - report( DEBUG, kModule ) << "\tdevice pointer for factors: " << m_pcDevAmpFact << endl; - report( DEBUG, kModule ) << "\thost pointer for amplitude: " << a.m_pdAmps << endl; - report( DEBUG, kModule ) << "\thost pointer for factors: " << a.m_pdAmpFactors << endl; - report( DEBUG, kModule ) << "\tsize of amplitude array: " << m_iAmpArrSize << " bytes" << endl; - report( DEBUG, kModule ) << "\tsize of factor array: " << m_iGDoubleFactArrSize << " bytes" << endl; - - cudaDeviceSynchronize(); gpuErrChk( cudaMemcpy( a.m_pdAmps, m_pfDevAmps, m_iAmpArrSize, @@ -444,7 +441,7 @@ void GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, const vector< vector< int > >* pvPermutations, unsigned int userVarsOffset, - unsigned int startEvent ) + unsigned int startEvent, unsigned int chunkSize ) { #ifdef SCOREP SCOREP_USER_REGION_DEFINE( calcAmplitudeAll_gpuMgr ) @@ -454,7 +451,9 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, dim3 dimBlock( m_iDimThreadX, m_iDimThreadY ); dim3 dimGrid( m_iDimGridXAmpFact, m_iDimGridYAmpFact ); - report( DEBUG, kModule ) << "Calculating amplitude factors for amplitude " << amp->name() << endl; + unsigned int nEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize ); + + report( DEBUG, kModule ) << "Calculating amplitude factors for amplitude " << amp->identifier() << endl; // do the computation for all events for each permutation in the // vector of permunations @@ -466,7 +465,7 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, unsigned int udLocalOffset = 0; unsigned int permOffset = 0; for( ; permItr != pvPermutations->end(); ++permItr ){ - + // copy the permutation to global memory gpuErrChk( cudaMemcpy( m_piDevPerm, &((*permItr)[0]), m_iNParticles * sizeof( int ), @@ -493,8 +492,8 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, // increment the offset so that we place the computation for the // next permutation after the previous in pcResAmp - permOffset += 2 * m_iNEvents; - udLocalOffset += amp->numUserVars() * m_iNEvents; + permOffset += 2 * nEvents; + udLocalOffset += amp->numUserVars() * nEvents; } #ifdef SCOREP SCOREP_USER_REGION_END( calcAmplitudeAll_gpuMgr ) @@ -517,8 +516,6 @@ GPUManager::assembleTerms( int iAmpInd, int nFact, int nPerm, unsigned int nEven report( DEBUG, kModule ) << "Assembling terms for amplitude index " << iAmpInd << " with " << nFact << " factors and " << nPerm << " permutations for " << nEvents << " events." << endl; - report( DEBUG, kModule ) << "\tdevice pointer for amplitude: " << &(m_pfDevAmps[2*nEvents*iAmpInd]) << endl; - report( DEBUG, kModule ) << "\tdevice pointer for factors: " << m_pcDevAmpFact << endl; GPU_ExecFactPermKernel( dimGrid, dimBlock, &(m_pfDevAmps[2*nEvents*iAmpInd]), m_pcDevAmpFact, nFact, nPerm, nEvents ); diff --git a/AmpTools/GPUManager/GPUManager.h b/AmpTools/GPUManager/GPUManager.h index ad84566d..36c29984 100644 --- a/AmpTools/GPUManager/GPUManager.h +++ b/AmpTools/GPUManager/GPUManager.h @@ -84,7 +84,7 @@ class GPUManager void calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, const vector< vector< int > >* pvPermutations, unsigned int userVarsOffset, - unsigned int startEvent ); + unsigned int startEvent = 0, unsigned int chunkSize = 0 ); void assembleTerms( int iAmpInd, int nFact, int nPerm, unsigned int nEvents ); diff --git a/AmpTools/IUAmpTools/Amplitude.cc b/AmpTools/IUAmpTools/Amplitude.cc index c4518cd3..8e6fe46f 100644 --- a/AmpTools/IUAmpTools/Amplitude.cc +++ b/AmpTools/IUAmpTools/Amplitude.cc @@ -78,7 +78,7 @@ Amplitude::calcUserVarsAll( GDouble* pdData, GDouble* pdUserVars, int iNEvents, SCOREP_USER_REGION_DEFINE( calcUserVarsAll ) #endif - report( DEBUG, kModule ) << "Caculating user data for " << name() << endl; + report( DEBUG, kModule ) << "Caculating user data for " << identifier() << endl; unsigned int numVars = numUserVars(); diff --git a/AmpTools/IUAmpTools/AmplitudeManager.cc b/AmpTools/IUAmpTools/AmplitudeManager.cc index 4d83657d..4c825367 100644 --- a/AmpTools/IUAmpTools/AmplitudeManager.cc +++ b/AmpTools/IUAmpTools/AmplitudeManager.cc @@ -123,8 +123,9 @@ IntensityManager( reaction, reactionName ) if( m_symmCombos.size() > 1 ){ - report( INFO, kModule ) << "The following " << numberOfCombos << " orderings of the particles are" << endl - << "indistinguishable and will be permuted when computing amplitudes." << endl; + report( INFO, kModule ) << "The following " << numberOfCombos + << " orderings of the particles are indistinguishable" << endl; + report( INFO, kModule ) << "and will be permuted when computing amplitudes." << endl; for( unsigned int i = 0; i < m_symmCombos.size(); ++i ){ @@ -543,9 +544,8 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON a.m_pdUserVars + userVarsOffset, startEvent, nEvents ); #else - a.m_gpuMan.calcAmplitudeAll( pCurrAmp, uAmpFactOffset, - &vvPermuations, - userVarsOffset, startEvent ); + a.m_gpuMan.calcAmplitudeAll( pCurrAmp, uAmpFactOffset, &vvPermuations, + userVarsOffset, startEvent, nEvents ); #endif//GPU_ACCELERATION } From 154386ae6702dba24c9c500d827728c276cd8219 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Tue, 30 Jun 2026 17:30:00 -0400 Subject: [PATCH 24/32] add additional debugging statement --- AmpTools/GPUManager/GPUManager.cc | 3 +++ 1 file changed, 3 insertions(+) diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc index 8e01ad5e..ad984d6a 100644 --- a/AmpTools/GPUManager/GPUManager.cc +++ b/AmpTools/GPUManager/GPUManager.cc @@ -459,6 +459,9 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, // vector of permunations vector< vector< int > >::const_iterator permItr = pvPermutations->begin(); + report( DEBUG, kModule ) << "\tSize of permutation elements: " << permItr->size() << endl; + report( DEBUG, kModule ) << "\tNumber of permutations: " << pvPermutations->size() << endl; + report( DEBUG, kModule ) << "\tNumber of particles: " << m_iNParticles << endl; // if this is not true, AmplitudeManager hasn't been setup properly assert( permItr->size() == m_iNParticles ); From 1ae7ef3514267812983a908117bb37b8f93c5490 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Tue, 30 Jun 2026 17:31:39 -0400 Subject: [PATCH 25/32] fix bugs related to sharing of repeated data sets and getting the data sets copied to the GPU; adjust initData call order to be compatible with GPU operations --- AmpTools/IUAmpTools/AmpVecs.cc | 30 +++++++++++++++++++----------- AmpTools/IUAmpTools/AmpVecs.h | 2 +- 2 files changed, 20 insertions(+), 12 deletions(-) diff --git a/AmpTools/IUAmpTools/AmpVecs.cc b/AmpTools/IUAmpTools/AmpVecs.cc index ecf349f3..448db27c 100644 --- a/AmpTools/IUAmpTools/AmpVecs.cc +++ b/AmpTools/IUAmpTools/AmpVecs.cc @@ -209,19 +209,19 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent, m_iNTrueEvents = iNTrueEvents; m_iNEvents = iNTrueEvents; -#ifdef GPU_ACCELERATION - m_iNEvents = GPUManager::calcNEventsGPU(iNTrueEvents); -#endif - m_iNParticles = pKinematics->particleList().size(); assert(m_iNParticles); - - m_pdData = new GDouble[4*m_iNParticles*m_iNEvents]; - m_pdWeights = new GDouble[m_iNEvents]; -#ifdef GPU_ACCELERATION +#ifdef GPU_ACCELERATION + m_iNEvents = GPUManager::calcNEventsGPU(iNTrueEvents); m_gpuMan.initData( *this, !needsUserVarsOnly ); #endif + + report( DEBUG, kModule ) << "Allocating CPU memory for " << m_iNEvents + << " events with " << m_iNParticles << " particles" << endl; + + m_pdData = new GDouble[4*m_iNParticles*m_iNEvents]; + m_pdWeights = new GDouble[m_iNEvents]; } // check to be sure we won't exceed the bounds of the array @@ -413,7 +413,7 @@ AmpVecs::getEvent( int iEvent ){ } void -AmpVecs::shareDataWith( AmpVecs* targetAmpVecs ){ +AmpVecs::shareDataWith( AmpVecs* targetAmpVecs, bool needsUserVarsOnly ){ // if this object is already using shared data // then it should rely on the host to distribute @@ -421,7 +421,7 @@ AmpVecs::shareDataWith( AmpVecs* targetAmpVecs ){ // the set of shared data friends complete if( m_usesSharedData ){ - m_sharedDataHost->shareDataWith( targetAmpVecs ); + m_sharedDataHost->shareDataWith( targetAmpVecs, needsUserVarsOnly ); return; } @@ -440,7 +440,7 @@ AmpVecs::shareDataWith( AmpVecs* targetAmpVecs ){ // while we are really going to share the four-vectors, // we will give the target a copy of the weights -- this // avoids complications if the four-vectors can later - // be flusehd from memory but the weights need to remain + // be flushed from memory but the weights need to remain // for the fit targetAmpVecs->m_pdWeights = new GDouble[m_iNEvents]; @@ -454,6 +454,14 @@ AmpVecs::shareDataWith( AmpVecs* targetAmpVecs ){ m_sharedDataFriends.insert( targetAmpVecs ); targetAmpVecs->m_usesSharedData = true; targetAmpVecs->m_sharedDataHost = this; + + // even when using shared data the class will have its + // own GPU manager so we need to initialize it and + // copy the data to the GPU for the new object +#ifdef GPU_ACCELERATION + targetAmpVecs->m_gpuMan.initData( *targetAmpVecs, !needsUserVarsOnly ); + targetAmpVecs->m_gpuMan.copyDataToGPU( *targetAmpVecs, !needsUserVarsOnly ); +#endif } void diff --git a/AmpTools/IUAmpTools/AmpVecs.h b/AmpTools/IUAmpTools/AmpVecs.h index 71195592..3a11a5d4 100644 --- a/AmpTools/IUAmpTools/AmpVecs.h +++ b/AmpTools/IUAmpTools/AmpVecs.h @@ -307,7 +307,7 @@ struct AmpVecs * This function will share this classes data four vectors with the * with the targetAmpVecs object specified by the argument. */ - void shareDataWith( AmpVecs* targetAmpVecs ); + void shareDataWith( AmpVecs* targetAmpVecs, bool needsUserVarsOnly = false ); /** * This function will allow this class to take ownership of the From ba09a101dad357b56a7e61ff2a6fff002901a0db Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Tue, 30 Jun 2026 17:36:57 -0400 Subject: [PATCH 26/32] fix memory consuming problem with MPI jobs where lead node loads all data into memory in addition to data being loaded into memory on the follower nodes; use new loadMC method to control single process and MPI loading; change some unsigned long int to unsigned int --- AmpTools/IUAmpTools/NormIntInterface.cc | 103 ++++++++++-------- AmpTools/IUAmpTools/NormIntInterface.h | 10 +- AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc | 22 ++-- 3 files changed, 74 insertions(+), 61 deletions(-) diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc index 7208314d..63eac3f1 100644 --- a/AmpTools/IUAmpTools/NormIntInterface.cc +++ b/AmpTools/IUAmpTools/NormIntInterface.cc @@ -101,41 +101,7 @@ m_termNames( intenManager.getTermNames() ) assert( ( m_accMCReader != NULL ) && ( m_genMCReader != NULL ) ); m_termNames = intenManager.getTermNames(); - - std::map::iterator genVecs = m_uniqueDataSets.find( m_genMCReader ); - if( genVecs == m_uniqueDataSets.end() ){ - - report( INFO, kModule ) << "Loading generated Monte Carlo from file..." << endl; - m_genMCVecs.loadData( m_genMCReader, intenManager.needsUserVarsOnly() ); - - m_uniqueDataSets[m_genMCReader] = &m_genMCVecs; - } - else{ - - report( NOTICE, kModule ) << "Duplicated Monte Carlo set detected, " - << "using previously loaded version" << endl; - - genVecs->second->shareDataWith( &m_genMCVecs ); - } - m_nGenEvents = m_genMCVecs.m_iNTrueEvents; - - std::map::iterator accVecs = m_uniqueDataSets.find( m_accMCReader ); - if( accVecs == m_uniqueDataSets.end() ){ - - report( INFO, kModule ) << "Loading accepted Monte Carlo from file..." << endl; - m_accMCVecs.loadData( m_accMCReader, intenManager.needsUserVarsOnly() ); - - m_uniqueDataSets[m_accMCReader] = &m_accMCVecs; - } - else{ - - report( NOTICE, kModule ) << "Duplicated Monte Carlo set detected, " - << "using previously loaded version" << endl; - - accVecs->second->shareDataWith( &m_accMCVecs ); - } - m_sumAccWeights = m_accMCVecs.m_dSumWeights; - + initializeCache(); } @@ -386,26 +352,29 @@ NormIntInterface::forceCacheUpdate( bool normIntOnly ) const { report( DEBUG, kModule ) << "Update of the NI cache -- normIntOnly = " - << normIntOnly << ", emptyNormIntCache = " << m_emptyNormIntCache - << ", emptyAmpIntCache = " << m_emptyAmpIntCache << endl; - - // if the accepted MC is not available, then the data have likely - // not been loaded into AmpVecs and we can't reclaculate the integrals - // below - assert( m_accMCVecs.m_dataLoaded ); - + << normIntOnly << ", emptyNormIntCache = " << m_emptyNormIntCache + << ", emptyAmpIntCache = " << m_emptyAmpIntCache << endl; + // do "lazy" allocation of memory here -- this is important for MPI jobs // where forceCacheUpdate is only called on follower nodes, as it // avoids big memory allocations on the lead nodes + + // this will load both the accepted and generated MC into memory + // (For MPI jobs, loadMC is called explicitly on the follower nodes + // during the steup step, so this will not be called again on those nodes + // and forceCacheUpdate is never called on the MPI lead node) + if( !m_accMCVecs.m_dataLoaded ) loadMC(); + + // allocate the space for calculating the amplitudes if( m_accMCVecs.m_iNTerms == 0 ) m_accMCVecs.allocateTerms( *m_pIntenManager ); - + // we won't enter here if the cache is empty, which can happen // on the first pass through the data -- for subsequent passes // (during fitting) the loop below will execute and return if( !m_emptyNormIntCache && normIntOnly ){ report( DEBUG, kModule ) << "Asking IntensityManager to update integrals " - << "using the accepted MC." << endl; + << "using the accepted MC." << endl; m_pIntenManager->calcIntegrals( m_accMCVecs, m_nGenEvents ); setNormIntMatrix( m_accMCVecs.m_pdIntegralMatrix ); @@ -558,6 +527,7 @@ NormIntInterface::setNormIntMatrix( const double* input ) const { m_emptyNormIntCache = false; } +#ifndef __ACLIC__ unsigned int NormIntInterface::genMCChunkSize() const { @@ -565,8 +535,8 @@ NormIntInterface::genMCChunkSize() const { // need to know the sizes of the data sets assert( m_accMCVecs.m_dataLoaded && m_genMCVecs.m_dataLoaded ); - unsigned long nGen = m_genMCVecs.m_iNEvents; - unsigned long nAcc = m_accMCVecs.m_iNEvents; + unsigned int nGen = m_genMCVecs.m_iNEvents; + unsigned int nAcc = m_accMCVecs.m_iNEvents; // get the chunk size by doing integer division by 2 // until the size of the generated MC is less than the @@ -581,7 +551,44 @@ NormIntInterface::genMCChunkSize() const { return( nGen >> iPow ); } -#ifndef __ACLIC__ +void +NormIntInterface::loadMC() const { + + std::map::iterator genVecs = m_uniqueDataSets.find( m_genMCReader ); + if( genVecs == m_uniqueDataSets.end() ){ + + report( INFO, kModule ) << "Loading generated Monte Carlo from file..." << endl; + m_genMCVecs.loadData( m_genMCReader, m_pIntenManager->needsUserVarsOnly() ); + + m_uniqueDataSets[m_genMCReader] = &m_genMCVecs; + } + else{ + + report( NOTICE, kModule ) << "Duplicated Monte Carlo set detected, " + << "using previously loaded version" << endl; + + genVecs->second->shareDataWith( &m_genMCVecs, m_pIntenManager->needsUserVarsOnly() ); + } + m_nGenEvents = m_genMCVecs.m_iNTrueEvents; + + std::map::iterator accVecs = m_uniqueDataSets.find( m_accMCReader ); + if( accVecs == m_uniqueDataSets.end() ){ + + report( INFO, kModule ) << "Loading accepted Monte Carlo from file..." << endl; + m_accMCVecs.loadData( m_accMCReader, m_pIntenManager->needsUserVarsOnly() ); + + m_uniqueDataSets[m_accMCReader] = &m_accMCVecs; + } + else{ + + report( NOTICE, kModule ) << "Duplicated Monte Carlo set detected, " + << "using previously loaded version" << endl; + + accVecs->second->shareDataWith( &m_accMCVecs, m_pIntenManager->needsUserVarsOnly() ); + } + m_sumAccWeights = m_accMCVecs.m_dSumWeights; +} + void NormIntInterface::invalidateTerms(){ diff --git a/AmpTools/IUAmpTools/NormIntInterface.h b/AmpTools/IUAmpTools/NormIntInterface.h index 6f58a2e3..dbe87743 100644 --- a/AmpTools/IUAmpTools/NormIntInterface.h +++ b/AmpTools/IUAmpTools/NormIntInterface.h @@ -69,7 +69,7 @@ class NormIntInterface istream& loadNormIntCache( istream& in ); void operator+=( const NormIntInterface& nii ); - unsigned long int numGenEvents() const { return m_nGenEvents; } + unsigned int numGenEvents() const { return m_nGenEvents; } double numAccEvents() const { return m_sumAccWeights; } // this integral folds in detector acceptance @@ -101,7 +101,7 @@ class NormIntInterface const double* ampIntMatrix() const { return m_ampIntCache; } const double* normIntMatrix() const { return m_normIntCache; } - void setGenEvents( unsigned long int events ) { m_nGenEvents = events; } + void setGenEvents( unsigned int events ) { m_nGenEvents = events; } void setAccEvents( double sumWeights ) { m_sumAccWeights = sumWeights; } protected: @@ -116,6 +116,8 @@ class NormIntInterface void setAmpIntMatrix( const double* input ) const; void setNormIntMatrix( const double* input ) const; + + void loadMC() const; private: @@ -133,8 +135,8 @@ class NormIntInterface mutable bool m_emptyNormIntCache; mutable bool m_emptyAmpIntCache; - unsigned long int m_nGenEvents; - double m_sumAccWeights; + mutable unsigned int m_nGenEvents; + mutable double m_sumAccWeights; #ifndef __ACLIC__ diff --git a/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc b/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc index 0a1af4e0..a9201d99 100644 --- a/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc +++ b/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc @@ -96,22 +96,22 @@ NormIntInterfaceMPI::setupMPI() m_isLeader = ( m_rank == 0 ); // this is unsigned elsewhere, but MPI - unsigned long int totalGenEvents = 0; + unsigned int totalGenEvents = 0; double totalAccWeights = 0; if( m_isLeader ){ for( int i = 1; i < m_numProc; ++i ){ - unsigned long int thisEvents; + unsigned int thisEvents; double thisWeights; // trigger sending of events from followers -- data is irrelevant - MPI_Send( &thisEvents, 1, MPI_UNSIGNED_LONG, i, MPITag::kAcknowledge, + MPI_Send( &thisEvents, 1, MPI_UNSIGNED, i, MPITag::kAcknowledge, MPI_COMM_WORLD ); // now receive actual data - MPI_Recv( &thisEvents, 1, MPI_UNSIGNED_LONG, i, MPITag::kLongIntSend, + MPI_Recv( &thisEvents, 1, MPI_UNSIGNED, i, MPITag::kIntSend, MPI_COMM_WORLD, &status ); totalGenEvents += thisEvents; @@ -120,7 +120,7 @@ NormIntInterfaceMPI::setupMPI() totalAccWeights += thisWeights; // send acknowledgment - MPI_Send( &thisEvents, 1, MPI_UNSIGNED_LONG, i, MPITag::kAcknowledge, + MPI_Send( &thisEvents, 1, MPI_UNSIGNED, i, MPITag::kAcknowledge, MPI_COMM_WORLD ); } @@ -129,7 +129,11 @@ NormIntInterfaceMPI::setupMPI() } else{ - unsigned long int thisEvents; + // load the MC into CPU memory on the follower nodes -- this is effectively + // a copy from the DataReaderMPI cache to the AmpVecs structure + loadMC(); + + unsigned int thisEvents; double thisWeights; // if we are not the leader, send generated and accepted events @@ -140,16 +144,16 @@ NormIntInterfaceMPI::setupMPI() // to signal that it is ready to accept numbers of events // data is irrelevant for this receive - MPI_Recv( &thisEvents, 1, MPI_UNSIGNED_LONG, 0, MPITag::kAcknowledge, MPI_COMM_WORLD, + MPI_Recv( &thisEvents, 1, MPI_UNSIGNED, 0, MPITag::kAcknowledge, MPI_COMM_WORLD, &status ); thisEvents = numGenEvents(); - MPI_Send( &thisEvents, 1, MPI_UNSIGNED_LONG, 0, MPITag::kLongIntSend, MPI_COMM_WORLD ); + MPI_Send( &thisEvents, 1, MPI_UNSIGNED, 0, MPITag::kIntSend, MPI_COMM_WORLD ); thisWeights = numAccEvents(); MPI_Send( &thisWeights, 1, MPI_DOUBLE, 0, MPITag::kDoubleSend, MPI_COMM_WORLD ); - MPI_Recv( &thisEvents, 1, MPI_UNSIGNED_LONG, 0, MPITag::kAcknowledge, MPI_COMM_WORLD, + MPI_Recv( &thisEvents, 1, MPI_UNSIGNED, 0, MPITag::kAcknowledge, MPI_COMM_WORLD, &status ); } } From 25f5d933f7f72f6c1c945d05a8c0a11635b88a19 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Tue, 30 Jun 2026 17:37:42 -0400 Subject: [PATCH 27/32] add unit test products to gitignore --- .gitignore | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 3c9079d7..e3a94763 100644 --- a/.gitignore +++ b/.gitignore @@ -30,4 +30,17 @@ Tutorials/Dalitz/doc/*.aux Tutorials/Dalitz/doc/*.log Tutorials/Dalitz/doc/*.toc Tutorials/Dalitz/doc/*.synctex.gz - +UnitTests/AmpToolsInterfaceTest +UnitTests/AmpToolsInterfaceTest_GPU +UnitTests/AmpToolsInterfaceTestMPI +UnitTests/AmpToolsInterfaceTestMPI_GPU +UnitTests/configurationInfoTest +UnitTests/configurationInfoTest_GPU +UnitTests/fitResultsTest +UnitTests/fitResultsTest_GPU +UnitTests/fitResultsTestMPI +UnitTests/fitResultsTestMPI_GPU +UnitTests/parserTest +UnitTests/parserTest_GPU +UnitTests/writeModels +UnitTests/writeModels_GPU From c9094fb490e61d4db5d86c678d953886054904e4 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Tue, 30 Jun 2026 18:41:46 -0400 Subject: [PATCH 28/32] more clear comment --- AmpTools/GPUManager/GPUManager.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc index ad984d6a..043a2eca 100644 --- a/AmpTools/GPUManager/GPUManager.cc +++ b/AmpTools/GPUManager/GPUManager.cc @@ -205,7 +205,7 @@ GPUManager::initData( const AmpVecs& a, bool use4Vectors ) totalMemory /= (1024*1024); report( INFO, kModule ) << "Attempting to allocate " << (int)totalMemory - << " MB of global GPU memory for data." << endl; + << " MB of global GPU memory for event data." << endl; report( DEBUG, kModule ) << "AmpVecs details (initData): " << endl; report( DEBUG, kModule ) << "\tiNEvents: " << m_iNEvents << endl; report( DEBUG, kModule ) << "\tiNTrueEvents: " << m_iNTrueEvents << endl; From f2619884d8a925ce4d796e353f7cba0545bdbebb Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Tue, 30 Jun 2026 18:42:01 -0400 Subject: [PATCH 29/32] improve message printing to screen --- AmpTools/IUAmpTools/LikelihoodCalculator.cc | 31 ++++++++++----------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/AmpTools/IUAmpTools/LikelihoodCalculator.cc b/AmpTools/IUAmpTools/LikelihoodCalculator.cc index ed05a91a..af659398 100644 --- a/AmpTools/IUAmpTools/LikelihoodCalculator.cc +++ b/AmpTools/IUAmpTools/LikelihoodCalculator.cc @@ -251,14 +251,13 @@ SCOREP_USER_REGION_BEGIN( dataTerm, "dataTerm", SCOREP_USER_REGION_TYPE_COMMON ) m_ampVecsBkgnd.allocateTerms( m_intenManager, true ); if( m_ampVecsBkgnd.m_hasMixedSignWeights ){ - report( NOTICE, kModule ) << "\n" - << "***************************************************************\n" - << "* NOTICE: Weights with both positive and negative signs were *\n" - << "* detected in the background file. This may be desirable for *\n" - << "* some applications. Older versions of AmpTools (v0.10.x and *\n" - << "* prior) will not properly handle this case and will also not *\n" - << "* print this notification to the screen. *\n" - << "***************************************************************\n" << endl; + report( NOTICE, kModule ) << "***************************************************************" << endl; + report( NOTICE, kModule ) << "* NOTICE: Weights with both positive and negative signs were *" << endl; + report( NOTICE, kModule ) << "* detected in the background file. This may be desirable for *" << endl; + report( NOTICE, kModule ) << "* some applications. Older versions of AmpTools (v0.10.x and *" << endl; + report( NOTICE, kModule ) << "* prior) will not properly handle this case and will also not *" << endl; + report( NOTICE, kModule ) << "* print this notification to the screen. *" << endl; + report( NOTICE, kModule ) << "***************************************************************" << endl; } m_sumBkgWeights = m_ampVecsBkgnd.m_dSumWeights; @@ -267,14 +266,14 @@ SCOREP_USER_REGION_BEGIN( dataTerm, "dataTerm", SCOREP_USER_REGION_TYPE_COMMON ) // may fail on one of the follower nodes if the background sample // is sparse if( m_sumBkgWeights < 0 && !suppressError ){ - report( ERROR, kModule ) << "\n" - << "****************************************************************\n" - << "* ERROR: The sum of all background weights is negative. This *\n" - << "* implies a negative background in the signal region, which *\n" - << "* unphysical. The weighted sum of the background events *\n" - << "* should represent the background contribution to the signal *\n" - << "* region. *\n" - << "****************************************************************\n" << endl; + report( ERROR, kModule ) << "****************************************************************" << endl; + report( ERROR, kModule ) << "* ERROR: The sum of all background weights is negative. This *" << endl; + report( ERROR, kModule ) << "* implies a negative background in the signal region, which *" << endl; + report( ERROR, kModule ) << "* unphysical. The weighted sum of the background events *" << endl; + report( ERROR, kModule ) << "* should represent the background contribution to the signal *" << endl; + report( ERROR, kModule ) << "* region. *" << endl; + report( ERROR, kModule ) << "****************************************************************" << endl; + assert( false ); } From 9e4bfb20d42319c48944643d0fa8d5aa696f6b9f Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Wed, 8 Jul 2026 15:08:43 -0400 Subject: [PATCH 30/32] add debugging statement --- AmpTools/IUAmpTools/AmplitudeManager.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/AmpTools/IUAmpTools/AmplitudeManager.cc b/AmpTools/IUAmpTools/AmplitudeManager.cc index 4c825367..23fd7ba2 100644 --- a/AmpTools/IUAmpTools/AmplitudeManager.cc +++ b/AmpTools/IUAmpTools/AmplitudeManager.cc @@ -975,7 +975,8 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T << ")" << endl; } #endif - + if( nChunk > 1 ) report( DEBUG, kModule ) << "Integral calculation: chunk " + << iChunk << " of " << nChunk << " complete." << endl; } // in a "chunked" calculation nCompute, iIndex, and jIndex get reset From 5475fefbadf10355aed2dcba8038a778458e7b59 Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Wed, 8 Jul 2026 15:11:16 -0400 Subject: [PATCH 31/32] functionality to improve memory use; GPU manager can share data on the GPU in parallel with AmpVecs data sharing; NormIntInterface cleans up memory allocated for terms after gen MC integral calculation --- AmpTools/GPUManager/GPUManager.cc | 28 ++++++++- AmpTools/GPUManager/GPUManager.h | 3 + AmpTools/IUAmpTools/AmpVecs.cc | 84 +++++++++++++++---------- AmpTools/IUAmpTools/AmpVecs.h | 7 +++ AmpTools/IUAmpTools/NormIntInterface.cc | 6 +- 5 files changed, 91 insertions(+), 37 deletions(-) diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc index 043a2eca..1b21ed55 100644 --- a/AmpTools/GPUManager/GPUManager.cc +++ b/AmpTools/GPUManager/GPUManager.cc @@ -63,7 +63,7 @@ bool GPUManager::m_cudaDisplay = false; template void reduce(int size, int threads, int blocks, T *d_idata, T *d_odata); -GPUManager::GPUManager() +GPUManager::GPUManager() : m_ownsData( true ) { m_iNEvents=0; m_iNTrueEvents=0; @@ -225,6 +225,28 @@ GPUManager::initData( const AmpVecs& a, bool use4Vectors ) calcCUDADims(); } +void +GPUManager::useDataFrom( const AmpVecs& a ){ + + clearData(); + + report( DEBUG, kModule ) << "Using shared GPU device arrays." << endl; + + m_ownsData = false; + + // copy over some info from the AmpVecs object for array dimensions + m_iNTrueEvents = a.m_iNTrueEvents; + m_iNEvents = a.m_iNEvents; + m_iNParticles = a.m_iNParticles; + + // the rest of the data are derived: + m_iGDoubleDataArrSize = sizeof(GDouble) * m_iNEvents; + + m_pfDevWeights = a.m_gpuMan.m_pfDevWeights; + m_pfDevData = a.m_gpuMan.m_pfDevData; + + calcCUDADims(); +} void GPUManager::initTerms( const AmpVecs& a, unsigned int chunkSize ) @@ -710,11 +732,11 @@ GPUManager::clearData(){ m_iGDoubleDataArrSize=0; if(m_pfDevData) - cudaFree(m_pfDevData); + if( m_ownsData )cudaFree(m_pfDevData); m_pfDevData=0; if(m_pfDevWeights) - cudaFree(m_pfDevWeights); + if( m_ownsData ) cudaFree(m_pfDevWeights); m_pfDevWeights=0; //CUDA Thread and Grid sizes diff --git a/AmpTools/GPUManager/GPUManager.h b/AmpTools/GPUManager/GPUManager.h index 36c29984..2a0d619a 100644 --- a/AmpTools/GPUManager/GPUManager.h +++ b/AmpTools/GPUManager/GPUManager.h @@ -76,6 +76,7 @@ class GPUManager void init( const AmpVecs& a, bool use4Vectors = true ); void initData( const AmpVecs& a, bool use4Vectors = true ); + void useDataFrom( const AmpVecs& a ); void initTerms( const AmpVecs& a, unsigned int chunkSize = 0 ); void copyDataToGPU( const AmpVecs& a, bool use4Vectors = true ); @@ -107,6 +108,8 @@ class GPUManager return( (1<m_gpuMan.initData( *targetAmpVecs, !needsUserVarsOnly ); - targetAmpVecs->m_gpuMan.copyDataToGPU( *targetAmpVecs, !needsUserVarsOnly ); + targetAmpVecs->m_gpuMan.useDataFrom( *this ); #endif } @@ -483,4 +496,9 @@ AmpVecs::claimDataOwnership( set< AmpVecs* > sharedFriends ){ (*avItr)->m_sharedDataHost = this; } + +#ifdef GPU_ACCELERATION + m_gpuMan.m_ownsData = true; +#endif + } diff --git a/AmpTools/IUAmpTools/AmpVecs.h b/AmpTools/IUAmpTools/AmpVecs.h index 3a11a5d4..c0eb10b6 100644 --- a/AmpTools/IUAmpTools/AmpVecs.h +++ b/AmpTools/IUAmpTools/AmpVecs.h @@ -234,6 +234,13 @@ struct AmpVecs void allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity = false, unsigned int chunkSize = 0 ); + + /** + * This routine deallocates the arrays that hold the calculated terms and + * (optionally) the calculated intensities. It should be called before + * the AmpVecs object is destroyed or before reallocation of the arrays. + */ + void deallocTerms(); #ifdef GPU_ACCELERATION /** diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc index 63eac3f1..3ba9aff4 100644 --- a/AmpTools/IUAmpTools/NormIntInterface.cc +++ b/AmpTools/IUAmpTools/NormIntInterface.cc @@ -413,7 +413,11 @@ NormIntInterface::forceCacheUpdate( bool normIntOnly ) const m_pIntenManager->calcIntegrals( m_genMCVecs, m_nGenEvents, chunkSize ); setAmpIntMatrix( m_genMCVecs.m_pdIntegralMatrix ); - + + // try to keep memory usage down by freeing the memory needed to calculate the + // integrals -- the data will remain in memory + m_genMCVecs.deallocTerms(); + m_emptyAmpIntCache = false; } From 2161ce5e5fd24176f9af91db48f8a0abb0b668da Mon Sep 17 00:00:00 2001 From: Matthew Shepherd Date: Wed, 8 Jul 2026 18:49:30 -0400 Subject: [PATCH 32/32] fix indexing bug in user variables array when doing chunked amplitude calculation --- AmpTools/GPUManager/GPUManager.cc | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc index 1b21ed55..68d5aa55 100644 --- a/AmpTools/GPUManager/GPUManager.cc +++ b/AmpTools/GPUManager/GPUManager.cc @@ -462,7 +462,7 @@ GPUManager::copyAmpsFromGPU( AmpVecs& a ) void GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, const vector< vector< int > >* pvPermutations, - unsigned int userVarsOffset, + unsigned int userVarsOffset, unsigned int startEvent, unsigned int chunkSize ) { #ifdef SCOREP @@ -518,7 +518,11 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, // increment the offset so that we place the computation for the // next permutation after the previous in pcResAmp permOffset += 2 * nEvents; - udLocalOffset += amp->numUserVars() * nEvents; + + // unlike the amplitude factor array the user variable array is not "chunked" + // when doing a chunked calculation, so we need to increment the offset by the + // the number of user variables times the total number of events, not the chunk size + udLocalOffset += amp->numUserVars() * m_iNEvents; } #ifdef SCOREP SCOREP_USER_REGION_END( calcAmplitudeAll_gpuMgr )