Odpowiedź na tytułowe pytanie brzmi: To zależy, zależy od tego co rozumiesz przez słowo „szybko”. W niniejszym artykule postaramy się rozwiązać ten problem na kilka różnych sposobów.
Jeżeli musimy wykonać obliczenia tylko jeden raz/dziennie/tygodniowo/miesięcznie to najlepszym rozwiązaniem będzie użycie narzędzia, które pozwoli na prostą implementację algorytmu, ale niekoniecznie wykona obliczania w najkrótszym możliwym czasie. Może to być sytuacja, w której dysponujemy zapisanymi na dysku danymi pomiarowymi z 10 czujników w postaci wektorów o długości 4194304 i po prostu chcemy je w jakiś sposób zbadać. Do tego celu idealnie sprawdzi się na przykład środowisko MATLAB lub Python – czyli jakiś interpreter i instrukcje wysokiego poziomu. My posłużymy się tym pierwszym (MATLAB R2022a), a w roli maszyny testowej wykorzystamy komputer Intel NUC 8 opisany tutaj, gdyż jest to doskonały reprezentant sprzętu za rozsądną cenę. Podczas wszystkich testów używaliśmy systemu Ubuntu 20.04.
Zaczniemy od zasymulowania wektorów danych pomiarowych:
clear all;
close all;
NELEMENTS = 4194304
var = zeros(10, NELEMENTS);
tic;
var(1,:) = sin([0:1:(NELEMENTS - 1)]) + cos([0:1:(NELEMENTS - 1)]);
var(2,:) = exp(var(1,:)) + exp(-1 * [0:1:(NELEMENTS - 1)]);
var(3,:) = sin(var(2,:)) .* cos(var(1,:)) + var(2,:);
var(4,:) = (var(1,:).^2 + var(3,:).^2).^.5;
var(5,:) = nthroot(var(1,:), 3);
var(6,:) = sin(var(2,:)) + cos(var(1,:));
var(7,:) = exp(var(3,:)) + exp(-1 * var(5,:));
var(8,:) = sin(var(2,:)) .* cos(var(1,:)) + cos(var(4,:)) .* sin(var(3,:));
var(9,:) = (var(3,:).^2 + var(2,:).^2).^.5;
var(10,:) = nthroot(var(4,:), 3);
toc
Następnie obliczymy macierz kowariancji.
tic
cov(var.');
toc
W ten sposób powstała macierz kowariancji o wymiarach 10×10. Czas obliczeń na naszym komputerze wynosi około 130 ms. Trzeba więc przyznać, że biorąc pod uwagę rozmiar danych wejściowych wynik otrzymaliśmy całkiem szybko. Stało się tak dlatego, że MATLAB jako środowisko matematycznie jest dość dobrze zoptymalizowany – wręcz stworzony do obliczania takiej właśnie kowariancji.
Co zrobić gdy to nie wystarczy?
Załóżmy radar wyposażony w 10 kanałów, który generuje bloki rzeczywistych próbek sygnału o długości 4194304, 20 razy na sekundę. Z jakiegoś powodu nasze wspaniałe algorytmy do cyfrowego przetwarzania sygnałów wymagają obliczenia takiej macierzy kowariancji. To oznacza, że do dyspozycji mamy raptem 50 ms, czyli mniej niż połowę czasu jakiego wymaga algorytm w środowisku MATLAB. Trzeba przygotować jakąś szybszą implementację. Optymalizację zaczynamy od przyjrzenia się samemu algorytmowi obliczania macierzy kowariancji:
- Od wszystkich wektorów odjąć średnią
- Obliczyć sumę iloczynów dla każdej pary wektorów
Łatwo zauważyć, że ze względu na przemienność operacji mnożenia, wystarczy obliczyć tylko połowę tej macierzy, gdyż druga będzie zawierała dokładnie te same wartości – taka macierz jest symetryczna. Tylko czy wykorzystanie tego faktu pozwoli wystarczająco przyspieszyć obliczenia? Tym razem wykorzystamy język C++ z opcjami kompilatora std=c++14, O3 oraz ffast-math. Dane testowe wygenerowaliśmy w identyczny sposób jak w środowisku MATLAB (to pozwala na łatwe porównanie wyników w poszukiwaniu błędów).
void generateVectors(float *var[10])
{
unsigned long long n;
b_t();
for (n = 0; n < NELEMENTS; ++n)
{
var[0][n] = sin(static_cast<float>(n)) + cos(static_cast<float>(n));
var[1][n] = exp(var[0][n]) + exp(-1* static_cast<float>(n));
var[2][n] = sin(var[1][n]) * cos(var[0][n]) + var[1][n];
var[3][n] = hypot(var[0][n], var[2][n]);
var[4][n] = cbrt(var[0][n]);
var[5][n] = sin(var[1][n]) + cos(var[0][n]);
var[6][n] = exp(var[2][n]) + exp(-1* var[4][n]);
var[7][n] = sin(var[1][n]) * cos(var[0][n]) + cos(var[3][n]) * sin(var[2][n]);
var[8][n] = hypot(var[2][n], var[1][n]);
var[9][n] = cbrt(var[3][n]);
}
double tgen = e_t(); // stop timing
printf("# GENERATION TIME: %f ms\n", tgen*1000);
}
Alokację pamięci realizujemy przy pomocy aligned_alloc(), przyczyna zostanie opisana w części poświęconej implementacji przy pomocy AVX2. Nie wpływało to na zmierzony czas obliczeń. Nasza najprostsza implementacja kowariancji w języku C++ znajduje się poniżej.
namespace ser{
template <typename T>
static T calcMeanGlobal(T* vec, unsigned long size)
{
double reg = 0;
unsigned i;
for (i = 0; i < size; ++i)
{
reg += vec[i];
}
return reg / (T)size;
}
template <typename T>
static void rmVal(T* vec, unsigned long size, T val)
{
unsigned i;
for (i = 0; i < size; ++i)
{
vec[i] -= val;
}
}
template <typename T>
static T mulReduceGlobal(T* vec1, T* vec2, unsigned long size)
{
double reg = 0;
unsigned i;
for (i = 0; i < size; ++i)
{
reg += vec1[i] * vec2[i];
}
return reg / (T)(size - 1);
}
template <typename T>
void cov(T* vec[], unsigned long size, T out[10][10])
{
int i, j;
for (i = 0; i < 10; ++i)
{
rmVal(vec[i], size, calcMeanGlobal(vec[i], size));
}
for (i = 0; i < 10; i++)
{
for (j = 0; j <= i; j++)
{
out[i][j] = mulReduceGlobal(vec[i], vec[j], size);
}
}
}
}
Użyliśmy szablonu, żeby umożliwić wykorzystanie danych wejściowych w postaci wektorów liczb zmiennoprzecinkowych pojedynczej lub podwójnej precyzji. Warto jednak zwrócić uwagę, że akumulatory (zmienna reg) w funkcjach mulReduceGlobal() i calcMeanGlobal() muszą być typu double, żeby algorytm zwrócił wynik zgodny z tym uzyskanym w środowisku MATLAB (błąd względny < 1e-6) – duże redukcje sprzyjają błędom numerycznym. Niestety pomimo obliczania jedynie połowy macierzy kowariancji czas obliczeń wyniósł 377.058572 ms (czas średni za 100 uruchomień), czyli o 247 ms wolniej niż prosty kod w środowisku MATLAB. Przyglądając się implementacji C++ zauważamy dwie rzeczy. Po pierwsze, odejmując wartość średnią od wektorów wykonujemy N (długość wektora) zbędnych dostępów do pamięci RAM. Można je połączyć z obliczaniem kowariancji, choćby nawet oznaczało to odejmowanie jej za każdym razem dla wszystkich par wektorów – odejmowanie jest o wiele szybsze niż grzebanie w pamięci RAM. Z drugiej strony odejmowanie stałej ma dużo mniejszą złożoność niż suma iloczynów kombinacji par wektorów – zajmuje ułamek czasu obliczeń i dlatego nie będzie to przedmiotem dalszej optymalizacji. Drugą rzeczą wartą uwagi jest to, że wykorzystaliśmy tylko mały kawałek mocy obliczeniowej tylko jednego rdzenia procesora, a dostępne są aż 4.
OpenMP
Do równoległej implementacji algorytmu obliczającego macierz kowariancji wykorzystaliśmy bibliotekę OpenMP. Jest ona bardzo prosta w użyciu i ma ogromne możliwości pod względem drobnoziarnistego zrównoleglenia obliczeń. Kod znajduje się poniżej:
#define CHUNK_SIZE 1024
namespace omp{
template <typename T>
static T calcMeanGlobal(T* vec, unsigned long size)
{
double reg = 0;
unsigned i;
#pragma omp parallel shared(vec) private(i) reduction(+:reg)
{
#pragma omp for schedule(dynamic, CHUNK_SIZE)
for (i = 0; i < size; ++i)
{
reg += vec[i];
}
}
return reg / (T)size;
}
template <typename T>
static void rmVal(T* vec, unsigned long size, T val)
{
unsigned i;
#pragma omp parallel shared(vec, val) private(i)
{
#pragma omp for schedule(dynamic, CHUNK_SIZE)
for (i = 0; i < size; ++i)
{
vec[i] -= val;
}
}
}
template <typename T>
static T mulReduceGlobal(T* vec1, T* vec2, unsigned long size)
{
double reg = 0;
unsigned i;
#pragma omp parallel shared(vec1, vec2) private(i) reduction(+:reg)
{
#pragma omp for schedule(dynamic, CHUNK_SIZE)
for (i = 0; i < size; ++i)
{
reg += vec1[i] * vec2[i];
}
}
return reg / (T)(size - 1);
}
template <typename T>
void cov(T* vec[], unsigned long size, T out[10][10])
{
int i, j;
for (i = 0; i < 10; ++i)
{
rmVal(vec[i], size, calcMeanGlobal(vec[i], size));
}
for (i = 0; i < 10; i++)
{
for (j = 0; j <= i; j++)
{
out[i][j] = mulReduceGlobal(vec[i], vec[j], size);
}
}
}
}
Jak widać praktycznie nic się nie zmieniło. Użycie OpenMP może być tak proste jak pętla parfor w MATLAB. O wykorzystaniu równoległych mechanizmów z biblioteki w ramach sekcji kodu decyduje pierwsze #pragma. Określa się tu obszary pamięci współdzielone pomiędzy wątkami (shared) – w tym przypadku wektory wejściowe. Poza tym są też zmienne prywatne (private) oraz akumulator (reduction(+:reg)). Druga sekcja #pragma definiuje plan realizacji pętli for() (schedule). Wątki, wykonują obliczenia blokami o długości zdefiniowanej przez CHUNK_SIZE i mogą to robić na kilka sposobów:
- Tablica podzielona na pomiędzy wątki, tak że najpierw są bloki obliczane przez pierwszy wątek, potem drugi itd.
- Bloki ułożone na przemian, tak że pierwszy blok oblicza wątek pierwszy, drugi blok oblicza wątek drugi itd.
- Każdy następny blok jest obliczany przez aktualnie swobodny wątek (taki, który aktualnie nic nie liczy)
My posłużyliśmy się schematem trzecim (dynamic), gdyż pozwala on na liniowe dostępy do pamięci. Zauważmy, że w schemacie pierwszym, każdy z wątków pracuje na odległych skrawkach pamięci RAM. W schemacie drugim, teoretycznie kolejne wątki są zawsze oddalone od siebie o CHUNK_SIZE, ale w praktyce mogą one wykonywać obliczenia z różną szybkością w wyniku czego z czasem się od siebie oddalą. Opcja trzecia gwarantuje to, że wszystkie wątki zawsze będą działały w obrębie sąsiadujących ze sobą bloków.
Ta prosta implementacja OpenMP pozwoliła na skrócenie czasu obliczania macierzy kowariancji do 75.721334 ms. To jest ponad pięciokrotnie szybciej niż w przypadku wersji szeregowej (więcej niż 4 może być skutkiem działania funkcji HT w procesorach Intel) i prawie dwa razy szybciej niż w środowisku MATLAB. Niestety nadal jest to więcej niż docelowe 50 ms. Zauważmy, że chociaż użyliśmy wszystkich rdzeni procesora to nie wykorzystaliśmy ich całkowicie.
AVX
Jeden rdzeń może przetwarzać więcej niż jedną operację zmiennoprzecinkową jednocześnie. Do tego mogą śłużyć instrukcje z zestawów takich jak SSE/AVX/AVX512. Nasz procesor wspiera maksymalnie standard AVX2. Oznacza to możliwość wykorzystania rejestrów o długości 256 bit, które pomieszczą 8 liczb zmiennoprzecinkowych pojedynczej precyzji (_m256) lub 4 podwójnej (_m256d). Wykonanie operacji zwykle ogranicza się do
- wczytania do rejestru zawartości tablicy _mm256_load…(),
- obliczeń np. c = _mm256_add_ps(a, b) (c = a + b, gdzie wszystkie zmienne są typu _m256 – dodanie do siebie 8 par liczb float jednocześnie, wynikiem również jest 8 liczb)
- zapisania wyniku _mm256_store…().
Ważne jest odpowiednie wyrównanie bufora w pamięci RAM, którego się używa zarówno wczytując jak i zapisując (stąd aligned_alloc()). Pełny zbiór dostępnych instrukcji można znaleźć tutaj.
Implementacja naszego algorytmu z użyciem AVX2 znajduje się poniżej. Niestety musieliśmy zmodyfikować również kod wykorzystujący OpenMP i ręcznie zaimplementowaliśmy schemat typu drugiego (jest mniej wydajny, ale użycie AVX2 to rekompensuje).
#define CHUNK_SIZE 1024
namespace simd{
double calcMeanGlobal(float* vec, unsigned long size)
{
unsigned i;
int tx, gl;
float out = 0;
#pragma omp parallel shared(vec, out) private(tx, gl)
{
tx = omp_get_thread_num();
gl = omp_get_num_threads() * CHUNK_SIZE;
__m256 __reg;
double reg = 0;
float tt[8] __attribute__((aligned (32)));
while (tx < size)
{
__reg = _mm256_set_ps(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
for (unsigned i = 0; i < CHUNK_SIZE; i += 8)
{
__reg = _mm256_add_ps(__reg, _mm256_load_ps(&vec[tx + i]));
}
//__reg = _mm256_add_ps(__reg, _mm256_permute2f128_ps(__reg , __reg , 1));
//__reg = _mm256_hadd_ps(__reg, __reg);
//__reg = _mm256_hadd_ps(__reg, __reg);
_mm256_store_ps(tt, __reg);
reg += tt[0] + tt[1] + tt[2] + tt[3] + tt[4] + tt[5] + tt[6] + tt[7];
tx += gl;
}
#pragma omp atomic update
out += reg;
}
return out / (double)size;
}
void rmVal(float* vec, unsigned long size, float val)
{
unsigned i;
#pragma omp parallel shared(vec, val) private(i)
{
#pragma omp for schedule(dynamic, CHUNK_SIZE)
for (i = 0; i < size; ++i)
{
vec[i] -= val;
}
}
}
double mulReduceGlobal(float* vec1, float* vec2, unsigned long size)
{
int tx, gl;
double out = 0;
#pragma omp parallel shared(vec1, vec2, out) private(tx, gl)
{
tx = omp_get_thread_num();
gl = omp_get_num_threads() * CHUNK_SIZE;
__m256 __reg;
__m256 v1;
double reg = 0;
float tt[8] __attribute__((aligned (16)));
while (tx < size)
{
__reg = _mm256_set_ps(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
for (unsigned i = 0; i < CHUNK_SIZE; i += 8)
{
v1 = _mm256_mul_ps(_mm256_load_ps(&vec1[tx + i]),
_mm256_loadu_ps(&vec2[tx + i]));
__reg = _mm256_add_ps(__reg, v1);
}
//__reg = _mm256_add_ps(__reg, _mm256_permute2f128_ps(__reg , __reg , 1));
//__reg = _mm256_hadd_ps(__reg, __reg);
//__reg = _mm256_hadd_ps(__reg, __reg);
_mm256_store_ps(tt, __reg);
reg += tt[0] + tt[1] + tt[2] + tt[3] + tt[4] + tt[5] + tt[6] + tt[7];
tx += gl;
}
#pragma omp atomic update
out += reg;
}
return out / (double)(size - 1);
}
void cov(float* vec[], unsigned long size, float out[10][10])
{
int i, j;
for (i = 0; i < 10; ++i)
{
rmVal(vec[i], size, calcMeanGlobal(vec[i], size));
}
for (i = 0; i < 10; i++)
{
for (j = 0; j <= i; j++)
{
out[i][j] = mulReduceGlobal(vec[i], vec[j], size);
}
}
}
}
Tym razem udało się obliczyć kowariancję w czasie 39.095279 ms, czyli prawie 10 razy szybciej niż w przypadku najprostszej implementacji C++. Jest to również szybciej niż potrzeba w założonym przez nas przypadku (50 ms). Moglibyśmy na tym poprzestać, ale… Można obliczyć taką macierz kowariancji w czasie, krótszym niż 50 ms bez zużywania tak wielu cennych zasobów procesora, które będą mogły zostać wykorzystane w innych celach.
OpenCL
Nasz Intel NUC8 posiada kartę graficzną zintegrowaną z procesorem, która wspiera standard OpenCL 3.0 – czemu z tego nie skorzystać? Jesteśmy zbyt leniwi, żeby używać czystego OpenCL dlatego wykorzystaliśmy nakładkę z bibliotek boost::compute. Z drugiej strony lubimy sobie utrudniać życie więc zamiast kombinacji funkcji lambda i compute::transform zdefiniowaliśmy własne kernele przy pomocy makra BOOST_COMPUTE_STRINGIZE_SOURCE.
const char kernelFunctionSource[] = BOOST_COMPUTE_STRINGIZE_SOURCE(
kernel void sum1(__global const float* var, __global float* results, int size, int i)
{
//Wymiary przestrzenie roboczych - indeksowanie tablic
int tx = get_global_id(0);
int gl = get_global_size(0);
int lx = get_local_id(0);
int ll = get_local_size(0);
int txD = get_group_id(0);
int ll2;
__local float buffer[256];
float reg = 0;
while(tx < size)
{
reg += var[tx + i * size];
tx += gl;
}
buffer[lx] = reg;
ll2=ll/2; // redukcja hierarchiczna
while (ll2 > 0)
{
barrier(CLK_LOCAL_MEM_FENCE);
if (lx < ll2)
{
buffer[lx] += buffer[lx + ll2];
}
ll2 = ll2 / 2;
}
// barrier(CLK_LOCAL_MEM_FENCE);
if (lx == 0)
{
results[txD] = buffer[0];
}
}
kernel void sum2(__global float* bufferGlobal, __global float* result, float div, int offset)
{
int tx = get_global_id(0);
int gl = get_global_size(0);
int lx = get_local_id(0);
int ll = get_local_size(0);
int txD = get_group_id(0);
int ll2;
__local float kot[256];
float reg = 0;
reg += bufferGlobal[tx];
kot[lx] = reg;
ll2=ll/2;
while (ll2 > 0)
{
barrier(CLK_LOCAL_MEM_FENCE);
if (lx < ll2)
{
kot[lx] += kot[lx + ll2];
}
ll2 = ll2 / 2;
}
//barrier(CLK_LOCAL_MEM_FENCE);
if (lx == 0)
{
result[txD + offset] = kot[0] / div;
}
}
kernel void cov1(global float* var, global const float* sumBuffer, global float* results, int size, int i)
{
int tx = get_global_id(0);
int gl = get_global_size(0);
int lx = get_local_id(0);
int ll = get_local_size(0);
int txD = get_group_id(0);
int gr = gl / ll;
int ll2;
local float buffer[256];
float reg[10];
float t;
float mean = sumBuffer[0];
// Szybciej na CPU
reg[i] = 0;
while(tx < size)
{
t = var[tx + i * size] - mean;
var[tx + i * size] = t;
reg[i] += t * t;
tx += gl;
}
for (int r = i - 1; r >= 0; --r)
{
tx = get_global_id(0);
reg[r] = 0;
while(tx < size)
{
reg[r] += var[tx + i * size] * var[tx + r * size];
tx += gl;
}
}
//Szybciej na GPU
/*for (int r = 0; r <= i; ++r)
{
reg[r] = 0;
}
while(tx < size)
{
t = var[tx + i * size] - mean;
var[tx + i * size] = t;
reg[i] += t * t;
for (int r = 0; r < i; ++r)
{
reg[r] += t * var[tx + r * size];
}
tx += gl;
}*/
for (int r = 0; r <= i; ++r)
{
barrier(CLK_LOCAL_MEM_FENCE);
buffer[lx] = reg[r];
ll2=ll/2;
while (ll2 > 0)
{
barrier(CLK_LOCAL_MEM_FENCE);
if (lx < ll2)
{
buffer[lx] += buffer[lx + ll2];
}
ll2 = ll2 / 2;
}
if (lx == 0)
{
results[txD + r * gr] = buffer[0];
}
}
}
);
Redukcje rozbiliśmy na dwa kernele. Pierwszy z nich (sum1()) redukuje rozmiar wektora do ilości grup wątków (GLOBAL_WORK_SIZE/LOCAL_WORK_SIZE). Kolejne wątki sumują elementy wektora do swojego prywatnego rejestru reg. Potem zapisują wynik do lokalnego (współdzielonego w ramach grupy wątków) bufora buffer. Następnie zachodzi redukcja hierarchiczna w obrębie wszystkich grup. Pierwszy wątek każdej grupy zwraca wynik w tablicy wyjściowej pod indeksem grupy. Drugi etap sumy (kernel sum2()) uruchamia jedną grupę wątków do zsumowania wyników z poprzedniego kernela (również metodą sumy hierarchicznej). Po zsumowaniu, pierwszy wątek w grupie dzieli wynik przez pierwotny rozmiar wektora i zapisuje do tablicy wyjściowej pod indeksem swojej grupy.
W trzecim kernelu (cov1()) każdy wątek posiada prywatną tablicę reg o rozmiarze 10. Służy ona do przechowywania sum iloczynów liczb z wektorów, z których obliczamy kowariancję. Najpierw przelicza się ostatni wektor w tablicy (o indeksie aktualnie obliczanego wiersza macierzy kowariancji) tak, żeby odjąć od niego średnią i jednocześnie policzyć wspomniany iloczyn. Następnie w pętli for() iloczyn elementów jest obliczany, dla każdej kombinacji tego wektora z poprzednimi – tak, kernele będą wywoływane 10 razy 
. Potem znowu pojawia się suma hierarchiczna kolejno dla wszystkich elementów tablicy reg. Uzyskanie kowariancji wymaga ponownego uruchomienia kernela sum2() z liczbą grup wątków taką jak indeks obliczanego wiersza macierzy kowariancji. Takie kernele należy odpowiednio wykorzystać z poziomu C++. Powstał więc następujący kod:
#define GLOBAL_WORK_SIZE 65536
#define LOCAL_WORK_SIZE 256
namespace ocl{
ComputingPlatformPtr platform;
std::shared_ptr<CLUnit> device;
std::shared_ptr<compute::vector<cl_float>> inputBuffer;
std::shared_ptr<compute::vector<cl_float>> reduceBuffer;
std::shared_ptr<compute::vector<cl_float>> reduceResult;
std::shared_ptr<compute::vector<cl_float>> resultBuffer;
std::shared_ptr<compute::vector<cl_float>> resultCL;
std::vector<float> result(100);
compute::kernel reduceKernel1;
compute::kernel reduceKernel2;
compute::kernel covKernel1;
void setup(unsigned id, unsigned size)
{
ComputingPlatform::Config conf(std::vector<unsigned>({id}));
platform = std::make_unique<ComputingPlatform>(conf);
device = platform->get(0);
std::cout << *device << std::endl;
inputBuffer = std::make_shared<compute::vector<cl_float>>(size * 10, *device->ctx());
reduceBuffer = std::make_shared<compute::vector<cl_float>>(GLOBAL_WORK_SIZE / LOCAL_WORK_SIZE * 10, *device->ctx());
reduceResult = std::make_shared<compute::vector<cl_float>>(10, *device->ctx());
resultCL = std::make_shared<compute::vector<cl_float>>(100, *device->ctx());
std::string options = "-cl-fast-relaxed-math";
compute::program program = compute::program::build_with_source(kernelFunctionSource, *device->ctx(), options);
reduceKernel1 = program.create_kernel("sum1");
reduceKernel2 = program.create_kernel("sum2");
covKernel1 = program.create_kernel("cov1");
}
void cov(float* vec[], unsigned size, float out[10][10])
{
for (unsigned i = 0; i < 10; ++i)
{
//device->defaultCmdQueue.enqueue_barrier();
compute::copy_async(vec[i],
vec[i] + size,
inputBuffer->begin() + i * size,
device->defaultCmdQueue);
//device->defaultCmdQueue.enqueue_barrier();
reduceKernel1.set_args(inputBuffer->get_buffer(),
reduceBuffer->get_buffer(),
static_cast<int>(size),
i);
device->defaultCmdQueue.enqueue_1d_range_kernel(reduceKernel1, 0, GLOBAL_WORK_SIZE, LOCAL_WORK_SIZE);
//device->defaultCmdQueue.enqueue_barrier();
reduceKernel2.set_args(reduceBuffer->get_buffer(),
reduceResult->get_buffer(),
static_cast<float>(size),
0);
device->defaultCmdQueue.enqueue_1d_range_kernel(reduceKernel2, 0, GLOBAL_WORK_SIZE/LOCAL_WORK_SIZE, GLOBAL_WORK_SIZE/LOCAL_WORK_SIZE);
//device->defaultCmdQueue.enqueue_barrier();
covKernel1.set_args(inputBuffer->get_buffer(),
reduceResult->get_buffer(),
reduceBuffer->get_buffer(),
static_cast<int>(size),
i);
device->defaultCmdQueue.enqueue_1d_range_kernel(covKernel1, 0, GLOBAL_WORK_SIZE, LOCAL_WORK_SIZE);
//device->defaultCmdQueue.enqueue_barrier();
reduceKernel2.set_args(reduceBuffer->get_buffer(),
resultCL->get_buffer(),
static_cast<float>(size - 1),
i * 10);
device->defaultCmdQueue.enqueue_1d_range_kernel(reduceKernel2, 0, (i + 1)*(GLOBAL_WORK_SIZE/LOCAL_WORK_SIZE), GLOBAL_WORK_SIZE/LOCAL_WORK_SIZE);
}
device->defaultCmdQueue.finish();
device->defaultCmdQueue.flush();
compute::copy(resultCL->begin(),
resultCL->end(),
result.begin(),
device->defaultCmdQueue);
for (unsigned i = 0; i < 10; ++i)
{
for (unsigned j = 0; j <= i; j++)
{
out[i][j] = result[j + i * 10];
}
}
}
}
Funkcja setup() pozwala na wybranie urządzenia do obliczeń, utworzenie buforów oraz kompilację kerneli. Do obsługi urządzenia użyliśmy kodu, dostępnego wyłącznie dla członków KN RDSP (nasza wewnętrzna biblioteka, upraszczająca korzystanie z OpenCL). Tutaj zaznaczamy, że obiekt klasy CLUnit, zawiera domyślną kolejkę oraz kontekst, z których korzystamy.
Funcja cov() to właściwe wywołanie kerneli. W pętli for() obliczane są kolejne wiersze macierzy kowariancji. Najpierw kopiowany jest wektor z danymi. Potem uruchamiane są kernele. Polega to na ustawieniu im odpowiednich argumentów i umieszczeniu w kolejce, która działa asynchronicznie. Oczywiście argumenty kernela mogą być ustawione jednorazowo. W naszym wypadku zmieniają się w pętli i konieczne jest ich aktualizowanie za każdym razem (zmienna i). Po wprowadzeniu wszystkich kerneli do kolejki czekamy na jej zakończenie (funkcja finish()/flush()) i kopiujemy wyniki do tablicy wyjściowej out[][].
Wykonanie tego kodu pozwala na obliczenie kowariancji w czasie 108.465610 ms na procesorze (przy użyciu ICD intela) oraz 57.415763 ms na zintegrowanej karcie graficznej. Wynik osiągnięty przy pomocy procesora interesuje nas o tyle, że możemy go porównać z wcześniejszymi. Plasuje się on pomiędzy czasem uzyskanym przy pomocy środowiska MATLAB a C++/OpenMP. Wynik uzyskany przy pomocy IGPU, również nie jest zadowalający (o 7 ms za wolno). Przyczyną jest na przykład to, że pomiar czasu obejmuje kopiowanie danych do bufora OpenCL (compute::vector). Mimo wszystko zależy nam na skróceniu czasu obliczania kowariancji poniżej 50 ms. Z tego powodu zmodyfikowaliśmy kod zmieniając sposób dostępów do pamięci. Modyfikacja ta polega na zakomentowaniu kodu od komentarza //Szybciej na CPU do komentarza //Szybciej na GPU oraz odkomentowaniu kodu pod komentarzem //Szybciej na GPU. Skutkiem jest spowolnienie obliczeń CPU (421.027041 ms) oraz przyspieszenie IGPU (44.449896 ms). W ten sposób pokonaliśmy granicę 50 ms, nawet nie dotykając procesora! Z drugiej strony uwidaczniają się tutaj dwa problemy. Po pierwsze OpenCL jest rozwiązaniem uniwersalnym i kompatybilnym z różnymi platformami sprzętowymi, ale jeżeli oczekujemy od programu szybkości to i tak konieczna jest optymalizacja kodu pod konkretne urządzenie. Po drugie jeżeli wykorzystujemy na przykład kartę graficzną może być konieczne wysłanie/pobranie jakiś danych do/z jej pamięci, a to może trwać długo i lepiej tego unikać.
Podsumowanie
Wszystkie zmierzone czasy znajdują się na wykresie poniżej.
Niezależnie od tego co rozumiesz przez słowo szybko, najważniejsze jest znalezienie funkcji celu i użycie właściwego narzędzia do jej rozwiązania. Można walczyć o każdą mikrosekundę i optymalizować jedną linię kodu przez resztę swoich dni. Można też po prostu zrobić to co trzeba i iść na kebab.