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 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__ 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 ); 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]; 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 ); } 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__ diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc index 2660dddd..68d5aa55 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"; @@ -62,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; @@ -70,8 +71,10 @@ GPUManager::GPUManager() m_iNAmps=0; m_iNUserVars=0; - m_iGDoubleEventArrSize=0; - m_iDoubleEventArrSize=0; + m_iGDoubleDataArrSize=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,27 +176,106 @@ GPUManager::~GPUManager() // Initialization routines: -void +void GPUManager::init( const AmpVecs& a, bool use4Vectors ) { - clearAll(); + initData( a, use4Vectors ); + initTerms( a ); +} +void +GPUManager::initData( const AmpVecs& a, bool use4Vectors ) +{ + 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; - m_iNAmps = a.m_iNTerms; - m_iNUserVars = a.m_userVarsPerEvent; + + // 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_iNParticles * sizeof( int ); + + totalMemory /= (1024*1024); + + report( INFO, kModule ) << "Attempting to allocate " << (int)totalMemory + << " 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; + 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 ) ) ; + + report( INFO, kModule ) << "GPU memory allocated for storage of " + << m_iNEvents << " events (" << m_iNTrueEvents + << " actual events)." << endl; + + //CUDA Dims + 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_iGDoubleEventArrSize = sizeof(GDouble) * m_iNEvents; + 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 ) +{ + 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 + unsigned int iNAmpEvents = ( chunkSize == 0 ? a.m_iNEvents : chunkSize ); + + report( DEBUG, kModule ) << "\tChunk size for amplitude calculations: " << iNAmpEvents << endl; + + m_iNUserVars = a.m_userVarsPerEvent; + m_iNAmps = a.m_iNTerms; - // double precision - m_iDoubleTrueEventArrSize = sizeof(double) * m_iNTrueEvents; - m_iDoubleEventArrSize = sizeof(double) * m_iNEvents; + // double precision intensity calculation + m_iDoubleIntenArrSize = sizeof(double) * a.m_iNEvents; + + // GDouble precision for factors and amplitudes + // ... and use the chunkSize for these (iNAmpEvents) when it is provided + m_iGDoubleFactArrSize = 2 * sizeof(GDouble) * a.m_maxFactPerEvent * iNAmpEvents; + // size needed to store amplitudes for each event - m_iAmpArrSize = 2 * sizeof(GDouble) * m_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; @@ -206,46 +288,56 @@ GPUManager::init( const AmpVecs& a, bool use4Vectors ) // 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 += 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 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 << " B ( " + << m_iVArrSize/sizeof(GDouble) << " elements )" << endl; + 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_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_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_pcDevAmpFact , 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 " + << iNAmpEvents << " events (" << m_iNTrueEvents << " total events)." << endl; - - //CUDA Dims - calcCUDADims(); + + calcCUDADimsAmpFact( iNAmpEvents ); } @@ -259,7 +351,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 +396,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; @@ -325,10 +417,15 @@ 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_iGDoubleEventArrSize, - cudaMemcpyHostToDevice ) ); + m_iNUserVars * m_iGDoubleDataArrSize, + cudaMemcpyHostToDevice ) ); #ifdef SCOREP SCOREP_USER_REGION_END( copyUserVarsToGPU ) #endif @@ -347,12 +444,15 @@ 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; + gpuErrChk( cudaMemcpy( a.m_pdAmps, m_pfDevAmps, m_iAmpArrSize, cudaMemcpyDeviceToHost ) ); - gpuErrChk( cudaMemcpy( a.m_pdAmpFactors, m_pcDevCalcAmp, - a.m_maxFactPerEvent * m_iGDoubleEventArrSize, + gpuErrChk( cudaMemcpy( a.m_pdAmpFactors, m_pcDevAmpFact, + m_iGDoubleFactArrSize, cudaMemcpyDeviceToHost ) ); #ifdef SCOREP SCOREP_USER_REGION_END( copyAmpsFromGPU ) @@ -360,9 +460,10 @@ GPUManager::copyAmpsFromGPU( AmpVecs& a ) } void -GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long long offset, +GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, const vector< vector< int > >* pvPermutations, - unsigned long long iUserVarsOffset ) + unsigned int userVarsOffset, + unsigned int startEvent, unsigned int chunkSize ) { #ifdef SCOREP SCOREP_USER_REGION_DEFINE( calcAmplitudeAll_gpuMgr ) @@ -370,19 +471,26 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long long offset, #endif dim3 dimBlock( m_iDimThreadX, m_iDimThreadY ); - dim3 dimGrid( m_iDimGridX, m_iDimGridY ); + dim3 dimGrid( m_iDimGridXAmpFact, m_iDimGridYAmpFact ); + + 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 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 ); - unsigned long long udLocalOffset = 0; - unsigned long long permOffset = 0; + 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 ), @@ -393,9 +501,9 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long 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 @@ -409,7 +517,11 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long long offset, // increment the offset so that we place the computation for the // next permutation after the previous in pcResAmp - permOffset += 2 * m_iNEvents; + permOffset += 2 * 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 @@ -418,24 +530,28 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned long 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; + + 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; @@ -449,7 +565,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 ) @@ -478,11 +594,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(); @@ -502,8 +623,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 ){ @@ -541,15 +662,16 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef, void GPUManager::calcIntegrals( double* result, int nElements, - const vector& iIndex, - const vector& jIndex ){ + const vector& iIndex, + 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]; @@ -574,38 +696,71 @@ 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: -void GPUManager::clearAll() -{ +void +GPUManager::clearAll() { + + clearData(); + clearTerms(); +} + +void +GPUManager::clearData(){ + m_iNParticles=0; m_iNEvents=0; m_iNTrueEvents=0; - m_iGDoubleEventArrSize=0; - m_iDoubleEventArrSize=0; - m_iDoubleTrueEventArrSize=0; - m_iNAmps=0; + m_iGDoubleDataArrSize=0; + + if(m_pfDevData) + if( m_ownsData )cudaFree(m_pfDevData); + m_pfDevData=0; + + if(m_pfDevWeights) + if( m_ownsData ) 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_iNUserVars=0; + + m_iDoubleIntenArrSize=0; + m_iGDoubleFactArrSize=0; m_iAmpArrSize=0; m_iVArrSize=0; + m_iNICalcSize=0; //Host Memory @@ -621,17 +776,13 @@ 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; + if(m_pcDevAmpFact) + cudaFree(m_pcDevAmpFact); + m_pcDevAmpFact=0; if(m_piDevPerm) cudaFree(m_piDevPerm); @@ -649,9 +800,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 +808,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: @@ -717,3 +859,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 82a8f2a1..2a0d619a 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; @@ -78,21 +70,24 @@ class GPUManager ~GPUManager(); void clearAll(); - void clearAmpCalc(); - void clearLikeCalc(); + void clearData(); + void clearTerms(); void init( const AmpVecs& a, bool use4Vectors = true ); - // Interface Utils - // First Amplitude calculation interface + 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 ); void copyUserVarsToGPU( const AmpVecs& a ); - void calcAmplitudeAll( const Amplitude* amp, unsigned long long offset, + void calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset, const vector< vector< int > >* pvPermutations, - unsigned long long userVarsOffset ); + unsigned int userVarsOffset, + unsigned int startEvent = 0, unsigned int chunkSize = 0 ); - void assembleTerms( int iAmpInd, int nFact, int nPerm ); + void assembleTerms( int iAmpInd, int nFact, int nPerm, unsigned int nEvents ); void copyAmpsFromGPU( AmpVecs& a ); @@ -100,8 +95,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 ){ @@ -112,25 +108,29 @@ class GPUManager return( (1<>> - ( nElements, pdDevNICalc, pfDevAmps, pfDevWeights, nEvents, nTrueEvents ); + ( nElements, pdDevNICalc, pfDevAmps, pfDevWeights, startEvent, nEvents, nTrueEvents ); } 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] << " "; } diff --git a/AmpTools/IUAmpTools/AmpVecs.cc b/AmpTools/IUAmpTools/AmpVecs.cc index 28423bf4..8252b9cb 100644 --- a/AmpTools/IUAmpTools/AmpVecs.cc +++ b/AmpTools/IUAmpTools/AmpVecs.cc @@ -88,12 +88,7 @@ AmpVecs::deallocAmpVecs() m_iNTrueEvents = 0 ; m_dSumWeights = 0 ; m_iNParticles = 0 ; - m_iNTerms = 0 ; - m_maxFactPerEvent = 0 ; - m_userVarsPerEvent = 0; - - m_termsValid = false ; - m_integralValid = false ; + m_dataLoaded = false ; m_hasNonUnityWeights = false; @@ -109,31 +104,9 @@ AmpVecs::deallocAmpVecs() delete[] m_pdWeights; m_pdWeights=0; - if(m_pdIntensity) - delete[] m_pdIntensity; - m_pdIntensity=0; - - if(m_pdIntegralMatrix) - delete[] m_pdIntegralMatrix; - m_pdIntegralMatrix=0; + deallocTerms(); - if(m_pdUserVars) - delete[] m_pdUserVars; - m_pdUserVars=0; - - m_userVarsOffset.clear(); - -#ifndef GPU_ACCELERATION - - if(m_pdAmps) - delete[] m_pdAmps; - m_pdAmps=0; - - if(m_pdAmpFactors) - delete[] m_pdAmpFactors; - m_pdAmpFactors=0; - -#else +#ifdef GPU_ACCELERATION //Deallocate "pinned memory" if(m_pdAmps) cudaFreeHost(m_pdAmps); @@ -198,26 +171,30 @@ 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, 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){ m_iNTrueEvents = iNTrueEvents; m_iNEvents = iNTrueEvents; + m_iNParticles = pKinematics->particleList().size(); + assert(m_iNParticles); + #ifdef GPU_ACCELERATION m_iNEvents = GPUManager::calcNEventsGPU(iNTrueEvents); + m_gpuMan.initData( *this, !needsUserVarsOnly ); #endif - - m_iNParticles = pKinematics->particleList().size(); - assert(m_iNParticles); - + + 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]; + m_pdWeights = new GDouble[m_iNEvents]; } // check to be sure we won't exceed the bounds of the array @@ -232,6 +209,7 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned long long iEvent, m_pdWeights[iEvent] = pKinematics->weight(); + m_termsValid = false; m_integralValid = false; m_dataLoaded = true; @@ -240,7 +218,7 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned long long iEvent, void -AmpVecs::loadData( DataReader* pDataReader ){ +AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly ){ // Make sure no data is already loaded @@ -272,9 +250,9 @@ 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 ); + loadEvent(pKinematics, iEvent, m_iNTrueEvents, needsUserVarsOnly ); float weight = pKinematics->weight(); @@ -292,13 +270,19 @@ 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++){ - loadEvent(pKinematics, iEvent, m_iNTrueEvents ); + for (unsigned int iEvent = m_iNTrueEvents; iEvent < m_iNEvents; iEvent++){ + loadEvent(pKinematics, iEvent, m_iNTrueEvents, needsUserVarsOnly ); } if( m_iNTrueEvents ) delete pKinematics; - + +#ifdef GPU_ACCELERATION + + m_gpuMan.copyDataToGPU( *this, !needsUserVarsOnly ); + +#endif + m_termsValid = false; m_integralValid = false; m_dataLoaded = true; @@ -307,8 +291,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 ampsEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize ); + m_iNTerms = intenMan.getTermNames().size(); m_maxFactPerEvent = intenMan.maxFactorStoragePerEvent(); m_userVarsPerEvent = intenMan.userVarsPerEvent(); @@ -346,13 +332,14 @@ 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]; + // 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.init( *this, !intenMan.needsUserVarsOnly() ); - m_gpuMan.copyDataToGPU( *this, !intenMan.needsUserVarsOnly() ); + m_gpuMan.initTerms( *this, chunkSize ); #endif // GPU_ACCELERATION @@ -360,6 +347,47 @@ AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity ) m_integralValid = false; } +void +AmpVecs::deallocTerms(){ + + m_iNTerms = 0; + m_maxFactPerEvent = 0; + m_userVarsPerEvent = 0; + + m_termsValid = false; + m_integralValid = false; + + if( m_pdIntegralMatrix ) + delete[] m_pdIntegralMatrix; + m_pdIntegralMatrix = 0; + + if( m_pdIntensity ) + delete[] m_pdIntensity; + m_pdIntensity = 0; + + if( m_pdUserVars ) + delete[] m_pdUserVars; + m_pdUserVars = 0; + + m_userVarsOffset.clear(); + +#ifndef GPU_ACCELERATION + +if( m_pdAmps ) + delete[] m_pdAmps; + m_pdAmps = 0; + if( m_pdAmpFactors ) + delete[] m_pdAmpFactors; + m_pdAmpFactors = 0; + +#else + + m_gpuMan.clearTerms(); + +#endif // GPU_ACCELERATION + +} + #ifdef GPU_ACCELERATION void AmpVecs::allocateCPUAmpStorage( const IntensityManager& intenMan ){ @@ -369,7 +397,7 @@ AmpVecs::allocateCPUAmpStorage( const IntensityManager& intenMan ){ // 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_pdAmpFactors, 2 * m_iNEvents * m_maxFactPerEvent * sizeof(GDouble)); cudaError_t cudaErr = cudaGetLastError(); if( cudaErr != cudaSuccess ){ @@ -378,7 +406,7 @@ AmpVecs::allocateCPUAmpStorage( const IntensityManager& intenMan ){ assert( false ); } } -#endif +#endif // GPU_ACCELERATION Kinematics* AmpVecs::getEvent( int iEvent ){ @@ -399,7 +427,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 @@ -407,7 +435,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; } @@ -426,7 +454,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]; @@ -440,6 +468,13 @@ 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.useDataFrom( *this ); +#endif } void @@ -461,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 c2f61a87..c0eb10b6 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,21 @@ 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 ); + + /** + * 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 /** @@ -240,6 +252,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 @@ -251,11 +264,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 @@ -271,8 +286,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, bool needsUserVarsOnly = false ); /** * A helper routine to get an event i from the array of data and weights. @@ -299,7 +314,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 diff --git a/AmpTools/IUAmpTools/Amplitude.cc b/AmpTools/IUAmpTools/Amplitude.cc index 053583e8..8e6fe46f 100644 --- a/AmpTools/IUAmpTools/Amplitude.cc +++ b/AmpTools/IUAmpTools/Amplitude.cc @@ -78,9 +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(); @@ -132,9 +130,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* pdAmpFact, unsigned int iNEvents, + const vector< vector< int > >* pvPermutations, + GDouble* pdUserVars, unsigned int startEvent, + unsigned int chunkSize ) const { #ifdef SCOREP @@ -154,19 +153,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 pdAmpFact 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; /** diff --git a/AmpTools/IUAmpTools/AmplitudeManager.cc b/AmpTools/IUAmpTools/AmplitudeManager.cc index 69e796d7..23fd7ba2 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 ){ @@ -273,7 +274,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 +301,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 +310,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 +338,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 +385,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 +418,35 @@ 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 + + return calcTerms( a, 0, 0 ); +} + +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 +461,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 +523,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 uAmpFactOffset = 0; for( iFactor=0; iFactor < iNFactors; - iFactor++, iLocalOffset += 2 * a.m_iNEvents * 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 long long uOffset = + unsigned int userVarsOffset = ( pCurrAmp->areUserVarsStatic() ? a.m_userVarsOffset[pCurrAmp->name()] : a.m_userVarsOffset[pCurrAmp->identifier()] ); @@ -524,17 +539,16 @@ 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, - &vvPermuations, - uOffset ); + a.m_gpuMan.calcAmplitudeAll( pCurrAmp, uAmpFactOffset, &vvPermuations, + userVarsOffset, startEvent, nEvents ); #endif//GPU_ACCELERATION } - - + // now assemble all the factors in an amplitude into a single // symmetrized amplitude for each event @@ -543,29 +557,26 @@ 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; - - // re-ordering of data will be useful to not fall out of (CPU) memory cache!!! - + unsigned int iOffsetA, iOffsetP, iOffsetF; + // 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; @@ -584,15 +595,18 @@ 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 - a.m_gpuMan.assembleTerms( iAmpIndex, iNFactors, iNPermutations ); + a.m_gpuMan.assembleTerms( iAmpIndex, iNFactors, iNPermutations, nEvents ); #endif } - + a.m_termsValid = true; // finally make a loop and record the iteration of parameters @@ -796,7 +810,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 +823,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 : a.m_iNEvents / chunkSize ); + if( ( chunkSize !=0 ) && ( a.m_iNEvents % 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 +840,151 @@ 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 = 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 = a.m_iNEvents % 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+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; } - } - -#ifndef GPU_ACCELERATION - for( int iTerm = 0; iTerm < nCompute; ++iTerm ){ - - int i = iIndex[iTerm]; - int j = jIndex[iTerm]; - for( int iEvent = 0; iEvent < a.m_iNTrueEvents; iEvent++ ) { +#else + + vector temp(maxNIElements*2); - //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] ); + // use the GPU manager to compute the result for the chunk + a.m_gpuMan.calcIntegrals( &(temp[0]), nCompute, iIndex, jIndex, + startEvent, nEvents ); - 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] ); + // 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 + if( nChunk > 1 ) report( DEBUG, kModule ) << "Integral calculation: chunk " + << iChunk << " of " << nChunk << " complete." << endl; } - -#else - // use the GPU manager to compute the result - a.m_gpuMan.calcIntegrals( &(result[0]), nCompute, iIndex, jIndex ); + // 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 -#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 +995,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 diff --git a/AmpTools/IUAmpTools/LikelihoodCalculator.cc b/AmpTools/IUAmpTools/LikelihoodCalculator.cc index c40420b1..af659398 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,18 +247,17 @@ 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 ){ - 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 ); } diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc index 7f02bf4c..3ba9aff4 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 ); - - 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 ); - - 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 ); @@ -428,23 +397,33 @@ 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, false, 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 ); - + + // 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; } // 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 @@ -552,8 +531,68 @@ NormIntInterface::setNormIntMatrix( const double* input ) const { m_emptyNormIntCache = false; } - #ifndef __ACLIC__ + +unsigned int +NormIntInterface::genMCChunkSize() const { + + // need to know the sizes of the data sets + assert( m_accMCVecs.m_dataLoaded && m_genMCVecs.m_dataLoaded ); + + 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 + // 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/4 ) ++iPow; + + report( DEBUG, kModule ) << "Chunk size for generated MC: " << ( nGen >> iPow ) << endl; + + return( nGen >> iPow ); +} + +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 976420d7..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,12 +116,16 @@ class NormIntInterface void setAmpIntMatrix( const double* input ) const; void setNormIntMatrix( const double* input ) const; + + void loadMC() const; private: void initializeCache(); int m_cacheSize; + unsigned int genMCChunkSize() const; + vector< string > m_termNames; map< string, int > m_termIndex; @@ -131,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 ); } }