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);
}
```