CUDA Dot Product

Question

I am trying to implement the classic dot-product core for double precision arrays with atomic calculation of the final sum over various blocks. I used atomAdd for double precision, as indicated on page 116 of the programming guide. Maybe I'm doing something wrong. The partial sums of the flows in each block are calculated correctly, but the afterword atomic operation does not seem to work properly since every time I start my kernel with the same data, I get different results. I would be grateful if anyone could find a mistake or provide an alternative solution! Here is my core:

__global__ void cuda_dot_kernel(int *n,double *a, double *b, double *dot_res) { __shared__ double cache[threadsPerBlock]; //thread shared memory int global_tid=threadIdx.x + blockIdx.x * blockDim.x; int i=0,cacheIndex=0; double temp = 0; cacheIndex = threadIdx.x; while (global_tid < (*n)) { temp += a[global_tid] * b[global_tid]; global_tid += blockDim.x * gridDim.x; } cache[cacheIndex] = temp; __syncthreads(); for (i=blockDim.x/2; i>0; i>>=1) { if (threadIdx.x < i) { cache[threadIdx.x] += cache[threadIdx.x + i]; } __syncthreads(); } __syncthreads(); if (cacheIndex==0) { *dot_res=cuda_atomicAdd(dot_res,cache[0]); } } 

And here is my device atomicAdd function:

 __device__ double cuda_atomicAdd(double *address, double val) { double assumed,old=*address; do { assumed=old; old= __longlong_as_double(atomicCAS((unsigned long long int*)address, __double_as_longlong(assumed), __double_as_longlong(val+assumed))); }while (assumed!=old); return old; } 

Answer 1

You are using the cuda_atomicAdd function cuda_atomicAdd . This section of your kernel:

 if (cacheIndex==0) { *dot_res=cuda_atomicAdd(dot_res,cache[0]); } 

is the culprit. Here you add atomically to dot_res . then non-atomically sets dot_res with the return of the result. The result of returning from this function is the previous location value, which was updated using the atom, and it was intended for "information" or local use only by the caller. You do not assign it to the fact that you update atomically, which completely strikes the goal of using access to atomic memory in the first place. Do something like this:

 if (cacheIndex==0) { double result=cuda_atomicAdd(dot_res,cache[0]); } 

Answer 2

Getting the right shorthand using the special CUDA code can be difficult, so here is an alternative solution that uses the Thrust algorithm, which is included in the CUDA Toolkit:

 #include <thrust/inner_product.h> #include <thrust/device_ptr.h> double do_dot_product(int n, double *a, double *b) { // wrap raw pointers to device memory with device_ptr thrust::device_ptr<double> d_a(a), d_b(b); // inner_product implements a mathematical dot product return thrust::inner_product(d_a, d_a + n, d_b, 0.0); } 

Answer 3

Did not check your code in depth, but here are some tips.
I would advise using Thrust if you only use your GPU for such general tasks, because if there is a difficult problem, people have no idea to effectively program parallel on gpu.

• Start a new parallel reduction core to summarize the point product.
  Since the data is already on the device, you will not see a decrease in performance, starting with the new kernel.

• The kernel does not seem to scale according to the maximum number of possible blocks on the latest GPU. If that were the case, and your kernel could calculate the point product from millions of values, productivity would drop dramatically due to serialized atomic operation.

• Entry-level error: have your input and shared memory access range been verified? Or are you sure that the input is always a multiple of your block size? Otherwise, you will read the trash. Most of my incorrect results were caused by this error.

• optimize parallel contraction. My Thesis or Optimization Mark Harris

Unconfirmed, I just wrote it in a notebook:

 /* * @param inCount_s unsigned long long int Length of both input arrays * @param inValues1_g double* First value array * @param inValues2_g double* Second value array * @param outDots_g double* Output dots of each block, length equals the number of blocks */ __global__ void dotProduct(const unsigned long long int inCount_s, const double* inValuesA_g, const double* inValuesB_g, double* outDots_g) { //get unique block index in a possible 3D Grid const unsigned long long int blockId = blockIdx.x //1D + blockIdx.y * gridDim.x //2D + gridDim.x * gridDim.y * blockIdx.z; //3D //block dimension uses only x-coordinate const unsigned long long int tId = blockId * blockDim.x + threadIdx.x; /* * shared value pair products array, where BLOCK_SIZE power of 2 * * To improve performance increase its size by multiple of BLOCK_SIZE, so that each threads loads more then 1 element! * (outDots_g length decreases by same factor, and you need to range check and initialize memory) * -> see harris gpu optimisations / parallel reduction slides for more informations. */ __shared__ double dots_s[BLOCK_SIZE]; /* * initialize shared memory array and calculate dot product of two values, * shared memory always needs to be initialized, its never 0 by default, else garbage is read later! */ if(tId < inCount_s) dots_s[threadIdx.x] = inValuesA_g[tId] * inValuesB_g[tId]; else dots_s[threadIdx.x] = 0; __syncthreads(); //do parallel reduction on shared memory array to sum up values reductionAdd(dots_s, dots_s[0]) //see my thesis link //output value if(threadIdx.x == 0) outDots_g[0] = dots_s[0]; //start new parallel reduction kernel to sum up outDots_g! } 

Edit: Removed unnecessary points.