TERCOM算法 - 在CUDA中从单线程切换到多线程

Question

我目前正致力于使用仅使用1个线程的端口TERCOM algorithm来使用多线程。简而言之，TERCOM算法接收5个测量值和航向，并将这些测量结果与预先存储的地图进行比较。该算法将选择最佳匹配，即最低平均绝对差（MAD），并返回该位置。

该代码完美工作与一个线程和for循环，但是当我尝试使用多个线程和块它返回错误的答案。看起来多线程版本不像“单线程”版本那样“运行”计算。有谁知道我做错了什么？

下面是使用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作为某种线索d本地暂存空间和输出变量。每个并发线程将尝试同时将值写入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工具包提供的示例来查找这些代码。将它们调整到您需要保留位置以及最小值的指定情况应该是微不足道的。