CUDA strange behavior accessing vector

Question

I have implemented a simple fft program in cuda. This is the kernel function:

__global__
void fftKernel(cuComplex* dev_samples,
              size_t length,
              size_t llog,
              Direction direction)
{
    int tid = threadIdx.x + blockDim.x * blockIdx.x;
    if (tid < length / 2) {

        // First step, sorts data with bit reversing and compute new coefficients;
        cuComplex c1 = dev_samples[bit_reverse(tid * 2) >> (32 - llog)];
        cuComplex c2 = dev_samples[bit_reverse(tid * 2 + 1) >> (32 - llog)];

        dev_samples[tid * 2] = cuCaddf(c1, c2);
        dev_samples[tid * 2 + 1] = cuCsubf(c1, c2);

        __syncthreads();

        int k = tid * 2;
        int block_start = tid * 2;

        // Butterfly propagation
        for (int i = 1, n = 4; i < llog; i++, n *= 2) {

            int new_block_start = (k/n) * n;
            if (new_block_start != block_start) {
                block_start = new_block_start;
                k -= n/4;
            }

            // We compute
            // X(k) = E(k) + TF(k, N) * O(k)
            // and
            // X(k + N/2) = E(k) - TF(k, N) * O(k)
            c1 = dev_samples[k];
            c2 = cuCmulf(twiddle_factor(k % n, n, direction), dev_samples[k + n/2]);

            dev_samples[k] = cuCaddf(c1, c2);
            dev_samples[k + n / 2] = cuCsubf(c1, c2);

            __syncthreads();
        }

        // Normalize if backward transforming
        if (direction == Backward) {
            dev_samples[tid].x /= length;
            dev_samples[tid].y /= length;
            dev_samples[tid + length/2].x /= length;
            dev_samples[tid + length/2].y /= length;
        }


        // Strange behavior if using this snipset 
        // instead of the previous one
        /*
        if (direction == Backward) {
            dev_samples[tid * 2].x /= length;
            dev_samples[tid * 2].y /= length;
            dev_samples[tid * 2 + 1].x /= length;
            dev_samples[tid * 2 + 1].y /= length;
        } */
    }
}

Now if I use this code I get the expected output, that is

 (1.5, 0)
 (-0.5, -0.5)
 (-0.5, 0)
 (-0.5, 0.5)

where length = 4. But if i use the commented out section I get the wrong result

 (1.5, 0)
 (-0.5, -0.5)
 (1, 0)  <--- wrong
 (-0.5, 0.5)

This is odd because what changes between the two versions is only the access pattern to the vector. In the first case I access two elements per thread that are length/2 apart. In the second case I access two adjacent elements. This makes no sense to me... Can it be a hw issue?

EDIT: running with cuda-memcheck signals no errors and... the error does not maifest!

EDIT: full (non)working copy:

#include <iostream>
#include <stdio.h>
#include <cuda.h>
#include <cuComplex.h>
#include <math_constants.h>

static void cuCheck( cudaError_t err,
                     const char *file,
                     int line ) {
    if (err != cudaSuccess) {
    printf( "%s in %s at line %d\n", cudaGetErrorString( err ),
            file, line );
    exit( EXIT_FAILURE );
    }
}
#define CU_CHECK( err ) (cuCheck( err, __FILE__, __LINE__ ))

#define BLOCK_SIZE 512

enum Direction {
    Forward,
    Backward
};

/* Computes the in place fft of dev_samples.
 * Retrun true in the fft has successfully computed, false otherwise.*/
bool fft(cuComplex* dev_samples, size_t length, Direction direction);

/* Returns true if n is a power of 2 and return the logarithm
 * of the greater number smaller than n that is a power of 2
 * in the variable pointed by exp */
__device__ __host__
bool isPowerOf2(int n, int* exp = NULL);

/* Computes the twiddle factor */
__device__
cuComplex twiddle_factor(int k, size_t n, Direction direction);
__device__
unsigned int bit_reverse(unsigned int i);

__global__
void fftKernel(cuComplex* dev_samples,
           size_t length,
           size_t llog,
           Direction direction);

__device__ __host__
bool isPowerOf2(int n, int* exp)
{
    int tmp_exp = 0;
    int i;
    for (i = 1; i < n; i *= 2, tmp_exp++);
    if (exp) {
    *exp = tmp_exp;
    }
    return (i == n);
}

bool fft(cuComplex* dev_samples, size_t length, Direction direction)
{
    int llog;
    if (!isPowerOf2(length, &llog))
    return false;

    int gridSize = max((int)length/2/BLOCK_SIZE, (int)1);
    fftKernel<<<BLOCK_SIZE, gridSize>>>(dev_samples,
                                    length,
                                    llog,
                                    direction);
    CU_CHECK(cudaDeviceSynchronize());


    return true;
}

__global__
void fftKernel(cuComplex* dev_samples,
           size_t length,
           size_t llog,
           Direction direction)
{
    int tid = threadIdx.x + blockDim.x * blockIdx.x;
    if (tid < length / 2) {

    // First step, sorts data with bit reversing and compute new coefficients;
    cuComplex c1 = dev_samples[bit_reverse(tid * 2) >> (32 - llog)];
    cuComplex c2 = dev_samples[bit_reverse(tid * 2 + 1) >> (32 - llog)];

    __syncthreads();

    dev_samples[tid * 2] = cuCaddf(c1, c2);
    dev_samples[tid * 2 + 1] = cuCsubf(c1, c2);

    __syncthreads();

    int k = tid * 2;
    int block_start = tid * 2;

    // Butterfly propagation
    for (int i = 1, n = 4; i < llog; i++, n *= 2) {

        int new_block_start = (k/n) * n;
        if (new_block_start != block_start) {
            block_start = new_block_start;
            k -= n/4;
        }

        // We compute
        // X(k) = E(k) + TF(k, N) * O(k)
        // and
        // X(k + N/2) = E(k) - TF(k, N) * O(k)
        c1 = dev_samples[k];
        c2 = cuCmulf(twiddle_factor(k % n, n, direction), dev_samples[k + n/2]);

        dev_samples[k] = cuCaddf(c1, c2);
        dev_samples[k + n / 2] = cuCsubf(c1, c2);

        __syncthreads();
     }

    //--------------------- ERROR HERE -----------------------
    // Normalize if backward transforming
    if (direction == Backward) {
        dev_samples[tid * 2].x /= length;
        dev_samples[tid * 2].y /= length;
        dev_samples[tid * 2+ 1].x /= length;
        dev_samples[tid * 2+ 1].y /= length;
    }

    // This works!!
    /*if (direction == Backward) {
        dev_samples[tid * 2].x /= length;
        dev_samples[tid * 2].y /= length;
        dev_samples[tid * 2+ 1].x /= length;
        dev_samples[tid * 2+ 1].y /= length;
    }*/
       //---------------------------------------------------------
    }
}

__device__
cuComplex twiddle_factor(int k, size_t n, Direction direction)
{
    if (direction == Forward)
    return make_cuFloatComplex(cos(-2.0*CUDART_PI_F*k/n),
                               sin(-2.0*CUDART_PI_F*k/n));
    else
    return make_cuFloatComplex(cos(2.0*CUDART_PI_F*k/n),
                               sin(2.0*CUDART_PI_F*k/n));
}

__device__
unsigned int bit_reverse(unsigned int i) {
  register unsigned int mask = 0x55555555; // 0101...
  i = ((i & mask) << 1) | ((i >> 1) & mask);
  mask = 0x33333333; // 0011...
  i = ((i & mask) << 2) | ((i >> 2) & mask);
  mask = 0x0f0f0f0f; // 00001111...
  i = ((i & mask) << 4) | ((i >> 4) & mask);
  mask = 0x00ff00ff; // 0000000011111111...
  i = ((i & mask) << 8) | ((i >> 8) & mask);
  // 00000000000000001111111111111111 no need for mask
  i = (i << 16) | (i >> 16);
  return i;
}

#define SIGNAL_LENGTH 4

using namespace std;

int main(int argc, char** argv)
{
    size_t mem_size = SIGNAL_LENGTH*sizeof(cuComplex);
    cuComplex* signal = (cuComplex*)malloc(mem_size);
    cuComplex* dev_signal;


    for (int i = 0; i < SIGNAL_LENGTH; i++) {
    signal[i].x = i;
    signal[i].y = 0;
    }

    CU_CHECK(cudaMalloc((void**)&dev_signal, mem_size));
    CU_CHECK(cudaMemcpy(dev_signal, signal, mem_size, cudaMemcpyHostToDevice));

    fft(dev_signal, SIGNAL_LENGTH, Backward);

    CU_CHECK(cudaMemcpy(signal, dev_signal, mem_size, cudaMemcpyDeviceToHost));

    cout << "polito_Z_out:" << endl;
    for (int i = 0; i < SIGNAL_LENGTH; i++) {
    cout << " (" << signal[i].x << ", " << signal[i].y << ")" << endl;
    }

    cudaFree(dev_signal);
    free(signal);
}

Answer 1

I'm pretty sure this is not correct:

fftKernel<<<BLOCK_SIZE, gridSize>>>(...);

The first kernel configuration parameter should be the grid dimensions. The second should be the block dimensions.

I was able to reproduce the error you indicate in the code you have posted. When I reverse those two parameters:

fftKernel<<<gridSize, BLOCK_SIZE>>>(...);

the error goes away for me.

FWIW, the synccheck tool does report a barrier error, even with the above "fix". synccheck is new in CUDA 7 and also requires a cc3.5 or higher device.

Answer 2

I think you may need another __syncthreads(), here:

    cuComplex c1 = dev_samples[bit_reverse(tid * 2) >> (32 - llog)];
    cuComplex c2 = dev_samples[bit_reverse(tid * 2 + 1) >> (32 - llog)];

    __syncthreads(); // <<< need this, otherwise possible race condition ?

    dev_samples[tid * 2] = cuCaddf(c1, c2);
    dev_samples[tid * 2 + 1] = cuCsubf(c1, c2);