title: Cuda Dot Product Failing for Non Multiples of 1024
body: <p>I'm just looking for some help here when it comes to calculating the dot product of two arrays.</p>
  
  <p>Let's say I set the array size to 2500 and the max thread count per block to 1024.</p>
  
  <p>In essence, I want to calculate the dot product of each block, and then sum the dot products in another kernel function. I calculate the number of blocks as such:</p>
  
  <pre><code>nblcks = (n + 1024 -1)/1024
  </code></pre>
  
  <p>So, <code>nblcks = 3</code></p>
  
  <p>This is my kernel function:</p>
  
  <pre><code>// calculate the dot product block by block
  __global__ void dotProduct(const float* a, const float* b, float* c, int n){
      // store the product of a[i] and b[i] in shared memory
      // sum the products in shared memory
      // store the sum in c[blockIdx.x]
  
      __shared__ float s[ntpb];
      int tIdx = threadIdx.x;
      int i = blockDim.x * blockIdx.x + threadIdx.x;
  
      //calc product
      if (i &lt; n)
          s[tIdx] = a[i] * b[i];
      __syncthreads();
  
      for (int stride = 1; stride &lt; blockDim.x; stride &lt;&lt;= 1) {
           if (tIdx % (2 * stride) == 0)
               s[tIdx] += s[tIdx + stride];
           __syncthreads();
       }
  
      if (threadIdx.x == 0){
          c[blockIdx.x] = s[0];
      }
  
  }
  </code></pre>
  
  <p>I call the kernel:</p>
  
  <pre><code>dotProduct&lt;&lt;&lt;nblocks, ntpb&gt;&gt;&gt;(d_a, d_b, d_c, n);
  </code></pre>
  
  <p>And everything works! Well, almost.</p>
  
  <p>d_c, which has 3 elements - each one the dot product of the block is thrown off on the last element.</p>
  
  <pre><code>d_c[0] = correct
  d_c[1] = correct
  d_c[2] = some massive number of 10^18
  </code></pre>
  
  <p>Can someone point out why this is occurring? It only seems to work for multiples of 1024. So... 2048, 3072, etc... Am I iterating over null values or stack overflowing? </p>
  
  <p>Thank you!</p>
  
  <p>Edit:</p>
  
  <pre><code> // host vectors
      float* h_a = new float[n];
      float* h_b = new float[n];
      init(h_a, n);
      init(h_b, n);
      // device vectors (d_a, d_b, d_c)
      float* d_a;
      float* d_b;
      float* d_c;
      cudaMalloc((void**)&amp;d_a, n * sizeof(float));
      cudaMalloc((void**)&amp;d_b, n * sizeof(float));
      cudaMalloc((void**)&amp;d_c, nblocks * sizeof(float));
  
      // copy from host to device h_a -&gt; d_a, h_b -&gt; d_b
      cudaMemcpy(d_a, h_a, n * sizeof(float), cudaMemcpyHostToDevice);
      cudaMemcpy(d_b, h_b, n * sizeof(float), cudaMemcpyHostToDevice);
  </code></pre>
  
  <p>Initialization of the array's are done in this function (n times):</p>
  
  <pre><code>void init(float* a, int n) {
      float f = 1.0f / RAND_MAX;
      for (int i = 0; i &lt; n; i++)
          a[i] = std::rand() * f; // [0.0f 1.0f]
  }
  </code></pre>