Karatsuba - 多项式乘法与 CUDA

Posted

技术标签:

【中文标题】Karatsuba - 多项式乘法与 CUDA【英文标题】:Karatsuba - polynomials multiplication with CUDA 【发布时间】:2018-10-06 10:13:40 【问题描述】:

我将 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);

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

【问题讨论】:

我从未在 CUDA 工作过。话虽如此,如果答案总是不同,则问题很可能在于数据依赖性。我可能会抓挠我的头发来弄清楚。你可能会有更好的运气。 在这段代码中理解CUDA很简单。您只需要替换 for 循环if 条件 并使用大量线程运行该代码。但是谢谢,我会试着弄清楚为什么会有数据依赖。 【参考方案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 应该是可用的。

其他几个cmets:

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

    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 会有所帮助)。

【讨论】:

现在我可以看到结果总是一样的,但仍然是错误的。你能配置为什么吗? TYPE 只是 int 我很确定在这种情况下这种类型是可以的。我相信问题出在我的网格尺寸和块数上。在我看来,我的代码现在对多项式的每个偶数系数执行了两次。 看看我编辑的答案中的注释 1 和 3。如果这些不能解决问题,那么我需要查看minimal reproducible example

以上是关于Karatsuba - 多项式乘法与 CUDA的主要内容,如果未能解决你的问题,请参考以下文章

最小二乘法拟合与多项式拟合的关系是啥?

c++编程 多项式的乘法

C语言求多项式乘法

[PAT] 一元多项式的乘法与加法运算 C语言实现

7-22 一元多项式的乘法与加法运算 (20 分)

02-线性结构2 一元多项式的乘法与加法运算