7. cvičenie – základy CUDA 🔢

Úvod

Zdrojové súbory na cvičenie si stiahnite spoločne ako archív.

1. Veľkosť gridu a veľkosť bloku

Dokončite program nižšie podľa TODO komentárov. Aká je najlepšia konfigurácia gridu? Akej konfigurácii sa vyvarovať vzhľadom na výkon programu?

import numpy as np
from numba import cuda
import math


# TODO: Urobte z funkcie kernel.
def math_heavy_kernel(data_in, data_out):
    # TODO: Zistite globalnu poziciu vlakna v gride
    pos = 0

    if pos < data_in.size:
        val = data_in[pos]
        for i in range(100):
            val = math.sqrt(val * 1.000001) + math.sin(val)
        data_out[pos] = val


def run_benchmark():
    device = cuda.get_current_device()
    cc = device.compute_capability
    print(f"Zariadenie  : {device.name.decode()}")
    print(f"Compute cap : {cc[0]}.{cc[1]}")

    # Priprava dat
    n = 10**7
    h_in = np.random.uniform(1, 10, n).astype(np.float32)
    h_out = np.empty(n, dtype=np.float32)

    # TODO: Nakopirujte polia `h_in` a `h_out` do globalnej pamati grafickej
    # karty.
    d_in = h_in
    d_out = h_out

    # TODO: Pripravte niekolko variant poctu vlakien na jeden blok. Kolko cisel
    # vlozite do zoznamu, tolkokrat sa vykona experiment s inym poctom vlakien.
    # 1. Ake su vhodne pocty vlakien na jeden blok?
    # 2. Ake su nevhodne pocty vlakien na jeden blok? Vyskusajte aj tie a
    #    pozorujte vysledky.
    tpb_options = []

    N_RUNS = 5

    # Prve spustenie kernelu, aby sa skompiloval a vedeli sme neskor odmerat
    # dlzku jeho vykonavania.
    _bpg_warmup = math.ceil(n / 256)
    # TODO: Zavolajte kernel `math_heavy_kernel` s argumentmi `d_in` a `d_out`
    # s `_bpg_warmup` blokmi a 256 vlaknami na jeden blok.
    math_heavy_kernel()
    cuda.synchronize()

    # Hlavicka tabulky vypisu casov
    print(f"\n{' '*62}")
    print(
        f"{'Threads/Block':<15}   {'Blokov/Grid':<13}   {'Priem. [ms]':<13}   Min [ms]"
    )
    print(f"{' '*62}")

    timing_results = {}

    # Pokus opakujeme pre niekolko konfiguracii
    for tpb in tpb_options:
        bpg = math.ceil(n / tpb)
        times = []

        # A pre kazdu konfiguraciu niekolkokrat, aby sme ziskali presnejsi cas
        # behu bez externych zavislosti
        for _ in range(N_RUNS):
            t_start = cuda.event(timing=True)
            t_end = cuda.event(timing=True)

            t_start.record()

            # TODO: Zavolajte kernel `math_heavy_kernel` s argumentami `d_in` a
            # `d_out` s `bpg` blokmi a `tbp` vlaknami na jeden blok.
            math_heavy_kernel()

            t_end.record()
            t_end.synchronize()

            times.append(cuda.event_elapsed_time(t_start, t_end))

        avg_ms = float(np.mean(times))
        min_ms = float(np.min(times))
        timing_results[tpb] = avg_ms
        print(f"{tpb:<15}   {bpg:<13}   {avg_ms:<13.3f}   {min_ms:.3f}")


if __name__ == "__main__":
    run_benchmark()

import numpy as np

from numba import cuda

import math

# TODO: Urobte z funkcie kernel.

def math_heavy_kernel(data_in, data_out):

# TODO: Zistite globalnu poziciu vlakna v gride

pos = 0

if pos < data_in.size:

val = data_in[pos]

for i in range(100):

val = math.sqrt(val * 1.000001) + math.sin(val)

data_out[pos] = val

def run_benchmark():

device = cuda.get_current_device()

cc = device.compute_capability

print(f"Zariadenie : {device.name.decode()}")

print(f"Compute cap : {cc[0]}.{cc[1]}")

# Priprava dat

n = 10**7

h_in = np.random.uniform(1, 10, n).astype(np.float32)

h_out = np.empty(n, dtype=np.float32)

# TODO: Nakopirujte polia `h_in` a `h_out` do globalnej pamati grafickej

# karty.

d_in = h_in

d_out = h_out

# TODO: Pripravte niekolko variant poctu vlakien na jeden blok. Kolko cisel

# vlozite do zoznamu, tolkokrat sa vykona experiment s inym poctom vlakien.

# 1. Ake su vhodne pocty vlakien na jeden blok?

# 2. Ake su nevhodne pocty vlakien na jeden blok? Vyskusajte aj tie a

# pozorujte vysledky.

tpb_options = []

N_RUNS = 5

# Prve spustenie kernelu, aby sa skompiloval a vedeli sme neskor odmerat

# dlzku jeho vykonavania.

_bpg_warmup = math.ceil(n / 256)

# TODO: Zavolajte kernel `math_heavy_kernel` s argumentmi `d_in` a `d_out`

# s `_bpg_warmup` blokmi a 256 vlaknami na jeden blok.

math_heavy_kernel()

cuda.synchronize()

# Hlavicka tabulky vypisu casov

print(f"\n{' '*62}")

print(

f"{'Threads/Block':<15} {'Blokov/Grid':<13} {'Priem. [ms]':<13} Min [ms]"

)

print(f"{' '*62}")

timing_results = {}

# Pokus opakujeme pre niekolko konfiguracii

for tpb in tpb_options:

bpg = math.ceil(n / tpb)

times = []

# A pre kazdu konfiguraciu niekolkokrat, aby sme ziskali presnejsi cas

# behu bez externych zavislosti

for _ in range(N_RUNS):

t_start = cuda.event(timing=True)

t_end = cuda.event(timing=True)

t_start.record()

# TODO: Zavolajte kernel `math_heavy_kernel` s argumentami `d_in` a

# `d_out` s `bpg` blokmi a `tbp` vlaknami na jeden blok.

math_heavy_kernel()

t_end.record()

t_end.synchronize()

times.append(cuda.event_elapsed_time(t_start, t_end))

avg_ms = float(np.mean(times))

min_ms = float(np.min(times))

timing_results[tpb] = avg_ms

print(f"{tpb:<15} {bpg:<13} {avg_ms:<13.3f} {min_ms:.3f}")

if __name__ == "__main__":

run_benchmark()

2. Prúdy

V programe nižšie sa robí výpočet na 200 poliach. Aj keď je tento výpočet vykonávaný na grafickej karte, jednotlivé výpočty a pamäťové presuny idú za sebou sekvenčne v predvolenom prúde (default stream). Upravte program tak, aby bolo každé pole spracované vo svojom samostatnom prúde.

"""Running multiple kernels on the GPU sequentially."""

from numba import cuda
import numpy as np
from time import perf_counter


NUM_ARRAYS = 200
ARRAY_LEN = 2**20


@cuda.jit
def kernel(array):
    """Do some dummy operations on an array.

    :param array: array where the operations will be done. Should be
    left with the same values after the kernel is done.
    """
    thd = cuda.grid(1)
    num_iters = array.size // cuda.blockDim.x
    for j in range(num_iters):
        i = j * cuda.blockDim.x + thd
        for k in range(50):
            array[i] *= 2
            array[i] /= 2


data = []
data_gpu = []
gpu_out = []

for _ in range(NUM_ARRAYS):
    data.append(np.random.randn(ARRAY_LEN).astype('float32'))
    # TODO 1: Vytvorte novy prud pre spracovanie pola vyssie.

t_start = perf_counter()

# TODO 2: Upravte kopirovanie poli z RAM na GPU, vypocet kernelu a kopirovanie
# poli z GPU do RAM tak, aby kazde pole bolo spracovane v SAMOSTATNOM prude.
#
# Porovnajte cas spracovania v default prude a cas spracovania v samostatnych
# prudoch.

for array in data:
    data_gpu.append(cuda.to_device(array))

for array in data_gpu:
    kernel[1, 64](array)

for array in data_gpu:
    gpu_out.append(array.copy_to_host())

t_end = perf_counter()

# Check if the arrays are the same (within a small floating point error)
for gpu_array, array in zip(gpu_out, data):
    assert np.allclose(gpu_array, array)

print(f'Total time: {t_end - t_start:.2f}')

"""Running multiple kernels on the GPU sequentially."""

from numba import cuda

import numpy as np

from time import perf_counter

NUM_ARRAYS = 200

ARRAY_LEN = 2**20

@cuda.jit

def kernel(array):

"""Do some dummy operations on an array.

:param array: array where the operations will be done. Should be

left with the same values after the kernel is done.

"""

thd = cuda.grid(1)

num_iters = array.size // cuda.blockDim.x

for j in range(num_iters):

i = j * cuda.blockDim.x + thd

for k in range(50):

array[i] *= 2

array[i] /= 2

data = []

data_gpu = []

gpu_out = []

for _ in range(NUM_ARRAYS):

data.append(np.random.randn(ARRAY_LEN).astype('float32'))

# TODO 1: Vytvorte novy prud pre spracovanie pola vyssie.

t_start = perf_counter()

# TODO 2: Upravte kopirovanie poli z RAM na GPU, vypocet kernelu a kopirovanie

# poli z GPU do RAM tak, aby kazde pole bolo spracovane v SAMOSTATNOM prude.

# Porovnajte cas spracovania v default prude a cas spracovania v samostatnych

# prudoch.

for array in data:

data_gpu.append(cuda.to_device(array))

for array in data_gpu:

kernel[1, 64](array)

for array in data_gpu:

gpu_out.append(array.copy_to_host())

t_end = perf_counter()

# Check if the arrays are the same (within a small floating point error)

for gpu_array, array in zip(gpu_out, data):

assert np.allclose(gpu_array, array)

print(f'Total time: {t_end - t_start:.2f}')

3. Súčin matíc so zdieľanou pamäťou

Upravte program súčinu matíc z prednášky tak, aby využíval zdieľanú pamäť. Cieľom je rozdeliť veľkú maticu na dlažice, skopírovať dáta z matíc do zdieľanej pamäte a medzivýsledok počítať zo zdieľanej pamäte v rámci jedného bloku vlákien.

"""Nasobenie matic.

Zdroj: https://nyu-cds.github.io/python-numba/05-cuda/
"""

from numba import cuda
import numpy
import math


# Velkost bloku pre matice bude 16x16
TPB = 16


@cuda.jit
def matmul(A, B, C):
    """Vypocitaj sucin matic na GPU.

    :param a: Lava matica.
    :param b: Prava matica.
    :param c: Vysledok.
    """
    # TODO 1: ALOKACIA ZDIELANEJ PAMATE
    # Alokujte dve polia sA a sB s velkostou TPB x TPB v ZDIELANEJ pamati typu
    # float32.
    sA = None
    sB = None

    tx = cuda.threadIdx.x
    ty = cuda.threadIdx.y
    row, col = cuda.grid(2)
    if row >= C.shape[0] and col >= C.shape[1]:
        return

    tmp = 0.0

    # TODO 2: Upravte cyklus tak, aby index i oznacoval "dlazdicu" matice.
    for i in range(A.shape[1]):
        # TODO 3: NACITANIE DAT DO DLAZDICE
        # Kazde vlakno nakopiruje jeden prvok z matice A do matice sA a prvok
        # z B do matice sB. Odkial sa kopiruje je urcene BLOKOM ("dlazdicou"),
        # v ktorom sa vlakno nachadza.

        # TODO 4: Bariera, tu sa vsetky vlakna pockaju, kym nebude zdielana pamat
        # naplnena. Pouzite metodu modulu cuda.

        # TODO 5: Upravte vypocet tak, aby ste vypocitali ciastocny sucin
        # sA a sB. Riadok/stlpec matice urcite podla pozicie v bloku.
        tmp += A[row, i] * B[i, col]

        # TODO 6: Bariera, tu sa vsetky vlakna pockaju, kym nebudu matice
        # v zdielanej pamati spracovane.

    C[row, col] = tmp


A = numpy.full((TPB * 2, TPB * 3), 3, numpy.float32)
B = numpy.full((TPB * 3, TPB * 1), 4, numpy.float32)

A_global_mem = cuda.to_device(A)
B_global_mem = cuda.to_device(B)

C_global_mem = cuda.device_array((TPB * 2, TPB * 1))

threadsperblock = (TPB, TPB)
blockspergrid_x = int(math.ceil(A.shape[0] / threadsperblock[0]))
blockspergrid_y = int(math.ceil(B.shape[1] / threadsperblock[1]))
blockspergrid = (blockspergrid_x, blockspergrid_y)

print(f"Spustam matmul na dlazdiciach {TPB}x{TPB}...")
matmul[blockspergrid, threadsperblock](
    A_global_mem, B_global_mem, C_global_mem
)

C = C_global_mem.copy_to_host()

print(C)

"""Nasobenie matic.

Zdroj: https://nyu-cds.github.io/python-numba/05-cuda/

"""

from numba import cuda

import numpy

import math

# Velkost bloku pre matice bude 16x16

TPB = 16

@cuda.jit

def matmul(A, B, C):

"""Vypocitaj sucin matic na GPU.

:param a: Lava matica.

:param b: Prava matica.

:param c: Vysledok.

"""

# TODO 1: ALOKACIA ZDIELANEJ PAMATE

# Alokujte dve polia sA a sB s velkostou TPB x TPB v ZDIELANEJ pamati typu

# float32.

sA = None

sB = None

tx = cuda.threadIdx.x

ty = cuda.threadIdx.y

row, col = cuda.grid(2)

if row >= C.shape[0] and col >= C.shape[1]:

return

tmp = 0.0

# TODO 2: Upravte cyklus tak, aby index i oznacoval "dlazdicu" matice.

for i in range(A.shape[1]):

# TODO 3: NACITANIE DAT DO DLAZDICE

# Kazde vlakno nakopiruje jeden prvok z matice A do matice sA a prvok

# z B do matice sB. Odkial sa kopiruje je urcene BLOKOM ("dlazdicou"),

# v ktorom sa vlakno nachadza.

# TODO 4: Bariera, tu sa vsetky vlakna pockaju, kym nebude zdielana pamat

# naplnena. Pouzite metodu modulu cuda.

# TODO 5: Upravte vypocet tak, aby ste vypocitali ciastocny sucin

# sA a sB. Riadok/stlpec matice urcite podla pozicie v bloku.

tmp += A[row, i] * B[i, col]

# TODO 6: Bariera, tu sa vsetky vlakna pockaju, kym nebudu matice

# v zdielanej pamati spracovane.

C[row, col] = tmp

A = numpy.full((TPB * 2, TPB * 3), 3, numpy.float32)

B = numpy.full((TPB * 3, TPB * 1), 4, numpy.float32)

A_global_mem = cuda.to_device(A)

B_global_mem = cuda.to_device(B)

C_global_mem = cuda.device_array((TPB * 2, TPB * 1))

threadsperblock = (TPB, TPB)

blockspergrid_x = int(math.ceil(A.shape[0] / threadsperblock[0]))

blockspergrid_y = int(math.ceil(B.shape[1] / threadsperblock[1]))

blockspergrid = (blockspergrid_x, blockspergrid_y)

print(f"Spustam matmul na dlazdiciach {TPB}x{TPB}...")

matmul[blockspergrid, threadsperblock](

A_global_mem, B_global_mem, C_global_mem

)

C = C_global_mem.copy_to_host()

print(C)