8. cvičenie – pokročilá CUDA ⚛️

Úvod

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

1. Atomické inštrukcie

Nižšie je program na výpočet histogramu. Zatiaľ obsahuje naivnú implementáciu, ktorá počíta histogram paralelne, ale bez žiadnych pamäťových ochrán a zámkov. Týmto je výsledný histogram nesprávny. Opravte výpočet pomocou atomických inštrukcií.

Po oprave je výpočet správny, ale použitie atomických inštrukcií spôsobilo veľké súperenie o výstupné pole s histogramom hist. Dokončite implementáciu tretieho kernelu, ktorý počíta histogram nezávisle v zdieľanej pamäti na každom bloku a po výpočte prenesie výsledky do globálneho histogramu.

Program implementujte podľa pokynov v TODO komentároch.

"""Paralelny vypocet histogramu."""

import numpy as np
from numba import cuda
import math
import os

NUM_BINS = 256
THREADS_PER_BLOCK = 256
BLOCKS_PER_GRID = 1024
DATA_SIZE = 1 << 24


@cuda.jit
def histogram_non_atomic(hist, data):
    """Vypocet histogramu bez atomickych instrukcii.

    Neupravujte a neopravujte. Vysledny histogram je nespravny.
    """
    idx = cuda.grid(1)
    if idx < data.shape[0]:
        val = data[idx]
        hist[val] += 1


# TODO 1: Opravte kernel tak, aby na navysenie poctu hodnot v histograme pouzil
#         ATOMICKU instrukciu. Vhodnu metodu najdite v dokumentacii:
#         https://nvidia.github.io/numba-cuda/reference/kernel.html
@cuda.jit
def histogram_global_atomic(hist, data):
    """Vypocet histogramu pomocou atomickych instrukcii.

    Vypocet je spravny, ale atomicky pristup k jednemu polu sposobi velke
    superenie.
    """
    idx = cuda.grid(1)
    if idx < data.shape[0]:
        val = data[idx]
        hist[val] += 1


@cuda.jit
def histogram_shared_atomic(hist, data):
    """Vypocet histogramu pomocou atomickych instrukcii.

    Superenie je nizsie, kazdy blok ma vlastny histogram.
    """
    # Alokacia zdielanej pamate
    # Kazdy blok bude mat vlastny lokalny histogram.
    s_hist = cuda.shared.array(NUM_BINS, dtype=np.int32)

    tx = cuda.threadIdx.x

    # Vynulovanie lokalneho histogramu
    for i in range(tx, NUM_BINS, cuda.blockDim.x):
        s_hist[i] = 0
    cuda.syncthreads()

    idx = cuda.grid(1)
    if idx < data.shape[0]:
        val = data[idx]
        # TODO 2: Pouzite ATOMICKU instrukciu na inkrementaciu hodnoty
        #         v histograme
        s_hist[val] += 1
        # TODO 3: Kernel bol spusteny s mensim poctom vlakien ako je pocet
        #         prvkov v poli. Ktore DALSIE PRVKY spracuje toto vlakno?
    cuda.syncthreads()

    # TODO 4: Preneste lokalne histogramy blokov do globalneho histogramu.


def measure_kernel(kernel, d_hist, d_data, blocks, threads):
    # Spustenie nanecisto, aby bol kernel skompilovany a mohli sme ho casovat
    kernel[blocks, threads](d_hist, d_data)
    cuda.synchronize()

    # Vynulovanie histogramu v pamati GPU po trenovacom spusteni
    d_hist[:] = 0
    cuda.synchronize()

    start_event = cuda.event()
    end_event = cuda.event()

    # Casovanie vykonania kernelu
    start_event.record()
    kernel[blocks, threads](d_hist, d_data)
    end_event.record()

    # Pockame na ukoncenie kernelu
    end_event.synchronize()
    elapsed_time_ms = cuda.event_elapsed_time(start_event, end_event)

    return elapsed_time_ms


if __name__ == "__main__":
    # Aby vo vypisoch neboli varovania o malej velkosti gridu
    os.environ['NUMBA_CUDA_LOW_OCCUPANCY_WARNINGS'] = '0'

    print("=== Paralelny vypocet histogramu ===\n")
    print(f"Data: {DATA_SIZE} prvkov ({DATA_SIZE / 1e6:.1f} MB)\n")

    h_data = np.random.randint(0, 256, DATA_SIZE, dtype=np.int32)
    d_data = cuda.to_device(h_data)
    d_hist = cuda.device_array(NUM_BINS, dtype=np.int32)

    # --- Experiment 1: kernel bez atomickych instrukcii ---
    d_hist[:] = 0
    t_non_atomic = measure_kernel(
        histogram_non_atomic,
        d_hist,
        d_data,
        math.ceil(DATA_SIZE / THREADS_PER_BLOCK),
        THREADS_PER_BLOCK,
    )
    histogram_sum = d_hist.copy_to_host().sum()
    result = 'ok' if histogram_sum == DATA_SIZE else 'error'
    print(f"1. non-atomic:      prvkov = {histogram_sum:} {result}")

    # --- Experiment 2: kernel s atomickymi instrukciami, vysoke superenie ---
    t_global = measure_kernel(
        histogram_global_atomic,
        d_hist,
        d_data,
        math.ceil(DATA_SIZE / THREADS_PER_BLOCK),
        THREADS_PER_BLOCK,
    )
    histogram_sum = d_hist.copy_to_host().sum()
    result = 'ok' if histogram_sum == DATA_SIZE else 'error'
    print(
        f"2. global atomic:   prvkov = {histogram_sum:} {result}  cas = {t_global:.3f} ms"
    )

    # --- Experiment 3: kernel s atomickymi instrukciami, nizke superenie ---
    t_shared = measure_kernel(
        histogram_shared_atomic,
        d_hist,
        d_data,
        BLOCKS_PER_GRID,
        THREADS_PER_BLOCK,
    )
    histogram_sum = d_hist.copy_to_host().sum()
    result = 'ok' if histogram_sum == DATA_SIZE else 'error'
    print(
        f"3. shared atomic:   prvkov = {histogram_sum:} {result}  cas = {t_shared:.3f} ms "
        f"({t_global / t_shared:.1f}x rychlejsie)"
    )

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

"""Paralelny vypocet histogramu."""

import numpy as np

from numba import cuda

import math

import os

NUM_BINS = 256

THREADS_PER_BLOCK = 256

BLOCKS_PER_GRID = 1024

DATA_SIZE = 1 << 24

@cuda.jit

def histogram_non_atomic(hist, data):

"""Vypocet histogramu bez atomickych instrukcii.

Neupravujte a neopravujte. Vysledny histogram je nespravny.

"""

idx = cuda.grid(1)

if idx < data.shape[0]:

val = data[idx]

hist[val] += 1

# TODO 1: Opravte kernel tak, aby na navysenie poctu hodnot v histograme pouzil

# ATOMICKU instrukciu. Vhodnu metodu najdite v dokumentacii:

# https://nvidia.github.io/numba-cuda/reference/kernel.html

@cuda.jit

def histogram_global_atomic(hist, data):

"""Vypocet histogramu pomocou atomickych instrukcii.

Vypocet je spravny, ale atomicky pristup k jednemu polu sposobi velke

superenie.

"""

idx = cuda.grid(1)

if idx < data.shape[0]:

val = data[idx]

hist[val] += 1

@cuda.jit

def histogram_shared_atomic(hist, data):

"""Vypocet histogramu pomocou atomickych instrukcii.

Superenie je nizsie, kazdy blok ma vlastny histogram.

"""

# Alokacia zdielanej pamate

# Kazdy blok bude mat vlastny lokalny histogram.

s_hist = cuda.shared.array(NUM_BINS, dtype=np.int32)

tx = cuda.threadIdx.x

# Vynulovanie lokalneho histogramu

for i in range(tx, NUM_BINS, cuda.blockDim.x):

s_hist[i] = 0

cuda.syncthreads()

idx = cuda.grid(1)

if idx < data.shape[0]:

val = data[idx]

# TODO 2: Pouzite ATOMICKU instrukciu na inkrementaciu hodnoty

# v histograme

s_hist[val] += 1

# TODO 3: Kernel bol spusteny s mensim poctom vlakien ako je pocet

# prvkov v poli. Ktore DALSIE PRVKY spracuje toto vlakno?

cuda.syncthreads()

# TODO 4: Preneste lokalne histogramy blokov do globalneho histogramu.

def measure_kernel(kernel, d_hist, d_data, blocks, threads):

# Spustenie nanecisto, aby bol kernel skompilovany a mohli sme ho casovat

kernel[blocks, threads](d_hist, d_data)

cuda.synchronize()

# Vynulovanie histogramu v pamati GPU po trenovacom spusteni

d_hist[:] = 0

cuda.synchronize()

start_event = cuda.event()

end_event = cuda.event()

# Casovanie vykonania kernelu

start_event.record()

kernel[blocks, threads](d_hist, d_data)

end_event.record()

# Pockame na ukoncenie kernelu

end_event.synchronize()

elapsed_time_ms = cuda.event_elapsed_time(start_event, end_event)

return elapsed_time_ms

if __name__ == "__main__":

# Aby vo vypisoch neboli varovania o malej velkosti gridu

os.environ['NUMBA_CUDA_LOW_OCCUPANCY_WARNINGS'] = '0'

print("=== Paralelny vypocet histogramu ===\n")

print(f"Data: {DATA_SIZE} prvkov ({DATA_SIZE / 1e6:.1f} MB)\n")

h_data = np.random.randint(0, 256, DATA_SIZE, dtype=np.int32)

d_data = cuda.to_device(h_data)

d_hist = cuda.device_array(NUM_BINS, dtype=np.int32)

# --- Experiment 1: kernel bez atomickych instrukcii ---

d_hist[:] = 0

t_non_atomic = measure_kernel(

histogram_non_atomic,

d_hist,

d_data,

math.ceil(DATA_SIZE / THREADS_PER_BLOCK),

THREADS_PER_BLOCK,

)

histogram_sum = d_hist.copy_to_host().sum()

result = 'ok' if histogram_sum == DATA_SIZE else 'error'

print(f"1. non-atomic: prvkov = {histogram_sum:} {result}")

# --- Experiment 2: kernel s atomickymi instrukciami, vysoke superenie ---

t_global = measure_kernel(

histogram_global_atomic,

d_hist,

d_data,

math.ceil(DATA_SIZE / THREADS_PER_BLOCK),

THREADS_PER_BLOCK,

)

histogram_sum = d_hist.copy_to_host().sum()

result = 'ok' if histogram_sum == DATA_SIZE else 'error'

print(

f"2. global atomic: prvkov = {histogram_sum:} {result} cas = {t_global:.3f} ms"

)

# --- Experiment 3: kernel s atomickymi instrukciami, nizke superenie ---

t_shared = measure_kernel(

histogram_shared_atomic,

d_hist,

d_data,

BLOCKS_PER_GRID,

THREADS_PER_BLOCK,

)

histogram_sum = d_hist.copy_to_host().sum()

result = 'ok' if histogram_sum == DATA_SIZE else 'error'

print(

f"3. shared atomic: prvkov = {histogram_sum:} {result} cas = {t_shared:.3f} ms "

f"({t_global / t_shared:.1f}x rychlejsie)"

)

2. Redukcia

Dokončite implementáciu jednoduchej redukcie. Uvedomte si, že záverečný súčet v kerneli znižuje efektivitu výpočtu, nakoľko kvôli konzistencii dát musí byť serializovaný. Ako by vyzerala efektívnejšia implementácia?

Vytvorte redukčný kernel pomocou dekorátora, ktorý poskytuje Numba. Porovnajte dobu výpočtu oboch verzií redukcie na GPU a výpočtu na CPU pomocou Numpy.

import numpy as np
from numba import cuda
import math
import time
import os

TPB = 256


@cuda.jit
def custom_reduce_kernel(d_in, d_out):
    # Alokacia zdielanej pamate pre blok
    sdata = cuda.shared.array(TPB, dtype=np.float32)

    tid = cuda.threadIdx.x
    pos = cuda.grid(1)

    # TODO 1: Nacitanie dat do zdielanej pamate:
    #         Skopirujte cislo zo vstupneho pola do zdielanej pamate.
    #         1. Ktory prvok skopirujete?
    #         2. Do akeho miesta v zdielanej pamati prvok ulozite?
    #         3. Ak je globalna pozicia vlakna mimo pola, nahrajte do zdielanej
    #            pamate neutralny prvok vzhladom na redukciu.





    # Pockame, kym vsetky vlakna v bloku nacitaju svoje data
    cuda.syncthreads()

    # TODO 2: Implementujte redukciu v ramci bloku.
    #         V prvej faze kazde vlakno z dolnej polovice vlakien v bloku
    #         k svojmu cislu pripocita cislo z vlakna o krok (`stride`) vyssie.
    #         V druhej faze to uz robi iba polovica z tychto vlakien, atd.
    #         1. Aka je pociatocna velkost kroku?
    #         2. Pred zaciatkom dalsej fazy sa musite uistit, ze vlakna
    #            z predchdzajucej fazy maju vypocet hotovy.
    stride = 0
    while stride > 0:
        pass




    # TODO 3: Sucet lokalnych suctov:
    #         Jedno vlakno z kazdeho bloku zapise lokalny sucet do vystupneho
    #         pola `d_out`.
    #         Dajte pozor na konzistenciu dat!




# TODO 4: Implementujte redukcny kernel s operaciou sucet.
#         Metoda `numba_reduce` implementuje binarnu operaciu sucet. Na
#         vytvorenie redukcneho kernelu pouzite dekorator `cuda.reduce`.
#         https://nvidia.github.io/numba-cuda/user/reduction.html
def numba_reduce():
    pass


def run_reduction_benchmark():
    # Aby vo vypisoch neboli varovania o malej velkosti gridu
    os.environ['NUMBA_CUDA_LOW_OCCUPANCY_WARNINGS'] = '0'

    # 32 Milionov prvkov (128 MB dat pri 32 bit. floate)
    N = 1024 * 1024 * 32

    print(f"Inicializacia pola s velkostou {N} prvkov...")
    h_in = np.random.rand(N).astype(np.float32)

    # Referencny vypocet na CPU (NumPy)
    start_cpu = time.perf_counter()
    expected_sum = np.sum(h_in)
    time_cpu = time.perf_counter() - start_cpu
    print(
        f"[  Redukcia na CPU   ] cas: {time_cpu:.5f} s, vysledok: {expected_sum}"
    )

    # Prenos dat na GPU
    d_in = cuda.to_device(h_in)

    # Vystupne pole o velkosti 1 (bude obsahovat vysledok redukcie)
    h_out = np.zeros(1, dtype=np.float32)
    d_out = cuda.to_device(h_out)

    blocks = math.ceil(d_in.size / TPB)

    print("\nSpustam GPU kernely...")

    # --- TEST 1: vlastna redukcia ---
    # Prvy beh nanecisto
    custom_reduce_kernel[blocks, TPB](d_in, d_out)
    cuda.synchronize()

    # Vystupne pole zresetujeme
    d_out.copy_to_device(h_out)

    start = time.perf_counter()
    custom_reduce_kernel[blocks, TPB](d_in, d_out)
    cuda.synchronize()
    time_custom = time.perf_counter() - start

    result_custom = d_out.copy_to_host()[0]
    print(
        f"[  Vlastna redukcia  ] cas: {time_custom:.5f} s, vysledok: {result_custom}"
    )

    # --- TEST 2: redukcia dekoratorom ---
    # Prvy beh nanecisto
    _ = numba_reduce(d_in)

    start = time.perf_counter()
    result_builtin = numba_reduce(d_in)
    time_builtin = time.perf_counter() - start
    print(
        f"[Redukcia dekoratorom] cas: {time_builtin:.5f} s, vysledok: {result_builtin}"
    )

    # --- VYHODNOTENIE ---
    print("\n--- POROVNANIE VYKONU ---")
    if time_builtin < time_custom:
        print(
            f"Redukcia dekoratorom je {time_custom / time_builtin:.2f}x RYCHLEJSIA ako vlastna implementacia."
        )
    else:
        print(
            f"Vlastna redukcia je {time_builtin / time_custom:.2f}x RYCHLEJSIA (toto sme necakali, asi bude v implementacii nejaka chyba)."
        )

    # Kontrola spravnosti
    is_custom_ok = np.isclose(result_custom, expected_sum, rtol=1e-4)
    is_builtin_ok = np.isclose(result_builtin, expected_sum, rtol=1e-4)

    print("\n--- KONTROLA SPRAVNOSTI ---")
    print(f"Vlastna redukcia     {'OK' if is_custom_ok else 'CHYBA'}")
    print(f"Redukcia dekoratorom {'OK' if is_builtin_ok else 'CHYBA'}")


if __name__ == "__main__":
    run_reduction_benchmark()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

import numpy as np

from numba import cuda

import math

import time

import os

TPB = 256

@cuda.jit

def custom_reduce_kernel(d_in, d_out):

# Alokacia zdielanej pamate pre blok

sdata = cuda.shared.array(TPB, dtype=np.float32)

tid = cuda.threadIdx.x

pos = cuda.grid(1)

# TODO 1: Nacitanie dat do zdielanej pamate:

# Skopirujte cislo zo vstupneho pola do zdielanej pamate.

# 1. Ktory prvok skopirujete?

# 2. Do akeho miesta v zdielanej pamati prvok ulozite?

# 3. Ak je globalna pozicia vlakna mimo pola, nahrajte do zdielanej

# pamate neutralny prvok vzhladom na redukciu.

# Pockame, kym vsetky vlakna v bloku nacitaju svoje data

cuda.syncthreads()

# TODO 2: Implementujte redukciu v ramci bloku.

# V prvej faze kazde vlakno z dolnej polovice vlakien v bloku

# k svojmu cislu pripocita cislo z vlakna o krok (`stride`) vyssie.

# V druhej faze to uz robi iba polovica z tychto vlakien, atd.

# 1. Aka je pociatocna velkost kroku?

# 2. Pred zaciatkom dalsej fazy sa musite uistit, ze vlakna

# z predchdzajucej fazy maju vypocet hotovy.

stride = 0

while stride > 0:

pass

# TODO 3: Sucet lokalnych suctov:

# Jedno vlakno z kazdeho bloku zapise lokalny sucet do vystupneho

# pola `d_out`.

# Dajte pozor na konzistenciu dat!

# TODO 4: Implementujte redukcny kernel s operaciou sucet.

# Metoda `numba_reduce` implementuje binarnu operaciu sucet. Na

# vytvorenie redukcneho kernelu pouzite dekorator `cuda.reduce`.

# https://nvidia.github.io/numba-cuda/user/reduction.html

def numba_reduce():

pass

def run_reduction_benchmark():

# Aby vo vypisoch neboli varovania o malej velkosti gridu

os.environ['NUMBA_CUDA_LOW_OCCUPANCY_WARNINGS'] = '0'

# 32 Milionov prvkov (128 MB dat pri 32 bit. floate)

N = 1024 * 1024 * 32

print(f"Inicializacia pola s velkostou {N} prvkov...")

h_in = np.random.rand(N).astype(np.float32)

# Referencny vypocet na CPU (NumPy)

start_cpu = time.perf_counter()

expected_sum = np.sum(h_in)

time_cpu = time.perf_counter() - start_cpu

print(

f"[ Redukcia na CPU ] cas: {time_cpu:.5f} s, vysledok: {expected_sum}"

)

# Prenos dat na GPU

d_in = cuda.to_device(h_in)

# Vystupne pole o velkosti 1 (bude obsahovat vysledok redukcie)

h_out = np.zeros(1, dtype=np.float32)

d_out = cuda.to_device(h_out)

blocks = math.ceil(d_in.size / TPB)

print("\nSpustam GPU kernely...")

# --- TEST 1: vlastna redukcia ---

# Prvy beh nanecisto

custom_reduce_kernel[blocks, TPB](d_in, d_out)

cuda.synchronize()

# Vystupne pole zresetujeme

d_out.copy_to_device(h_out)

start = time.perf_counter()

custom_reduce_kernel[blocks, TPB](d_in, d_out)

cuda.synchronize()

time_custom = time.perf_counter() - start

result_custom = d_out.copy_to_host()[0]

print(

f"[ Vlastna redukcia ] cas: {time_custom:.5f} s, vysledok: {result_custom}"

)

# --- TEST 2: redukcia dekoratorom ---

# Prvy beh nanecisto

_ = numba_reduce(d_in)

start = time.perf_counter()

result_builtin = numba_reduce(d_in)

time_builtin = time.perf_counter() - start

print(

f"[Redukcia dekoratorom] cas: {time_builtin:.5f} s, vysledok: {result_builtin}"

)

# --- VYHODNOTENIE ---

print("\n--- POROVNANIE VYKONU ---")

if time_builtin < time_custom:

print(

f"Redukcia dekoratorom je {time_custom / time_builtin:.2f}x RYCHLEJSIA ako vlastna implementacia."

)

else:

print(

f"Vlastna redukcia je {time_builtin / time_custom:.2f}x RYCHLEJSIA (toto sme necakali, asi bude v implementacii nejaka chyba)."

)

# Kontrola spravnosti

is_custom_ok = np.isclose(result_custom, expected_sum, rtol=1e-4)

is_builtin_ok = np.isclose(result_builtin, expected_sum, rtol=1e-4)

print("\n--- KONTROLA SPRAVNOSTI ---")

print(f"Vlastna redukcia {'OK' if is_custom_ok else 'CHYBA'}")

print(f"Redukcia dekoratorom {'OK' if is_builtin_ok else 'CHYBA'}")

if __name__ == "__main__":

run_reduction_benchmark()

3. Redukcia vo warpe pomocou shuffle inštrukcie

Nižšie je jednoduchá implementácia násobenia matice vektorom. Výpočet paralelizujeme po riadkoch, pričom každý blok vlákien má na starosti vynásobenie jedného riadka matice. Počet vlákien na blok bol zvolený schválne – je rovný počtu vlákien vo warpe (32). Využite špeciálnu inštrukciu shfl_down_sync, ktorá umožňuje vykonať operáciu súčtu párovo na niekoľkých vláknach warpu paralelne. Implementujte pomocou nej redukciu čiastočných súčtov val každého vlákna vo warpe. Po redukcii by mal byť konečný súčet uložený vo vlákne 0 každého bloku.

"""Vypocet matica ⨯ vektor.

Zalozene na: Maharshi Pandya, https://github.com/Maharshi-Pandya/cudacodes
Do Numba CUDA prelozil Google, Gemini 3.
"""

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


@cuda.jit(device=True)
def warp_reduce_sum(val):
    """TODO: Kazde vlakno vo warpe ma svoju hodnotu premennej `val`. Urobte na
    tejto premennej redukciu suctom pomocou instrukcie shfl_down_sync.

    Redukciu je potrebne robit v cykle niekokymi priebehmi. Kolko ich bude?
    """
    return val


@cuda.jit
def coalesced_warp_sgmev_kernel(matd, vecd, resd):
    # Kazdy blok spracuje jeden riadok matice
    row = cuda.blockIdx.x

    if row < matd.shape[0]:
        tid = cuda.threadIdx.x
        N = matd.shape[1]

        # Kazde vlakno vypocita lokalny ciastocny sucet
        partial_sum = 0.0
        for col in range(tid, N, cuda.blockDim.x):
            # Nasobenie riadka matice vektorom
            partial_sum += matd[row, col] * vecd[col]

        # Redukcia napriec warpom
        total_sum = warp_reduce_sum(partial_sum)

        # Iba jedno vlakno z bloku/warpu (v tomto pripade bol kernel spusteny
        # tak, aby warp = blok) zapise vysledok
        if tid == 0:
            resd[row] = total_sum


def run_benchmark(M=4096, N=8192):
    print(f"------------ Spustam test pre M = {M}, N = {N} ---------------")

    h_mat = np.random.normal(0, 1, (M, N)).astype(np.float32)
    h_vec = np.random.normal(0, 1, N).astype(np.float32)

    d_mat = cuda.to_device(h_mat)
    d_vec = cuda.to_device(h_vec)
    d_res = cuda.device_array(M, dtype=np.float32)
    cuda.synchronize()

    # Tento kernel vyzaduje, aby bol blok velkosti warpu (32)
    threads_per_block = 32
    blocks_per_grid = M

    # Spustenie kernelu nanecisto kvoli meraniu casu
    coalesced_warp_sgmev_kernel[blocks_per_grid, threads_per_block](
        d_mat, d_vec, d_res
    )
    cuda.synchronize()

    start_evt = cuda.event()
    end_evt = cuda.event()

    start_evt.record()
    coalesced_warp_sgmev_kernel[blocks_per_grid, threads_per_block](
        d_mat, d_vec, d_res
    )
    end_evt.record()
    end_evt.synchronize()

    ms = cuda.event_elapsed_time(start_evt, end_evt)

    h_res_gpu = d_res.copy_to_host()
    h_res_cpu = h_mat @ h_vec

    # Citlivost na chybu 0.001
    if np.allclose(h_res_gpu, h_res_cpu, atol=1e-3):
        print("Kotrola spravnosti: OK")
    else:
        print("Kotrola spravnosti: CHYBA")
        diff = np.abs(h_res_gpu - h_res_cpu)
        print(f"Maximalny rozdiel: {np.max(diff)}")


if __name__ == "__main__":
    run_benchmark()

"""Vypocet matica ⨯ vektor.

Zalozene na: Maharshi Pandya, https://github.com/Maharshi-Pandya/cudacodes

Do Numba CUDA prelozil Google, Gemini 3.

"""

import numpy as np

from numba import cuda

import time

import math

@cuda.jit(device=True)

def warp_reduce_sum(val):

"""TODO: Kazde vlakno vo warpe ma svoju hodnotu premennej `val`. Urobte na

tejto premennej redukciu suctom pomocou instrukcie shfl_down_sync.

Redukciu je potrebne robit v cykle niekokymi priebehmi. Kolko ich bude?

"""

return val

@cuda.jit

def coalesced_warp_sgmev_kernel(matd, vecd, resd):

# Kazdy blok spracuje jeden riadok matice

row = cuda.blockIdx.x

if row < matd.shape[0]:

tid = cuda.threadIdx.x

N = matd.shape[1]

# Kazde vlakno vypocita lokalny ciastocny sucet

partial_sum = 0.0

for col in range(tid, N, cuda.blockDim.x):

# Nasobenie riadka matice vektorom

partial_sum += matd[row, col] * vecd[col]

# Redukcia napriec warpom

total_sum = warp_reduce_sum(partial_sum)

# Iba jedno vlakno z bloku/warpu (v tomto pripade bol kernel spusteny

# tak, aby warp = blok) zapise vysledok

if tid == 0:

resd[row] = total_sum

def run_benchmark(M=4096, N=8192):

print(f"------------ Spustam test pre M = {M}, N = {N} ---------------")

h_mat = np.random.normal(0, 1, (M, N)).astype(np.float32)

h_vec = np.random.normal(0, 1, N).astype(np.float32)

d_mat = cuda.to_device(h_mat)

d_vec = cuda.to_device(h_vec)

d_res = cuda.device_array(M, dtype=np.float32)

cuda.synchronize()

# Tento kernel vyzaduje, aby bol blok velkosti warpu (32)

threads_per_block = 32

blocks_per_grid = M

# Spustenie kernelu nanecisto kvoli meraniu casu

coalesced_warp_sgmev_kernel[blocks_per_grid, threads_per_block](

d_mat, d_vec, d_res

)

cuda.synchronize()

start_evt = cuda.event()

end_evt = cuda.event()

start_evt.record()

coalesced_warp_sgmev_kernel[blocks_per_grid, threads_per_block](

d_mat, d_vec, d_res

)

end_evt.record()

end_evt.synchronize()

ms = cuda.event_elapsed_time(start_evt, end_evt)

h_res_gpu = d_res.copy_to_host()

h_res_cpu = h_mat @ h_vec

# Citlivost na chybu 0.001

if np.allclose(h_res_gpu, h_res_cpu, atol=1e-3):

print("Kotrola spravnosti: OK")

else:

print("Kotrola spravnosti: CHYBA")

diff = np.abs(h_res_gpu - h_res_cpu)

print(f"Maximalny rozdiel: {np.max(diff)}")

if __name__ == "__main__":

run_benchmark()