Hello,

I am trying to implement a matrix vector kernel with reduction tree as below code snippet. This implementation is based on spmv_csr_vector_kernel explained in this paper. [http://mgarland.org/files/papers/nvr-2008-004.pdf]

But it does NOT work without adding lot of syncthreads after each if condition of the reduction tree.

If I remove the syncthreads after each if condition in the kernel, the CPU and GPU results don’t match and the test reports failure.

Based on my understanding these synchthreads are not necessary. So wanted some help to figure out if this implementation has any bug.

Thank you,

Ram

```
#include <stdio.h>
#define THREADS_PER_ROW 32
#define BLOCK_SIZE 256
#define NUM_ROWS 256
#define NUM_COLS 256
#define cudaErrorCheck(err) \
if(err != cudaSuccess) { \
printf("##### error in cuda call at %s: %d. %s: %s #####\n", __FILE__, __LINE__, cudaGetErrorName(err), cudaGetErrorString(err)); \
exit(1); \
}
__global__ void matrix_vector_kernel ( const int num_rows , int num_cols, const double * data , const double * x , double * y)
{
__shared__ double vals [BLOCK_SIZE];
int thread_id = blockDim.x * blockIdx.x + threadIdx.x ; // global thread index
int warp_id = thread_id / THREADS_PER_ROW; // global warp index
int lane = thread_id & (THREADS_PER_ROW-1); // thread index within the warp
int row = warp_id;
if(row < num_rows){
vals [ threadIdx.x ] = 0;
for ( int jj = lane ; jj < num_cols ; jj += THREADS_PER_ROW)
vals [ threadIdx.x ] += data [ (row*num_cols)+jj ] * x [jj];
__syncthreads();
// parallel reduction in shared memory
if (lane < 16) vals[threadIdx.x] += vals[threadIdx.x + 16];
__syncthreads();
if (lane < 8) vals[threadIdx.x] += vals[threadIdx.x + 8];
__syncthreads();
if (lane < 4) vals[threadIdx.x] += vals[threadIdx.x + 4];
__syncthreads();
if (lane < 2) vals[threadIdx.x] += vals[threadIdx.x + 2];
__syncthreads();
if (lane < 1) vals[threadIdx.x] += vals[threadIdx.x + 1];
__syncthreads();
// first thread OF EACH WARP ACCUMULATES the result
if ( lane == 0)
y[row] = vals [ threadIdx.x ];
__syncthreads();
}
}
int main(){
double M[NUM_ROWS*NUM_COLS], N[NUM_COLS], P[NUM_ROWS];
//init
for(int i=0; i<NUM_ROWS; i++){
for(int j=0; j<NUM_COLS; j++)
M[i*NUM_COLS + j] = i*NUM_COLS + j;
P[i] = 0;
}
for(int j=0; j<NUM_COLS; j++)
N[j] = j;
//GPU calc
double *m, *n, *p;
cudaErrorCheck(cudaMalloc (( void **)&m , NUM_ROWS*NUM_COLS*sizeof(double)));
cudaErrorCheck(cudaMalloc (( void **)&n , NUM_COLS*sizeof(double)));
cudaErrorCheck(cudaMalloc (( void **)&p , NUM_ROWS*sizeof(double)));
cudaErrorCheck(cudaMemcpy(m, M, NUM_ROWS*NUM_COLS*sizeof(double), cudaMemcpyHostToDevice));
cudaErrorCheck(cudaMemcpy(n, N, NUM_COLS*sizeof(double), cudaMemcpyHostToDevice));
cudaErrorCheck(cudaMemcpy(p, P, NUM_ROWS*sizeof(double), cudaMemcpyHostToDevice));
int rows_per_block = BLOCK_SIZE/THREADS_PER_ROW;
int num_of_blocks = (NUM_ROWS + rows_per_block - 1) / rows_per_block;
matrix_vector_kernel<<<num_of_blocks, BLOCK_SIZE>>>(NUM_ROWS, NUM_COLS, m, n, p);
cudaErrorCheck(cudaMemcpy(P, p, NUM_ROWS*sizeof(double), cudaMemcpyDeviceToHost));
cudaErrorCheck(cudaDeviceSynchronize());
//cpu check
for(int i=0; i<NUM_ROWS; i++){
double temp = 0;
for(int j=0; j<NUM_COLS; j++)
temp += M[i*NUM_COLS + j] * N[j];
if(temp != P[i]){
printf("error for element %d. CPU: %f, GPU: %f\n", i, temp, P[i]);
return 0;
}
}
printf("success!\n");
return 0;
}
```