algorithm - TERCOM 算法 - 在 CUDA 中从单线程变为多线程

Question

我目前正在移植 TERCOM algorithm从仅使用 1 个线程到使用多个线程。简而言之，TERCOM 算法接收 5 个测量值和航向，并将该测量值与预存 map 进行比较。该算法将选择最佳匹配，即最低平均绝对差 (MAD)，并返回位置。

该代码在一个线程和 for 循环中工作得很好，但是当我尝试使用多个线程和 block 时，它返回了错误的答案。似乎多线程版本不像单线程版本那样“运行”计算。有谁知道我做错了什么？

这是使用for循环的代码

__global__ void kernel (int m, int n, int h, int N, float *f, float heading, float *measurements) 
{
    //Without threads
    float pos[2]={0};
    float theta=heading*(PI/180);
    float MAD=0;

    // Calculate how much to move in x and y direction
    float offset_x = h*cos(theta);
    float offset_y = -h*sin(theta); 

    float min=100000; //Some High value

    //Calculate Mean Absolute Difference
    for(float row=0;row<m;row++)
    {
        for(float col=0;col<n;col++)
        {
            for(float g=0; g<N; g++)
            {
                f[(int)g] = tex2D (tex, col+(g-2)*offset_x+0.5f, row+(g-2)*offset_y+0.5f);
                MAD += abs(measurements[(int)g]-f[(int)g]);
            }
            if(MAD<min) 
            {
                min=MAD;
                pos[0]=col;
                pos[1]=row;
            }
            MAD=0;                  //Reset MAD
        }
    }

    f[0]=min;
    f[1]=pos[0];
    f[2]=pos[1];
}

这是我尝试使用多线程

__global__ void kernel (int m, int n, int h, int N, float *f, float heading, float *measurements) 
{
    // With threads
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int idy = blockIdx.y * blockDim.y + threadIdx.y;
    float pos[2]={0};
    float theta=heading*(PI/180);
    float MAD=0;

    // Calculate how much to move in x and y direction
    float offset_x = h*cos(theta);
    float offset_y = -h*sin(theta); 

    float min=100000; //Some High value

    if(idx < n && idy < m)
    {
        for(float g=0; g<N; g++)
        {
            f[(int)g] = tex2D (tex, idx+(g-2)*offset_x+0.5f, idy+(g-2)*offset_y+0.5f);
            MAD += abs(measurements[(int)g]-f[(int)g]); 
        }

        if(MAD<min) 
        {
            min=MAD;
            pos[0]=idx;
            pos[1]=idy;
        }
        MAD=0;                  //Reset MAD
    }
    f[0]=min;
    f[1]=pos[0];
    f[2]=pos[1];
}

启动内核

dim3 dimBlock( 16,16 );
dim3 dimGrid;
dimGrid.x = (n + dimBlock.x - 1)/dimBlock.x;
dimGrid.y = (m + dimBlock.y - 1)/dimBlock.y;

kernel <<< dimGrid,dimBlock >>> (m, n, h, N, dev_results, heading, dev_measurements);

Answer 1

这里的基本问题是您在代码中存在内存竞争，主要围绕使用 f 作为某种线程局部暂存空间和输出变量。每个并发线程都将尝试同时将值写入 f 中的相同位置，这将产生未定义的行为。

据我所知，根本不需要使用 f 作为暂存空间，内核的主要计算部分可以写成这样:

if(idx < n && idy < m)
{
    for(float g=0; g<N; g++)
    {
        float fval = tex2D (tex, idx+(g-2)*offset_x+0.5f, idy+(g-2)*offset_y+0.5f);
        MAD += abs(measurements[(int)g]-fval); 
    }
    min=MAD;
    pos[0]=idx;
    pos[1]=idy;
}

[免责声明:用浏览器编写，使用风险自负]

在该计算结束时，每个线程 都有自己的min 和pos 值。至少这些必须存储在唯一的全局内存中(即输出必须为每个线程结果提供足够的空间)。然后，您将需要执行某种缩减操作，以从线程局部值集中获得全局最小值。这可能在主机中，或在设备代码中，或两者的某种组合。已经有很多代码可用于 CUDA 并行缩减，您应该能够通过搜索和/或查看 CUDA 工具包提供的示例找到这些代码。使它们适应您需要保留位置和最小值的特定情况应该是微不足道的。