title: cuda shared memory - inconsistent results
body: <p>I'm trying to do a parallel reduction to sum an array in CUDA. Currently i pass an array in which to store the sum of the elements in each block. This is my code:</p>
  
  <pre><code>#include &lt;cstdlib&gt;
  #include &lt;iostream&gt;
  #include &lt;cuda.h&gt;
  #include &lt;cuda_runtime_api.h&gt;
  #include &lt;helper_cuda.h&gt;
  #include &lt;host_config.h&gt;
  #define THREADS_PER_BLOCK 256
  #define CUDA_ERROR_CHECK(ans) { gpuAssert((ans), __FILE__, __LINE__); }
  
  using namespace std;
  
  inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
  {
     if (code != cudaSuccess) 
     {
        fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
        if (abort) exit(code);
     }
  }
  
  struct double3c {
      double x; 
      double y;
      double z;
      __host__ __device__ double3c() : x(0), y(0), z(0) {}
      __host__ __device__ double3c(int x_, int y_, int z_) : x(x_), y(y_), z(z_) {}
      __host__ __device__ double3c&amp; operator+=(const double3c&amp; rhs) { x += rhs.x; y += rhs.y; z += rhs.z;}
      __host__ __device__ double3c&amp; operator/=(const double&amp; rhs) { x /= rhs; y /= rhs; z /= rhs;}
  
  };
  
  class VectorField {
  public:
      double3c *data;
      int size_x, size_y, size_z;
      bool is_copy;  
  
      __host__ VectorField () {}
  
      __host__ VectorField (int x, int y, int z) {
          size_x = x; size_y = y; size_z = z;
          is_copy = false;
          CUDA_ERROR_CHECK (cudaMalloc(&amp;data, x * y * z * sizeof(double3c))); 
      }
  
      __host__ VectorField (const VectorField&amp; other) {
          size_x = other.size_x; size_y = other.size_y; size_z = other.size_z;
          this-&gt;data = other.data;
          is_copy = true;
      }
  
      __host__ ~VectorField() {     
          if (!is_copy) CUDA_ERROR_CHECK (cudaFree(data));
      }
  };
  
  __global__ void KernelCalculateMeanFieldBlock (VectorField m, double3c* result) {
      __shared__ double3c blockmean[THREADS_PER_BLOCK];    
      int index = threadIdx.x + blockIdx.x * blockDim.x;
      if (index &lt; m.size_x * m.size_y * m.size_z) blockmean[threadIdx.x] = m.data[index] = double3c(0, 1, 0);
      else blockmean[threadIdx.x] = double3c(0,0,0);
      __syncthreads();
      for(int s = THREADS_PER_BLOCK / 2; s &gt; 0; s /= 2) {
          if (threadIdx.x &lt; s) blockmean[threadIdx.x] += blockmean[threadIdx.x + s];
          __syncthreads();
      }
  
  
      if(threadIdx.x == 0) result[blockIdx.x] = blockmean[0];   
  }
  
  double3c CalculateMeanField (VectorField&amp; m) { 
      int blocknum = (m.size_x * m.size_y * m.size_z - 1) / THREADS_PER_BLOCK + 1;
      double3c *mean = new double3c[blocknum]();
      double3c *cu_mean;
      CUDA_ERROR_CHECK (cudaMalloc(&amp;cu_mean, sizeof(double3c) * blocknum));
      CUDA_ERROR_CHECK (cudaMemset (cu_mean, 0, sizeof(double3c) * blocknum));
  
          KernelCalculateMeanFieldBlock &lt;&lt;&lt;blocknum, THREADS_PER_BLOCK&gt;&gt;&gt; (m, cu_mean);
          CUDA_ERROR_CHECK (cudaPeekAtLastError());
          CUDA_ERROR_CHECK (cudaDeviceSynchronize());
          CUDA_ERROR_CHECK (cudaMemcpy(mean, cu_mean, sizeof(double3c) * blocknum, cudaMemcpyDeviceToHost));
  
      CUDA_ERROR_CHECK (cudaFree(cu_mean));
      for (int i = 1; i &lt; blocknum; i++) {mean[0] += mean[i];}
      mean[0] /= m.size_x * m.size_y * m.size_z;
      double3c aux = mean[0];
      delete[] mean;
      return aux;
  }
  
  
  
  int main() {
      VectorField m(100,100,100);
      double3c sum = CalculateMeanField (m);
      cout &lt;&lt;  sum.x &lt;&lt; '\t' &lt;&lt; sum.y &lt;&lt; '\t' &lt;&lt;sum.z;  
  
  
      return 0;
  }
  </code></pre>
  
  <p><strong>EDIT</strong> </p>
  
  <p>Posted a functional code. Constructing a <code>VectorField</code> with 10x10x10 elements works fine and gives mean 1, but constructing it with 100x100x100 elements gives mean ~0.97 (it varies from run to run). Is this a right way to do a parallel reduction, or should I stick to one kernel launch per block?</p>