Cuda распараллелить ядро - PullRequest
Новая
       19

Cuda распараллелить ядро

0 голосов
/ 01 сентября 2018

Я пытаюсь распараллелить простой цикл обновления симуляции на GPU. По сути, существует группа «существ», представленных кружками, которые в каждом цикле обновления будут двигаться, а затем будет проверяться, пересекаются ли какие-либо из них. радиус - это радиус различных типов существ. 

import numpy as np
import math
from numba import cuda


@cuda.jit('void(float32[:], float32[:], float32[:], uint8[:], float32[:], float32[:], float32, uint32, uint32)')
def update(p_x, p_y, radii, types, velocities, max_velocities, acceleration, num_creatures, cycles):
    for c in range(cycles):
        for i in range(num_creatures):
            velocities[i] = velocities[i] + acceleration
            if velocities[i] > max_velocities[i]:
                velocities[i] = max_velocities[i]
            p_x[i] = p_x[i] + (math.cos(1.0) * velocities[i])
            p_y[i] = p_y[i] + (math.sin(1.0) * velocities[i])
        for i in range(num_creatures):
            for j in range(i, num_creatures):
                delta_x = p_x[j] - p_x[i]
                delta_y = p_y[j] - p_y[i]
                distance_squared = (delta_x * delta_x) + (delta_y * delta_y)
                sum_of_radii = radii[types[i]] + radii[types[i]]
                if distance_squared < sum_of_radii * sum_of_radii:
                    pass


acceleration = .1
creature_radius = 10
spacing = 20
food_radius = 3

max_num_creatures = 1500
num_creatures = 0
max_num_food = 500
num_food = 0
max_num_entities = max_num_creatures + max_num_food
num_entities = 0
cycles = 1


p_x = np.zeros(max_num_entities, dtype=np.float32)
p_y = np.zeros(max_num_entities, dtype=np.float32)
radii = np.array([creature_radius, creature_radius, food_radius], dtype=np.float32)
types = np.zeros(max_num_entities, dtype=np.uint8)

velocities = np.zeros(max_num_creatures, dtype=np.float32)
max_velocities = np.zeros(max_num_creatures, dtype=np.float32)
# types:
# male - 0
# female - 1
# food - 2
for x in range(1, 800 // spacing):
    for y in range(1, 600 // spacing):
        if num_creatures % 2 == 0:
            types[num_creatures] = 0
        else:
            types[num_creatures] = 1
        p_x[num_creatures] = x * spacing
        p_y[num_creatures] = y * spacing
        max_velocities[num_creatures] = 5
        num_creatures += 1


device_p_x = cuda.to_device(p_x)
device_p_y = cuda.to_device(p_y)
device_radii = cuda.to_device(radii)
device_types = cuda.to_device(types)
device_velocities = cuda.to_device(velocities)
device_max_velocities = cuda.to_device(max_velocities)
threadsperblock = 64
blockspergrid = 16
update[blockspergrid, threadsperblock](device_p_x, device_p_y, device_radii, device_types, device_velocities, device_max_velocities,
        acceleration, num_creatures, cycles)
print(device_p_x.copy_to_host())

1.0 в math.cos и math.sin - просто заполнитель для указаний отдельных существ.

Как сейчас, есть несколько потоков, но они выполняют один и тот же код. Какие изменения нужно внести в ядро, чтобы распараллелить его?

1 Ответ

0 голосов
/ 01 сентября 2018

Самым очевидным измерением для распараллеливания мне кажется цикл в i в вашем ядре, который перебирает num_creatures. Поэтому я опишу, как это сделать.

  1. Наша цель будет удалить цикл на num_creatures и вместо этого позволить каждой итерации цикла обрабатываться отдельным потоком CUDA. Это возможно, потому что работа, выполняемая в каждой итерации цикла, (в основном) независима - она ​​не зависит от результатов других итераций цикла (но см. 2 ниже).

  2. Проблема, с которой мы столкнемся, состоит в том, что 2-й цикл i for в num_creatures предположительно зависит от результатов первого. Если мы оставляем все как последовательный код, работающий в одном потоке, то об этой зависимости заботится природа выполнения последовательного кода. Однако мы хотим распараллелить это. Поэтому нам нужна глобальная синхронизация между первым циклом for в num_creatures и вторым. Простая и удобная глобальная синхронизация в CUDA - это запуск ядра, поэтому мы разделим код ядра на две функции ядра. Мы назовем их update1 и update2

  3. Затем возникает проблема с тем, что делать с циклом перегиба в cycles. Мы не можем просто повторить этот цикл в обоих ядрах, потому что это изменило бы функциональное поведение - мы бы вычислили cycles обновлений до p_x, прежде чем вычислять, например, один delta_x. Это, по-видимому, не то, что нужно. Итак, для простоты мы выведем этот цикл из кода ядра обратно в код хоста. Затем код хоста будет вызывать ядра update1 и update2 для cycles итераций.

  4. Мы также хотим сделать обработку ядра адаптируемой к различным размерам num_creatures. Таким образом, мы выберем жестко заданный размер для threadsperblock, но сделаем количество запускаемых блоков переменным в зависимости от размера num_creatures. Для этого нам понадобится проверка потока (начальный оператор if) в каждом из наших ядер, чтобы «лишние» потоки ничего не делали.

С этим описанием мы получим что-то вроде этого: 

$ cat t11.py
import numpy as np
import math
from numba import cuda


@cuda.jit('void(float32[:], float32[:], float32[:], float32[:], float32, uint32)')
def update1(p_x, p_y, velocities, max_velocities, acceleration, num_creatures):
    i = cuda.grid(1)
    if i < num_creatures:
            velocities[i] = velocities[i] + acceleration
            if velocities[i] > max_velocities[i]:
                velocities[i] = max_velocities[i]
            p_x[i] = p_x[i] + (math.cos(1.0) * velocities[i])
            p_y[i] = p_y[i] + (math.sin(1.0) * velocities[i])

@cuda.jit('void(float32[:], float32[:], float32[:], uint8[:], uint32)')
def update2(p_x, p_y, radii, types, num_creatures):
    i = cuda.grid(1)
    if i < num_creatures:
            for j in range(i, num_creatures):
                delta_x = p_x[j] - p_x[i]
                delta_y = p_y[j] - p_y[i]
                distance_squared = (delta_x * delta_x) + (delta_y * delta_y)
                sum_of_radii = radii[types[i]] + radii[types[i]]
                if distance_squared < sum_of_radii * sum_of_radii:
                    pass


acceleration = .1
creature_radius = 10
spacing = 20
food_radius = 3

max_num_creatures = 1500
num_creatures = 0
max_num_food = 500
num_food = 0
max_num_entities = max_num_creatures + max_num_food
num_entities = 0
cycles = 1


p_x = np.zeros(max_num_entities, dtype=np.float32)
p_y = np.zeros(max_num_entities, dtype=np.float32)
radii = np.array([creature_radius, creature_radius, food_radius], dtype=np.float32)
types = np.zeros(max_num_entities, dtype=np.uint8)

velocities = np.zeros(max_num_creatures, dtype=np.float32)
max_velocities = np.zeros(max_num_creatures, dtype=np.float32)
# types:
# male - 0
# female - 1
# food - 2
for x in range(1, 800 // spacing):
    for y in range(1, 600 // spacing):
        if num_creatures % 2 == 0:
            types[num_creatures] = 0
        else:
            types[num_creatures] = 1
        p_x[num_creatures] = x * spacing
        p_y[num_creatures] = y * spacing
        max_velocities[num_creatures] = 5
        num_creatures += 1


device_p_x = cuda.to_device(p_x)
device_p_y = cuda.to_device(p_y)
device_radii = cuda.to_device(radii)
device_types = cuda.to_device(types)
device_velocities = cuda.to_device(velocities)
device_max_velocities = cuda.to_device(max_velocities)
threadsperblock = 64
blockspergrid = (num_creatures // threadsperblock) + 1
for i in range(cycles):
    update1[blockspergrid, threadsperblock](device_p_x, device_p_y, device_velocities, device_max_velocities, acceleration, num_creatures)
    update2[blockspergrid, threadsperblock](device_p_x, device_p_y, device_radii, device_types, num_creatures)
print(device_p_x.copy_to_host())
$ python t11.py
[ 20.05402946  20.05402946  20.05402946 ...,   0.           0.           0.        ]
$

, который выдает тот же результат, что и исходная опубликованная версия (в исходном коде запущено 16 блоков из 64 потоков, выполняющих одну и ту же функцию и наступающих друг на друга при записи в одни и те же данные. Поэтому я имею в виду оригинальная опубликованная версия, но с одним блоком из одного потока).

Обратите внимание, что, безусловно, есть и другие способы распараллеливания, и, возможно, возможна другая оптимизация, но это должно дать вам разумную отправную точку для работы CUDA.

Как уже упоминалось в в вашем предыдущем вопросе , второе ядро ​​здесь действительно ничего полезного не дает, но я предполагаю, что это просто заполнитель для будущей работы. Я также не уверен, что вы получите то, что хотите, используя radii здесь, но это также не является предметом этого вопроса.

Так каков эффект от всего этого исполнения? Опять же, мы сравним исходную опубликованную версию (t12.py, ниже), выполняющую всего один блок одного потока (не 16 блоков из 64 потоков, что в любом случае было бы хуже), с этой версией, которая работает с 18 блоками из 64 темы (t11.py, ниже): 

$ nvprof --print-gpu-trace python t11.py
==3551== NVPROF is profiling process 3551, command: python t11.py
[ 20.05402946  20.05402946  20.05402946 ...,   0.           0.           0.        ]
==3551== Profiling application: python t11.py
==3551== Profiling result:
   Start  Duration            Grid Size      Block Size     Regs*    SSMem*    DSMem*      Size  Throughput  SrcMemType  DstMemType           Device   Context    Stream  Name
446.77ms  1.8240us                    -               -         -         -         -  7.8125KB  4.0847GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
446.97ms  1.7600us                    -               -         -         -         -  7.8125KB  4.2333GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
447.35ms  1.2160us                    -               -         -         -         -       12B  9.4113MB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
447.74ms  1.3440us                    -               -         -         -         -  1.9531KB  1.3859GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
447.93ms  1.5040us                    -               -         -         -         -  5.8594KB  3.7154GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
448.13ms  1.5360us                    -               -         -         -         -  5.8594KB  3.6380GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
448.57ms  5.4720us             (18 1 1)        (64 1 1)        36        0B        0B         -           -           -           -  Quadro K2000 (0         1         7  cudapy::__main__::update1$241(Array<float, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, float, unsigned int) [49]
448.82ms  1.1200us             (18 1 1)        (64 1 1)         8        0B        0B         -           -           -           -  Quadro K2000 (0         1         7  cudapy::__main__::update2$242(Array<float, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, Array<unsigned char, int=1, A, mutable, aligned>, unsigned int) [50]
448.90ms  2.1120us                    -               -         -         -         -  7.8125KB  3.5277GB/s      Device    Pageable  Quadro K2000 (0         1         7  [CUDA memcpy DtoH]

Regs: Number of registers used per CUDA thread. This number includes registers used internally by the CUDA driver and/or tools and can be more than what the compiler shows.
SSMem: Static shared memory allocated per CUDA block.
DSMem: Dynamic shared memory allocated per CUDA block.
SrcMemType: The type of source memory accessed by memory operation/copy
DstMemType: The type of destination memory accessed by memory operation/copy

$ python t12.py
[ 20.05402946  20.05402946  20.05402946 ...,   0.           0.           0.        ]
$ nvprof --print-gpu-trace python t12.py
==3604== NVPROF is profiling process 3604, command: python t12.py
[ 20.05402946  20.05402946  20.05402946 ...,   0.           0.           0.        ]
==3604== Profiling application: python t12.py
==3604== Profiling result:
   Start  Duration            Grid Size      Block Size     Regs*    SSMem*    DSMem*      Size  Throughput  SrcMemType  DstMemType           Device   Context    Stream  Name
296.22ms  1.8240us                    -               -         -         -         -  7.8125KB  4.0847GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
296.41ms  1.7920us                    -               -         -         -         -  7.8125KB  4.1577GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
296.79ms  1.2160us                    -               -         -         -         -       12B  9.4113MB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
297.21ms  1.3440us                    -               -         -         -         -  1.9531KB  1.3859GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
297.40ms  1.5040us                    -               -         -         -         -  5.8594KB  3.7154GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
297.60ms  1.5360us                    -               -         -         -         -  5.8594KB  3.6380GB/s    Pageable      Device  Quadro K2000 (0         1         7  [CUDA memcpy HtoD]
298.05ms  1.8453ms              (1 1 1)         (1 1 1)        36        0B        0B         -           -           -           -  Quadro K2000 (0         1         7  cudapy::__main__::update$241(Array<float, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, Array<unsigned char, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, Array<float, int=1, A, mutable, aligned>, float, unsigned int, unsigned int) [38]
299.91ms  2.1120us                    -               -         -         -         -  7.8125KB  3.5277GB/s      Device    Pageable  Quadro K2000 (0         1         7  [CUDA memcpy DtoH]

Regs: Number of registers used per CUDA thread. This number includes registers used internally by the CUDA driver and/or tools and can be more than what the compiler shows.
SSMem: Static shared memory allocated per CUDA block.
DSMem: Dynamic shared memory allocated per CUDA block.
SrcMemType: The type of source memory accessed by memory operation/copy
DstMemType: The type of destination memory accessed by memory operation/copy
$

Мы видим, что профилировщик сообщает для исходной версии t12.py, что работает одно ядро ​​update с 1 блоком и 1 потоком, и оно занимает 1,8453 миллисекунды. Для модифицированной версии t11.py, опубликованной в этом ответе, профилировщик сообщает о 18 блоках по 64 потока в каждом для ядер update1 и update2, а общее время выполнения этих двух ядер составляет приблизительно 5,47 + 1,12 = 6,59. микросекунд.

EDIT: Исходя из некоторого обсуждения в комментариях, должно быть возможно объединить оба ядра в одно ядро, используя схему двойной буферизации на p_x и p_y: 

$ cat t11.py
import numpy as np
import math
from numba import cuda


@cuda.jit('void(float32[:], float32[:], float32[:], float32[:], float32[:], uint8[:], float32[:], float32[:], float32, uint32)')
def update(p_x, p_y, p_x_new, p_y_new, radii, types, velocities, max_velocities, acceleration, num_creatures):
    i = cuda.grid(1)
    if i < num_creatures:
            velocities[i] = velocities[i] + acceleration
            if velocities[i] > max_velocities[i]:
                velocities[i] = max_velocities[i]
            p_x_new[i] = p_x[i] + (math.cos(1.0) * velocities[i])
            p_y_new[i] = p_y[i] + (math.sin(1.0) * velocities[i])
            for j in range(i, num_creatures):
                delta_x = p_x[j] - p_x[i]
                delta_y = p_y[j] - p_y[i]
                distance_squared = (delta_x * delta_x) + (delta_y * delta_y)
                sum_of_radii = radii[types[i]] + radii[types[i]]
                if distance_squared < sum_of_radii * sum_of_radii:
                    pass


acceleration = .1
creature_radius = 10
spacing = 20
food_radius = 3

max_num_creatures = 1500000
num_creatures = 0
max_num_food = 500
num_food = 0
max_num_entities = max_num_creatures + max_num_food
num_entities = 0
cycles = 2


p_x = np.zeros(max_num_entities, dtype=np.float32)
p_y = np.zeros(max_num_entities, dtype=np.float32)
radii = np.array([creature_radius, creature_radius, food_radius], dtype=np.float32)
types = np.zeros(max_num_entities, dtype=np.uint8)

velocities = np.zeros(max_num_creatures, dtype=np.float32)
max_velocities = np.zeros(max_num_creatures, dtype=np.float32)
# types:
# male - 0
# female - 1
# food - 2
for x in range(1, 80000 // spacing):
    for y in range(1, 6000 // spacing):
        if num_creatures % 2 == 0:
            types[num_creatures] = 0
        else:
            types[num_creatures] = 1
        p_x[num_creatures] = x * spacing
        p_y[num_creatures] = y * spacing
        max_velocities[num_creatures] = 5
        num_creatures += 1


device_p_x = cuda.to_device(p_x)
device_p_y = cuda.to_device(p_y)
device_p_x_new = cuda.to_device(p_x)
device_p_y_new = cuda.to_device(p_y)
device_radii = cuda.to_device(radii)
device_types = cuda.to_device(types)
device_velocities = cuda.to_device(velocities)
device_max_velocities = cuda.to_device(max_velocities)
threadsperblock = 64
blockspergrid = (num_creatures // threadsperblock) + 1
for i in range(cycles):
    if i % 2 == 0:
        update[blockspergrid, threadsperblock](device_p_x, device_p_y, device_p_x_new, device_p_y_new, device_radii, device_types,  device_velocities, device_max_velocities, acceleration, num_creatures)
    else:
        update[blockspergrid, threadsperblock](device_p_x_new, device_p_y_new, device_p_x, device_p_y, device_radii, device_types,  device_velocities, device_max_velocities, acceleration, num_creatures)

print(device_p_x_new.copy_to_host())
print(device_p_x.copy_to_host())
$ python t11.py
[ 20.05402946  20.05402946  20.05402946 ...,   0.           0.           0.        ]
[ 20.1620903  20.1620903  20.1620903 ...,   0.          0.          0.       ]
$

По-прежнему необходимо сохранить цикл вызова ядра в cycles в коде хоста, поскольку нам все еще требуется глобальная синхронизация, предоставляемая вызовом ядра. Но для данного числа cycles это уменьшит вклад накладных расходов вызова ядра.

При использовании этой техники необходимо соблюдать осторожность при выборе cycles, а также при использовании данных из буфера p_x или p_x_new для получения согласованных результатов.

...