使用共享内存计算点之间的距离
Calculating distances between points using shared memory
我正在尝试计算所有点之间的距离（公制加权）。为了加快速度，我在 gpu 上并通过 cuda 和 numba 执行此操作，因为我认为它更具可读性和更易于使用。
我有两个一维点数组，我想计算同一数组中所有点之间的距离以及两个数组之间所有点之间的距离。我写了两个 cuda 内核，一个只使用全局内存，我已经使用 cpu 代码验证了它给出了正确的答案。就是这个。
@cuda.jit
def gpuSameSample(A,arrSum):
tx = cuda.blockDim.x*cuda.blockIdx.x + cuda.threadIdx.x
temp = A[tx]
tempSum = 0.0
for i in range(tx+1,A.size):
distance = (temp - A[i])**2
tempSum += math.exp(-distance/sigma**2)
arrSum[tx] = tempSum
我现在正尝试通过使用共享内存进一步优化它。这是我目前所拥有的。
@cuda.jit
def gpuSharedSameSample(A,arrSum):
#my block size is equal to 32
sA = cuda.shared.array(shape=(tpb),dtype=float32)
bpg = cuda.gridDim.x
tx = cuda.threadIdx.x + cuda.blockIdx.x *cuda.blockDim.x
count = len(A)
#loop through block by block
tempSum = 0.0
#myPoint = A[tx]
if(tx < count):
myPoint = A[tx]
for currentBlock in range(bpg):
#load in a block to shared memory
copyIdx = (cuda.threadIdx.x + currentBlock*cuda.blockDim.x)
if(copyIdx < count):
sA[cuda.threadIdx.x] = A[copyIdx]
#syncthreads to ensure copying finishes first
cuda.syncthreads()
if((tx < count)):
for i in range(cuda.threadIdx.x,cuda.blockDim.x):
if(copyIdx != tx):
distance = (myPoint - sA[i])**2
tempSum += math.exp(-distance/sigma**2)
#syncthreads here to avoid race conditions if a thread finishes earlier
#arrSum[tx] += tempSum
cuda.syncthreads()
arrSum[tx] += tempSum
我相信我在同步线程方面一直很小心，但这个答案给出的答案总是太大（大约 5%）。我猜一定存在一些竞争条件，但据我了解，每个线程都写入一个唯一索引，并且 tempSum 变量是每个线程的本地变量，因此不应该存在任何竞争条件。我很确定我的 for 循环条件是正确的。任何建议将不胜感激。
谢谢
最好能提供完整的代码。通过对您所展示的内容进行微不足道的添加来做到这一点应该很简单 - 就像我在下面所做的那样。但是，即使使用一组限制性假设，您的两种实现之间也存在差异。
我假设：
- 您的整体数据集大小是线程块大小的整数倍。
- 您启动的线程总数与数据集的大小完全一样。
我也不打算尝试评论您的共享实现是否有意义，即应该比非共享实现表现得更好。这似乎不是您问题的症结所在，这就是为什么您在 2 个实现之间得到数值差异。
主要问题是您在每种情况下 select 计算成对 "distance" 的元素的方法不匹配。在非共享实现中，对于输入数据集中的每个元素 i，您正在计算 i 与每个大于 i:
的元素之间的距离总和
for i in range(tx+1,A.size):
^^^^^^^^^^^
要求和的这 selection 项目与共享实现不匹配：
for i in range(cuda.threadIdx.x,cuda.blockDim.x):
if(copyIdx != tx):
这里有几个问题，但应该清楚地表明，对于每个复制的块，位置 threadIdx.x 的给定元素仅在块（数据）中的目标元素时才更新其总和大于那个指数。这意味着当您逐块浏览整个数据集时，您将跳过每个块中的元素。这不可能与非共享实现匹配。如果这不明显，则只需 select for 循环范围的实际值。假设 cuda.threadIdx.x 为 5，而 cuda.blockDim.x 为 32。那么该特定元素将只计算整个数组中每个数据块中第 6-31 项的总和。
这个问题的解决方案是强制共享实现与非共享实现对齐，就select元素对运行总和的贡献而言。
此外，在非共享实现中，您只更新了一次输出点，并且您正在进行直接赋值：
arrSum[tx] = tempSum
在共享实现中，您仍然只更新一次输出点，但是您没有进行直接赋值。我更改了它以匹配非共享：
arrSum[tx] += tempSum
这是解决了这些问题的完整代码：
$ cat t49.py
from numba import cuda
import numpy as np
import math
import time
from numba import float32
sigma = np.float32(1.0)
tpb = 32
@cuda.jit
def gpuSharedSameSample(A,arrSum):
#my block size is equal to 32
sA = cuda.shared.array(shape=(tpb),dtype=float32)
bpg = cuda.gridDim.x
tx = cuda.threadIdx.x + cuda.blockIdx.x *cuda.blockDim.x
count = len(A)
#loop through block by block
tempSum = 0.0
#myPoint = A[tx]
if(tx < count):
myPoint = A[tx]
for currentBlock in range(bpg):
#load in a block to shared memory
copyIdx = (cuda.threadIdx.x + currentBlock*cuda.blockDim.x)
if(copyIdx < count): #this should always be true
sA[cuda.threadIdx.x] = A[copyIdx]
#syncthreads to ensure copying finishes first
cuda.syncthreads()
if((tx < count)): #this should always be true
for i in range(cuda.blockDim.x):
if(copyIdx-cuda.threadIdx.x+i > tx):
distance = (myPoint - sA[i])**2
tempSum += math.exp(-distance/sigma**2)
#syncthreads here to avoid race conditions if a thread finishes earlier
#arrSum[tx] += tempSum
cuda.syncthreads()
arrSum[tx] = tempSum
@cuda.jit
def gpuSameSample(A,arrSum):
tx = cuda.blockDim.x*cuda.blockIdx.x + cuda.threadIdx.x
temp = A[tx]
tempSum = 0.0
for i in range(tx+1,A.size):
distance = (temp - A[i])**2
tempSum += math.exp(-distance/sigma**2)
arrSum[tx] = tempSum
size = 128
threads_per_block = tpb
blocks = (size + (threads_per_block - 1)) // threads_per_block
my_in = np.ones( size, dtype=np.float32)
my_out = np.zeros(size, dtype=np.float32)
gpuSameSample[blocks, threads_per_block](my_in, my_out)
print(my_out)
gpuSharedSameSample[blocks, threads_per_block](my_in, my_out)
print(my_out)
$ python t49.py
[ 127. 126. 125. 124. 123. 122. 121. 120. 119. 118. 117. 116.
115. 114. 113. 112. 111. 110. 109. 108. 107. 106. 105. 104.
103. 102. 101. 100. 99. 98. 97. 96. 95. 94. 93. 92.
91. 90. 89. 88. 87. 86. 85. 84. 83. 82. 81. 80.
79. 78. 77. 76. 75. 74. 73. 72. 71. 70. 69. 68.
67. 66. 65. 64. 63. 62. 61. 60. 59. 58. 57. 56.
55. 54. 53. 52. 51. 50. 49. 48. 47. 46. 45. 44.
43. 42. 41. 40. 39. 38. 37. 36. 35. 34. 33. 32.
31. 30. 29. 28. 27. 26. 25. 24. 23. 22. 21. 20.
19. 18. 17. 16. 15. 14. 13. 12. 11. 10. 9. 8.
7. 6. 5. 4. 3. 2. 1. 0.]
[ 127. 126. 125. 124. 123. 122. 121. 120. 119. 118. 117. 116.
115. 114. 113. 112. 111. 110. 109. 108. 107. 106. 105. 104.
103. 102. 101. 100. 99. 98. 97. 96. 95. 94. 93. 92.
91. 90. 89. 88. 87. 86. 85. 84. 83. 82. 81. 80.
79. 78. 77. 76. 75. 74. 73. 72. 71. 70. 69. 68.
67. 66. 65. 64. 63. 62. 61. 60. 59. 58. 57. 56.
55. 54. 53. 52. 51. 50. 49. 48. 47. 46. 45. 44.
43. 42. 41. 40. 39. 38. 37. 36. 35. 34. 33. 32.
31. 30. 29. 28. 27. 26. 25. 24. 23. 22. 21. 20.
19. 18. 17. 16. 15. 14. 13. 12. 11. 10. 9. 8.
7. 6. 5. 4. 3. 2. 1. 0.]
$
请注意，如果违反了我的两个假设中的任何一个，则您的代码还有其他问题。
以后，我鼓励您提供一个简短、完整的代码，如上所示。对于这样的问题，它应该不会有太多额外的工作。如果没有其他原因（还有其他原因），在你已经有了代码的情况下，强迫别人从头开始编写这段代码是很乏味的，以证明所提供答案的敏感性。
我正在尝试计算所有点之间的距离（公制加权）。为了加快速度，我在 gpu 上并通过 cuda 和 numba 执行此操作，因为我认为它更具可读性和更易于使用。
我有两个一维点数组，我想计算同一数组中所有点之间的距离以及两个数组之间所有点之间的距离。我写了两个 cuda 内核，一个只使用全局内存，我已经使用 cpu 代码验证了它给出了正确的答案。就是这个。
@cuda.jit
def gpuSameSample(A,arrSum):
tx = cuda.blockDim.x*cuda.blockIdx.x + cuda.threadIdx.x
temp = A[tx]
tempSum = 0.0
for i in range(tx+1,A.size):
distance = (temp - A[i])**2
tempSum += math.exp(-distance/sigma**2)
arrSum[tx] = tempSum
我现在正尝试通过使用共享内存进一步优化它。这是我目前所拥有的。
@cuda.jit
def gpuSharedSameSample(A,arrSum):
#my block size is equal to 32
sA = cuda.shared.array(shape=(tpb),dtype=float32)
bpg = cuda.gridDim.x
tx = cuda.threadIdx.x + cuda.blockIdx.x *cuda.blockDim.x
count = len(A)
#loop through block by block
tempSum = 0.0
#myPoint = A[tx]
if(tx < count):
myPoint = A[tx]
for currentBlock in range(bpg):
#load in a block to shared memory
copyIdx = (cuda.threadIdx.x + currentBlock*cuda.blockDim.x)
if(copyIdx < count):
sA[cuda.threadIdx.x] = A[copyIdx]
#syncthreads to ensure copying finishes first
cuda.syncthreads()
if((tx < count)):
for i in range(cuda.threadIdx.x,cuda.blockDim.x):
if(copyIdx != tx):
distance = (myPoint - sA[i])**2
tempSum += math.exp(-distance/sigma**2)
#syncthreads here to avoid race conditions if a thread finishes earlier
#arrSum[tx] += tempSum
cuda.syncthreads()
arrSum[tx] += tempSum
我相信我在同步线程方面一直很小心，但这个答案给出的答案总是太大（大约 5%）。我猜一定存在一些竞争条件，但据我了解，每个线程都写入一个唯一索引，并且 tempSum 变量是每个线程的本地变量，因此不应该存在任何竞争条件。我很确定我的 for 循环条件是正确的。任何建议将不胜感激。
谢谢
最好能提供完整的代码。通过对您所展示的内容进行微不足道的添加来做到这一点应该很简单 - 就像我在下面所做的那样。但是，即使使用一组限制性假设，您的两种实现之间也存在差异。
我假设：
- 您的整体数据集大小是线程块大小的整数倍。
- 您启动的线程总数与数据集的大小完全一样。
我也不打算尝试评论您的共享实现是否有意义，即应该比非共享实现表现得更好。这似乎不是您问题的症结所在，这就是为什么您在 2 个实现之间得到数值差异。
主要问题是您在每种情况下 select 计算成对 "distance" 的元素的方法不匹配。在非共享实现中，对于输入数据集中的每个元素 i，您正在计算 i 与每个大于 i:
for i in range(tx+1,A.size):
^^^^^^^^^^^
要求和的这 selection 项目与共享实现不匹配：
for i in range(cuda.threadIdx.x,cuda.blockDim.x):
if(copyIdx != tx):
这里有几个问题，但应该清楚地表明，对于每个复制的块，位置 threadIdx.x 的给定元素仅在块（数据）中的目标元素时才更新其总和大于那个指数。这意味着当您逐块浏览整个数据集时，您将跳过每个块中的元素。这不可能与非共享实现匹配。如果这不明显，则只需 select for 循环范围的实际值。假设 cuda.threadIdx.x 为 5，而 cuda.blockDim.x 为 32。那么该特定元素将只计算整个数组中每个数据块中第 6-31 项的总和。
这个问题的解决方案是强制共享实现与非共享实现对齐，就select元素对运行总和的贡献而言。
此外，在非共享实现中，您只更新了一次输出点，并且您正在进行直接赋值：
arrSum[tx] = tempSum
在共享实现中，您仍然只更新一次输出点，但是您没有进行直接赋值。我更改了它以匹配非共享：
arrSum[tx] += tempSum
这是解决了这些问题的完整代码：
$ cat t49.py
from numba import cuda
import numpy as np
import math
import time
from numba import float32
sigma = np.float32(1.0)
tpb = 32
@cuda.jit
def gpuSharedSameSample(A,arrSum):
#my block size is equal to 32
sA = cuda.shared.array(shape=(tpb),dtype=float32)
bpg = cuda.gridDim.x
tx = cuda.threadIdx.x + cuda.blockIdx.x *cuda.blockDim.x
count = len(A)
#loop through block by block
tempSum = 0.0
#myPoint = A[tx]
if(tx < count):
myPoint = A[tx]
for currentBlock in range(bpg):
#load in a block to shared memory
copyIdx = (cuda.threadIdx.x + currentBlock*cuda.blockDim.x)
if(copyIdx < count): #this should always be true
sA[cuda.threadIdx.x] = A[copyIdx]
#syncthreads to ensure copying finishes first
cuda.syncthreads()
if((tx < count)): #this should always be true
for i in range(cuda.blockDim.x):
if(copyIdx-cuda.threadIdx.x+i > tx):
distance = (myPoint - sA[i])**2
tempSum += math.exp(-distance/sigma**2)
#syncthreads here to avoid race conditions if a thread finishes earlier
#arrSum[tx] += tempSum
cuda.syncthreads()
arrSum[tx] = tempSum
@cuda.jit
def gpuSameSample(A,arrSum):
tx = cuda.blockDim.x*cuda.blockIdx.x + cuda.threadIdx.x
temp = A[tx]
tempSum = 0.0
for i in range(tx+1,A.size):
distance = (temp - A[i])**2
tempSum += math.exp(-distance/sigma**2)
arrSum[tx] = tempSum
size = 128
threads_per_block = tpb
blocks = (size + (threads_per_block - 1)) // threads_per_block
my_in = np.ones( size, dtype=np.float32)
my_out = np.zeros(size, dtype=np.float32)
gpuSameSample[blocks, threads_per_block](my_in, my_out)
print(my_out)
gpuSharedSameSample[blocks, threads_per_block](my_in, my_out)
print(my_out)
$ python t49.py
[ 127. 126. 125. 124. 123. 122. 121. 120. 119. 118. 117. 116.
115. 114. 113. 112. 111. 110. 109. 108. 107. 106. 105. 104.
103. 102. 101. 100. 99. 98. 97. 96. 95. 94. 93. 92.
91. 90. 89. 88. 87. 86. 85. 84. 83. 82. 81. 80.
79. 78. 77. 76. 75. 74. 73. 72. 71. 70. 69. 68.
67. 66. 65. 64. 63. 62. 61. 60. 59. 58. 57. 56.
55. 54. 53. 52. 51. 50. 49. 48. 47. 46. 45. 44.
43. 42. 41. 40. 39. 38. 37. 36. 35. 34. 33. 32.
31. 30. 29. 28. 27. 26. 25. 24. 23. 22. 21. 20.
19. 18. 17. 16. 15. 14. 13. 12. 11. 10. 9. 8.
7. 6. 5. 4. 3. 2. 1. 0.]
[ 127. 126. 125. 124. 123. 122. 121. 120. 119. 118. 117. 116.
115. 114. 113. 112. 111. 110. 109. 108. 107. 106. 105. 104.
103. 102. 101. 100. 99. 98. 97. 96. 95. 94. 93. 92.
91. 90. 89. 88. 87. 86. 85. 84. 83. 82. 81. 80.
79. 78. 77. 76. 75. 74. 73. 72. 71. 70. 69. 68.
67. 66. 65. 64. 63. 62. 61. 60. 59. 58. 57. 56.
55. 54. 53. 52. 51. 50. 49. 48. 47. 46. 45. 44.
43. 42. 41. 40. 39. 38. 37. 36. 35. 34. 33. 32.
31. 30. 29. 28. 27. 26. 25. 24. 23. 22. 21. 20.
19. 18. 17. 16. 15. 14. 13. 12. 11. 10. 9. 8.
7. 6. 5. 4. 3. 2. 1. 0.]
$
请注意，如果违反了我的两个假设中的任何一个，则您的代码还有其他问题。
以后，我鼓励您提供一个简短、完整的代码，如上所示。对于这样的问题，它应该不会有太多额外的工作。如果没有其他原因（还有其他原因），在你已经有了代码的情况下，强迫别人从头开始编写这段代码是很乏味的，以证明所提供答案的敏感性。