title: Sum not being stored correctly
body: <p>This is my second try at a CUDA program and parallel development.</p>
  
  <p>I'm trying to find why the dot product two arrays.</p>
  
  <p>Example:</p>
  
  <p>Given  </p>
  
  <pre><code>A = [1 2 3] 
  B = [4 5 6]
  
  C[0] = (1)(4) + (1)(5) + (1)(6)  
  C[1] = (2)(4) + (2)(5) + (2)(6)  
  C[2] = (3)(4) + (3)(5) + (3)(6)  
  </code></pre>
  
  <p>I intialize 2 arrays, A and B , randomly fill the elements within a range, then I multiply each element in <code>A</code> with each in <code>B</code> and store the sum of the products in a third array identified as <code>C</code>. I've set the array's size of A, B and C to 100.</p>
  
  <p>That gives me 10'000 multiplications which are the ones that I parallelise by using 100 blocks and 128 threads (due to warp size).</p>
  
  <p>Here is my kernel function:</p>
  
  <pre><code>__global__ void kernel(float *a, float *b, float *c, const int N) {
      if( threadIdx.x &lt; N ) 
          c[blockIdx.x] += a[blockIdx.x] * b[threadIdx.x];
  }
  </code></pre>
  
  <p>This is my reasoning, since the aggregation has to be accumulated in <code>C</code> which has the same index as the pivot index in <code>A</code> so I can reuse the <code>blockidx.x</code> and that 'should' work correctly; but it doesn't.</p>
  
  <p>My suspicion is that the <code>C</code> index gets cleared or is not shared when the thread changes but I'm really not sure which is why I seek advice.</p>
  
  <p>Here's is the complete code, I've explicitly avoided the <code>HANDLE_ERROR</code> function wrappers for shortness</p>
  
  <pre><code>#include &lt;stdio.h&gt;
  #include &lt;cuda.h&gt;
  #include &lt;time.h&gt;
  
  #define M 100
  
  __global__ void kernel(float *a, float *b, float *c, const int N) {
      if(threadIdx.x &lt; N) 
          c[blockIdx.x] += a[blockIdx.x] * b[threadIdx.x];
  }
  
  void init_array(float*, const int);
  void fill_array(float*, const int, const float); 
  void print_array(float*, const int, *char);
  
  int main (void) {
      srand( time(NULL) );
  
      float a[M], b[M], c[M] = { 0.0 };
      float *dev_a, *dev_b, *dev_c;
      const int S = sizeof(float) * M;
  
      init_array(a, M);
      init_array(b, M);
  
      print_array(a, M, "a");
      print_array(b, M, "b");
      print_array(c, M, "c");
  
      cudaMalloc((void**)&amp;dev_a, S);
      cudaMalloc((void**)&amp;dev_b, S);
      cudaMalloc((void**)&amp;dev_c, S);
  
      cudaMemcpy(dev_a, a, S, cudaMemcpyHostToDevice);
      cudaMemcpy(dev_b, b, S, cudaMemcpyHostToDevice);
      cudaMemcpy(dev_c, c, S, cudaMemcpyHostToDevice);
  
      kernel&lt;&lt;&lt;M, M + 28&gt;&gt;&gt;(dev_a, dev_b, dev_c, M);
  
      cudaMemcpy(c, dev_c, S, cudaMemcpyDeviceToHost);
  
      cudaFree(dev_a);
      cudaFree(dev_b);
      cudaFree(dev_c);
  
      print_array(c, M, "c");
  
      return 0;
  }
  
  void init_array(float *a, const int N) {
     int i;  
     for(i=0; i&lt;N; i++)
         a[i] = rand() % M + 1;
  }
  
  void fill_array(float *a, const int N, const float v) {
     int i;  
     for(i=0; i&lt;N; i++)
         a[i] = v;
  }
  
  void print_array(float *a, const int N, char *d) {
     int i;  
     for(i=0; i&lt;N; i++)
         printf("\n%s[%d]: %f",d, i, a[i]);
  }
  </code></pre>