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

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

1. __kernel void collideAndUpdate(
2. __global Ball* b1,
3. __global Ball* b2,
4. const int ballCnt,
5. const int scrW,
6. const int scrH,
7. const double frameTimeMs)
8. {
9. // Get work-item identifiers.
10. int i = get_global_id(0);
11. Ball tmp = b1[i];
12. for (int j = 0; j < ballCnt; j++)
13. {
14. // check collision between one ball to others,
15. // but don't check collision to itself
16. if (j != i)
17. {
18. Ball tmp2 = b1[j];
19. tmp = checkCollision(tmp, tmp2);
20. }
21. }
22.  
23. // check borders
24. if (tmp.pos.x <= 0 || tmp.pos.x >= scrW)
25. {
26. tmp.f.x *= -1;
27. }
28. if (tmp.pos.y <= 0 || tmp.pos.y >= scrH)
29. {
30. tmp.f.y *= -1;
31. }
32.  
33. // update positions
34. tmp.pos.x += tmp.f.x * frameTimeMs;
35. tmp.pos.y += tmp.f.y * frameTimeMs;
36.  
37. b2[i] = tmp;
38. }

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

1. Ball checkCollision(Ball b1, const Ball b2)
2. {
3. float2 p1 = (float2)(b1.pos.x, b1.pos.y);
4. float2 p2 = (float2)(b2.pos.x, b2.pos.y);
5. const float dist = getDistanceBetween(p1, p2);
6. if (dist < b1.r + b2.r)
7. {
8. // direction to other ball
9. float2 to2 = p2 - p1;
10. float2 f = (float2)(b1.f.x, b1.f.y);
11. // calculate dot product
12. float dotProd = dot(f, to2);
13. // if dot product is negative then force directed away from B ball
14. // and we do nothing
15. if (dotProd > 0)
16. {
17. // angle between normal and force (moving) vectors
18. float angle = getAngleTo(f, to2);
19. // the angle of incidence is equal to the angle of reflection
20. f = rotVect(f, angle * 2);
21. f = f * (-1);
22. b1.f.x = f.x;
23. b1.f.y = f.y;
24. }
25. }
26. return b1;
27. }

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

1. float getDistanceBetween(float2 p1, float2 p2)
2. {
3. float2 p1p2 = p2 - p1;
4. float dist = native_sqrt(p1p2[0] * p1p2[0] + p1p2[1] * p1p2[1]);
5. return dist;
6. }
7.  
8. float getLen(float2 p)
9. {
10. float len = getDistanceBetween((float2)(0, 0), p);
11. return len;
12. }
13.  
14. float getCrossProd(float2 p1, float2 p2)
15. {
16. float crossProd = p1[0] * p2[1] - p1[1] * p2[0];
17. return crossProd;
18. }
19.  
20. float getAngleTo(float2 p1, float2 p2)
21. {
22. return asin(getCrossProd(p1, p2) / (getLen(p1) * getLen(p2)));
23. }
24.  
25. float2 rotVect(float2 v, float angle)
26. {
27. // first of all create rotation matrix
28. float c = cos(angle);
29. float s = sin(angle);
30. float2 mr0 = { c, -s }; // first row
31. float2 mr1 = { s, c }; // second row
32. // get rotated vector by multiply matrix to vector
33. float2 tmp;
34. tmp[0] = (v[0] * mr0[0] + v[1] * mr0[1]);
35. tmp[1] = (v[0] * mr1[0] + v[1] * mr1[1]);
36. return tmp;
37. }

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

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

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

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

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

OpenCL: продолжение работы с изображениями. Оптимизация kernel-кода и неожиданные результаты тестов. Использование printf при отладке kernel-кода

 Продолжение. Начало здесь.

Итак, теперь мы попробуем оптимизировать код обработки изображения при помощи локальной памяти и рабочих групп. Будем, как и в случае работы с перемножением матриц, кешировать значения в локальный буфер. Код выглядит так:

1. // This kernel function convolves an image input_image[imgWidth, imgHeight]
2. // with a mask of size maskSize by caching submatrices from the input image
3. // in the device local memory.
4. __kernel void filterImageCached(
5. __global unsigned char* inImg,
6. __global unsigned char* outImg,
7. const int bytesPerPix,
8. const unsigned int maskSize,
9. __constant float* mask)
10. {
11. // Get work-item identifiers.
12. int x = get_global_id(0);
13. int y = get_global_id(1);
14. int lx = get_local_id(0);
15. int ly = get_local_id(1);
16. int imgW = get_global_size(0);
17. int imgH = get_global_size(1);
18. int offset = ((y * imgW) + x) * bytesPerPix;
19.  
20. // Declare submatrix used to cache the input image on local memory.
21. __local unsigned char sub[SUB_SIZE][SUB_SIZE][4];
22.  
23. sub[ly][lx][0] = inImg[offset];
24. sub[ly][lx][1] = inImg[offset + 1];
25. sub[ly][lx][2] = inImg[offset + 2];
26. if (bytesPerPix == 4) // if we have alfa-channel
27. sub[ly][lx][3] = inImg[offset + 3];
28.  
29. // Synchronize all work-items in this work-group.
30. barrier(CLK_LOCAL_MEM_FENCE);
31.  
32. // Check if the mask cannot be applied to the current pixel
33. if (x < maskSize / 2
34. || y < maskSize / 2
35. || x >= imgW - maskSize / 2
36. || y >= imgH - maskSize / 2)
37. {
38. outImg[offset] = 0;
39. outImg[offset + 1] = 0;
40. outImg[offset + 2] = 0;
41. if (bytesPerPix == 4) // if we have alfa-channel
42. outImg[offset + 3] = inImg[offset + 3];
43. return;
44. }
45.  
46. // Apply mask based on the neighborhood of pixel inputImg.
47. int outSumB = 0;
48. int outSumG = 0;
49. int outSumR = 0;
50. for (size_t k = 0; k < maskSize; k++)
51. {
52. for (size_t l = 0; l < maskSize; l++)
53. {
54. // Calculate the current mask index.
55. size_t maskIdx = (maskSize - 1 - k) + (maskSize - 1 - l) * maskSize;
56. // Compute output pixel.
57. size_t maskLX = lx - maskSize / 2 + k;
58. size_t maskLY = ly - maskSize / 2 + l;
59. // Сheck if the current input pixel is in the local memory
60. if (maskLX >= 0 && maskLX < SUB_SIZE && maskLY >= 0 && maskLY < SUB_SIZE)
61. {
62. outSumB += sub[maskLY][maskLX][0] * mask[maskIdx];
63. outSumG += sub[maskLY][maskLX][1] * mask[maskIdx];
64. outSumR += sub[maskLY][maskLX][2] * mask[maskIdx];
65. }
66. else
67. {
68. // Read the current input pixel from the global memory
69. size_t maskX = x - maskSize / 2 + k;
70. size_t maskY = y - maskSize / 2 + l;
71. int offsetM = ((maskY * imgW) + maskX) * bytesPerPix;
72. outSumB += inImg[offsetM] * mask[maskIdx];
73. outSumG += inImg[offsetM + 1] * mask[maskIdx];
74. outSumR += inImg[offsetM + 2] * mask[maskIdx];
75. }
76. }
77. }
78.  
79. // Write output pixel.
80. outImg[offset] = MinMaxVal(0, 255, outSumB);
81. outImg[offset + 1] = MinMaxVal(0, 255, outSumG);
82. outImg[offset + 2] = MinMaxVal(0, 255, outSumR);
83. if (bytesPerPix == 4) // if we have alfa-channel
84. outImg[offset + 3] = inImg[offset + 3];
85. }

Вкратце работает так: создаем кеш размером SUB_SIZE на SUB_SIZE (в нашем случае 16 на 16) и записываем в него пиксели (строки 116-122). На краях рабочей группы озможны случаи, когда нам требуются пиксели, которые не закешированы и тогда мы извлекаем значения из глобальной памяти (строки 60-75). А дальше все как обычно: вычисляем сумму согласно коэффициентам фильтра и записываем ее в глобальный массив (строки 80-84).

Теперь самое интересное - результаты, которые довольно неожиданные. Использование кеша должно было улучшить производительность, по аналогии с прошлым разом, но получилось, что неоптимизированный метод (функция filterImage) работает быстрее, чем оптимизированный метод (функция filterImageCached), причем довольно заметно: 2.2-2.3мсек "неоптимизированного" против 2.7-2.8мсек "оптимизированного" кода для картинки 2048 на 1536. Я ожидал хотя бы незначительного, но прироста, а тут получилось значительное падение. Причины точно не ясны, но можно все же сделать несколько предположений. Например, аппаратный кеш видеокарты (у RTX 3060 он 3 мегабайта) вместе с быстрой памятью (192bit GDDR6) делает излишним дополнительные оптимизации в этом случае, но в то же время накладные расходы (как минимум дополнительные 4 сравнения в строке 60 и синхронизация тредов в строке 30) никуда не деваются. Тут, конечно, было бы интересно посмотреть на результаты с какой-то более слабой видеокартой c 64-битной шиной памяти - например GTX 1030. Так что не все так просто с оптимизациями и всегда надо очень тщательно проверять что и на каком оборудовании запускается.

Еще один нюанс с отладкой kernel-кода: оказывается, там можно использовать самый обычный printf(). Но, конечно, делать это надо с осторожностью, потому что обработка даже маленькой картинки 400 на 400 пикселей даст на выходе 16 тысяч трейсов в терминале, в которых можно потеряться.

Еще один важный нюанс при работе с OpenCL состоит в том, что при использовании рабочих групп непременно надо, чтобы размер входного массива был обязательно кратен размерам рабочей группы. Потому что я довольно долго не мог понять, почему неоптимизированная функция работает, а оптимизированная (и обработка разбита на группы тредов 16 на 16) не работает вообще. Оказалось, что enqueueNDRangeKernel возвращала мне значиние ошибки -54, что есть значение константы CL_INVALID_WORK_GROUP_SIZE. Об этом написали добрые люди на форуме: 

Error -54 is CL_INVALID_WORK_GROUP_SIZE - values in your “globalWorkSize” array are not divisible with values in your “localWorkSize” array, and that’s it (remember that global work size is really total number of threads along each dimension, and not the size of the “block” of threads along corresponding dimension).

Так что неудивительно, что OpenCL отказался обрабатывать мою картинку размером 1012 на 737 пикселей. После того, как я сменил разрешение на 1024 на 768, все заработало. Проверяйте коды ошибок - это может быть очень полезно.

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

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