# 3D vector CUDA kernel

I designed this CUDA kernel to compute a function on a 3D domain:

p and Ap are 3D vectors that are actually implemented as a single long array:

__global__ void update(P_REAL* data, P_REAL* tmp, P_REAL* f, P_REAL* reduced,
const unsigned int sizeX,const unsigned int sizeY,const unsigned int sizeZ,
const P_REAL aX, const P_REAL aY, const P_REAL aZ, const P_REAL aF,
const P_REAL cX, const P_REAL cY, const P_REAL cZ){

extern __shared__ P_REAL s_pAp[];

// Get grid size
unsigned int gY((sizeY-2)/blockDim.y),gZ((sizeZ-2)/blockDim.z);

unsigned int sX(sizeY*sizeZ);
unsigned int i((idX*sizeY+idY)*sizeZ+idZ);

Ap[i]=cX*(p[i-sX]+p[i+sX]-2.0f*p[i])
+cY*(p[i-sizeZ]+p[i+sizeZ]-2.0f*p[i])
+cZ*(p[i-1]+p[i+1]-2.0f*p[i]);
s_pAp[bi] = Ap[i]*p[i];



Reduction of s_pAp (to be summed):

// Take advantage of the geometry of the block
// At each iteration, there are 7 points to consider
// Use byte-shift to ease division by 2

for(int offset = blockDim.x/2;offset>0;offset>>=1){
s_pAp[bi]+=s_pAp[bi+offset]+s_pAp[bi+offset*blockDim.z]+s_pAp[bi+offset*blockDim.z*blockDim.y]
+s_pAp[bi+offset*(1+blockDim.z)]
+s_pAp[bi+offset*(1+blockDim.z*blockDim.y)]
+s_pAp[bi+offset*blockDim.z*(1+blockDim.y)]
+s_pAp[bi+offset*(1+blockDim.z*(1+blockDim.y))];
}

}

// Write final result with (0,0,0) thread
reduced[(blockIdx.x*gY+blockIdx.y)*gZ+blockIdx.z] = s_pAp[0];

}
}


This works like a charm, but when profiling I get awful performance. Could I get some insight on how to improve this code?

• Ummm... shouldn't the __syncthreads() in the for-loop, be outside the for-loop? – rolfl Jun 6 '14 at 15:34
• This synchronization is needed to prevent wrong data overwriting. Basically, the loop halves each dimension of the cube at each iteration, we need to do this synchronously. – repptilia Jun 6 '14 at 15:44
• I have some ideas, but I need more information before I could make concrete suggestions. For starters, what is the original algorithm you are implementing? (Hopefully this will answer why the reduction must consider 7 points) – t. fochtman Sep 10 '14 at 18:04

• The indentation within update() is misleading. Most of all, everything inside of the function should be indented. Otherwise, it'll be hard to tell what exactly is inside of the function.

I do, however, like that you've indented the (lengthy) parameter list by eight spaces instead of the standard four. It should be easier to tell them apart from the code inside of the function.

This line, for instance, is hard to read:

Ap[i]=cX*(p[i-sX]+p[i+sX]-2.0f*p[i])


You can add more whitespace as such:

Ap[i] = cX * (p[i-sX] + p[i+sX] - 2.0f * p[i])


Now you can really tell what everything is at a glance.

You can also add a space before the opening curly braces for functions. It may not be a big deal with functions with smaller signatures, but yours are pretty lengthy, and it's harder to see the opening curly brace with them.

You also have inconsistent whitespace, such as in here:

for(int offset = blockDim.x/2;offset>0;offset>>=1){


There is spacing between the =, but nowhere else. This should look clearer:

for (int offset = blockDim.x / 2; offset > 0; offset >>= 1) {

• Considering that CUDA's blockIdx and blockDim are primarily signed, you should make your thread ID variables of the same type. If they won't be too large, make them int. If they can be larger than an int, make them either long or long long and cast blockIdx to the same.