CUDA 奇怪的行为访问向量

Question

我在 cuda 中实现了一个简单的 fft 程序。这是核函数：

__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;
        } */
    }
}

现在，如果我使用此代码，我会得到预期的输出，即

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

其中长度 = 4。 但是如果我使用注释掉的部分，我会得到错误的结果

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

这很奇怪，因为两个版本之间的变化只是对向量的访问模式。 在第一种情况下，我访问每个线程的两个元素，它们的长度/2 分开。 在第二种情况下，我访问两个相邻的元素。 这对我来说毫无意义......这可能是硬件问题吗？

编辑： 使用 cuda-memcheck 运行表示没有错误，并且...错误不会出现！

编辑：完整（非）工作副本：

#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

我很确定这是不正确的：

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

第一个内核配置参数应该是网格尺寸。第二个应该是块尺寸。

我能够重现您在发布的代码中指出的错误。当我反转这两个参数时：

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

错误对我来说消失了。

FWIW，synccheck 工具确实会报告屏障错误，即使使用上述“修复”也是如此。 synccheck 是 CUDA 7 中的新功能，requires a cc3.5 or higher device。

Answer 2

我认为你可能需要另一个__syncthreads()，在这里：

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