Karatsuba

Question

我正在使用 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);

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

Answer 1

像这样的问题，你应该提供一个MCVE。（参见第 1 项 here）例如，我不知道 TYPE 表示什么类型，这对我提出的解决方案的正确性很重要。

在你的第一个内核中，整个网格中只有一个线程在读写位置result[i]。但是在你的第二个内核中，你现在有多个线程写入 result[i] 位置。他们互相冲突。 CUDA 没有指定线程运行的顺序，有些线程可能在其他线程之前、之后或同时运行。在这种情况下，某些线程可能会与其他线程同时读取 result[i]。然后，当线程写入它们的结果时，它们将不一致。它可能从运行到运行不等。你在那里有一个竞争条件（执行顺序依赖性，而不是数据依赖性）。

解决这个问题的规范方法是采用 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 类型和您使用的 TYPE 的实际类型，这可能无法按原样工作（可能无法编译）。但是由于您之前使用 TYPE 作为循环变量，我假设它是一个整数类型，并且应该提供必要的 atomicAdd。

其他一些评论：

这可能不是您期望的网格大小：

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

我想通常的计算应该是：

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

我总是建议 proper CUDA error checking 和运行将您的代码与 cuda-memcheck 结合使用，每当您在使用 CUDA 代码时遇到问题，以确保有没有运行时间错误。
在我看来也是这样的：
```
if(j >= start && j <= end ){
```
应该是这样的：
```
if(j >= start && j < end ){
```
以匹配您的 for 循环范围。我还假设 size 小于 resultSize（同样，MCVE 会有所帮助）。

Karatsuba - 使用 CUDA 进行多项式乘法

Karatsuba - polynomials multiplication with CUDA

c++

parallel-processing

cuda

nvidia