Шпаргалка: cpu

Улучшаем предыдущий пример. Теперь реализуем этот 2д-коллайдер на OpenCL. Учитывая, что нам надо будет параллельно вычислять и менять текущую позицию и вектор движения одновременно для всех обьектов, то использовать один массив обьектов будет проблематично, т.к. нужно будет использовать различные примитивы синхронизации, что безусловно пагубно скажется на производительности. Поэтому обход/модификация всех обьектов будет производиться следующим образом: один массив обьектов будет только для чтения, а второй только для записи и каждый тред на GPU будет записывать изменения в свою ячейку, не пересекаясь с другими.

__kernel void collideAndUpdate(
__global Ball* b1,
__global Ball* b2,
const int ballCnt,
const int scrW,
const int scrH,
const double frameTimeMs)
{
// Get work-item identifiers.
int i = get_global_id(0);
Ball tmp = b1[i];
for (int j = 0; j < ballCnt; j++)
{
// check collision between one ball to others,
// but don't check collision to itself
if (j != i)
{
Ball tmp2 = b1[j];
tmp = checkCollision(tmp, tmp2);
}
}
// check borders
if (tmp.pos.x <= 0 || tmp.pos.x >= scrW)
{
tmp.f.x *= -1;
}
if (tmp.pos.y <= 0 || tmp.pos.y >= scrH)
{
tmp.f.y *= -1;
}
// update positions
tmp.pos.x += tmp.f.x * frameTimeMs;
tmp.pos.y += tmp.f.y * frameTimeMs;
b2[i] = tmp;
}

Поскольку не получилось подключить C++ заголовочный файл в код ядра, то пришлось по сути переписывать все заново для OpenCL:

Ball checkCollision(Ball b1, const Ball b2)
{
float2 p1 = (float2)(b1.pos.x, b1.pos.y);
float2 p2 = (float2)(b2.pos.x, b2.pos.y);
const float dist = getDistanceBetween(p1, p2);
if (dist < b1.r + b2.r)
{
// direction to other ball
float2 to2 = p2 - p1;
float2 f = (float2)(b1.f.x, b1.f.y);
// calculate dot product
float dotProd = dot(f, to2);
// if dot product is negative then force directed away from B ball
// and we do nothing
if (dotProd > 0)
{
// angle between normal and force (moving) vectors
float angle = getAngleTo(f, to2);
// the angle of incidence is equal to the angle of reflection
f = rotVect(f, angle * 2);
f = f * (-1);
b1.f.x = f.x;
b1.f.y = f.y;
}
}
return b1;
}

И структуру Ball и все векторные операции тоже пришлось переписывать. Потому что OpenCL компилятор и С++ компилятор отличаются очень сильно. Фактически настолько, что и там и там можно использовать только какие-то простые структуры и константы с дефайнами. Вот как выглядят векторные операции:

float getDistanceBetween(float2 p1, float2 p2)
{
float2 p1p2 = p2 - p1;
float dist = native_sqrt(p1p2[0] * p1p2[0] + p1p2[1] * p1p2[1]);
return dist;
}
float getLen(float2 p)
{
float len = getDistanceBetween((float2)(0, 0), p);
return len;
}
float getCrossProd(float2 p1, float2 p2)
{
float crossProd = p1[0] * p2[1] - p1[1] * p2[0];
return crossProd;
}
float getAngleTo(float2 p1, float2 p2)
{
return asin(getCrossProd(p1, p2) / (getLen(p1) * getLen(p2)));
}
float2 rotVect(float2 v, float angle)
{
// first of all create rotation matrix
float c = cos(angle);
float s = sin(angle);
float2 mr0 = { c, -s }; // first row
float2 mr1 = { s, c }; // second row
// get rotated vector by multiply matrix to vector
float2 tmp;
tmp[0] = (v[0] * mr0[0] + v[1] * mr0[1]);
tmp[1] = (v[0] * mr1[0] + v[1] * mr1[1]);
return tmp;
}

Я не нашел как работать с матрицами, поэтому в качестве матрицы 2х2 я использовал просто два вектора float2, каждый из которых играет роль строки в матрице (строки 30-31). То есть довольно много заново написанного кода и если кто-то задумал перенести что-то с помощью OpenCL на видеокарту, то пусть имеют ввиду, что такой даже не копипасты, а полной переработки кода с учетом кучи нюансов будет очень много.

Чтобы увидеть реальную разницу в производительности, понадобилось увеличить количество шариков с 500 до 20000 и уменьшить изх диаметр до 1, чтобы они все поместились. Результат:

Производительность на GPU при 20 тысячах шариков получилась 29.3 кадра в секунду, а на CPU всего лишь 1.6 кадра в секунду. Разница, как говорится, налицо!

Тестовая платформа: Ryzen 3700X, 16GB RAM, RTX 3060 12GB

Весь код здесь.

Продолжение. Начало тут.

В предыдущий раз мы разбивали исходные матрицы на субматрицы меньшего размера, которые перемножали и потом складывали результат. Работает это так:

Почему это работает? Каждый элемент результирующей матрицы это результат перемножения соотв. элементов строки и столбца двух меремножаемых матриц. То есть мы по сути просто разбиваем строки на подстроки и столбцы на подстолбцы, которые соответсвенно содержатся в субматрицах. Вот так это работает:

Ну, в прошлый раз мы реализовали этот алгоритм, а теперь интересно посмотреть, насколько эффективно он работает по сравнению с "наивной" реализацией, когда мы просто напрямую из памяти устройства берем элементы строк и столбцов, не кешируя их в локальную память. Вот код "наивного" метода:

__kernel void multiplyMatrices(__global int* a,
__global int* b,
__global int* c,
const int M,
const int N,
const int K)
{
//Get work - item identifiers.
int colIndex = get_global_id(0);
int rowIndex = get_global_id(1);
//Compute element c[rowIndex, colIndex].
int sum = 0;
for (int k = 0; k < K; k++) {
sum += a[rowIndex*K + k] * b[k*N + colIndex];
}
// save result
c[(rowIndex * N) + colIndex] = sum;
}

При размере матриц 1024 на 1024, метод с кешированием работает приблизительно в 1.5 раза быстрее "наивного". Неплохо, но можно лучше! Нашел еще один способ оптимизации. Суть его в том, чтобы каждый тред обрабатывал не один выходной элемент, а несколько - таким образом уменьшится число обращений в память и должна увеличиться производительность. Вот код:

__kernel void multiplyMatrices(__global int* a,
__global int* b,
__global int* c,
const int M,
const int N,
const int K)
{
const int RTS = SUB_SIZE / WPT;
// Get work-item identifiers.
const int row = get_local_id(0); // Local row ID (max: SUB_SIZE)
const int col = get_local_id(1); // Local col ID (max: SUB_SIZE/WPT == RTS)
const int globalRow = SUB_SIZE * get_group_id(0) + row; // Row ID of C (0..M)
const int globalCol = SUB_SIZE * get_group_id(1) + col; // Col ID of C (0..N)
// Create submatrices that will cache the matrices A and B in local memory.
__local int aSub[SUB_SIZE][SUB_SIZE];
__local int bSub[SUB_SIZE][SUB_SIZE];
// Initialize accumulator
float acc[WPT];
for (int w = 0; w < WPT; w++)
{
acc[w] = 0.0f;
}
// Loop over all submatrices.
const int nSub = K / SUB_SIZE;
for(int s = 0; s < nSub; s++)
{
// Load one tile of A and B into local memory
for (int w = 0; w < WPT; w++)
{
const int tiledRow = SUB_SIZE * s + row;
const int tiledCol = SUB_SIZE * s + col;
aSub[col + w * RTS][row] = a[(tiledCol + w * RTS)*M + globalRow];
bSub[col + w * RTS][row] = b[(globalCol + w * RTS)*K + tiledRow];
}
// Synchronize all work-items in this work-group.
barrier(CLK_LOCAL_MEM_FENCE);
// Perform the computation for a single tile
for (int k = 0; k < SUB_SIZE; k++)
{
for (int w = 0; w < WPT; w++)
{
acc[w] += aSub[k][row] * bSub[col + w * RTS][k];
}
}
// Synchronize all work-items in this work-group.
barrier(CLK_LOCAL_MEM_FENCE);
}
// Store the final results in C
for (int w = 0; w < WPT; w++)
{
c[(globalCol + w * RTS)*M + globalRow] = acc[w];
}
}

В этом коде SUB_SIZE и WPT - это дефайны, переданные в kernel-код. Первый это размер субматрицы, а второй это Work Per Thread или количество вычисляемых элементов для одного треда. У нас он равен 8. Для того, чтобы это сработало, надо уменьшить количество тредов в рабочей группе. Т.е. нужно сообщить коду ядра, как именно мы хотим разбить работу над массивами. Вот примеры для всех трех методов:

// naive way: each thread calculate one element of out matrix 
cl::CommandQueue q(context, device);
q.enqueueNDRangeKernel(kernel, cl::NullRange, cl::NDRange(N, M));
q.enqueueReadBuffer(cBuf, CL_TRUE, 0, M * N * sizeof(int), c);
 
//...
 
// submatrix way: group threads into work groups SUB_SIZExSUB_SIZE (8x8)
// and each group caches one tile 8 x 8 submatrix (one thread for one out item)
cl::CommandQueue qe(context, device);
q.enqueueNDRangeKernel(kernel, cl::NullRange, cl::NDRange(N,M), cl::NDRange(SUB_SIZE,SUB_SIZE));
q.enqueueReadBuffer(cBuf, CL_TRUE, 0, M * N * sizeof(int), c);
 
//...
 
// submatrix way + WPT: group threads into work groups SUB_SIZExSUB_SIZE (8x8)
// and each group caches one tile 8 x 8 and calculate 8 out matrice elements.
// Also we need to reduce threads amount in group
cl::CommandQueue qe(context, device);
q.enqueueNDRangeKernel(kernel, cl::NullRange, cl::NDRange(N,M/WPT), cl::NDRange(SUB_SIZE,SUB_SIZE/WPT));
q.enqueueReadBuffer(cBuf, CL_TRUE, 0, M * N * sizeof(int), c);

В последнем случае мы уменьшаем количество тредов в рабочей группе и так же уменьшаем количество столбцов, которые обратаывает каждая группа (поскольку теперь она обрабатывает их в 8 раз больше). То есть вот такая как бы разбивка матрицы между рабочими группами.

Теперь самое интересное - результаты. А они, оказывается, очень зависят от размера матрицы. Вот, например, результаты для матриц 1024 на 1024:

Как видно, тут "наивный" способ дает 9.3 секунды, субматричный 6.9 сек, а субматрица+WPT 5.7 секунд. То есть прогресс есть во всех случаях, хотя между submatrix и submatrix+WPT он минимален, хотя и заметен.

Вот результат на матрицах размером 2048 на 2048:

Тут отличия уже гораздо заметнее: наивный способ в 3 раза медленнее, чем submatrix и в 4 раза медленне, чем submatrix+WPT. Это уже очень заметно!

Теперь посмотрим на 4096 на 4096 элементов:

Тут у нас разница еще немного увеличивается, хотя и не сильно отличается от предыдущего замера (2k на 2k элементов). Итого в результате получается, что самый быстрый способ перемножения матриц на GPU превосходит самый быстрый способ на CPU в 3040ms/133ms = 23 раза! Это уже очень значительно.

Весь код можно посмотреть здесь.

Дополнительные пояснения к коду.

Я все три метода вычислений на GPU (naive, submatrix, submatrix+WPT) поместил в один файл, хотя можно было и в разные, но я хотел посмотреть как работает передача дефайнов в компилятор OpenCL-программы и оказалось, что работает. С помощью дирепктив #if define и #endif я отключил те или иные участки в OpenCL-коде, что может оказаться очень удобно в ряде случаев. А еще использовать дефайны очень удобно, если хочется передать какой-то дополнительный параметр, но не хочется менять сигнатуру kernel-функции.

Шпаргалка

пятница, 29 декабря 2023 г.

2d-коллайдер и его реализация на OpenCL

вторник, 21 ноября 2023 г.

OpenCL: продолжение. Оптимизация кода и тестирование на бОльших объемах данных

пятница, 29 декабря 2023 г.

2d-коллайдер и его реализация на OpenCL

вторник, 21 ноября 2023 г.

OpenCL: продолжение. Оптимизация кода и тестирование на бОльших объемах данных

пятница, 29 декабря 2023 г.

вторник, 21 ноября 2023 г.