title: CUDA: sum-reduction --- data lost in call to device function
body: <p>I am writing a CUDA sum reduction code taking the sum of the absolute values of an array starting on element begin_index through end_index (I am using one block with a variable number of threads). However, when I pass the array <code>A</code> to <code>reduce_fabs()</code>, everything from <code>A[4]</code> and higher indexes end up somehow being set to zero. Here is the code and the function calls and more explanation of what I have tried following.</p>
  
  <p>This the kernel that calls the sum-reduction kernel:</p>
  
  <pre><code>__device__ void tridiag(float *A,int *dim, float *diag,float *offdiag) {
  
      A[0] = 1.0f; A[1] = 2.0f; A[2] = 3.0f;
      A[3] = 4.0f; A[4] = 5.0f; A[5] = 6.0f;
      diag[0] = reduce_fabs(A,0,3);
      __syncthreads();
      diag[1] = reduce_fabs(A,0,4);
  
      return;
  }
  </code></pre>
  
  <p>This is the sum-reduction kernel:</p>
  
  <pre><code>__device__ float reduce_fabs(float *v, int begin_index, int end_index) {
      extern __shared__ float sum_array[];
      int tid = threadIdx.x;
      if(tid &gt;= begin_index &amp;&amp; tid &lt;= end_index) {
          sum_array[tid-begin_index] = fabs(v[tid]);
          sum_array[tid+end_index-begin_index+1] = 0;
      }
      __syncthreads();
      for(int j=1;j&lt;=(end_index-begin_index);j*=2) {
          if((tie-begin_index)%(2*j) == 0 &amp;&amp; tie &gt;= begin_index &amp;&amp; tid &lt;= end_index) {
              sum_array[(tie-begin_index)] += sum_array[(tie-begin_index)+j];
          }
          __syncthreads();
      }
  
      return sum_array[0];
  }
  </code></pre>
  
  <p>Having seen the code, a more specific description of the problem is that when I pass the array <code>A</code> to <code>reduce_fabs()</code>, the value of element 4 becomes <code>A[4] = 0</code> which is wrong. <code>A[0]</code>,<code>A[1]</code>,<code>A[2]</code>, and <code>A[3]</code> are fine for any <code>end_index</code>, but every <code>end_index</code> larger than 3 results in the elements of <code>A</code> beyond 3 being set to zero upon passage to the reduction kernel.</p>
  
  <p>Here is what I tried already:</p>
  
  <ul>
  <li>I tried <code>diag[0] = A[4]</code> to make sure the initial assignment was working. It was.</li>
  <li>I eliminated the summation portion of <code>reduce_fabs()</code>, stopping after the first <code>__syncthreads()</code> and the problem still persisted.</li>
  <li>I eliminated the zeroing of elements beyond those of interest for the sum; i.e. I commented <code>//sum_array[tid+end_index-begin_index+1] = 0</code> (the summation was commented out for this part too). No success.</li>
  <li>I eliminated <code>if(tie &gt;= begin_index &amp;&amp; tid &lt;= end_index)</code> so that <code>sum_array</code> was assigned for all threads (summation also commented out) to no avail.</li>
  </ul>
  
  <p>Here is a main function to run it:</p>
  
  <pre><code>#include&lt;stdio.h&gt;
  #include&lt;stdlib.h&gt;
  #include&lt;cuda.h&gt;
  
  int main() {
      int n = 10;
      float *A = (float *)malloc(n*sizeof(*A));
      float *diag = (float *)malloc(n*sizeof(*A));
      float *offdiag = (float *)malloc(n*sizeof(*A));
  
      int *p_n;
      float *p_A, *p_diag, *p_offdiag;
  
      cudaMalloc((void**) &amp;p_A,n*sizeof(float));
      cudaMalloc((void**) &amp;p_diag,n*sizeof(float));
      cudaMalloc((void**) &amp;p_offdiag,n*sizeof(float));
      cudaMalloc((void**) &amp;p_n,sizeof(int));
  
      cudaMemcpy(p_n,&amp;n,n*sizeof(int),cudaMemcpyHostToDevice);
  
      tridiag&lt;&lt;&lt;1,n&gt;&gt;&gt;(p_A,p_n,p_diag,p_offdiag);
  
      cudaMemcpy(A,p_A,n*sizeof(float),cudaMemcpyDeviceToHost);
      cudaMemcpy(diag,p_diag,n*sizeof(float),cudaMemcpyDeviceToHost);
  
      printf("A[0] = %f  A[1] = %f  A[2] = %f\n",A[0],A[1],A[2]);
      printf("A[3] = %f  A[4] = %f  A[5] = %f\n",A[3],A[4],A[4]);
      printf("diag[0] = %f  diag[1] = %f\m",diag[0],diag[1]);
  
      cudaFree(p_A);
      cudaFree(p_diag);
      cudaFree(p_offdiag);
      free(A);
      free(diag);
      free(offdiag);
  
      return 0;
  }
  </code></pre>