cufftSetStream 导致垃圾输出。难道我做错了什么？答案

【问题标题】：cufftSetStream causes garbage output. Am I doing something wrong?cufftSetStream 导致垃圾输出。难道我做错了什么？
【发布时间】：2016-08-23 11:14:16
【问题描述】：

根据文档，cufftSetStream() 函数

将 CUDA 流与 cuFFT 计划相关联。现在，在计划执行期间进行的所有内核启动都通过关联的流完成 [...直到...] 流被另一个 cufftSetStream() 调用更改。

不幸的是，结果变成了垃圾。这是一个示例，通过两种方式执行一堆转换来演示这一点：一次是每个流都有自己的专用计划，一次是重复使用单个计划，如上面的文档所示。前者的行为符合预期，重用/cufftSetStream 方法在大多数转换中都有错误。这在我在 CentOS 7 linux 上尝试过的两张卡（GTX 750 ti，Titan X）上观察到 Cuda编译工具，发布7.0、V7.0.27；并发布 7.5、V7.5.17。

编辑：请参阅下面的“FIX”cmets，了解一种解决问题的方法。

#include <cufft.h>
#include <stdexcept>
#include <iostream>
#include <numeric>
#include <vector>

#define ck(cmd) if ( cmd) { std::cerr << "error at line " << __LINE__ << std::endl;exit(1);}


__global__
void fill_input(cufftComplex * buf, int batch,int nbins,int stride,int seed)
{
    for (int i = blockDim.y * blockIdx.y + threadIdx.y; i< batch;i += gridDim.y*blockDim.y)
        for (int j = blockDim.x * blockIdx.x + threadIdx.x; j< nbins;j += gridDim.x*blockDim.x)
            buf[i*stride + j] = make_cuFloatComplex( (i+seed)%101 - 50,(j+seed)%41-20);
}

__global__
void check_output(const float * buf1,const float * buf2,int batch, int nfft, int stride, int * errors)
{
    for (int i = blockDim.y * blockIdx.y + threadIdx.y; i< batch;i += gridDim.y*blockDim.y) {
        for (int j = blockDim.x * blockIdx.x + threadIdx.x; j< nfft;j += gridDim.x*blockDim.x) {
            float e=buf1[i*stride+j] - buf2[i*stride+j];
            if (e*e > 1) // gross error
                atomicAdd(errors,1);
        }
    }
}

void demo(bool reuse_plan)
{
    if (reuse_plan)
        std::cout << "Reusing the same fft plan with multiple stream via cufftSetStream ... ";
    else
        std::cout << "Giving each stream its own dedicated fft plan ... ";
    int nfft = 1024;
    int batch = 1024;
    int nstreams = 8;
    int nbins = nfft/2+1;
    int nit=100;
    size_t inpitch,outpitch;

    std::vector<cufftComplex*> inbufs(nstreams);
    std::vector<float*> outbufs(nstreams);
    std::vector<float*> checkbufs(nstreams);
    std::vector<cudaStream_t> streams(nstreams);
    std::vector<cufftHandle> plans(nstreams);
    for (int i=0;i<nstreams;++i) {
        ck( cudaStreamCreate(&streams[i]));
        ck( cudaMallocPitch((void**)&inbufs[i],&inpitch,nbins*sizeof(cufftComplex),batch) );
        ck( cudaMallocPitch((void**)&outbufs[i],&outpitch,nfft*sizeof(float),batch));
        ck( cudaMallocPitch((void**)&checkbufs[i],&outpitch,nfft*sizeof(float),batch) );
        if (i==0 || reuse_plan==false)
            ck ( cufftPlanMany(&plans[i],1,&nfft,&nbins,1,inpitch/sizeof(cufftComplex),&nfft,1,outpitch/sizeof(float),CUFFT_C2R,batch) );
    }

    // fill the input buffers and FFT them to get a baseline for comparison
    for (int i=0;i<nstreams;++i) {
        fill_input<<<20,dim3(32,32)>>>(inbufs[i],batch,nbins,inpitch/sizeof(cufftComplex),i);
        ck (cudaGetLastError());
        if (reuse_plan) {
            ck (cufftExecC2R(plans[0],inbufs[i],checkbufs[i]));
        }else{
            ck (cufftExecC2R(plans[i],inbufs[i],checkbufs[i]));
            ck( cufftSetStream(plans[i],streams[i]) ); // only need to set the stream once
        }
        ck( cudaDeviceSynchronize());
    }
    // allocate a buffer for the error count
    int * errors;
    cudaMallocHost((void**)&errors,sizeof(int)*nit);
    memset(errors,0,sizeof(int)*nit);

    /* FIX: an event can protect the plan internal buffers 
    by serializing access to the plan
    cudaEvent_t ev;
    cudaEventCreateWithFlags(&ev,cudaEventDisableTiming);
    */

    // perform the FFTs and check the outputs on streams
    for (int it=0;it<nit;++it) {
        int k = it % nstreams;
        ck( cudaStreamSynchronize(streams[k]) ); // make sure any prior kernels have completed
        if (reuse_plan) {
            // FIX: ck(cudaStreamWaitEvent(streams[k],ev,0 ) );
            ck(cufftSetStream(plans[0],streams[k]));
            ck(cufftExecC2R(plans[0],inbufs[k],outbufs[k]));
            // FIX: ck(cudaEventRecord(ev,streams[k] ) );
        }else{
            ck(cufftExecC2R(plans[k],inbufs[k],outbufs[k]));
        }
        check_output<<<100,dim3(32,32),0,streams[k]>>>(outbufs[k],checkbufs[k],batch,nfft,outpitch/sizeof(float),&errors[it]);
        ck (cudaGetLastError());
    }
    ck(cudaDeviceSynchronize());

    // report number of errors
    int errcount=0;
    for (int it=0;it<nit;++it)
        if (errors[it])
            ++errcount;
    std::cout << errcount << " of " << nit << " transforms had errors\n";

    for (int i=0;i<nstreams;++i) {
        cudaFree(inbufs[i]);
        cudaFree(outbufs[i]);
        cudaStreamDestroy(streams[i]);
        if (i==0 || reuse_plan==false)
            cufftDestroy(plans[i]);
    }
}

int main(int argc,char ** argv)
{
    demo(false);
    demo(true);
    return 0;
}

典型输出

为每个流提供自己专用的 fft 计划 ... 100 次转换中有 0 次出现错误
通过 cufftSetStream 重用具有多个流的相同 fft 计划 ... 100 次转换中有 87 次出现错误

【问题讨论】：

当我在一个非常普通的移动 GPU 上编译并运行您使用 CUDA 7.0 发布的代码时，两种情况下都出现 0 个错误。
@talonmies，感谢您提供的数据点。我刚刚用 cuda 7.0 尝试过——也失败了（见编辑）也许你卡的“谦虚”使它不会失败（即更少的资源 == 更少的竞争条件）。你在什么操作系统上？
带有计算 2.1 设备的 Windows 10
@RobertCrovella 我认为您关于追溯应用 SetStream 的理论是有道理的。可能某些内核在 exec 函数返回（错误的流）之后启动。虽然我同意文档的承诺有点含糊，但“直到 [] 再次调用 cufftSetStream”表明我正在尝试做的事情应该得到支持。
@MarkBorgerding 我不认为这是在您更改流后启动了一些内核（实际上，这只有在 cuFFT 有自己的线程时才会发生）。我假设 cuFFT 计划通常有一些相关的暂存空间。如果您尝试同时对两个流使用相同的计划，它们都会从相同的暂存空间读取和写入，从而导致值不正确。

标签： cuda cufft

【解决方案1】：

为了以您希望的方式重用计划，您需要手动管理 cuFFT 工作区。

每个计划都有一个中间计算结果的空间。如果您想同时将计划句柄用于两个或多个不同的计划执行，您需要为每个并发 cufftExec* 调用提供临时缓冲区。

您可以使用 cufftSetWorkArea 来执行此操作 - 请查看 cuFFT 文档中的第 3.7 节。第 2.2 节也有助于理解其工作原理。

这是一个工作示例，显示了为此对您的代码所做的更改：

$ cat t1241.cu
#include <cufft.h>
#include <stdexcept>
#include <iostream>
#include <numeric>
#include <vector>

#define ck(cmd) if ( cmd) { std::cerr << "error at line " << __LINE__ << std::endl;exit(1);}


__global__
void fill_input(cufftComplex * buf, int batch,int nbins,int stride,int seed)
{
    for (int i = blockDim.y * blockIdx.y + threadIdx.y; i< batch;i += gridDim.y*blockDim.y)
        for (int j = blockDim.x * blockIdx.x + threadIdx.x; j< nbins;j += gridDim.x*blockDim.x)
            buf[i*stride + j] = make_cuFloatComplex( (i+seed)%101 - 50,(j+seed)%41-20);
}

__global__
void check_output(const float * buf1,const float * buf2,int batch, int nfft, int stride, int * errors)
{
    for (int i = blockDim.y * blockIdx.y + threadIdx.y; i< batch;i += gridDim.y*blockDim.y) {
        for (int j = blockDim.x * blockIdx.x + threadIdx.x; j< nfft;j += gridDim.x*blockDim.x) {
            float e=buf1[i*stride+j] - buf2[i*stride+j];
            if (e*e > 1) // gross error
                atomicAdd(errors,1);
        }
    }
}

void demo(bool reuse_plan)
{
    if (reuse_plan)
        std::cout << "Reusing the same fft plan with multiple stream via cufftSetStream ... ";
    else
        std::cout << "Giving each stream its own dedicated fft plan ... ";
    int nfft = 1024;
    int batch = 1024;
    int nstreams = 8;
    int nbins = nfft/2+1;
    int nit=100;
    size_t inpitch,outpitch;

    std::vector<cufftComplex*> inbufs(nstreams);
    std::vector<float*> outbufs(nstreams);
    std::vector<float*> checkbufs(nstreams);
    std::vector<cudaStream_t> streams(nstreams);
    std::vector<cufftHandle> plans(nstreams);
    // if plan reuse, set up independent work areas
    std::vector<char *> wk_areas(nstreams);
    for (int i=0;i<nstreams;++i) {
        ck( cudaStreamCreate(&streams[i]));
        ck( cudaMallocPitch((void**)&inbufs[i],&inpitch,nbins*sizeof(cufftComplex),batch) );
        ck( cudaMallocPitch((void**)&outbufs[i],&outpitch,nfft*sizeof(float),batch));
        ck( cudaMallocPitch((void**)&checkbufs[i],&outpitch,nfft*sizeof(float),batch) );
        if (i==0 || reuse_plan==false)
            ck ( cufftPlanMany(&plans[i],1,&nfft,&nbins,1,inpitch/sizeof(cufftComplex),&nfft,1,outpitch/sizeof(float),CUFFT_C2R,batch) );
    }
    if (reuse_plan){
      size_t ws;
      ck(cufftGetSize(plans[0], &ws));
      for (int i = 0; i < nstreams; i++)
        ck(cudaMalloc(&(wk_areas[i]), ws));
      ck(cufftSetAutoAllocation(plans[0], 0));
      ck(cufftSetWorkArea(plans[0], wk_areas[0]));
      }
    // fill the input buffers and FFT them to get a baseline for comparison
    for (int i=0;i<nstreams;++i) {
        fill_input<<<20,dim3(32,32)>>>(inbufs[i],batch,nbins,inpitch/sizeof(cufftComplex),i);
        ck (cudaGetLastError());
        if (reuse_plan) {
            ck (cufftExecC2R(plans[0],inbufs[i],checkbufs[i]));
        }else{
            ck (cufftExecC2R(plans[i],inbufs[i],checkbufs[i]));
            ck( cufftSetStream(plans[i],streams[i]) ); // only need to set the stream once
        }
        ck( cudaDeviceSynchronize());
    }
    // allocate a buffer for the error count
    int * errors;
    cudaMallocHost((void**)&errors,sizeof(int)*nit);
    memset(errors,0,sizeof(int)*nit);

    // perform the FFTs and check the outputs on streams
    for (int it=0;it<nit;++it) {
        int k = it % nstreams;
        ck( cudaStreamSynchronize(streams[k]) ); // make sure any prior kernels have completed
        if (reuse_plan) {
            ck(cufftSetStream(plans[0],streams[k]));
            ck(cufftSetWorkArea(plans[0], wk_areas[k])); // update work area pointer in plan
            ck(cufftExecC2R(plans[0],inbufs[k],outbufs[k]));
        }else{
            ck(cufftExecC2R(plans[k],inbufs[k],outbufs[k]));
        }
        check_output<<<100,dim3(32,32),0,streams[k]>>>(outbufs[k],checkbufs[k],batch,nfft,outpitch/sizeof(float),&errors[it]);
        ck (cudaGetLastError());
    }
    ck(cudaDeviceSynchronize());

    // report number of errors
    int errcount=0;
    for (int it=0;it<nit;++it)
        if (errors[it])
            ++errcount;
    std::cout << errcount << " of " << nit << " transforms had errors\n";

    for (int i=0;i<nstreams;++i) {
        cudaFree(inbufs[i]);
        cudaFree(outbufs[i]);
        cudaFree(wk_areas[i]);
        cudaStreamDestroy(streams[i]);
        if (i==0 || reuse_plan==false)
            cufftDestroy(plans[i]);
    }
}

int main(int argc,char ** argv)
{
    demo(false);
    demo(true);
    return 0;
}
$ nvcc -o t1241 t1241.cu -lcufft
$ ./t1241
Giving each stream its own dedicated fft plan ... 0 of 100 transforms had errors
Reusing the same fft plan with multiple stream via cufftSetStream ... 0 of 100 transforms had errors
$

【讨论】：

我接受了这个答案，因为它触及了问题的核心，并且提供了 a 的解决方法。问题是 cufft 计划有内部工作缓冲区，在并行流中使用是不安全的。在我看来，将访问与事件同步是一个更好的解决方案。首先，一个合适的袖带计划执行使设备非常忙碌。因此，串行计划执行本身并没有那么昂贵。额外的同步成本似乎可以忽略不计。这保留了没有垃圾输出的单计划内存占用。请参阅我更新的代码中的“FIX”cmets。
额外的内存占用对你来说是不是很大的负担？
是的。我的应用程序中可能有数百个“通道”，每个通道可能需要不同的 FFT 大小。我宁愿不要将内存需求乘以流的数量。
如果我们不使用流，使用所有 cufftExec 的最大大小的工作区是否安全？
您的意思是为不同的计划重复使用单个工作区？是的，那会奏效。