PyCUDA Kernel 针对特定计算返回不一致的除法结果答案

【问题标题】：PyCUDA Kernel returns inconsistent division result for specific calculationsPyCUDA Kernel 针对特定计算返回不一致的除法结果
【发布时间】：2020-09-23 01:17:41
【问题描述】：

我正在尝试实现一个计算百分比平均值的内核。

示例- 取 3D 数组（在下面的代码中）块 [[2,4],[3,6],[4,8]] 并计算 (4+6+8)/((4+6+8)+(2+3+4))

这是一个快速运行以下代码的 colab notebook：https://colab.research.google.com/drive/1k_XfOVOYWOTnNQFA9Vo_H93D9l-xWO8K?usp=sharing

# -*- coding: utf-8 -*-

import numpy as np
import pycuda.autoinit
import pycuda.driver as cuda
import pycuda.gpuarray as gpuarray
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule

# set dimentions
ROWS = 3
COLS = 2

h_perms = np.array([[
    [ 1,1],
    [ 1,1],
    [ 1,1]
],[
    [ 2,7],
    [ 3,11],
    [ 4,13]
],[
    [ 2,4],
    [ 3,6],
    [ 4,8]
],[
    [ 2,7],
    [ 3,11],
    [ 4,13]
],[
    [ 2,4],
    [ 3,6],
    [ 4,8]
],[
    [ 1,1],
    [ 1,1],
    [ 1,1]
]
], dtype=np.float32).flatten()

# send to device
d_perms = gpuarray.to_gpu(h_perms)

kernel = SourceModule("""
  __global__ 
  void calc(float *permutations, int *permutationShape, float *results)
  {
    __shared__ float c;
    __shared__ float b;
    int bIdx = blockIdx.y * gridDim.x + blockIdx.x;
    int tIdx = threadIdx.y * blockDim.x + threadIdx.x;
    int rowCount = permutationShape[0];
    int colCount = permutationShape[1];
    int i = (bIdx * rowCount * colCount) + (tIdx * colCount);
    c += permutations[i];
    b += permutations[i+1];
    __syncthreads();
    results[bIdx] = b / (b + c);
  }
  """)

calc = kernel.get_function('calc')

# prepare results array
d_results = gpuarray.zeros((6,1), np.float32)
d_results = gpuarray.to_gpu(d_results)

h_perms_shape = np.array([ROWS,COLS], np.int32);
d_perms_shape = gpuarray.to_gpu(h_perms_shape);

start = cuda.Event()
end = cuda.Event()
start.record()
calc(d_perms, d_perms_shape, d_results, block=(ROWS,1,1), grid=(ROWS*COLS,1,1))
end.record()
secs = start.time_till(end)*1e-3
print(secs)

print(d_results)

我希望得到这个-

array([[0.5      ],
       [0.775],
       [0.6666667],
       [0.775],
       [0.6666667],
       [0.5      ]], dtype=float32)

但我明白了-

array([[0.5      ],
       [0.7777778],
       [0.6666667],
       [0.7777778],
       [0.6666667],
       [0.5      ]], dtype=float32)

我试图理解为什么 (7+11+13)/((7+11+13)+(2+3+4)) 的特定计算结果不是 0.775

【问题讨论】：

你在那个内核中有几个内存竞争。让多个线程同时写入内存位置是未定义的行为
@talonmies 感谢您的评论，希望您能帮助解决后续问题-（1）__syncthreads() 不应该解决这个问题吗？ (2) 你认为我应该在每个线程计算中解决它并让每个线程的for 循环来处理它吗？
线程同步在这里没有任何帮助。您有多个线程读取和写入（并可能在寄存器中缓存）b 和 c，然后在输出计算中使用它们。除非您以某种方式对它们进行序列化（每个线程一个计算、原子操作、并行缩减），否则这将永远不会起作用。这是一个玩具示例，我不知道什么是最好的——也许atomicAdd。
请注意您对事件的使用不正确。在执行时间计算之前，您需要在end 上进行同步。同样，这是一个简短的运行玩具示例，因此它可能会意外运行，但如果内核执行时间落后于 python 解释器，我猜你会在 secs 计算时遇到运行时错误
@talonmies 谢谢！我会进行修改并检查新结果并更新:)

标签： cuda pycuda

【解决方案1】：

您发布的代码在此处包含内存竞争：

int i = (bIdx * rowCount * colCount) + (tIdx * colCount);
c += permutations[i];
b += permutations[i+1];

因为b 和c 在共享内存中，您将有多个线程同时尝试从/向相同的内存位置读取和写入，这在 CUDA 中是未定义的行为（除非在极其特殊的条件下不'不适用于此处）。

如果我把它写成一个玩具例子，我可能会这样做：

  __global__ 
  void calc(float *permutations, int *permutationShape, float *results)
  {
    __shared__ float c;
    __shared__ float b;
    int bIdx = blockIdx.y * gridDim.x + blockIdx.x;
    int tIdx = threadIdx.y * blockDim.x + threadIdx.x;
    int rowCount = permutationShape[0];
    int colCount = permutationShape[1];
    int i = (bIdx * rowCount * colCount) + (tIdx * colCount);
    atomicAdd(&c, permutations[i]);
    atomicAdd(&b, permutations[i+1]);
    __syncthreads();
    if (tIdx == 0) {
        results[bIdx] =  b / (b + c);;
    }
  }

在此代码中，atomicAdd 确保添加和内存事务按顺序发生，从而避免内存竞争。从性能角度来看，这对于不那么琐碎的示例来说不是一个好的解决方案（请查看共享内存减少技术），但它应该可以按预期工作。

【讨论】：