diff --git a/AmpTools/GPUManager/GPUAmpProductKernel.cu b/AmpTools/GPUManager/GPUAmpProductKernel.cu
index 2a0716a..c719e7c 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, unsigned int nEvents, double* pdDevRes )
+            int nAmps, size_t 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, unsigned int nEvents, double* pdDevRes )
+     int nAmps, size_t 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 d78f50e..b6b291c 100755
--- a/AmpTools/GPUManager/GPUCustomTypes.h
+++ b/AmpTools/GPUManager/GPUCustomTypes.h
@@ -52,8 +52,8 @@ using namespace std;
 
 // standard arguments to ampliutde kernel launches
 #define GPU_AMP_PROTO GDouble* pfDevData, GDouble* pfDevUserVars, \
-                      WCUComplex* pcDevAmp, int* piDevPerm, unsigned int iNParticles, \
-                      unsigned int iNEvents, unsigned int startEvent
+                      WCUComplex* pcDevAmp, int* piDevPerm, size_t iNParticles, \
+                      size_t iNEvents, size_t startEvent
 #define GPU_AMP_ARGS  pfDevData, pfDevUserVars, pcDevAmp, piDevPerm, iNParticles, \
                       iNEvents, startEvent                   
 
diff --git a/AmpTools/GPUManager/GPUFactPermKernel.cu b/AmpTools/GPUManager/GPUFactPermKernel.cu
index faefbc0..3ec2511 100644
--- a/AmpTools/GPUManager/GPUFactPermKernel.cu
+++ b/AmpTools/GPUManager/GPUFactPermKernel.cu
@@ -38,7 +38,7 @@
 
 __global__ void
 fact_perm_kernel( GDouble* pfDevAmps, GDouble* pcDevAmpFact, int nFact,
-                  int nPerm, unsigned int nEvents )
+                  int nPerm, size_t nEvents )
 {
 
 	int i = threadIdx.x + GPU_BLOCK_SIZE_X * threadIdx.y + 
@@ -46,14 +46,14 @@ fact_perm_kernel( GDouble* pfDevAmps, GDouble* pcDevAmpFact, int nFact,
 
   for( int iPerm = 0; iPerm < nPerm; ++iPerm ){
   
-    unsigned int offsetP = 2*nEvents*iPerm + 2*i;
+    size_t offsetP = 2*nEvents*iPerm + 2*i;
     
     GDouble ampRe = pcDevAmpFact[offsetP];
     GDouble ampIm = pcDevAmpFact[offsetP+1];
 
     for( int iFactor = 1; iFactor < nFact; ++iFactor ){
     
-       unsigned int offsetF = 2*nEvents*nPerm*iFactor + offsetP;
+       size_t offsetF = 2*nEvents*nPerm*iFactor + offsetP;
 
        GDouble re = ampRe;
        GDouble im = ampIm;
@@ -71,7 +71,7 @@ fact_perm_kernel( GDouble* pfDevAmps, GDouble* pcDevAmpFact, int nFact,
 }
 
 extern "C" void GPU_ExecFactPermKernel( dim3 dimGrid, dim3 dimBlock,
-    GDouble* pfDevAmps, GDouble* pcDevAmpFact, int nFact, int nPerm, unsigned int nEvents )
+    GDouble* pfDevAmps, GDouble* pcDevAmpFact, int nFact, int nPerm, size_t nEvents )
 {
   fact_perm_kernel<<< dimGrid, dimBlock >>>( pfDevAmps, pcDevAmpFact, 
                                              nFact, nPerm, nEvents );
diff --git a/AmpTools/GPUManager/GPUKernel.h b/AmpTools/GPUManager/GPUKernel.h
index 44116c1..5217d03 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, unsigned int nEvents, double* pfDevRes);
+                                  int nAmps, size_t nEvents, double* pfDevRes);
 
 extern "C" void GPU_ExecFactPermKernel( dim3 dimGrid, dim3 dimBlock,
                                         GDouble* pfDevAmps, GDouble* pcDevAmpFact,
-                                        int nFact, int nPerm, unsigned int nEvents );
+                                        int nFact, int nPerm, size_t nEvents );
 
-extern "C" void GPU_ExecNICalcKernel( dim3 dimGrid, dim3 dimBlock, unsigned int sharedSize,
+extern "C" void GPU_ExecNICalcKernel( dim3 dimGrid, dim3 dimBlock, size_t sharedSize,
                                       int nElements, double* pdDevNICalc,
                                       GDouble* pfDevAmps, GDouble* pfDevWeights,
-                                      unsigned int startEvent, unsigned int nEvents, 
-                                      unsigned int nTrueEvents, unsigned int maxSize );
+                                      size_t startEvent, size_t nEvents, 
+                                      size_t nTrueEvents, size_t maxSize );
 
 #endif //__GPU_KERNEL__H__
diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc
index 68d5aa5..2e858be 100644
--- a/AmpTools/GPUManager/GPUManager.cc
+++ b/AmpTools/GPUManager/GPUManager.cc
@@ -103,9 +103,6 @@ GPUManager::GPUManager() : m_ownsData( true )
   
   int thisDevice = 0;
   
-  if( !m_cudaDisplay )
-    report( INFO, kModule ) << "################### CUDA DEVICE ##################" << endl;
-  
 #ifdef USE_MPI
   
   // Note that a better algorithm would be to utilize the jobs "local
@@ -119,15 +116,11 @@ GPUManager::GPUManager() : m_ownsData( true )
   MPI_Comm_rank( MPI_COMM_WORLD, &rank );
   cudaGetDeviceCount(&devs);
   thisDevice = rank % devs;
+
+  #endif
   
-  if( !m_cudaDisplay ) {
-    report( INFO, kModule ) << "Parallel GPU configuration requested." << endl;
-    report( INFO, kModule ) << "Number of CUDA devices available on this node:  " << devs << endl;
-    report( INFO, kModule ) << "MPI process " << rank << " is using device " << thisDevice << endl;
-  }
-#endif
-  
-  ///////CUDA INITIALIZATION
+  // CUDA INITIALIZATION
+
   gpuErrChk( cudaSetDevice( thisDevice ) );
   cudaDeviceProp devProp;
   gpuErrChk( cudaGetDeviceProperties( &devProp, thisDevice ) );
@@ -136,13 +129,18 @@ GPUManager::GPUManager() : m_ownsData( true )
   m_devProp_minor = devProp.minor;
 
   if( ! m_cudaDisplay ){
-    
+    report( INFO, kModule ) << "################### CUDA DEVICE ##################" << endl;
+#ifdef USE_MPI
+    report( INFO, kModule ) << "Parallel GPU configuration requested." << endl;
+    report( INFO, kModule ) << "Number of CUDA devices available on this node:  " << devs << endl;
+    report( INFO, kModule ) << "MPI process " << rank << " is using device " << thisDevice << endl;
+#endif
     report( INFO, kModule ) << "Current GPU Properites:\n";
     report( INFO, kModule ) << "\t Name: "<<devProp.name<<endl;
     report( INFO, kModule ) << "\t Total global memory: "<<devProp.totalGlobalMem/((float)1024*1024)<<" MB"<<endl;
     report( INFO, kModule ) << "\t Rev.: "<<devProp.major<<"."<<devProp.minor<<endl;
     report( INFO, kModule ) << "##################################################" << endl;
-    ///////END OF CUDA INITIALIZATION
+
     m_cudaDisplay = true;
   }
   
@@ -162,18 +160,11 @@ GPUManager::GPUManager() : m_ownsData( true )
      ( m_devProp_major == 7 && m_devProp_minor == 0 ) ) m_maxShared_bytes = 98304;
 }
 
-GPUManager::GPUManager( const AmpVecs& a )
-{
-  GPUManager();
-  init( a );
-}
-
 GPUManager::~GPUManager()
 {
   clearAll();
 }
 
-
 // Initialization routines:
 
 void
@@ -249,7 +240,7 @@ GPUManager::useDataFrom( const AmpVecs& a ){
 }
 
 void
-GPUManager::initTerms( const AmpVecs& a, unsigned int chunkSize )
+GPUManager::initTerms( const AmpVecs& a, size_t chunkSize )
 {
   clearTerms();
 
@@ -259,20 +250,23 @@ GPUManager::initTerms( const AmpVecs& a, unsigned int chunkSize )
   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 );
+  size_t 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 intensity calculation
-  m_iDoubleIntenArrSize = sizeof(double) * a.m_iNEvents;
+  // double precision intensity calculation on the device
+  m_iDoubleIntenArrSize = sizeof(double) * m_iNEvents;
+
+  // array size for reduction of intensity sum on the device
+  m_iDoubleReduceArrSize = sizeof(double) * ( m_iNTrueEvents <= m_iNBlocks ? m_iNTrueEvents : m_iNBlocks );
 
   // 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;
+  // m_maxFactPerEvent is in units of doubles, and includes a factor of 2 for complex numbers
+  m_iGDoubleFactArrSize = sizeof(GDouble) * a.m_maxFactPerEvent * iNAmpEvents;
 
   // size needed to store amplitudes for each event
   m_iAmpArrSize = 2 * sizeof(GDouble) * iNAmpEvents * m_iNAmps;
@@ -288,14 +282,15 @@ GPUManager::initTerms( const AmpVecs& a, unsigned int chunkSize )
   
   // host memory needed for intensity or integral calculation
   cudaMallocHost( (void**)&m_pfVVStar , m_iVArrSize     );
-  cudaMallocHost( (void**)&m_pdRes    , m_iDoubleIntenArrSize );
+  cudaMallocHost( (void**)&m_pdRes    , m_iDoubleReduceArrSize );
   
   double totalMemory = 0;
   
   totalMemory += m_iNUserVars * m_iGDoubleDataArrSize;
   totalMemory += m_iVArrSize;
   totalMemory += m_iNICalcSize;
-  totalMemory += 2*m_iDoubleIntenArrSize;
+  totalMemory += m_iDoubleIntenArrSize;
+  totalMemory += m_iDoubleReduceArrSize;
   totalMemory += m_iGDoubleFactArrSize;
   totalMemory += m_iAmpArrSize;
   totalMemory += m_iNParticles * sizeof( int );
@@ -327,7 +322,7 @@ GPUManager::initTerms( const AmpVecs& a, unsigned int chunkSize )
   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_pdDevREDUCE  , m_iDoubleReduceArrSize ) ) ;
   gpuErrChk( cudaMalloc( (void**)&m_pcDevAmpFact , m_iGDoubleFactArrSize ) ) ;
   gpuErrChk( cudaMalloc( (void**)&m_pfDevAmps    , m_iAmpArrSize         ) ) ;
 
@@ -460,10 +455,10 @@ GPUManager::copyAmpsFromGPU( AmpVecs& a )
 }
 
 void 
-GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset,
+GPUManager::calcAmplitudeAll( const Amplitude* amp, size_t uAmpFactOffset,
                               const vector< vector< int > >* pvPermutations,
-                              unsigned int userVarsOffset, 
-                              unsigned int startEvent, unsigned int chunkSize )
+                              size_t userVarsOffset, 
+                              size_t startEvent, size_t chunkSize )
 {
   #ifdef SCOREP
   SCOREP_USER_REGION_DEFINE( calcAmplitudeAll_gpuMgr )                                                                                    
@@ -473,9 +468,13 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset,
   dim3 dimBlock( m_iDimThreadX, m_iDimThreadY );
   dim3 dimGrid( m_iDimGridXAmpFact, m_iDimGridYAmpFact );
 
-  unsigned int nEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize );
+  size_t nEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize );
 
   report( DEBUG, kModule ) << "Calculating amplitude factors for amplitude " << amp->identifier() << endl;
+  report( DEBUG, kModule ) << "\tChunk size for amplitude calculations:  " << nEvents << endl;
+  report( DEBUG, kModule ) << "\tUser variable offset: " << userVarsOffset << endl;
+  report( DEBUG, kModule ) << "\tAmplitude factor offset: " << uAmpFactOffset << endl;
+  report( DEBUG, kModule ) << "\tStart event: " << startEvent << endl;
 
   // do the computation for all events for each permutation in the
   // vector of permunations
@@ -487,8 +486,8 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset,
   // if this is not true, AmplitudeManager hasn't been setup properly
   assert( permItr->size() == m_iNParticles );
   
-  unsigned int udLocalOffset = 0;
-  unsigned int permOffset = 0;
+  size_t udLocalOffset = 0;
+  size_t permOffset = 0;
   for( ; permItr != pvPermutations->end(); ++permItr ){
 
     // copy the permutation to global memory
@@ -498,13 +497,18 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset,
     
     // calculate amplitude factor for all events --
     // casting amp array to WCUComplex for 8 or 16 bit write 
-    // operation of both real and complex parts at once
+    // operation of both real and complex parts at once;
+    // note that here we need to pass the total number of events
+    // not the chunk size so that the data and user varaiable
+    // arrays can be properly indexed in the kernel
 
     amp->calcAmplitudeGPU( dimGrid, dimBlock, m_pfDevData,
                            &m_pfDevUserVars[userVarsOffset+udLocalOffset],
                           (WCUComplex*)&m_pcDevAmpFact[uAmpFactOffset+permOffset],
                            m_piDevPerm, m_iNParticles, m_iNEvents, startEvent,
                            *permItr );
+
+    cudaDeviceSynchronize();
 			   
     // check to be sure kernel execution was OK
     cudaError_t cerrKernel=cudaGetLastError();
@@ -530,7 +534,7 @@ GPUManager::calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset,
 }
 
 void
-GPUManager::assembleTerms( int iAmpInd, int nFact, int nPerm, unsigned int nEvents ){
+GPUManager::assembleTerms( int iAmpInd, int nFact, int nPerm, size_t nEvents ){
   #ifdef SCOREP
   SCOREP_USER_REGION_DEFINE( assembleTerms )                                                                                    
   SCOREP_USER_REGION_BEGIN( assembleTerms, "assembleTerms", SCOREP_USER_REGION_TYPE_COMMON )                           
@@ -572,7 +576,7 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef,
   SCOREP_USER_REGION_BEGIN( calcSumLogIntensity_gpuMgr, "calcSumLogIntensity_gpuMgr", SCOREP_USER_REGION_TYPE_COMMON )                           
   #endif
 
-  unsigned int i,j;
+  size_t i,j;
   
   // precompute the real and imaginary parts of ViVj* and copy to 
   // GPU global memory
@@ -623,7 +627,7 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef,
   if( m_iNTrueEvents <= m_iNBlocks )
   {
     gpuErrChk( cudaMemcpy( m_pdRes, m_pdDevRes,
-                           m_iDoubleIntenArrSize,cudaMemcpyDeviceToHost ) );
+                           m_iDoubleReduceArrSize, cudaMemcpyDeviceToHost ) );
     for( i=0; i < m_iNTrueEvents; i++ )
       dGPUResult += m_pdRes[i];
   }
@@ -637,6 +641,8 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef,
     // execute the kernel to sum partial sums from each block on CPU
     reduce<double>( m_iNEvents, m_iNThreads, m_iNBlocks, m_pdDevRes, m_pdDevREDUCE );
 
+    cudaDeviceSynchronize();
+
     cerrKernel = cudaGetLastError();
     if( cerrKernel!= cudaSuccess  ){
       
@@ -646,11 +652,10 @@ GPUManager::calcSumLogIntensity( const vector< complex< double > >& prodCoef,
     }
   
     // Copy result from device to host
-    gpuErrChk( cudaMemcpy( m_pdRes, m_pdDevREDUCE, m_iNBlocks*sizeof(double),
+    gpuErrChk( cudaMemcpy( m_pdRes, m_pdDevREDUCE, m_iDoubleReduceArrSize,
                            cudaMemcpyDeviceToHost) );
     for(i=0; i<m_iNBlocks; i++)
       dGPUResult += m_pdRes[i];
-
   }
     
   #ifdef SCOREP
@@ -664,11 +669,11 @@ void
 GPUManager::calcIntegrals( double* result, int nElements,
                            const vector<int>& iIndex,
                            const vector<int>& jIndex,
-                           unsigned int startEvent, unsigned int nEvents ){
+                           size_t startEvent, size_t nEvents ){
 
-  unsigned int resultSize = 2*sizeof(double)*nElements;
-  unsigned int indexSize = sizeof(int)*nElements;
-  unsigned int totalSize = resultSize + 2*indexSize;
+  size_t resultSize = 2*sizeof(double)*nElements;
+  size_t indexSize = sizeof(int)*nElements;
+  size_t totalSize = resultSize + 2*indexSize;
 
   dim3 dimBlock( m_iDimThreadX, m_iDimThreadY );
   dim3 dimGrid( m_iDimGridXAmpFact, m_iDimGridYAmpFact );
@@ -695,25 +700,23 @@ GPUManager::calcIntegrals( double* result, int nElements,
     << "Unable to continue -- reduce the number of amplitudes to perform this fit on the GPU." << endl;
     exit( 1 );
   }
-  
-  // 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 );
-    }
-  }
+  GPU_ExecNICalcKernel( dimGrid, dimBlock, totalSize, nElements,
+                        m_pdDevNICalc, m_pfDevAmps, m_pfDevWeights,
+                        startEvent, nEvents, m_iNTrueEvents,
+                        m_devProp_major >= 7 ? m_maxShared_bytes : 0 );
+
+  cudaDeviceSynchronize();
 
+  // 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 ) );
 }
 
@@ -819,8 +822,8 @@ void GPUManager::calcCUDADims()
   m_iDimThreadX=GPU_BLOCK_SIZE_X;
   m_iDimThreadY=GPU_BLOCK_SIZE_Y; 
   
-  unsigned int iBlockSizeSq=GPU_BLOCK_SIZE_SQ;
-  unsigned int iNBlocks=m_iNEvents/iBlockSizeSq;
+  size_t iBlockSizeSq=GPU_BLOCK_SIZE_SQ;
+  size_t iNBlocks=m_iNEvents/iBlockSizeSq;
   if(iNBlocks<=1)
   {
     m_iDimGridX=1;
@@ -828,7 +831,7 @@ void GPUManager::calcCUDADims()
   }
   else
   {
-    unsigned int iDivLo=1,iDivHi=iNBlocks;
+    size_t iDivLo=1,iDivHi=iNBlocks;
     for(iDivLo=static_cast<int>(sqrt(iNBlocks));iDivLo>=1;iDivLo--)
     {
       iDivHi=iNBlocks/iDivLo;
@@ -843,8 +846,8 @@ void GPUManager::calcCUDADims()
   report( DEBUG, kModule ) << "\tGrid dimensions:  ("<<m_iDimGridX<<","<<m_iDimGridY<<")\n";
   
   //Reduction Parameters
-  unsigned int maxThreads = ( m_devProp_major >= 2 ? 1024 : 512 );  // number of threads per block
-  unsigned int maxBlocks = 1024;  
+  size_t maxThreads = ( m_devProp_major >= 2 ? 1024 : 512 );  // number of threads per block
+  size_t maxBlocks = 1024;  
   
   if (m_iNEvents == 1) 
     m_iNThreads = 1;
@@ -859,10 +862,10 @@ void GPUManager::calcCUDADims()
   report( DEBUG, kModule ) << "\tNumber of blocks:   "<<m_iNBlocks<<"\n\n\n"<<flush;
 }
 
-void GPUManager::calcCUDADimsAmpFact( unsigned int chunkSize )
+void GPUManager::calcCUDADimsAmpFact( size_t chunkSize )
 {
-  unsigned int iBlockSizeSq = GPU_BLOCK_SIZE_SQ;
-  unsigned int iNBlocksAmpFact = chunkSize / iBlockSizeSq;
+  size_t iBlockSizeSq = GPU_BLOCK_SIZE_SQ;
+  size_t iNBlocksAmpFact = chunkSize / iBlockSizeSq;
   if( iNBlocksAmpFact <= 1 )
   {
     m_iDimGridXAmpFact = 1;
@@ -870,7 +873,7 @@ void GPUManager::calcCUDADimsAmpFact( unsigned int chunkSize )
   }
   else
   {
-    unsigned int iDivLo=1,iDivHi=iNBlocksAmpFact;
+    size_t iDivLo=1,iDivHi=iNBlocksAmpFact;
     for(iDivLo=static_cast<int>(sqrt(iNBlocksAmpFact));iDivLo>=1;iDivLo--)
     {
       iDivHi=iNBlocksAmpFact/iDivLo;
diff --git a/AmpTools/GPUManager/GPUManager.h b/AmpTools/GPUManager/GPUManager.h
index 2a0d619..fa7505a 100644
--- a/AmpTools/GPUManager/GPUManager.h
+++ b/AmpTools/GPUManager/GPUManager.h
@@ -66,7 +66,6 @@ class GPUManager
 public:
   
   GPUManager();
-  GPUManager( const AmpVecs& a );
   ~GPUManager();
   
   void clearAll();
@@ -77,17 +76,17 @@ class GPUManager
   
   void initData( const AmpVecs& a, bool use4Vectors = true );
   void useDataFrom( const AmpVecs& a );
-  void initTerms( const AmpVecs& a, unsigned int chunkSize = 0 );
+  void initTerms( const AmpVecs& a, size_t chunkSize = 0 );
   
   void copyDataToGPU( const AmpVecs& a, bool use4Vectors = true  );
   void copyUserVarsToGPU( const AmpVecs& a );
   
-  void calcAmplitudeAll( const Amplitude* amp, unsigned int uAmpFactOffset,
+  void calcAmplitudeAll( const Amplitude* amp, size_t uAmpFactOffset,
                          const vector< vector< int > >* pvPermutations,
-                         unsigned int userVarsOffset,
-                         unsigned int startEvent = 0, unsigned int chunkSize = 0 );
+                         size_t userVarsOffset,
+                         size_t startEvent = 0, size_t chunkSize = 0 );
   
-  void assembleTerms( int iAmpInd, int nFact, int nPerm, unsigned int nEvents );
+  void assembleTerms( int iAmpInd, int nFact, int nPerm, size_t nEvents );
   
   void copyAmpsFromGPU( AmpVecs& a );
 
@@ -97,15 +96,15 @@ class GPUManager
 
   void calcIntegrals( double* result, int nElements, 
                       const vector<int>& iIndex, const vector<int>& jIndex, 
-                      unsigned int startEvent, unsigned int nEvents );
+                      size_t startEvent, size_t nEvents );
   
   // General utils:
-  static int calcNEventsGPU( int iNEvents ){
+  static size_t calcNEventsGPU( size_t iNEvents ){
     
     //Should be a power of 2 for reduction to work, also multiple of GPU_BLOCK_SIZE_SQ    
     int iPow = 0;
-    while( ( 1 << iPow ) < iNEvents ) iPow++;
-    return(  (1<<iPow) < GPU_BLOCK_SIZE_SQ ? GPU_BLOCK_SIZE_SQ : (1<<iPow) );
+    while( ( 1L << iPow ) < iNEvents ) iPow++;
+    return(  ( 1L << iPow ) < GPU_BLOCK_SIZE_SQ ? GPU_BLOCK_SIZE_SQ : ( 1L << iPow ) );
   }
   
   bool m_ownsData;
@@ -115,21 +114,22 @@ class GPUManager
   static bool m_cudaDisplay;
   
   // array dimensions
-  unsigned int m_iNParticles;
-  unsigned int m_iNEvents;
-  unsigned int m_iNTrueEvents;
-  unsigned int m_iNAmps;
-  unsigned int m_iNUserVars;
+  size_t m_iNParticles;
+  size_t m_iNEvents;
+  size_t m_iNTrueEvents;
+  size_t m_iNAmps;
+  size_t m_iNUserVars;
   
   // array sizes (data)
-  unsigned int m_iGDoubleDataArrSize;
+  size_t 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;
+  size_t m_iDoubleIntenArrSize;
+  size_t m_iDoubleReduceArrSize;
+  size_t m_iGDoubleFactArrSize;
+  size_t m_iAmpArrSize;
+  size_t m_iVArrSize;
+  size_t m_iNICalcSize;
 
   //Host Arrays
   GDouble* m_pfVVStar;
@@ -150,24 +150,24 @@ class GPUManager
   double* m_pdDevREDUCE;
   
   // 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;
+  size_t m_iDimGridX;
+  size_t m_iDimGridY;
+  size_t m_iDimGridXAmpFact;
+  size_t m_iDimGridYAmpFact;
+  size_t m_iDimThreadX;
+  size_t m_iDimThreadY;
   
-  unsigned int m_iNBlocks;
-  unsigned int m_iNThreads;
+  size_t m_iNBlocks;
+  size_t m_iNThreads;
   
   // Internal Utils
   
   unsigned int m_devProp_major;
   unsigned int m_devProp_minor;
-  unsigned int m_maxShared_bytes;
+  size_t m_maxShared_bytes;
 
   void calcCUDADims();
-  void calcCUDADimsAmpFact( unsigned int chunkSize = 0 );
+  void calcCUDADimsAmpFact( size_t chunkSize = 0 );
 
   static const char* kModule;
 };
diff --git a/AmpTools/GPUManager/GPUNICalcKernel.cu b/AmpTools/GPUManager/GPUNICalcKernel.cu
index 244ab43..fc083c1 100644
--- a/AmpTools/GPUManager/GPUNICalcKernel.cu
+++ b/AmpTools/GPUManager/GPUNICalcKernel.cu
@@ -40,8 +40,8 @@
 __global__ void
 ni_calc_kernel( int nElements, double* pdDevNICalc,
                 GDouble* pfDevAmps, GDouble* pfDevWeights,
-                unsigned int startEvent, unsigned int nEvents, 
-                unsigned int nTrueEvents )
+                size_t startEvent, size_t nEvents, 
+                size_t nTrueEvents )
 {
 
   // used shared memory block for amplitude indices and results
@@ -75,11 +75,11 @@ ni_calc_kernel( int nElements, double* pdDevNICalc,
   __syncthreads();
   
   // this is the event index in the chunk
-  unsigned int iEvtChunk = iThread +
+  size_t 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;  
+  size_t iEvt = iEvtChunk + startEvent;  
 
   if( iEvt < nTrueEvents ){ // do not compute for the padding events
 
@@ -128,11 +128,11 @@ ni_calc_kernel( int nElements, double* pdDevNICalc,
 
 
 extern "C" void GPU_ExecNICalcKernel( dim3 dimGrid, dim3 dimBlock,
-       	   			                      unsigned int sharedSize,
+       	   			                      size_t sharedSize,
                                       int nElements, double* pdDevNICalc,
                                       GDouble* pfDevAmps, GDouble* pfDevWeights,
-                                      unsigned int startEvent, unsigned int nEvents, 
-                                      unsigned int nTrueEvents, unsigned int maxSize )
+                                      size_t startEvent, size_t nEvents, 
+                                      size_t nTrueEvents, size_t maxSize )
 {
 
    if( maxSize )
diff --git a/AmpTools/GPUManager/GPUReductionKernel.cu b/AmpTools/GPUManager/GPUReductionKernel.cu
index 3ee8b02..8e4f39f 100644
--- a/AmpTools/GPUManager/GPUReductionKernel.cu
+++ b/AmpTools/GPUManager/GPUReductionKernel.cu
@@ -142,17 +142,17 @@ struct SharedMemory<double>
     In other words if blockSize <= 32, allocate 64*sizeof(T) bytes.  
     If blockSize > 32, allocate blockSize*sizeof(T) bytes.
 */
-template <class T, unsigned int blockSize, bool nIsPow2>
+template <class T, size_t blockSize, bool nIsPow2>
 __global__ void
-reduce6(T *g_idata, T *g_odata, unsigned int n)
+reduce6(T *g_idata, T *g_odata, size_t n)
 {
     T *sdata = SharedMemory<T>();
 
     // perform first level of reduction,
     // reading from global memory, writing to shared memory
-    unsigned int tid = threadIdx.x;
-    unsigned int i = blockIdx.x*blockSize*2 + threadIdx.x;
-    unsigned int gridSize = blockSize*2*gridDim.x;
+    size_t tid = threadIdx.x;
+    size_t i = blockIdx.x*blockSize*2 + threadIdx.x;
+    size_t gridSize = blockSize*2*gridDim.x;
     
     T mySum = 0;
 
diff --git a/AmpTools/IUAmpTools/AmpVecs.cc b/AmpTools/IUAmpTools/AmpVecs.cc
index 8252b9c..c5a797a 100644
--- a/AmpTools/IUAmpTools/AmpVecs.cc
+++ b/AmpTools/IUAmpTools/AmpVecs.cc
@@ -171,8 +171,8 @@ AmpVecs::clearFourVecs(){
 }
 
 void
-AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent,
-                    unsigned int iNTrueEvents, bool needsUserVarsOnly ){
+AmpVecs::loadEvent( const Kinematics* pKinematics, size_t iEvent,
+                    size_t iNTrueEvents, bool needsUserVarsOnly ){
   
   // allocate memory and set variables
   // if this is the first call to this method
@@ -196,19 +196,8 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent,
     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
-  assert( iEvent < m_iNEvents );
     
-  for (int iParticle = 0; iParticle < m_iNParticles; iParticle++){
-    m_pdData[4*iEvent*m_iNParticles+4*iParticle+0]=pKinematics->particle(iParticle).E();
-    m_pdData[4*iEvent*m_iNParticles+4*iParticle+1]=pKinematics->particle(iParticle).Px();
-    m_pdData[4*iEvent*m_iNParticles+4*iParticle+2]=pKinematics->particle(iParticle).Py();
-    m_pdData[4*iEvent*m_iNParticles+4*iParticle+3]=pKinematics->particle(iParticle).Pz();
-  }
-
-  m_pdWeights[iEvent] = pKinematics->weight();
-  
+  loadDataArrayElement( pKinematics, iEvent );
 
   m_termsValid = false;
   m_integralValid = false;
@@ -216,13 +205,13 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, unsigned int iEvent,
   m_userVarsOffset.clear();
 }
 
-
 void
-AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly ){
+AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly, size_t chunkSize ){
   
   //  Make sure no data is already loaded
   
-  if( m_pdData!=0 || m_pdWeights!=0 ){
+  if( m_pdData != NULL || m_pdWeights != NULL ){
+
     report( ERROR, kModule ) << "Trying to load data into a non-empty AmpVecs object\n"<<flush;
     assert(false);
   }
@@ -231,6 +220,20 @@ AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly ){
   
   pDataReader->resetSource();
   m_iNTrueEvents = pDataReader->numEvents();
+  m_iNEvents = m_iNTrueEvents;
+
+  #ifdef GPU_ACCELERATION
+
+  // this pads the number of events up to the next power of 2 for GPU calculations
+  // if a chunk size is specified then the chunk size is already a power of 2
+  // and we need only a multiple of the chunk size that exceeds the number of true events
+  m_iNEvents = GPUManager::calcNEventsGPU(m_iNTrueEvents);
+  if( chunkSize != 0 ){
+    while( m_iNEvents > m_iNTrueEvents ) m_iNEvents -= chunkSize;
+    m_iNEvents += chunkSize;
+  }
+
+  #endif
 
   // try to print an informative message -- this can be normal behavior in
   // an MPI job with a sparse background source and many concurrent processess
@@ -247,35 +250,53 @@ AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly ){
     report( NOTICE, kModule ) << "\n does not contain any events." << endl;
   }
   
-  // Loop over events and load each one individually
-  
+
+  // loop over events and load each one individually
   Kinematics* pKinematics;
-  for(unsigned int iEvent = 0; iEvent < m_iNTrueEvents; iEvent++){
+  for(size_t iEvent = 0; iEvent < m_iNTrueEvents; iEvent++){
+
     pKinematics = pDataReader->getEvent();
-    loadEvent(pKinematics, iEvent, m_iNTrueEvents, needsUserVarsOnly );
 
-    float weight = pKinematics->weight();
+    if( iEvent == 0 ){  
+
+      // allocate memory -- pointers are guaranteed to be NULL
+
+      m_iNParticles = pKinematics->particleList().size();
+      m_pdData = new GDouble[4*m_iNParticles*m_iNEvents];
+      m_pdWeights = new GDouble[m_iNEvents]; 
+
+      #ifdef GPU_ACCELERATION
+
+      m_gpuMan.initData( *this, !needsUserVarsOnly );
+      
+      #endif
+    }
+
+    loadDataArrayElement( pKinematics, iEvent );
+
+    float weight = m_pdWeights[iEvent];
     
     // fill some booleans that contain collective information about the weights
-    if( weight != 1 ) m_hasNonUnityWeights = true;
-    if( m_lastWeightSign == 0 ) m_lastWeightSign = weight;
-    int thisWeightSign = ( weight > 0 ? 1 : 0 );
-    thisWeightSign = ( weight < 0 ? -1 : thisWeightSign );
-    if( thisWeightSign * m_lastWeightSign < 0 ) m_hasMixedSignWeights = true;
-    m_lastWeightSign = thisWeightSign;
+    // short circuit logic for speed once the booleans are true
+    if( !m_hasNonUnityWeights && weight != 1 ) m_hasNonUnityWeights = true;
+    if( !m_hasMixedSignWeights ){
+
+      if( m_lastWeightSign == 0 ) m_lastWeightSign = weight;
+      int thisWeightSign = ( weight > 0 ? 1 : 0 );
+      thisWeightSign = ( weight < 0 ? -1 : thisWeightSign );
+      if( thisWeightSign * m_lastWeightSign < 0 ) m_hasMixedSignWeights = true;
+      m_lastWeightSign = thisWeightSign;
+    }
 
     m_dSumWeights += pKinematics->weight();
-    if (iEvent < (m_iNTrueEvents - 1)) delete pKinematics;
+    if( iEvent < ( m_iNTrueEvents - 1) ) delete pKinematics;
   }
   
-  // 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, needsUserVarsOnly );
-  }
+  // fill any remaining space in the data array with the last event's kinematics
+  for (size_t iEvent = m_iNTrueEvents; iEvent < m_iNEvents; iEvent++)
+    loadDataArrayElement(pKinematics, iEvent );
 
-  if( m_iNTrueEvents )
-  	delete pKinematics;
+  if( m_iNTrueEvents ) delete pKinematics;
    
 #ifdef GPU_ACCELERATION
 
@@ -289,14 +310,29 @@ AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly ){
   m_userVarsOffset.clear();
 }
 
+void
+AmpVecs::loadDataArrayElement( const Kinematics* pKinematics, size_t iEvent ){
+
+  // check to be sure we won't exceed the bounds of the array
+  assert( iEvent < m_iNEvents );
+
+  for( int iParticle = 0; iParticle < m_iNParticles; iParticle++ ){
+    m_pdData[4*iEvent*m_iNParticles+4*iParticle+0]=pKinematics->particle(iParticle).E();
+    m_pdData[4*iEvent*m_iNParticles+4*iParticle+1]=pKinematics->particle(iParticle).Px();
+    m_pdData[4*iEvent*m_iNParticles+4*iParticle+2]=pKinematics->particle(iParticle).Py();
+    m_pdData[4*iEvent*m_iNParticles+4*iParticle+3]=pKinematics->particle(iParticle).Pz();
+  }
+
+  m_pdWeights[iEvent] = pKinematics->weight();
+}
 
 void
-AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity, unsigned int chunkSize ){
+AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity, size_t chunkSize ){
 
-  unsigned int ampsEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize );
+  size_t ampsEvents = ( chunkSize == 0 ? m_iNEvents : chunkSize );
   
   m_iNTerms           = intenMan.getTermNames().size();
-  m_maxFactPerEvent   = intenMan.maxFactorStoragePerEvent();
+  m_maxFactPerEvent   = intenMan.maxFactorStoragePerEvent(); // in units of doubles; includes factor of 2 for complex numbers
   m_userVarsPerEvent  = intenMan.userVarsPerEvent();
   
   if( m_pdAmps!=0 || m_pdAmpFactors!=0 || m_pdUserVars!=0 || m_pdIntensity!=0 )
@@ -397,7 +433,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, 2 * m_iNEvents * m_maxFactPerEvent * sizeof(GDouble));
+  cudaMallocHost( (void**)&m_pdAmpFactors, 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 c0eb10b..b91950d 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 int m_iNEvents;
+  size_t 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 int m_iNTrueEvents;
+  size_t m_iNTrueEvents;
   
   /**
    * A double that stores the absolute value of the sum of the weights.  (For
@@ -91,25 +91,25 @@ struct AmpVecs
   /**
    * An integer that stores the number of particles in the final state.
    */
-  unsigned int m_iNParticles;
+  size_t m_iNParticles;
   
   /**
    * An integer that stores the number of amplitudes for a particular
    * configuration of the AmplitudeManager.
    */
-  unsigned int m_iNTerms;
+  size_t m_iNTerms;
   
   /**
    * An integer that is number of doubles required to store all factors
    * and permutations for any term for an event.
    */
-  unsigned int m_maxFactPerEvent;
+  size_t m_maxFactPerEvent;
   
   /**
    * An integer that is the number of doubles required to store all
    * (optionally) user-calculated data per event.
    */
-  unsigned int m_userVarsPerEvent;
+  size_t m_userVarsPerEvent;
 
   /**
    * An array of length 4 * iNEvents * iNParticles that stores the four-vectors
@@ -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 int > m_userVarsOffset;
+  map< string, size_t > m_userVarsOffset;
 
   /**
    * These booleans track features of the set of weights are are adjusted
@@ -233,7 +233,7 @@ struct AmpVecs
    */
   void allocateTerms( const IntensityManager& intenMan,
                       bool bAllocIntensity = false,
-                      unsigned int chunkSize = 0 );
+                      size_t chunkSize = 0 );
 
   /** 
    * This routine deallocates the arrays that hold the calculated terms and
@@ -270,7 +270,7 @@ struct AmpVecs
    * \see DataReader::resetSource
    * \see DataReader::getEvent
    */
-  void loadData( DataReader* pDataReader, bool needsUserVarsOnly = false );
+  void loadData( DataReader* pDataReader, bool needsUserVarsOnly = false, size_t chunkSize = 0 );
   
   /**
    * This routine fills the arrays of data and weights event by event
@@ -286,8 +286,8 @@ struct AmpVecs
    *
    * \see loadData
    */
-  void loadEvent( const Kinematics* pKinematics, unsigned int iEvent = 0,
-                  unsigned int iNTrueEvents = 1, bool needsUserVarsOnly = false );
+  void loadEvent( const Kinematics* pKinematics, size_t iEvent = 0,
+                  size_t iNTrueEvents = 1, bool needsUserVarsOnly = false );
   
   /**
    * A helper routine to get an event i from the array of data and weights.
@@ -333,6 +333,8 @@ struct AmpVecs
   AmpVecs* m_sharedDataHost;
   
 private:
+
+  void loadDataArrayElement( const Kinematics* pKinematics, size_t iEvent );
   
   int m_lastWeightSign;
   set< AmpVecs* > m_sharedDataFriends;
diff --git a/AmpTools/IUAmpTools/Amplitude.cc b/AmpTools/IUAmpTools/Amplitude.cc
index 8e6fe46..a0b2b03 100644
--- a/AmpTools/IUAmpTools/Amplitude.cc
+++ b/AmpTools/IUAmpTools/Amplitude.cc
@@ -109,7 +109,7 @@ SCOREP_USER_REGION_BEGIN( calcUserVarsAll, "calcUserVarsAll", SCOREP_USER_REGION
       // doesn't need to deal with permutations in their calcAmplitude
       // routine
 
-      unsigned int eventOffset = 4*iNParticles*iEvent;
+      size_t eventOffset = 4*iNParticles*iEvent;
       
       for( i = 0; i < iNParticles; i++ ){
 
@@ -117,7 +117,7 @@ SCOREP_USER_REGION_BEGIN( calcUserVarsAll, "calcUserVarsAll", SCOREP_USER_REGION
         pKin[i] = &(pdData[eventOffset+4*j]);
       }
 
-      unsigned int userIndex = iNEvents*iPermutation*numVars + iEvent*numVars;
+      size_t userIndex = iNEvents*iPermutation*numVars + iEvent*numVars;
       calcUserVars( pKin, &(pdUserVars[userIndex]) );
     }
   }
@@ -130,10 +130,10 @@ SCOREP_USER_REGION_END( calcUserVarsAll )
 }
 
 void
-Amplitude::calcAmplitudeAll( GDouble* pdData, GDouble* pdAmpFact, unsigned int iNEvents,
+Amplitude::calcAmplitudeAll( GDouble* pdData, GDouble* pdAmpFact, size_t iNEvents,
                              const vector< vector< int > >* pvPermutations,
-                             GDouble* pdUserVars, unsigned int startEvent,
-                             unsigned int chunkSize ) const
+                             GDouble* pdUserVars, size_t startEvent,
+                             size_t chunkSize ) const
 {
 
 #ifdef SCOREP
@@ -155,7 +155,7 @@ SCOREP_USER_REGION_BEGIN( calcAmplitudeAll, "calcAmplitudeAll", SCOREP_USER_REGI
   
   // 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 );
+  size_t nEvents = ( chunkSize == 0 ? iNEvents : chunkSize );
   
   int i, j, iEvent;
   for( iPermutation = 0; iPermutation < iNPermutations; iPermutation++ ){
@@ -169,7 +169,7 @@ SCOREP_USER_REGION_BEGIN( calcAmplitudeAll, "calcAmplitudeAll", SCOREP_USER_REGI
       // doesn't need to deal with permutations in their calcAmplitude
       // routine
 
-      unsigned int eventOffset = 4*iNParticles*(iEvent + startEvent);
+      size_t eventOffset = 4*iNParticles*(iEvent + startEvent);
 
       for( i = 0; i < iNParticles; i++ ){
 
@@ -177,7 +177,7 @@ SCOREP_USER_REGION_BEGIN( calcAmplitudeAll, "calcAmplitudeAll", SCOREP_USER_REGI
         pKin[i] = &(pdData[eventOffset+4*j]);
       }
 
-      unsigned int userIndex = iNEvents*iPermutation*numVars + (iEvent + startEvent)*numVars;
+      size_t userIndex = iNEvents*iPermutation*numVars + (iEvent + startEvent)*numVars;
 
       if( numVars != 0 ){
       
diff --git a/AmpTools/IUAmpTools/Amplitude.h b/AmpTools/IUAmpTools/Amplitude.h
index f5d90e9..2cf993d 100644
--- a/AmpTools/IUAmpTools/Amplitude.h
+++ b/AmpTools/IUAmpTools/Amplitude.h
@@ -292,10 +292,10 @@ class Amplitude
    * \see calcAmplitudeAll
    * \see AmplitudeManager::addAmpPermutation
    */
-  virtual void calcAmplitudeAll( GDouble* pdData, GDouble* pdAmpFact, unsigned int iNEvents,
+  virtual void calcAmplitudeAll( GDouble* pdData, GDouble* pdAmpFact, size_t iNEvents,
                                 const vector< vector< int > >* pvPermutations,
-                                GDouble* pdUserVars = 0, unsigned int startEvent = 0,
-                                unsigned int chunkSize = 0 ) const;
+                                GDouble* pdUserVars = 0, size_t startEvent = 0,
+                                size_t chunkSize = 0 ) const;
   
   
   /**
diff --git a/AmpTools/IUAmpTools/AmplitudeManager.cc b/AmpTools/IUAmpTools/AmplitudeManager.cc
index 83c0fc7..1f9ca72 100644
--- a/AmpTools/IUAmpTools/AmplitudeManager.cc
+++ b/AmpTools/IUAmpTools/AmplitudeManager.cc
@@ -62,11 +62,11 @@ IntensityManager( reaction, reactionName )
   // group the switches by particle type
   // dump out some information
   map< string, vector< pair< int, int > > > swapsByType;
-  for( unsigned int i = 0; i < reaction.size(); ++i ){
+  for( size_t i = 0; i < reaction.size(); ++i ){
     
     report( INFO, kModule ) << "\t particle index assignment:  " << reaction[i] << " -> " << i << endl;
     
-    for( unsigned int j = i + 1; j < reaction.size(); ++j ){
+    for( size_t j = i + 1; j < reaction.size(); ++j ){
       
       if( reaction[i] == reaction[j] ){
         
@@ -95,7 +95,7 @@ IntensityManager( reaction, reactionName )
   // it with numberOfCombos copies of the default ordering
   // then go in and make the swaps
   vector< int > defaultOrder( reaction.size() );
-  for( unsigned int i = 0; i < reaction.size(); ++i ){
+  for( size_t i = 0; i < reaction.size(); ++i ){
     
     defaultOrder[i] = i;
   }
@@ -127,9 +127,9 @@ IntensityManager( reaction, reactionName )
                             << " 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 ){
+    for( size_t i = 0; i < m_symmCombos.size(); ++i ){
       
-      for( unsigned int j = 0; j < reaction.size(); ++j ){
+      for( size_t j = 0; j < reaction.size(); ++j ){
         
         report( INFO, kModule ) << "\t" << m_symmCombos[i][j];
       }
@@ -148,17 +148,17 @@ AmplitudeManager::~AmplitudeManager() {
   }
 }
 
-unsigned int
+size_t
 AmplitudeManager::maxFactorStoragePerEvent() const {
   
   vector< string > ampNames = getTermNames();
   
-  unsigned int nAmpFactorsAndPerms = 0;
+  size_t nAmpFactorsAndPerms = 0;
   
   for( int i = 0; i < getTermNames().size(); i++ ) {
     
-    unsigned int iNPermutations = getPermutations( ampNames[i] ).size();
-    unsigned int iNFactors = getFactors( ampNames[i] ).size();
+    size_t iNPermutations = getPermutations( ampNames[i] ).size();
+    size_t iNFactors = getFactors( ampNames[i] ).size();
     
     assert( iNPermutations*iNFactors );
     
@@ -172,7 +172,7 @@ AmplitudeManager::maxFactorStoragePerEvent() const {
   return 2 * nAmpFactorsAndPerms;
 }
 
-unsigned int
+size_t
 AmplitudeManager::termStoragePerEvent() const {
   
   // for each amplitude we need to store a complex
@@ -181,7 +181,7 @@ AmplitudeManager::termStoragePerEvent() const {
   return 2 * getTermNames().size();
 }
 
-unsigned int
+size_t
 AmplitudeManager::userVarsPerEvent() const {
   
   set< string > countedStaticAmps;
@@ -189,11 +189,11 @@ AmplitudeManager::userVarsPerEvent() const {
   
   vector< string > ampNames = getTermNames();
   
-  unsigned int userStorage = 0;
+  size_t userStorage = 0;
   
   for( int i = 0; i < getTermNames().size(); i++ ) {
 
-    unsigned int iNPermutations = getPermutations( ampNames[i] ).size();
+    size_t iNPermutations = getPermutations( ampNames[i] ).size();
     vector< const Amplitude* > factorVec =
       m_mapNameToAmps.find( ampNames[i] )->second;
 
@@ -228,13 +228,13 @@ AmplitudeManager::userVarsPerEvent() const {
   return userStorage;
 }
 
-unsigned int
+size_t
 AmplitudeManager::uniqueNIElements() const {
 
   // find the number of nontrival amplitude products
   // in the expression for the intensity
   
-  unsigned int num = 0;
+  size_t num = 0;
   int N = m_sumCoherently.size();
   for( int i = 0; i < N; ++i ){
     for( int j = 0; j <= i; ++j ){
@@ -274,7 +274,7 @@ SCOREP_USER_REGION_BEGIN( calcUserVars, "calcUserVars", SCOREP_USER_REGION_TYPE_
   int iNAmps = ampNames.size();
 
   int iAmpIndex;
-  unsigned int iUserVarsOffset = 0;
+  size_t iUserVarsOffset = 0;
   for( iAmpIndex = 0; iAmpIndex < iNAmps; iAmpIndex++ )
   {
     
@@ -301,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 int thisOffset = 0;
+      size_t thisOffset = 0;
 
       if( pCurrAmp->areUserVarsStatic() ){
         
@@ -310,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 int >::const_iterator offsetItr =
+        map< string, size_t >::const_iterator offsetItr =
         a.m_userVarsOffset.find( pCurrAmp->name() );
         
         if( offsetItr == a.m_userVarsOffset.end() ){
@@ -338,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 int >::const_iterator offsetItr =
+        map< string, size_t >::const_iterator offsetItr =
         a.m_userVarsOffset.find( pCurrAmp->identifier() );
         
         if( offsetItr == a.m_userVarsOffset.end() ){
@@ -385,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 int cpuIndex =
+            size_t cpuIndex =
               thisOffset + iPerm*a.m_iNEvents*iNVars + iEvt*iNVars + iVar;
-            unsigned int gpuIndex =
+            size_t gpuIndex =
               iPerm*a.m_iNEvents*iNVars + iVar*a.m_iNEvents + iEvt;
             
             tmpVarStorage[gpuIndex] = a.m_pdUserVars[cpuIndex];
@@ -429,7 +429,7 @@ AmplitudeManager::calcTerms( AmpVecs& a ) const
 }
 
 vector<bool>
-AmplitudeManager::calcTerms( AmpVecs& a, unsigned int startEvent, unsigned int chunkSize ) const
+AmplitudeManager::calcTerms( AmpVecs& a, size_t startEvent, size_t chunkSize ) const
 {
   
 #ifdef SCOREP
@@ -437,7 +437,7 @@ 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 );
+  size_t nEvents = ( chunkSize == 0 ? a.m_iNEvents : chunkSize );
   
   report( DEBUG, kModule ) << "Calculating terms...     termsValid = "
   << a.m_termsValid << endl;
@@ -523,7 +523,7 @@ 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 uAmpFactOffset = 0;
+    size_t uAmpFactOffset = 0;
     for( iFactor=0; iFactor < iNFactors;
          iFactor++, uAmpFactOffset += 2 * nEvents * iNPermutations ){
       
@@ -531,7 +531,7 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON
       
       // if we have static user data, look up the location in the data array
       // if not, then look up by identifier
-      unsigned int userVarsOffset =
+      size_t userVarsOffset =
       ( pCurrAmp->areUserVarsStatic() ?
         a.m_userVarsOffset[pCurrAmp->name()] :
         a.m_userVarsOffset[pCurrAmp->identifier()] );
@@ -557,9 +557,10 @@ 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 int iOffsetA, iOffsetP, iOffsetF;
+    size_t iOffsetA, iOffsetP, iOffsetF;
         
-    // zeroing out the entire range
+    // zeroing out the entire chunk where we will accumulate
+    // compute symmetrized amplitude for this particular amplitude index
     memset( (void*)( a.m_pdAmps + 2 * nEvents * iAmpIndex ), 0,
            2 * nEvents * sizeof(GDouble) );
     
@@ -594,7 +595,7 @@ SCOREP_USER_REGION_BEGIN( calcTerms, "calcTerms", SCOREP_USER_REGION_TYPE_COMMON
       a.m_pdAmps[iOffsetA]   *= dSymmFactor;
       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;
     
@@ -810,8 +811,8 @@ SCOREP_USER_REGION_END( calcSumLogIntensity )
 
 
 void
-AmplitudeManager::calcIntegrals( AmpVecs& a, unsigned int iNGenEvents,
-                                 unsigned int chunkSize ) const
+AmplitudeManager::calcIntegrals( AmpVecs& a, size_t iNGenEvents,
+                                 size_t chunkSize ) const
 {
 #ifdef SCOREP
 SCOREP_USER_REGION_DEFINE( calcIntegralsA )                                                                                    
@@ -826,7 +827,7 @@ 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 : a.m_iNEvents / chunkSize );
+  size_t nChunk = ( chunkSize == 0 ? 1 : a.m_iNEvents / chunkSize );
   if( ( chunkSize !=0 ) && ( a.m_iNEvents % chunkSize != 0 ) ) ++nChunk;
   
   int iNAmps = a.m_iNTerms;
@@ -844,20 +845,18 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T
   // a variable to track how many elements need computing
   int nCompute;
 
-  for( unsigned int iChunk = 0; iChunk < nChunk; ++iChunk ){
+  for( size_t iChunk = 0; iChunk < nChunk; ++iChunk ){
         
-    unsigned int startEvent = iChunk * chunkSize;
-    unsigned int nEvents = a.m_iNEvents;
+    size_t startEvent = iChunk * chunkSize;
+    size_t nEvents = a.m_iNEvents;
+
     // for  chunked 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;
-    }
+    if( chunkSize != 0 ) nEvents = min( chunkSize, a.m_iNEvents - startEvent );      
     
     // this returns a vector indicating which terms have changed
-    const vector<bool>& termChanged = calcTerms( a, startEvent, chunkSize );
+    const vector<bool>& termChanged = calcTerms( a, startEvent, nEvents );
     
     // if we are doing a "chunked" calculation, we should invalidate the
     // termsValid boolean since the terms array will only hold a fraction
@@ -957,7 +956,8 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T
     }
     
 #else
-      
+    
+    // C++ standard sets initial value to zero
     vector<double> temp(maxNIElements*2);
 
     // use the GPU manager to compute the result for the chunk
@@ -976,7 +976,7 @@ SCOREP_USER_REGION_BEGIN( calcIntegralsA, "calcIntegralsA", SCOREP_USER_REGION_T
     }
 #endif
     if( nChunk > 1 ) report( DEBUG, kModule ) << "Integral calculation:  chunk " 
-                     << iChunk << " of " << nChunk << " complete." << endl;
+                     << iChunk+1 << " of " << nChunk << " complete." << endl;
   }
   
   // in a "chunked" calculation nCompute, iIndex, and jIndex get reset
@@ -1177,7 +1177,7 @@ AmplitudeManager::setupFromConfigurationInfo( const ConfigurationInfo* configInf
   
   // loop over amplitudes in the ConfigurationInfo
   vector<AmplitudeInfo*> ampInfoVector = configInfo->amplitudeList(reactionName());
-  for (unsigned int i = 0; i < ampInfoVector.size(); i++){
+  for (size_t i = 0; i < ampInfoVector.size(); i++){
     
     string ampName = ampInfoVector[i]->fullName();
     string sumName = ampInfoVector[i]->sumName();
@@ -1185,7 +1185,7 @@ AmplitudeManager::setupFromConfigurationInfo( const ConfigurationInfo* configInf
     
     // add amplitudes
     vector< vector<string> > ampFactors = ampInfoVector[i]->factors();
-    for (unsigned int j = 0; j < ampFactors.size(); j++){
+    for (size_t j = 0; j < ampFactors.size(); j++){
       string factorName = ampFactors[j][0];
       vector<string> ampParameters = ampFactors[j];
       ampParameters.erase(ampParameters.begin());
@@ -1194,7 +1194,7 @@ AmplitudeManager::setupFromConfigurationInfo( const ConfigurationInfo* configInf
     
     // add permutations
     vector< vector<int> > permutations = ampInfoVector[i]->permutations();
-    for (unsigned int j = 0; j < permutations.size(); j++){
+    for (size_t j = 0; j < permutations.size(); j++){
       addAmpPermutation( ampName, permutations[j] );
     }
     
diff --git a/AmpTools/IUAmpTools/AmplitudeManager.h b/AmpTools/IUAmpTools/AmplitudeManager.h
index c9915be..e281129 100644
--- a/AmpTools/IUAmpTools/AmplitudeManager.h
+++ b/AmpTools/IUAmpTools/AmplitudeManager.h
@@ -118,7 +118,7 @@ class AmplitudeManager : public IntensityManager
    * of number of factors * number of permutations for any term.
    */
   
-  unsigned int maxFactorStoragePerEvent() const;
+  size_t maxFactorStoragePerEvent() const;
   
   /**
    * This function returns the number of doubles needed to store all complete
@@ -126,14 +126,14 @@ class AmplitudeManager : public IntensityManager
    * the size of a double.
    */
   
-  unsigned int termStoragePerEvent() const;
+  size_t termStoragePerEvent() const;
   
   /**
    * This function returns the number of doubles needed to store optional
    * user data for all factors and permutations.
    */
   
-  unsigned int userVarsPerEvent() const;
+  size_t userVarsPerEvent() const;
   
   /**
    * This function returns the total number of unique (up to complex
@@ -141,7 +141,7 @@ class AmplitudeManager : public IntensityManager
    *  to index functions useds for NI calculation.
    */
   
-  unsigned int uniqueNIElements() const;
+  size_t uniqueNIElements() const;
   
   /**
    * This function triggers the calculation of optional user data that
@@ -230,8 +230,8 @@ class AmplitudeManager : public IntensityManager
    * \see calcAmplitudes
    * \see calcIntensities
    */
-  void calcIntegrals( AmpVecs& ampVecs, unsigned int iNGenEvents,
-                     unsigned int chunkSize = 0 ) const;
+  void calcIntegrals( AmpVecs& ampVecs, size_t iNGenEvents,
+                     size_t chunkSize = 0 ) const;
 
   /**
    * The function returns a list of permutations that will be performed on
@@ -456,8 +456,8 @@ class AmplitudeManager : public IntensityManager
    * \see calcSumLogIntensity
    * \see calcIntegrals
    */
-  vector<bool> calcTerms( AmpVecs& ampVecs, unsigned int startEvent,
-                          unsigned int chunkSize ) const;
+  vector<bool> calcTerms( AmpVecs& ampVecs, size_t startEvent,
+                          size_t chunkSize ) const;
   
   
   // recursive routine to symmetrize final state
@@ -494,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 int > m_staticUserVarsOffset;
+  mutable map< string, size_t > m_staticUserVarsOffset;
   
   static const char* kModule;
 };
diff --git a/AmpTools/IUAmpTools/IntensityManager.h b/AmpTools/IUAmpTools/IntensityManager.h
index 49f5fb8..3dabfcb 100644
--- a/AmpTools/IUAmpTools/IntensityManager.h
+++ b/AmpTools/IUAmpTools/IntensityManager.h
@@ -76,7 +76,7 @@ class IntensityManager {
    * memory allocation.
    */
   
-  virtual unsigned int maxFactorStoragePerEvent() const = 0;
+  virtual size_t maxFactorStoragePerEvent() const = 0;
   
   /**
    * This function should return the number of doubles required to store
@@ -89,7 +89,7 @@ class IntensityManager {
    * factor) remains constant throughout the fit.
    */
   
-  virtual unsigned int termStoragePerEvent() const = 0;
+  virtual size_t termStoragePerEvent() const = 0;
   
   /**
    * This function should return the number of doubles required to
@@ -97,7 +97,7 @@ class IntensityManager {
    * for a single event.
    */
   
-  virtual unsigned int userVarsPerEvent() const = 0;
+  virtual size_t userVarsPerEvent() const = 0;
 
   /*
    * These functions perform computations based on the current state
@@ -169,8 +169,8 @@ class IntensityManager {
    * \see calcAmplitudes
    * \see calcIntensities
    */
-  virtual void calcIntegrals( AmpVecs& ampVecs, unsigned int iNGenEvents,
-                              unsigned int chunkSize = 0 ) const = 0;
+  virtual void calcIntegrals( AmpVecs& ampVecs, size_t iNGenEvents,
+                              size_t chunkSize = 0 ) const = 0;
  
   /**
    * This function calculates the intensity for one event using a Kinematics
diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc
index 3ba9aff..c04ba1a 100644
--- a/AmpTools/IUAmpTools/NormIntInterface.cc
+++ b/AmpTools/IUAmpTools/NormIntInterface.cc
@@ -396,21 +396,16 @@ NormIntInterface::forceCacheUpdate( bool normIntOnly ) const
     // now we need to have the generated MC in addition to the accepted
     // 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, false, chunkSize );
+    if( m_genMCVecs.m_iNTerms == 0 ) m_genMCVecs.allocateTerms( *m_pIntenManager, false, m_chunkSize );
     
     report( DEBUG, kModule ) << "Asking IntensityManager to calculate integrals "
     << "using the generated MC." << endl;
     
-    m_pIntenManager->calcIntegrals( m_genMCVecs, m_nGenEvents, chunkSize );
+    m_pIntenManager->calcIntegrals( m_genMCVecs, m_nGenEvents, m_chunkSize );
     
     setAmpIntMatrix( m_genMCVecs.m_pdIntegralMatrix );
 
@@ -533,18 +528,21 @@ NormIntInterface::setNormIntMatrix( const double* input ) const {
 
 #ifndef __ACLIC__
 
-unsigned int
-NormIntInterface::genMCChunkSize() const {
+size_t
+NormIntInterface::genMCChunkSize( size_t nGen ) 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;
-  
+  assert( m_accMCVecs.m_dataLoaded );
+  size_t nAcc = m_accMCVecs.m_iNEvents;
+
+  if( nGen == 0 ){
+    assert( m_genMCVecs.m_dataLoaded );
+    nGen = m_genMCVecs.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
+  // 1/4 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;
@@ -557,40 +555,49 @@ NormIntInterface::genMCChunkSize() const {
 
 void
 NormIntInterface::loadMC() const {
+
+  std::map<DataReader*,AmpVecs*>::iterator accVecs = m_uniqueDataSets.find( m_accMCReader );
+  if( accVecs == m_uniqueDataSets.end() ){
   
-  std::map<DataReader*,AmpVecs*>::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() );
+    report( INFO, kModule ) << "Loading accepted Monte Carlo from file..." << endl;
+    m_accMCVecs.loadData( m_accMCReader, m_pIntenManager->needsUserVarsOnly() );
     
-    m_uniqueDataSets[m_genMCReader] = &m_genMCVecs;
+    m_uniqueDataSets[m_accMCReader] = &m_accMCVecs;
   }
   else{
     
     report( NOTICE, kModule ) << "Duplicated Monte Carlo set detected, "
          << "using previously loaded version" << endl;
     
-    genVecs->second->shareDataWith( &m_genMCVecs, m_pIntenManager->needsUserVarsOnly() );
+    accVecs->second->shareDataWith( &m_accMCVecs, m_pIntenManager->needsUserVarsOnly() );
   }
-  m_nGenEvents = m_genMCVecs.m_iNTrueEvents;
+  m_sumAccWeights = m_accMCVecs.m_dSumWeights;
 
-  std::map<DataReader*,AmpVecs*>::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() );
+   // 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
+#ifdef GPU_ACCELERATION
+  m_chunkSize = genMCChunkSize( GPUManager::calcNEventsGPU( m_genMCReader->numEvents() ) );
+#else
+  m_chunkSize = genMCChunkSize( m_genMCReader->numEvents() );
+#endif 
+
+  std::map<DataReader*,AmpVecs*>::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_chunkSize );
     
-    m_uniqueDataSets[m_accMCReader] = &m_accMCVecs;
+    m_uniqueDataSets[m_genMCReader] = &m_genMCVecs;
   }
   else{
     
     report( NOTICE, kModule ) << "Duplicated Monte Carlo set detected, "
          << "using previously loaded version" << endl;
     
-    accVecs->second->shareDataWith( &m_accMCVecs, m_pIntenManager->needsUserVarsOnly() );
+    genVecs->second->shareDataWith( &m_genMCVecs, m_pIntenManager->needsUserVarsOnly() );
   }
-  m_sumAccWeights = m_accMCVecs.m_dSumWeights;
+  m_nGenEvents = m_genMCVecs.m_iNTrueEvents;
 }
 
 void
diff --git a/AmpTools/IUAmpTools/NormIntInterface.h b/AmpTools/IUAmpTools/NormIntInterface.h
index dbe8774..624ce4a 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 int numGenEvents() const { return m_nGenEvents; }
+  size_t 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 int events ) { m_nGenEvents = events; }
+  void setGenEvents( size_t events ) { m_nGenEvents = events; }
   void setAccEvents( double sumWeights ) { m_sumAccWeights = sumWeights; }
   
 protected:
@@ -124,7 +124,7 @@ class NormIntInterface
   void initializeCache();
   int m_cacheSize;
   
-  unsigned int genMCChunkSize() const;
+  size_t genMCChunkSize( size_t nGen = 0 ) const;
   
   vector< string > m_termNames;
   map< string, int > m_termIndex;
@@ -135,7 +135,7 @@ class NormIntInterface
   mutable bool m_emptyNormIntCache;
   mutable bool m_emptyAmpIntCache;
   
-  mutable unsigned int m_nGenEvents;
+  mutable size_t m_nGenEvents;
   mutable double m_sumAccWeights;
   
 #ifndef __ACLIC__
@@ -148,6 +148,7 @@ class NormIntInterface
   // caches for MC
   mutable AmpVecs m_accMCVecs;
   mutable AmpVecs m_genMCVecs;
+  mutable size_t m_chunkSize;
   
   static map< DataReader*, AmpVecs* > m_uniqueDataSets;
 #endif
diff --git a/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc b/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc
index a9201d9..01ef82d 100644
--- a/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc
+++ b/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc
@@ -95,23 +95,22 @@ NormIntInterfaceMPI::setupMPI()
   
   m_isLeader = ( m_rank == 0 );
   
-  // this is unsigned elsewhere, but MPI
-  unsigned int totalGenEvents = 0;
+  long int totalGenEvents = 0;
   double totalAccWeights = 0;
     
   if( m_isLeader ){
     
     for( int i = 1; i < m_numProc; ++i ){
       
-      unsigned int thisEvents;
+      size_t thisEvents;
       double thisWeights;
 
       // trigger sending of events from followers -- data is irrelevant
-      MPI_Send( &thisEvents, 1, MPI_UNSIGNED, i, MPITag::kAcknowledge,
+      MPI_Send( &thisEvents, 1, MPI_LONG, i, MPITag::kAcknowledge,
                MPI_COMM_WORLD );
       
       // now receive actual data
-      MPI_Recv( &thisEvents, 1, MPI_UNSIGNED, i, MPITag::kIntSend,
+      MPI_Recv( &thisEvents, 1, MPI_LONG, i, MPITag::kIntSend,
                MPI_COMM_WORLD, &status );
       totalGenEvents += thisEvents;
       
@@ -120,7 +119,7 @@ NormIntInterfaceMPI::setupMPI()
       totalAccWeights += thisWeights;
       
       // send acknowledgment 
-      MPI_Send( &thisEvents, 1, MPI_UNSIGNED, i, MPITag::kAcknowledge,
+      MPI_Send( &thisEvents, 1, MPI_LONG, i, MPITag::kAcknowledge,
                MPI_COMM_WORLD );
     }
     
@@ -133,7 +132,7 @@ NormIntInterfaceMPI::setupMPI()
     // a copy from the DataReaderMPI cache to the AmpVecs structure
     loadMC();
 
-    unsigned int thisEvents;
+    long int thisEvents;
     double thisWeights;
 
     // if we are not the leader, send generated and accepted events
@@ -144,16 +143,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, 0, MPITag::kAcknowledge, MPI_COMM_WORLD,
+    MPI_Recv( &thisEvents, 1, MPI_LONG, 0, MPITag::kAcknowledge, MPI_COMM_WORLD,
               &status );
 
     thisEvents = numGenEvents();
-    MPI_Send( &thisEvents, 1, MPI_UNSIGNED, 0, MPITag::kIntSend, MPI_COMM_WORLD );
+    MPI_Send( &thisEvents, 1, MPI_LONG, 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, 0, MPITag::kAcknowledge, MPI_COMM_WORLD,
+    MPI_Recv( &thisEvents, 1, MPI_LONG, 0, MPITag::kAcknowledge, MPI_COMM_WORLD,
              &status );
   }
 }
@@ -184,8 +183,8 @@ NormIntInterfaceMPI::sumIntegrals( IntType type ) const
 
   // now broadcast the total number of events from the leader to the
   // followers so that they may renormalize the sum properly
-  int totalEvents = numGenEvents();
-  MPI_Bcast( &totalEvents, 1, MPI_INT, 0, MPI_COMM_WORLD );
+  long int totalEvents = numGenEvents();
+  MPI_Bcast( &totalEvents, 1, MPI_LONG, 0, MPI_COMM_WORLD );
   
   // and renormalize the sum
   for( int i = 0; i < cacheSize(); ++i ) result[i] /= totalEvents;