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 idX(blockDim.x*blockIdx.x+threadIdx.x+1);
unsigned int idY(blockDim.y*blockIdx.y+threadIdx.y+1);
unsigned int idZ(blockDim.z*blockIdx.z+threadIdx.z+1);

unsigned int sX(sizeY*sizeZ);
unsigned int i((idX*sizeY+idY)*sizeZ+idZ);
unsigned int bi((threadIdx.x*blockDim.y+threadIdx.y)*blockDim.z+threadIdx.z);

s_pAp[bi] = Ap[i]*p[i];

// Wait for all threads

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){
    if(threadIdx.x < offset && threadIdx.y < offset && threadIdx.z < offset){

    // Wait for all threads

// Write final result with (0,0,0) thread
if(threadIdx.x==0 && threadIdx.y==0 && threadIdx.z == 0){
    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?

  • 1
    \$\begingroup\$ Ummm... shouldn't the __syncthreads() in the for-loop, be outside the for-loop? \$\endgroup\$
    – rolfl
    Commented Jun 6, 2014 at 15:34
  • \$\begingroup\$ 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. \$\endgroup\$
    – repptilia
    Commented Jun 6, 2014 at 15:44
  • \$\begingroup\$ 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) \$\endgroup\$ Commented Sep 10, 2014 at 18:04

1 Answer 1

  • 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.

  • More whitespace within statements can be added for more readability.

    This line, for instance, is hard to read:


    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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.