代码之家 › 专栏 › 技术社区 › user8005765

Karatsuba-多项式与CUDA相乘

karatsuba nvidia cuda parallel-processing c++

user8005765 · 技术社区 · 6 年前

我正在使用CUDA进行迭代Karatsuba算法,我想问一下,为什么一条直线的计算总是不同的。

首先,我实现了这个函数,它总是正确地计算结果:

__global__ void kernel_res_main(TYPE *A, TYPE *B, TYPE *D, TYPE *result, TYPE size, TYPE resultSize){
    int i = blockDim.x * blockIdx.x + threadIdx.x;

    if( i > 0 && i < resultSize - 1){

        TYPE start = (i >= size) ? (i % size ) + 1 : 0;


        TYPE end = (i + 1) / 2;


        for(TYPE inner = start; inner < end; inner++){
            result[i] += ( A[inner] + A[i - inner] ) * ( B[inner] + B[i - inner] );
            result[i] -= ( D[inner] + D[i-inner] );
        }
    }
}

现在我想使用2D网格并将CUDA用于for循环,因此我将函数改为:

__global__ void kernel_res_nested(TYPE *A, TYPE *B, TYPE *D, TYPE *result, TYPE size, TYPE resultSize){

    int i = blockDim.x * blockIdx.x + threadIdx.x;
    int j = blockDim.y * blockIdx.y + threadIdx.y;

    TYPE rtmp = result[i];

    if( i > 0 && i < resultSize - 1){

        TYPE start = (i >= size) ? (i % size ) + 1 : 0;
        TYPE end = (i + 1) >> 1;

        if(j >= start && j <= end ){

           // WRONG 
           rtmp += ( A[j] + A[i - j] ) * ( B[j] + B[i - j] ) - ( D[j] + D[i - j] );
        }
    }

    result[i] = rtmp;
}

我这样调用此函数:

dim3 block( 32, 8 );
dim3 grid( (resultSize+1/32) , (resultSize+7/8) );
kernel_res_nested <<<grid, block>>> (devA, devB, devD, devResult, size, resultSize);

结果总是错误的,总是不同的。我不明白为什么第二个实现是错误的,并且总是计算错误的结果。我看不出有任何与数据依赖性相关的逻辑问题。有人知道我如何解决这个问题吗?

1 回复 | 直到 5 年前

Robert Crovella 6 年前

对于这样的问题,您应该提供MCVE。(见第1项 here )例如,我不知道表示什么类型 TYPE ,这确实关系到我将提出的解决方案的正确性。

在第一个内核中,整个网格中只有一个线程在读写位置 result[i] 。但是在第二个内核中,现在有多个线程写入 结果[一] 地方他们相互冲突。CUDA没有指定线程的运行顺序,有些线程可以在其他线程之前、之后或与其他线程同时运行。在这种情况下,某些线程可能会读取 结果[一] 和其他人一样。然后,当线程写入其结果时,它们将不一致。而且每次跑步都可能有所不同。你有一个 比赛条件 存在(执行顺序依赖关系,而不是数据依赖关系)。

解决这个问题的规范方法是使用 reduction 技巧

然而,为了简单起见,我建议 atomics 可以帮你解决。根据您所展示的内容,这更容易实现,并有助于确认比赛条件。在那之后,如果你想尝试一种简化方法,这里有大量的教程(上面链接了一个)和大量的问题 cuda 标记它。

您可以将内核修改为如下内容,以整理竞争条件:

__global__ void kernel_res_nested(TYPE *A, TYPE *B, TYPE *D, TYPE *result, TYPE size, TYPE resultSize){

    int i = blockDim.x * blockIdx.x + threadIdx.x;
    int j = blockDim.y * blockIdx.y + threadIdx.y;

    if( i > 0 && i < resultSize - 1){

        TYPE start = (i >= size) ? (i % size ) + 1 : 0;
        TYPE end = (i + 1) >> 1;

        if(j >= start && j < end ){ // see note below

           atomicAdd(result+i, (( A[j] + A[i - j] ) * ( B[j] + B[i - j] ) - ( D[j] + D[i - j] )));
        }
    }

}

请注意,根据您的GPU类型和 类型 您正在使用,这可能无法正常工作(可能无法编译)。但是因为你以前用过 类型 作为循环变量,我假设它是整数类型,并且 atomicAdd 对于那些应该是可用的。

其他一些评论:

这可能无法提供您期望的网格大小:

dim3 grid( (resultSize+1/32) , (resultSize+7/8) );

我认为通常的计算是:

dim3 grid( (resultSize+31)/32, (resultSize+7)/8 );

我总是建议 proper CUDA error checking 并使用 cuda-memcheck ,以确保没有运行时错误。
在我看来也是这样:
```
if(j >= start && j <= end ){
```
应该是这样的:
```
if(j >= start && j < end ){
```
以匹配for循环范围。我还假设 size 小于 resultSize (同样,MCVE会有所帮助)。