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的主要内容,如果未能解决你的问题,请参考以下文章