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 A to reduce_fabs(), everything from A[4] 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.
This the kernel that calls the sum-reduction kernel:
__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;
}
This is the sum-reduction kernel:
__device__ float reduce_fabs(float *v, int begin_index, int end_index) {
    extern __shared__ float sum_array[];
    int tid = threadIdx.x;
    if(tid >= begin_index && tid <= 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<=(end_index-begin_index);j*=2) {
        if((tie-begin_index)%(2*j) == 0 && tie >= begin_index && tid <= end_index) {
            sum_array[(tie-begin_index)] += sum_array[(tie-begin_index)+j];
        }
        __syncthreads();
    }
    return sum_array[0];
}
Having seen the code, a more specific description of the problem is that when I pass the array A to reduce_fabs(), the value of element 4 becomes A[4] = 0 which is wrong. A[0],A[1],A[2], and A[3] are fine for any end_index, but every end_index larger than 3 results in the elements of A beyond 3 being set to zero upon passage to the reduction kernel.
Here is what I tried already:
- I tried diag[0] = A[4]to make sure the initial assignment was working. It was.
- I eliminated the summation portion of reduce_fabs(), stopping after the first__syncthreads()and the problem still persisted.
- I eliminated the zeroing of elements beyond those of interest for the sum; i.e. I commented //sum_array[tid+end_index-begin_index+1] = 0(the summation was commented out for this part too). No success.
- I eliminated if(tie >= begin_index && tid <= end_index)so thatsum_arraywas assigned for all threads (summation also commented out) to no avail.
Here is a main function to run it:
#include<stdio.h>
#include<stdlib.h>
#include<cuda.h>
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**) &p_A,n*sizeof(float));
    cudaMalloc((void**) &p_diag,n*sizeof(float));
    cudaMalloc((void**) &p_offdiag,n*sizeof(float));
    cudaMalloc((void**) &p_n,sizeof(int));
    cudaMemcpy(p_n,&n,n*sizeof(int),cudaMemcpyHostToDevice);
    tridiag<<<1,n>>>(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;
}
