reviewed uses of atomicAdd(),;

temporarily disallowed Scatter() reductions

Author: frankseide

Source/Math/Convolution.cuh

@@ -269,4 +269,4 @@ __global__ void kAveragePoolingBackward(int batchSize, const int* mpRowCol, cons
     }
 }
 
-} } }
+}}}

Source/Math/GPUMatrix.cu

@@ -885,11 +885,11 @@ __global__ void _doGatherColumnsOf(ElemType* us, size_t usStride, const ElemType
     CUDA_LONG jOut = id / usStride; // col index into 'us' and 'idx'
 
     auto jInF = idx[jOut * idxStride]; // this is the column we need to get
-    if (jInF < 0)                      // negative index means gap
+    if (isnan(jInF) || jInF < 0)       // negative index means gap
         return;
     size_t jIn = (size_t)jInF;
-    if (jIn >= aCols)
-        return; // actually a failure
+    //if (jIn >= aCols)
+    //    return; // actually a failure
 
     const ElemType&  ra = a[    i + jIn  *  aStride  ];
     ElemType&       rus = us[id/*i + jOut * usStride*/];
@@ -928,6 +928,21 @@ GPUMatrix<ElemType>& GPUMatrix<ElemType>::DoGatherColumnsOf(ElemType beta, const
     return *this;
 }
 
+// little helper for debugging
+template <class ElemType>
+static void Peek(const GPUMatrix<ElemType>& m, const char* which)
+{
+    size_t rows = m.GetNumRows();
+    size_t cols = m.GetNumCols();
+    ElemType buf[10000] = { 0 };
+    size_t n = min(rows * cols, _countof(buf));
+    CUDA_CALL(cudaMemcpy(buf, m.Data(), sizeof(ElemType) * n, cudaMemcpyDeviceToHost));
+    UNUSED(which); UNUSED(rows); UNUSED(cols); sin(1.0f); // set breakpoint here
+    //CUDA_CALL(cudaMemcpy(const_cast<ElemType*>(m.Data()), buf, sizeof(ElemType) * n, cudaMemcpyHostToDevice));
+}
+
+#undef ALLOW_ATOMIC_SCATTER // allow to disable this, until we know atomicAdd() works properly here
+
 template <class ElemType>
 __global__ void _doScatterColumnsOf(ElemType* us, size_t usStride, size_t usCols, const ElemType* idx, size_t idxStride, const ElemType* a, size_t aStride, const ElemType alpha, CUDA_LONG numElements)
 {
@@ -941,34 +956,25 @@ __global__ void _doScatterColumnsOf(ElemType* us, size_t usStride, size_t usCols
     CUDA_LONG jIn = id / aStride; // col index into 'a' and 'idx'
 
     auto jOutF = idx[jIn * idxStride]; // this is the column we copy/add into
-    if (jOutF < 0)                     // negative index means gap
+    if (isnan(jOutF) || jOutF < 0)     // negative index means gap
         return;
     size_t jOut = (size_t)jOutF;
-    if (jOut >= usCols)
-        return; // actually a failure  --TODO: This should not be necessary. Why is it?
+    //if (jOut >= usCols)
+    //    return; // actually a failure  --TODO: This should not be necessary. Why is it?
 
     const ElemType&  ra =  a[id/*i + jIn  *  aStride*/];
     ElemType&       rus = us[    i + jOut * usStride  ];
 
     ElemType res = ra * alpha;
+#ifdef ALLOW_ATOMIC_SCATTER
     if (res != 0)             // avoid memory conflict if e.g. an entire column has no gradient
         atomicAdd(&rus, res); // rus += res;
+#else
+    rus += res;
+#endif
     // Note: atomicAdd() is supposed to be fast in case of no conflict (the simple case of Scatter())
 }
 
-// little helper for debugging
-template <class ElemType>
-static void Peek(const GPUMatrix<ElemType>& m, const char* which)
-{
-    size_t rows = m.GetNumRows();
-    size_t cols = m.GetNumCols();
-    ElemType buf[10000] = { 0 };
-    size_t n = min(rows * cols, _countof(buf));
-    CUDA_CALL(cudaMemcpy(buf, m.Data(), sizeof(ElemType) * n, cudaMemcpyDeviceToHost));
-    UNUSED(which); UNUSED(rows); UNUSED(cols); sin(1.0f); // set breakpoint here
-    //CUDA_CALL(cudaMemcpy(const_cast<ElemType*>(m.Data()), buf, sizeof(ElemType) * n, cudaMemcpyHostToDevice));
-}
-
 // *this[:,idx[j]] = a[:,j] * alpha + *this[:,idx[j]] * beta
 template <class ElemType>
 GPUMatrix<ElemType>& GPUMatrix<ElemType>::DoScatterColumnsOf(ElemType beta, const GPUMatrix<ElemType>& idx, const GPUMatrix<ElemType>& a, ElemType alpha)
@@ -986,6 +992,27 @@ GPUMatrix<ElemType>& GPUMatrix<ElemType>::DoScatterColumnsOf(ElemType beta, cons
 
     auto& us = *this;
 
+#ifndef ALLOW_ATOMIC_SCATTER // verify that atomicAdd is not needed  --this is not efficient
+    {
+        vector<ElemType> buf(idx.GetNumRows() * idx.GetNumCols()); // idx(,)are the column(s) we copy/add into
+        CUDA_CALL(cudaMemcpy(buf.data(), idx.Data(), sizeof(ElemType) * buf.size(), cudaMemcpyDeviceToHost));
+        vector<bool> writtenTo(GetNumCols(), false); // remember whether an output column is in fact a target
+        for (size_t i = 0; i < buf.size(); i++)
+        {
+            auto colF = buf[i];
+            if (isnan(colF) || colF < 0)
+                continue;
+            size_t col = (size_t)colF;
+            if (col >= GetNumCols())
+                LogicError("DoScatterColumnsOf: Index value out of bounds.");
+            if (writtenTo[col])
+                LogicError("DoScatterColumnsOf: #ifndef ALLOW_ATOMIC_SCATTER then columns must be unique. Column idx(%d,%d)=%d is used twice.", (int)(i % idx.GetNumCols()), (int)(i / idx.GetNumCols()), (int)col);
+            else
+                writtenTo[col] = true;
+        }
+    }
+#endif
+
     // pre-scale with beta upfront
     // Scatter may add more than one source column to the same target, so we must pre-scale with beta, and then just keep adding.
     Scale(beta, us); // if beta is 0, then this will be a memset()

Source/Math/GPUMatrixCUDAKernels.cuh

@@ -43,6 +43,7 @@
 #define IDX2C(i, j, ld) (((j) * (ld)) + (i)) // 0 based indexing
 
 // CUDA atomicAdd() only exists for 'float'. This is the 'double' version.
+// TODO: This may need to be guarded by CUDA version; newer devices may support this.
 static __inline__ __device__ double atomicAdd(double* address, double val)
 {
     unsigned long long int* address_as_ull = (unsigned long long int*) address;
@@ -3152,7 +3153,8 @@ __global__ void _scaleSparseBlockAndAddToDense(
     rhs[IDX2C(row, col, numRows)] += alpha * lhsValues[index];
 }
 
-// compute predictions in cross entory node
+#if 0
+// compute predictions in cross entropy node
 template <class ElemType>
 __global__ void _computePrediction(
     int nv,
@@ -3335,6 +3337,7 @@ __global__ void _computeGradientOfInput(
 
     atomicAdd(&grd[IDX2C(h, j, numrows)], sum);
 }
+#endif
 
 template <class ElemType>
 __global__ void computeNCEForwardProp(
@@ -3713,6 +3716,8 @@ __global__ void _assignNceDerivativeNew(
             atomicAdd(&c[wid], -er);
     }
 }
+
+#if 0
 // compute gradients of weights in cross entropy node
 template <class ElemType>
 __global__ void _computeGradientOfWeight(
@@ -3774,6 +3779,7 @@ __global__ void _computeGradientOfWeight(
         blockIds[ii] = i;
     }
 }
+#endif
 
 // used in clipping gradients
 template <class ElemType>

Source/Math/GPUTensor.cu

@@ -393,6 +393,8 @@ struct TensorOpElement<ElemType, N, M, K, /*parallelReduce=*/false, /*k=*/-1>
     }
 };
 
+#define ALLOW_ATOMIC_REDUCTION // undefine to disable use of atomicAdd() below, for testing it
+
 // specialization for k = -1 terminates the template recursion, and computes reductions in parallel
 template <class ElemType, C_size_t N, C_int M, C_int K>
 struct TensorOpElement<ElemType, N, M, K, /*parallelReduce=*/true, /*k=*/-1>
@@ -403,8 +405,8 @@ struct TensorOpElement<ElemType, N, M, K, /*parallelReduce=*/true, /*k=*/-1>
                                    const FixedArray<C_unsigned_int, K>& /*regularOpStrides*/, const FixedMatrix<C_int, N, K>& /*regularStrides*/,
                                    const FixedArray<C_unsigned_int, M>& reducingOpDims, const FixedMatrix<C_int, N, M>& reducingStrides, CUDA_LONG reductionBegin, CUDA_LONG reductionChunkSize)
     {
-        CUDA_LONG reductionBlock = blockIdx.z; // block index  --larger reductions are split into blocks
-        CUDA_LONG reductionBlocks = gridDim.z; // number of blocks
+        CUDA_LONG reductionBlock = blockIdx.z; // reduction-block index  --larger reductions are split into blocks
+        CUDA_LONG reductionBlocks = gridDim.z; // number of reduction blocks. If >1 we need atomicAdd
         CUDA_LONG tid = threadIdx.x;           // thread index
         CUDA_LONG tids = blockDim.x;           // out of how many threads  --note: last block is partial
 
@@ -427,7 +429,7 @@ struct TensorOpElement<ElemType, N, M, K, /*parallelReduce=*/true, /*k=*/-1>
         }
 
         // reduce    --cf https://docs.nvidia.com/cuda/samples/6_Advanced/reduction/doc/reduction.pdf
-        __shared__ ReduceElemType accumulators[GridDim::maxThreadsPerBlock /*tids*/];
+        __shared__ ReduceElemType volatile accumulators[GridDim::maxThreadsPerBlock /*tids*/];
         accumulators[tid] = sum;
         __syncthreads();
         static_assert(GridDim::maxThreadsPerBlock <= 512, "GridDim::maxThreadsPerBlock too large, need to add manually unrolled steps");
@@ -450,8 +452,12 @@ struct TensorOpElement<ElemType, N, M, K, /*parallelReduce=*/true, /*k=*/-1>
             auto* pout = pointers[pointers.size() - 1];
             if (reductionBlocks > 1) // multiple blocks: need to use atomicAdd()
             {
+#ifdef ALLOW_ATOMIC_REDUCTION
                 // in this case, outer calling code must pass beta = 1
                 atomicAdd(pout, val);
+#else
+                *pout = 1000000.0f;  // something that can't be missed? How to crash it?
+#endif
             }
             else
             {
@@ -560,8 +566,9 @@ static void LaunchTensorOpWithReduction(ElemType beta, array<ElemType*, N> point
     C_size_t reductionDim = 1; // number of elements to reduce over
     for (C_size_t k = 0; k < reducingOpDimVector.size(); k++)
         reductionDim *= (C_size_t) reducingOpDimVector[k];
-    let& props = GridDim::GetDeviceProps();
     GridDim grid(NN);
+#ifdef ALLOW_ATOMIC_REDUCTION // temporarily disabled to ensure it is not causing the non-reproducability
+    let& props = GridDim::GetDeviceProps();
     if (reductionDim > 1 && grid.m_blocksPerGrid < props.multiProcessorCount) // TODO: <= multiProcessorCount?
     {
         // we are reducing and are underutilizing the multiprocs we have: get more parallelism by doing reduction in parallel
@@ -603,6 +610,7 @@ static void LaunchTensorOpWithReduction(ElemType beta, array<ElemType*, N> point
         }
     }
     else
+#endif
     {
         // we got enough elements to generate: do one element per thread, and reduction inside
         _launchTensorOp<ElemType, N, M, K><<<grid.m_blocksPerGrid, grid.m_threadsPerBlock, 0, t_stream>>>(beta, pointers, alpha, op, regularOpStrides, regularStrides, grid.m_N, reducingOpDims, reducingStrides);

Source/Math/latticefunctionskernels.h

@@ -16,11 +16,7 @@
 #include "latticestorage.h"
 #include <limits>
 
-namespace msra { namespace cuda {
-
-class passtextureref;
-}
-}
+namespace msra { namespace cuda { class passtextureref; } }
 
 #ifdef CPUONLY
 #define __kernel_emulation__
@@ -34,7 +30,8 @@ using namespace std;
 #define __device__
 #endif
 #define CUDART_MIN_DENORM_F numeric_limits<float>::denorm_min()
-#define atomicAdd(address, value) (*(address) += (value)) // don't forget to #undef (#praga pop_macro)! Otherwise CUDA might compile with this...
+// renamed to x- so we make sure to not accidentally use these; rename back if ever needed again
+#define xatomicAdd(address, value) (*(address) += (value)) // don't forget to #undef (#praga pop_macro)! Otherwise CUDA might compile with this...
 #define atomicCAS(address, compare, val) \
     *address;                            \
     *address = *address == compare ? val : *address;
@@ -47,8 +44,8 @@ using namespace std;
 #if __CUDA_ARCH__ < 200
 //#warning Sequence training not supported on 1.x CUDA machines.
 #define force_crash() (*((int *) -1) = 0)         // TODO: this does not in fact seem to crash it...
-#define atomicAdd(a, v) (force_crash(), *(a) = v) // force a crash if used with 1.x devices
-#define atomicCAS(address, compare, val) (*(address) = compare + val, *((int *) -1) = 0)
+#define xatomicAdd(a, v) (force_crash(), *(a) = v) // force a crash if used with 1.x devices
+#define xatomicCAS(address, compare, val) (*(address) = compare + val, *((int *) -1) = 0)
 #define __double_as_longlong(in) (force_crash(), in)
 #define __longlong_as_double(in) (force_crash(), in)
 #define __float_as_int(in) (force_crash(), in)
@@ -956,8 +953,8 @@ struct latticefunctionskernels
         }
     }
 };
-};
-};
+
+}};
 
 #pragma pop_macro("atomicCAS")
 #pragma pop_macro("atomicAdd")