2017-08-07 39 views
0

我有个问题。我正在尝试学习OpenCl,所以我一直试图用OpenCl来实现FFT算法。我试图重新创建:OpenCl算法执行结果不同

void FFT (cmplx* data, int dataSize){ 
    if(dataSize == 1){ 
     return; 
    }else{ 
     cmplx* even = (cmplx*)malloc(dataSize/2*sizeof(cmplx)); 
     cmplx* odd = (cmplx*)malloc(dataSize/2*sizeof(cmplx)); 
     for (int i = 0;i<dataSize;i+=2){ 
      even[i/2] = data[i]; 
      odd[i/2] = data[i+1]; 
     } 

     FFT(even,dataSize/2); 
     FFT(odd,dataSize/2); 

     for (int i = 0;i<dataSize;i++){ 
      cmplx C = cmplx(-2*M_PI/dataSize*i); 
      data[i].real = even[i].real + C.real*odd[i].real - C.imag*odd[i].imag; 
      data[i].imag = even[i].imag + C.real*odd[i].imag + C.imag*odd[i].real; 
     } 
    } 
} 

CMPLX仅仅是持有复数两个浮体实部和虚部并具有与欧拉方程建立复数一个构造函数的类。和其他一切是相当简单

我可能不知道一些小的细微差别,在我的理解,我可以做的周期计算在独立线程,以便循环是这样的:

for (int i = 0;i<dataSize;i++){ 
    cmplx C = cmplx(-2*M_PI/dataSize*i); 
    data[i].real = even[i].real + C.real*odd[i].real - C.imag*odd[i].imag; 
    data[i].imag = even[i].imag + C.real*odd[i].imag + C.imag*odd[i].real; 
} 

随着OpenCL的这样的代码:

__kernel void FFTComplexSum(__global float *evenReal,__global float *evenImag, 
         __global float *oddReal,__global float *oddImag, 
         __global float *real,__global float *imag, 
         __global float *C){ 
    int gid = get_global_id(0); 
    real[gid] = evenReal[gid] + cos(C[gid])*oddReal[gid] - sin(C[gid])*oddImag[gid]; 
    imag[gid] = evenImag[gid] + cos(C[gid])*oddImag[gid] + sin(C[gid])*oddReal[gid]; 
} 

但如果运行这个命令:

.... // instantiating stuff like platform, device_id, kernel, program... 
    size_t buffer_size; 
    cl_mem evenReal_mem, evenImag_mem, oddReal_mem, oddImag_mem, real_mem, imag_mem, c_mem; 

    float evenReal[dataSize]; 
    float evenImag[dataSize]; 
    float tReal[dataSize]; 
    float tImag[dataSize]; 
    float oddReal[dataSize]; 
    float oddImag[dataSize]; 
    float C[dataSize]; 

    for (int i = 0;i<dataSize;i+=2){ 
     evenReal[i/2] = real[i]; 
     evenImag[i/2] = imag[i]; 
     oddReal[i/2] = real[i+1]; 
     oddImag[i/2] = imag[i+1]; 
     C[i] = -2*M_PI/dataSize*i; 
     C[i+1] = -2*M_PI/dataSize*(i+1); 
    } 

    doubleArray(evenReal,dataSize); // doubleArray function just makes array to loop 
    doubleArray(evenImag,dataSize); 

    doubleArray(oddReal,dataSize); 
    doubleArray(oddImag,dataSize); 

    buffer_size = sizeof(float) * dataSize; 

    evenReal_mem = clCreateBuffer(context, CL_MEM_READ_ONLY, buffer_size, NULL, NULL); 
    err = clEnqueueWriteBuffer(cmd_queue, evenReal_mem, CL_TRUE, 0, buffer_size,(void*)evenReal, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); // Fail check 

    evenImag_mem = clCreateBuffer(context, CL_MEM_READ_ONLY, buffer_size, NULL, NULL); 
    err = clEnqueueWriteBuffer(cmd_queue, evenImag_mem, CL_TRUE, 0, buffer_size,(void*)evenImag, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); // Fail check 

    oddReal_mem = clCreateBuffer(context, CL_MEM_READ_ONLY, buffer_size, NULL, NULL); 
    err = clEnqueueWriteBuffer(cmd_queue, oddReal_mem, CL_TRUE, 0, buffer_size,(void*)oddReal, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); // Fail check 

    oddImag_mem = clCreateBuffer(context, CL_MEM_READ_ONLY, buffer_size, NULL, NULL); 
    err = clEnqueueWriteBuffer(cmd_queue, oddImag_mem, CL_TRUE, 0, buffer_size,(void*)oddImag, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); // Fail check 

    real_mem = clCreateBuffer(context, CL_MEM_WRITE_ONLY, buffer_size, NULL, NULL); 
    err = clEnqueueWriteBuffer(cmd_queue, real_mem, CL_TRUE, 0, buffer_size,(void*)real, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); // Fail check 

    imag_mem = clCreateBuffer(context, CL_MEM_WRITE_ONLY, buffer_size, NULL, NULL); 
    err = clEnqueueWriteBuffer(cmd_queue, imag_mem, CL_TRUE, 0, buffer_size,(void*)imag, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); // Fail check 

    c_mem = clCreateBuffer(context, CL_MEM_READ_ONLY, sizeof(float), NULL, NULL); 
    err = clEnqueueWriteBuffer(cmd_queue, c_mem, CL_TRUE, 0, sizeof(float),(void*)C, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); // Fail check 

    clFinish(cmd_queue); 

    err = clSetKernelArg(kernel[0], 0, sizeof(cl_mem), &evenReal_mem); 
    err = clSetKernelArg(kernel[0], 1, sizeof(cl_mem), &evenImag_mem); 
    err = clSetKernelArg(kernel[0], 2, sizeof(cl_mem), &oddReal_mem); 
    err = clSetKernelArg(kernel[0], 3, sizeof(cl_mem), &oddImag_mem); 
    err = clSetKernelArg(kernel[0], 4, sizeof(cl_mem), &real_mem); 
    err = clSetKernelArg(kernel[0], 5, sizeof(cl_mem), &imag_mem); 
    err = clSetKernelArg(kernel[0], 6, sizeof(cl_mem), &c_mem); 
    assert(err == CL_SUCCESS); // Fail check 

    size_t global_work_size = dataSize; 
    err = clEnqueueNDRangeKernel(cmd_queue, kernel[0], 1, NULL, &global_work_size, NULL, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); 
    clFinish(cmd_queue); 

    printf("test data:\n"); 

    for (int i = 0;i<dataSize;i++){ 
     float r,I; 
     r = evenReal[i] + cos(C[i])*oddReal[i] - sin(C[i])*oddImag[i]; 
     I = evenImag[i] + cos(C[i])*oddImag[i] + sin(C[i])*oddReal[i]; 
     printf("%f + %f\n",r,I); 
    } 

    err = clEnqueueReadBuffer(cmd_queue, real_mem, CL_TRUE, 0, buffer_size, tReal, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); 
    clFinish(cmd_queue); 
    err = clEnqueueReadBuffer(cmd_queue, imag_mem, CL_TRUE, 0, buffer_size, tImag, 0, NULL, NULL); 
    assert(err == CL_SUCCESS); 
    clFinish(cmd_queue); 

    clReleaseMemObject(evenReal_mem); 
    clReleaseMemObject(evenImag_mem); 
    clReleaseMemObject(oddReal_mem); 
    clReleaseMemObject(oddImag_mem); 
    clReleaseMemObject(real_mem); 
    clReleaseMemObject(imag_mem); 
    clReleaseMemObject(c_mem); 

    prinf("data:"); 

    for (int i = 0;i<dataSize;i++){ 
     printf("%f + %f\n",tReal[i],tImag[i]); 
    } 

它[R eturns这:

test data: 
1.000000 + 0.000000 
0.000000 + 1.000000 
-1.000000 + 0.000000 
-0.000000 + -1.000000 
data: 
1.000000 + 0.000000 
-1.000000 + 0.000000 
1.000000 + 0.000000 
-1.000000 + 0.000000 

我真的很困惑,为什么我得到错误的答案。我错过了一些非常明显的东西?

对不起,长期以来的问题。

回答

0

c_mem有问题。

您在内核中访问C,如C[gid],但是仅在sizeof(float)的范围内创建了它。所以主数据不能适应那个(4字节)的存储空间,而你只写了4个字节。将其创建大小和写入大小与data_size相乘应该足够了。

这就是为什么实数值得到1和-1而虚数零(sin(0))。如果你有幸,溢出的C会给垃圾并给出真实和虚构的元素垃圾结果,这会立即显示错误的根源。

+0

Omg ......我觉得不好意思看不到这个...谢谢。 –