Adding error message for large patches on nbond kernels 82/4782/2
authorJulio Maia <jmaia@ks.uiuc.edu>
Wed, 7 Nov 2018 17:29:14 +0000 (11:29 -0600)
committerJulio Maia <jmaia@ks.uiuc.edu>
Wed, 7 Nov 2018 19:39:43 +0000 (13:39 -0600)
We had an 'invalid argument' error when large patches (>11k atoms) were
accounted for on CudaTileList::buildTileLists. This happened because
the amount of shMem allocated inside the tilelist kernels
is dependant on the number of atoms in the largest patch.

I've added a query to get the maximum shared memory available per block.
We can compare that against the necessary shared memory to launch the
tilelist kernels, and if those launches are about to go wrong, NAMD prints
a nice message to warn the user about the patches sizes.

Change-Id: Ie35b40d83f87b06ec33cd60b0a1936cec18b6114

src/CudaComputeNonbonded.C
src/CudaComputeNonbonded.h
src/CudaTileListKernel.cu
src/CudaTileListKernel.h

index 5c0bee8..63cd619 100644 (file)
@@ -34,7 +34,7 @@ extern "C" void CcdCallBacksReset(void *ignored, double curWallTime);  // fix Ch
 // Class constructor
 //
 CudaComputeNonbonded::CudaComputeNonbonded(ComputeID c, int deviceID,
-  CudaNonbondedTables& cudaNonbondedTables, bool doStreaming) : 
+  CudaNonbondedTables& cudaNonbondedTables, bool doStreaming) :
 Compute(c), deviceID(deviceID), doStreaming(doStreaming), nonbondedKernel(deviceID, cudaNonbondedTables, doStreaming),
 tileListKernel(deviceID, doStreaming), GBISKernel(deviceID) {
 
@@ -81,7 +81,7 @@ tileListKernel(deviceID, doStreaming), GBISKernel(deviceID) {
   dEdaSumHSize = 0;
   dHdrPrefixH = NULL;
   dHdrPrefixHSize = 0;
-
+  maxShmemPerBlock = 0;
   cudaPatches = NULL;
 
   atomsChangedIn = true;
@@ -542,7 +542,7 @@ void CudaComputeNonbonded::assignPatches(ComputeMgr* computeMgrIn) {
   // Set counters of avoidable PEs to high numbers
   for (int i=0;i < pesToAvoid.size();i++) {
     int pe = pesToAvoid[i];
-    peProxyPatchCounter[CkRankOf(pe)] = (1 << 20);    
+    peProxyPatchCounter[CkRankOf(pe)] = (1 << 20);
   }
 #endif
   // Avoid master Pe somewhat
@@ -613,6 +613,12 @@ void CudaComputeNonbonded::initialize() {
     allocate_host<VirialEnergy>(&h_virialEnergy, 1);
     allocate_device<VirialEnergy>(&d_virialEnergy, 1);
 
+  /* JM: Queries for maximum sharedMemoryPerBlock on deviceID
+   */
+   cudaDeviceProp props;
+   cudaCheck(cudaGetDeviceProperties(&props, deviceID)); //Gets properties of 'deviceID device'
+   maxShmemPerBlock = props.sharedMemPerBlock;
+
 #if CUDA_VERSION >= 5050
     int leastPriority, greatestPriority;
     cudaCheck(cudaDeviceGetStreamPriorityRange(&leastPriority, &greatestPriority));
@@ -629,7 +635,7 @@ void CudaComputeNonbonded::initialize() {
     lock = CmiCreateLock();
 
     reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
-  }  
+  }
 }
 
 //
@@ -899,7 +905,7 @@ void CudaComputeNonbonded::reallocateArrays() {
     reallocate_host<float4>(&h_forcesSlow, &h_forcesSlowSize, atomStorageSize, 1.4f);
   }
   reallocate_device<float4>(&d_forces, &d_forcesSize, atomStorageSize, 1.4f);
-  reallocate_device<float4>(&d_forcesSlow, &d_forcesSlowSize, atomStorageSize, 1.4f);  
+  reallocate_device<float4>(&d_forcesSlow, &d_forcesSlowSize, atomStorageSize, 1.4f);
 
   if (simParams->GBISOn) {
     reallocate_host<float>(&intRad0H, &intRad0HSize, atomStorageSize, 1.2f);
@@ -1050,9 +1056,9 @@ void CudaComputeNonbonded::launchWork() {
     if (gbisPhase == 1) {
       doGBISphase1();
     } else if (gbisPhase == 2) {
-      doGBISphase2(); 
+      doGBISphase2();
     } else if (gbisPhase == 3) {
-      doGBISphase3(); 
+      doGBISphase3();
     }
   }
 
@@ -1081,7 +1087,7 @@ void CudaComputeNonbonded::launchWork() {
 //
 void CudaComputeNonbonded::doGBISphase1() {
   cudaCheck(cudaSetDevice(deviceID));
-  
+
   if (atomsChanged) {
     GBISKernel.updateIntRad(atomStorageSize, intRad0H, intRadSH, stream);
   }
@@ -1159,7 +1165,7 @@ void CudaComputeNonbonded::doForce() {
     // Build initial tile lists and sort
     tileListKernel.buildTileLists(numTileLists, patches.size(), atomStorageSize,
       maxTileListLen, lata, latb, latc,
-      cudaPatches, (const float4*)atoms, plcutoff2, stream);
+      cudaPatches, (const float4*)atoms, plcutoff2, maxShmemPerBlock, stream);
     // Prepare tile list for atom-based refinement
     tileListKernel.prepareTileList(stream);
   }
@@ -1205,7 +1211,7 @@ void CudaComputeNonbonded::finishReductions() {
 
   if (CkMyPe() != masterPe)
     NAMD_bug("CudaComputeNonbonded::finishReductions() called on non masterPe");
-  
+
   // fprintf(stderr, "%d finishReductions doSkip %d doVirial %d doEnergy %d\n", CkMyPe(), doSkip, doVirial, doEnergy);
 
   if (!doSkip) {
@@ -1415,7 +1421,7 @@ void CudaComputeNonbonded::finishTimers() {
       traceUserBracketEvent(CUDA_GBIS3_KERNEL_EVENT, beforeForceCompute, CkWallTimer());
   } else {
     traceUserBracketEvent(CUDA_NONBONDED_KERNEL_EVENT, beforeForceCompute, CkWallTimer());
-  }  
+  }
 }
 
 //
@@ -1509,7 +1515,7 @@ void CudaComputeNonbonded::forceDoneCheck(void *arg, double walltime) {
     cudaDie(errmsg,cudaSuccess);
   }
 
-  // Call again 
+  // Call again
   CcdCallBacksReset(0, walltime);
   CcdCallFnAfter(forceDoneCheck, arg, 0.1);
 }
@@ -1547,7 +1553,7 @@ struct cr_sortop_distance {
 };
 
 static inline bool sortop_bitreverse(int a, int b) {
-  if ( a == b ) return 0; 
+  if ( a == b ) return 0;
   for ( int bit = 1; bit; bit *= 2 ) {
     if ( (a&bit) != (b&bit) ) return ((a&bit) < (b&bit));
   }
@@ -1637,7 +1643,7 @@ void CudaComputeNonbonded::buildExclusions() {
 #ifdef MEM_OPT_VERSION
   int natoms = mol->exclSigPoolSize;
 #else
-  int natoms = mol->numAtoms; 
+  int natoms = mol->numAtoms;
 #endif
 
        if (exclusionsByAtom != NULL) delete [] exclusionsByAtom;
index 5a3f078..f38dca1 100644 (file)
@@ -86,6 +86,7 @@ private:
   bool computesChanged;
 
   const int deviceID;
+  size_t maxShmemPerBlock;
   cudaStream_t stream;
 
   // PME and VdW CUDA kernels
index 810f348..b8a7c94 100644 (file)
@@ -334,7 +334,7 @@ __global__ void calcTileListPosKernel(const int numComputes,
     BlockScan(tempStorage).ExclusiveSum(numTiles1, tilePosVal);
 
     // Store into global memory
-    if (k < numComputes) {      
+    if (k < numComputes) {
       tilePos[k] = shTilePos0 + tilePosVal;
     }
 
@@ -483,7 +483,7 @@ buildTileListsBBKernel(const int numTileLists,
     numJtiles  = cub::ShuffleIndex(numJtiles,  WARPSIZE-1, WARPSIZE, WARP_FULL_MASK);
     jtileStart = cub::ShuffleIndex(jtileStart, WARPSIZE-1, WARPSIZE, WARP_FULL_MASK);
     if (jtileStart + numJtiles > tileJatomStartSize) {
-      // tileJatomStart out of memory, exit 
+      // tileJatomStart out of memory, exit
       if (wid == 0) tileListStat->tilesSizeExceeded = true;
       return;
     }
@@ -540,7 +540,7 @@ repackTileListsKernel(const int numTileLists, const int begin_bit, const int* __
   const int wid = threadIdx.x % WARPSIZE;
 
   // One warp does one tile list
-  for (int i = threadIdx.x/WARPSIZE + blockDim.x/WARPSIZE*blockIdx.x;i < numTileLists;i+=blockDim.x/WARPSIZE*gridDim.x) 
+  for (int i = threadIdx.x/WARPSIZE + blockDim.x/WARPSIZE*blockIdx.x;i < numTileLists;i+=blockDim.x/WARPSIZE*gridDim.x)
   {
     int j = tileListOrder[i];
     int start = tileListPos[i];
@@ -669,7 +669,7 @@ __global__ void reOrderTileListDepth(const int numTileLists, const int* __restri
     tileListDepthDst[i] = tileListDepthSrc[j];
   }
 
-} 
+}
 
 //
 // Bit shift tileListDepth so that only lower 16 bits are used
@@ -874,7 +874,7 @@ void CudaTileListKernel::updateComputes(const int numComputesIn,
 
 void CudaTileListKernel::writeTileList(const char* filename, const int numTileLists,
   const TileList* d_tileLists, cudaStream_t stream) {
-  
+
   TileList* h_tileLists = new TileList[numTileLists];
   copy_DtoH<TileList>(d_tileLists, h_tileLists, numTileLists, stream);
   cudaCheck(cudaStreamSynchronize(stream));
@@ -981,7 +981,8 @@ void CudaTileListKernel::markJtileOverlap(const int width, const int numTileList
 void CudaTileListKernel::buildTileLists(const int numTileListsPrev,
   const int numPatchesIn, const int atomStorageSizeIn, const int maxTileListLenIn,
   const float3 lata, const float3 latb, const float3 latc,
-  const CudaPatchRecord* h_cudaPatches, const float4* h_xyzq, const float plcutoff2In,
+  const CudaPatchRecord* h_cudaPatches, const float4* h_xyzq,
+  const float plcutoff2In, const size_t maxShmemPerBlock,
   cudaStream_t stream) {
 
   numPatches = numPatchesIn;
@@ -1059,6 +1060,9 @@ void CudaTileListKernel::buildTileLists(const int numTileListsPrev,
     int nblock = min(deviceCUDA->getMaxNumBlocks(), (numTileListsPrev-1)/nthread+1);
 
     int shmem_size = buildTileListsBBKernel_shmem_sizePerThread(maxTileListLen)*nthread;
+    if(shmem_size > maxShmemPerBlock){
+      NAMD_die("CudaTileListKernel::buildTileLists, maximum shared memory allocation exceeded. Too many atoms in a patch");
+    }
 
     // NOTE: In the first call numJtiles = 1. buildTileListsBBKernel will return and
     //       tell the required size in h_tileListStat->numJtiles. In subsequent calls,
@@ -1172,13 +1176,13 @@ void CudaTileListKernel::sortTileLists(
   // Check that the array sizes are adequate
   if (numTileListsSrc > tileListsSrc.size || numJtilesSrc > tileJatomStartSrc.size ||
     numTileListsSrc > tileListDepthSrc.size || numTileListsSrc > tileListOrderSrc.size ||
-    (patchPairsSrc.ptr != NULL && numTileListsSrc > patchPairsSrc.size) || 
+    (patchPairsSrc.ptr != NULL && numTileListsSrc > patchPairsSrc.size) ||
     (tileExclsSrc.ptr != NULL && numJtilesSrc > tileExclsSrc.size))
     NAMD_die("CudaTileListKernel::sortTileLists, Src allocated too small");
 
   if (numTileListsDst > tileListsDst.size || numJtilesDst > tileJatomStartDst.size ||
     numTileListsSrc > tileListDepthDst.size || numTileListsSrc > tileListOrderDst.size ||
-    (patchPairsDst.ptr != NULL && numTileListsDst > patchPairsDst.size) || 
+    (patchPairsDst.ptr != NULL && numTileListsDst > patchPairsDst.size) ||
     (tileExclsDst.ptr != NULL && numJtilesDst > tileExclsDst.size))
     NAMD_die("CudaTileListKernel::sortTileLists, Dst allocated too small");
 
@@ -1301,7 +1305,7 @@ void CudaTileListKernel::sortTileLists(
 
         int nthread = 1024;
         int nblock = min(deviceCUDA->getMaxNumBlocks(), (numTileListsSrc-1)/nthread+1);
-        buildRemoveZerosSortKey <<< nblock, nthread, 0, stream >>> 
+        buildRemoveZerosSortKey <<< nblock, nthread, 0, stream >>>
         (numTileListsSrc, tileListDepthSrc.ptr, begin_bit, sortKeySrc);
         cudaCheck(cudaGetLastError());
       }
@@ -1338,7 +1342,7 @@ void CudaTileListKernel::sortTileLists(
       {
         int nthread = 1024;
         int nblock = min(deviceCUDA->getMaxNumBlocks(), (numTileListsDst-1)/nthread+1);
-        reOrderTileListDepth <<< nblock, nthread, 0, stream >>> 
+        reOrderTileListDepth <<< nblock, nthread, 0, stream >>>
         (numTileListsDst, tileListOrderDst.ptr,
           tileListDepthSrc.ptr, tileListDepthDst.ptr);
         cudaCheck(cudaGetLastError());
@@ -1357,7 +1361,7 @@ void CudaTileListKernel::sortTileLists(
 
         int nthread = 1024;
         int nblock = min(deviceCUDA->getMaxNumBlocks(), (numTileListsSrc-1)/nthread+1);
-        setupSortKey <<< nblock, nthread, 0, stream >>> 
+        setupSortKey <<< nblock, nthread, 0, stream >>>
         (numTileListsSrc, maxTileListLen_sortKeys, tileListsSrc.ptr, tileListDepthSrc.ptr, begin_bit, sortKeys, sortKeySrc);
         cudaCheck(cudaGetLastError());
 
@@ -1397,7 +1401,7 @@ void CudaTileListKernel::sortTileLists(
       {
         int nthread = 1024;
         int nblock = min(deviceCUDA->getMaxNumBlocks(), (numTileListsDst-1)/nthread+1);
-        reOrderTileListDepth <<< nblock, nthread, 0, stream >>> 
+        reOrderTileListDepth <<< nblock, nthread, 0, stream >>>
         (numTileListsDst, tileListOrderDst.ptr,
           tileListDepthSrc.ptr, tileListDepthDst.ptr);
         cudaCheck(cudaGetLastError());
@@ -1407,7 +1411,7 @@ void CudaTileListKernel::sortTileLists(
       {
         int nthread = 32;
         int nblock = min(deviceCUDA->getMaxNumBlocks(), (numTileListsDst-1)/nthread+1);
-        localSort<32> <<< nblock, nthread, 0, stream >>> 
+        localSort<32> <<< nblock, nthread, 0, stream >>>
         (numTileListsDst, begin_bit, num_bit, tileListDepthDst.ptr, tileListOrderDst.ptr);
         cudaCheck(cudaGetLastError());
 
@@ -1521,7 +1525,7 @@ void CudaTileListKernel::sortTileLists(
     // Fill in patchNumLists[0...numPatches-1]
     int nthread = 512;
     int nblock = min(deviceCUDA->getMaxNumBlocks(), (numTileListsDst-1)/nthread+1);
-    calcPatchNumLists <<< nblock, nthread, 0, stream >>> 
+    calcPatchNumLists <<< nblock, nthread, 0, stream >>>
     (numTileListsDst, numPatches, tileListsDst.ptr, patchNumLists);
     cudaCheck(cudaGetLastError());
 
@@ -1530,7 +1534,7 @@ void CudaTileListKernel::sortTileLists(
 
     // Fill in tileListsDst[0...numTileListsDst-1].patchNumLists
     // and find empty patches into emptyPatches[0 ... numEmptyPatches - 1]
-    setPatchNumLists_findEmptyPatches <<< nblock, nthread, 0, stream >>> 
+    setPatchNumLists_findEmptyPatches <<< nblock, nthread, 0, stream >>>
     (numTileListsDst, tileListsDst.ptr, patchNumLists,
       numPatches, &emptyPatches[numPatches], emptyPatches);
     cudaCheck(cudaGetLastError());
index 9191592..5b8d9e9 100644 (file)
@@ -13,7 +13,7 @@ struct TileList {
   int jtileEnd;
   float3 offsetXYZ;
   int2 patchInd;        // Patch indices for this list
-  union {  
+  union {
     int2 patchNumList;    // Number of lists contributing to each patch
     // int icompute;
   };
@@ -144,7 +144,7 @@ private:
   int sortKeyDstSize;
 
   int maxTileListLen_sortKeys;
-  
+
   unsigned int* sortKeys;
   int sortKeysSize;
 
@@ -310,7 +310,8 @@ public:
   void buildTileLists(const int numTileListsPrev,
     const int numPatchesIn, const int atomStorageSizeIn, const int maxTileListLenIn,
     const float3 lata, const float3 latb, const float3 latc,
-    const CudaPatchRecord* h_cudaPatches, const float4* h_xyzq, const float plcutoff2In, cudaStream_t stream);
+    const CudaPatchRecord* h_cudaPatches, const float4* h_xyzq, const float plcutoff2In,
+    const size_t maxShmemPerBlock, cudaStream_t stream);
 
   void reSortTileLists(const bool doGBIS, cudaStream_t stream);
   // void applyOutputOrder(cudaStream_t stream);