Моя совокупная сумма в numba cuda дает неверные результаты при использовании 1024 потоков
Я пытаюсь реализовать накопительную сумму в numba cuda с параллельным расчетом. Эта функция возьмет одномерный массив ( a = [1, 1, 1, 1]) и вычислит его так, чтобы каждый элемент массива был суммой всех предыдущих элементов ( ans = [1, 2, 3, 4]). Вот функция, как я ее реализовал:
from numba import cuda
import math
import numpy as np
@cuda.jit
def c_sum(array):
tidx = cuda.threadIdx.x
nb_it = math.log2(cuda.blockDim.x)
if tidx < array.size:
for l in range(nb_it):
if tidx > 2**l-1:
array[tidx] += array[tidx-(2**l)]
cuda.syncthreads()
input_ar = np.ones([1024])
c_sum[1,1024](input_ar)
Для работы этой функции требуется, чтобы количество потоков было больше или равно размеру массива. Я тестировал его для массивов размером от 16,32,64,128 до 512, и он работает нормально! По какой-то причине он не работает с массивом, достигающим 1024, а это, к сожалению, то, что мне нужно для моего приложения. 1024 должно быть максимально допустимым количеством потоков в блоке потоков. Последние 5 элементов массива должны быть [1020,1021,1022,1023,1024]. Вот что я получаю:
print(input_ar[-5:])
[1352. 1348. 1344. 1348. 1344.]
На самом деле оно каждый раз меняется. У кого-нибудь есть идея?
Ответы 1
У вас глобальная гонка памяти. Не все потоки выполняются синхронно. Поэтому эта строка приводит к неоднозначному поведению:
array[tidx] += array[tidx-(2**l)]
Например, что, если поток 32 (т. е. tidx == 32) полностью выполнил свою работу, а затем запустился поток 31 (если это необходимо для ясности, предположим, что мы находимся на итерации цикла, когда l равно нулю). И наоборот, что, если поток 31 полностью выполнил свою работу, а затем запустился поток 32. Для этих двух разных случаев можно было бы ожидать двух разных результатов.
Исправить это можно, отделив грузы от магазинов соответствующими барьерами:
$ cat t3.py
from numba import cuda
import math
import numpy as np
@cuda.jit
def c_sum(array):
tidx = cuda.threadIdx.x
nb_it = math.log2(cuda.blockDim.x)
if tidx < array.size:
for l in range(nb_it):
if tidx > 2**l-1:
c = array[tidx-(2**l)]
cuda.syncthreads()
if tidx > 2**l-1:
array[tidx] += c
cuda.syncthreads()
input_ar = np.ones([1024])
c_sum[1,1024](input_ar)
print(input_ar[-5:])
$ python3 t3.py
/home/bob/.local/lib/python3.10/site-packages/numba/cuda/dispatcher.py:536: NumbaPerformanceWarning: Grid size 1 will likely result in GPU under-utilization due to low occupancy.
warn(NumbaPerformanceWarning(msg))
/home/bob/.local/lib/python3.10/site-packages/numba/cuda/cudadrv/devicearray.py:888: NumbaPerformanceWarning: Host array used in CUDA kernel will incur copy overhead to/from device.
warn(NumbaPerformanceWarning(msg))
[1020. 1021. 1022. 1023. 1024.]
$
Конечно, это не единственный возможный способ исправить ситуацию. Вы можете использовать временный массив и играть в пинг-понг взад и вперед, чтобы операция не работала полностью «на месте». И, вероятно, есть другие возможные исправления. Из соображений производительности вам, вероятно, также посоветуют не выполнять такого рода работу исключительно в глобальной памяти, а также использовать и общую память. Но проблема явно игрушечная (действительно, с некоторыми шаблонными предупреждениями, представленными в последних версиях numba cuda на этот счет, как можно видеть в выводе выше). Достаточно сказать, что для такого рода работы существуют лучшие, быстрые и эффективные методы.
Эту методологию нельзя напрямую расширить за пределы одного блока потоков, поскольку cuda.syncthreads() — это всего лишь барьер выполнения для потоков в блоке. Чтобы выйти за пределы одного блока, вам нужно будет использовать барьер всей сетки (всего устройства) или же использовать альтернативную формулировку для получения суммы префикса.
Вау, спасибо. Я почему-то не осознавал, что ссылаюсь на другие темы, и думал, что cuda.syncthreads() решит любые проблемы. Большое спасибо. По практическим соображениям я думаю попытаться реализовать это с помощью временного массива.