Wprowadzenie do problemu Redukcji Wektorów

Problem Redukcji Wektorów (ang. Vector Reduction) polega na obliczeniu sumy wszystkich elementów wektora (tablicy jednowymiarowej) i redukcji jego postaci do wartości skalarnej. Formalnie, dla danego wektora wejściowego A o długości N, zadanie sprowadza się do obliczenia wartości:

$$ S = \sum_{i = 0}^{N-1} A[i] $$

gdzie A[i] to i-ty element wektora wejściowego.

Metody Redukcji Wektorów

Redukcja sąsiednich elementów

Binarna redukcja wektora poprzez sumowanie elementów na sąsiednich pozycjach (Pairwise Summation Reduction) to technika pozwalająca na równoległe obliczanie sumy elementów wektora poprzez stopniowe redukowanie liczby operacji w każdej iteracji. W tej metodzie elementy wektora są sumowane parami, co stopniowo zmniejsza liczbę elementów do zsumowania o połowę w każdej iteracji.

image.png

Idea algorytmu

Dla wektora wejściowego A o długości N następują kolejne kroki:

  1. Pierwsza iteracja redukująca liczbę elementów do N / 2
    1. A[0] = A[0]+A[1],A[2] = A[2]+A[3], …, A[N-2] = A[N−2]+A[N−1]
    2. N = N / 2, przy zachowaniu pierwotnej długości wektora w pamięci
  2. Druga iteracja redukująca liczbę elementów do N / 4
    1. A[0] = A[0]+A[1],A[2] = A[2]+A[3], …, A[N-2] = A[N−2]+A[N−1]
    2. N = N / 2, przy zachowaniu pierwotnej długości wektora w pamięci
  3. Proces powtarzany jest do momentu aż wektor pozostanie wartością skalarną.

Implementacja

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <iostream>

// kernel wywola sie dla roznych (8) poziomow redukcji wektora, 
// kernel dla kazdego elementu pierwotnego wektora, 
// w nastepnych krokach czesc watkow jest IDLE
// mniej operacji redukcji jest w kazdym kolejnym zagniezdzeniu kernel'a
__global__ void ReduceInPlace(float* input, int n) {
    int tid = threadIdx.x;
		// index pierwszego elementu tablicy do zsumowania
    int index = blockIdx.x * blockDim.x + threadIdx.x;

		// petla ktora najpierw wykona sume dla sasiednich elementow, 
		// potem dla co 2, co 4 itd
    for (int stride = 1; stride < blockDim.x; stride *= 2) {
				// czeka az poprzednie watki zakoncza prace - race condition
        __syncthreads();

				// jezeli jest na pozycji co 2, 4, 8 czyli zgodnej ze stride 
				// i nie sa poza zakresem dlugosci
        if (tid % (2 * stride) == 0 && index + stride < n) { 
						// to w pierwszej pozycji wektora wejsciowego wykonaj sume
            input[index] += input[index + stride];
        }
    }

    if (tid == 0) { // jezeli to pierwszy watek w bloku
        // na pozycji wektora wejsciwego rownej indksowi bloku przypisz pierwsza 
        // pozycje danych przetwarzana przez blok
        // czyli wtedy na pierwszych 4096 pozycjach tego wektora beda sumy czastkowe
        input[blockIdx.x] = input[blockIdx.x * blockDim.x];
    }
}

float cpuReduce(float* input, int size) {
    float sum = 0.0f;
    for (int i = 0; i < size; ++i) {
        sum += input[i];
    }
    return sum;
}

int main() {
    int n = 1024 * 1024;
    size_t bytes = n * sizeof(float);

    float* h_input = new float[n];
    float* d_input;

    for (int i = 0; i < n; i++) {
        h_input[i] = static_cast<float>(i + 1);
    }

    cudaMalloc(&d_input, bytes);

    cudaMemcpy(d_input, h_input, bytes, cudaMemcpyHostToDevice);

    int block_size = 256;
    int grid_size = (n + block_size - 1) / block_size;
    float total_sum = cpuReduce(h_input, 1024 * 1024);
    std::cout << "Final sum (CPU): " << total_sum << std::endl;
    
    ReduceInPlace <<<grid_size, block_size>>> (d_input, n);
    cudaDeviceSynchronize();
    ReduceInPlace <<<16, block_size>>> (d_input, 4096);
    cudaDeviceSynchronize();

    ReduceInPlace <<<1, block_size >>> (d_input, 16);
    cudaDeviceSynchronize();
    cudaMemcpy(h_input, d_input, sizeof(float), cudaMemcpyDeviceToHost);

    std::cout << "Final sum (GPU): " << h_input[0] << std::endl;

    cudaFree(d_input);
    delete[] h_input;

    return 0;
}

Metoda sumowania elementów znajdujących się po przeciwnych stronach bloku

Metoda sumowania elementów znajdujących się po przeciwnych stronach bloku (In-Place Opposite-Side Reduction), to metoda równoległego sumowania elementów wektora, w której w każdej iteracji sumowane są elementy znajdujące się po przeciwnych stronach bloku wątków. W przeciwieństwie do klasycznej redukcji sąsiedniej, gdzie sumowane są elementy leżące obok siebie, tutaj operacje odbywają się na elementach oddalonych o rosnący współczynnik przesunięcia (stride). Dzięki takiemu podejściu zmniejsza się liczba operacji wymagających synchronizacji wątków, a sumowanie odbywa się, podobnie jak w metodzie binarnej - w miejscu. Metoda ta eliminuje konieczność kopiowania danych między różnymi tablicami. W kolejnych iteracjach rozmiar zbioru aktywnych elementów zmniejsza się dwukrotnie, aż do uzyskania końcowego wyniku. Dzięki zastosowaniu tego podejścia znacząco zmniejsza się liczba dostępów do pamięci globalnej.

image.png