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

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

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

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

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

1. __kernel void multiplyMatrices(__global int* a,
2. __global int* b,
3. __global int* c,
4. const int M,
5. const int N,
6. const int K)
7. {
8. //Get work - item identifiers.
9. int colIndex = get_global_id(0);
10. int rowIndex = get_global_id(1);
11.  
12. //Compute element c[rowIndex, colIndex].
13. int sum = 0;
14. for (int k = 0; k < K; k++) {
15. sum += a[rowIndex*K + k] * b[k*N + colIndex];
16. }
17. // save result
18. c[(rowIndex * N) + colIndex] = sum;
19. }

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

1. __kernel void multiplyMatrices(__global int* a,
2. __global int* b,
3. __global int* c,
4. const int M,
5. const int N,
6. const int K)
7. {
8. const int RTS = SUB_SIZE / WPT;
9. // Get work-item identifiers.
10. const int row = get_local_id(0); // Local row ID (max: SUB_SIZE)
11. const int col = get_local_id(1); // Local col ID (max: SUB_SIZE/WPT == RTS)
12. const int globalRow = SUB_SIZE * get_group_id(0) + row; // Row ID of C (0..M)
13. const int globalCol = SUB_SIZE * get_group_id(1) + col; // Col ID of C (0..N)
14.  
15. // Create submatrices that will cache the matrices A and B in local memory.
16. __local int aSub[SUB_SIZE][SUB_SIZE];
17. __local int bSub[SUB_SIZE][SUB_SIZE];
18.  
19. // Initialize accumulator
20. float acc[WPT];
21. for (int w = 0; w < WPT; w++)
22. {
23. acc[w] = 0.0f;
24. }
25.  
26. // Loop over all submatrices.
27. const int nSub = K / SUB_SIZE;
28. for(int s = 0; s < nSub; s++)
29. {
30. // Load one tile of A and B into local memory
31. for (int w = 0; w < WPT; w++)
32. {
33. const int tiledRow = SUB_SIZE * s + row;
34. const int tiledCol = SUB_SIZE * s + col;
35. aSub[col + w * RTS][row] = a[(tiledCol + w * RTS)*M + globalRow];
36. bSub[col + w * RTS][row] = b[(globalCol + w * RTS)*K + tiledRow];
37. }
38.  
39. // Synchronize all work-items in this work-group.
40. barrier(CLK_LOCAL_MEM_FENCE);
41.  
42. // Perform the computation for a single tile
43. for (int k = 0; k < SUB_SIZE; k++)
44. {
45. for (int w = 0; w < WPT; w++)
46. {
47. acc[w] += aSub[k][row] * bSub[col + w * RTS][k];
48. }
49. }
50.  
51. // Synchronize all work-items in this work-group.
52. barrier(CLK_LOCAL_MEM_FENCE);
53. }
54.  
55. // Store the final results in C
56. for (int w = 0; w < WPT; w++)
57. {
58. c[(globalCol + w * RTS)*M + globalRow] = acc[w];
59. }
60. }
61.  

В этом коде 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-функции. 

Тестирование скорости перемножения матриц OpenCL vs CPU. Полезность кеширования данных.

В предыдущем посте мы настроили OpenCL, а теперь мы попробуем провести тест производительности. За основу возьмем пример cached_matrix_multiplication, где перемножаются две матрицы и результат записывается в третью. Kernel-функция там имеет такой вид:

1. /**
2.  * This kernel function efficiently multiplies two matrices a[M,K] and b[K,N]
3.  * by caching submatrices from those input matrices in the device local memory.
4.  */
5.  
6. __kernel void multiplyMatricesWithCache(__global int* a,
7. __global int* b,
8. __global int* c,
9. const int M,
10. const int N,
11. const int K){
12.  
13. /**
14.   * Declare the size of each submatrix (it must be
15.   * the same work-group size declared in the host code).
16.   */
17.  
18. const int SUB_SIZE = 16;
19.  
20. /**
21.   * Get work-item identifiers.
22.   */
23.  
24. int colIndex = get_local_id(0);
25. int rowIndex = get_local_id(1);
26. int globalColIndex = get_global_id(0);
27. int globalRowIndex = get_global_id(1);
28. int index = (globalRowIndex * N) + globalColIndex;
29.  
30. /**
31.   * Create submatrices that will cache the matrices A and B in local memory.
32.   */
33.  
34. __local int aSub[SUB_SIZE][SUB_SIZE];
35. __local int bSub[SUB_SIZE][SUB_SIZE];
36.  
37. /**
38.   * Initialize accumulator register.
39.   */
40.  
41. int sum = 0;
42.  
43. /**
44.   * Loop over all submatrices.
45.   */
46.  
47. const int nSub = K / SUB_SIZE;
48. for(int s = 0; s < nSub; s++){
49.  
50. /**
51.   * Load submatrices into local memory.
52.   */
53.  
54. const int sCol = SUB_SIZE * s + colIndex;
55. const int sRow = SUB_SIZE * s + rowIndex;
56. aSub[rowIndex][colIndex] = a[globalRowIndex * K + sCol];
57. bSub[rowIndex][colIndex] = b[sRow * N + globalColIndex];
58.  
59. /**
60.   * Synchronize all work-items in this work-group.
61.   */
62.  
63. barrier(CLK_LOCAL_MEM_FENCE);
64.  
65. /**
66.   * Perform the computation for a single submatrix.
67.   */
68.  
69. for(int k = 0; k < SUB_SIZE; k++){
70. sum += aSub[rowIndex][k] * bSub[k][colIndex];
71. }
72.  
73. /**
74.   * Synchronize all work-items in this work-group.
75.   */
76.  
77. barrier(CLK_LOCAL_MEM_FENCE);
78. }
79.  
80. /**
81.   * Store the final result in the matrix C.
82.   */
83.  
84. c[index] = sum;
85. }
86.  

Попробуем разобраться, что же тут происходит. Функция вычисляет один элемент результирующей матрицы. В более наивной реализации используется значение непосредственно из входных матриц, но такой подход очень медлительный, поскольку каждый раз приходится обращаться в глобальную память видеокарты, а она гораздо медленнее, чем локальная память группы. Мы формируем из отдельных тредов на GPU группы по 16 на 16 тредов и они последовательно в цикле обрабатывают блоки матрицы такого же размера (16 на 16), затем все эти полученные матрицы складываются и получается результирующая матрица (более подробно разобрано тут). Но перед тем, как вычислять сумму перемноженных строк и столбцов, сначала надо их закешировать и именно это и делается в строках 54-63. Делается это следующим образом: сначала мы из глобальной памяти (переменные помечены словом __global) устройства заполняем массив в локальной памяти (переменные помечены __local), причем чтобы получить элемент из глобального массива, мы используем глобальный идентификатор (global id) треда, а чтобы поместить в локальную память группы используем локальный идентификатор внутри группы (local id). Каждый тред внутри группы копирует в локальную память только свой один элемент из матрицы a и один элемент из матрицы b, но для вычисления результата ему нужна целиком строка и целиком столбец, поэтому, чтобы убедиться что локальный массив заполнен (т.е. каждый тред из группы параллельно работающих тредов заполнил свой элемент), мы синхронизируем работу тредов  специальной функцией  barrier(CLK_LOCAL_MEM_FENCE). 

Количество тредов внутри группы определяется вот тут:

cl::CommandQueue queue(context, device);
queue.enqueueNDRangeKernel(kernel, cl::NullRange, cl::NDRange(N, M), cl::NDRange(WG_SIZE[0], WG_SIZE[1]));
queue.enqueueReadBuffer(cBuf, CL_TRUE, 0, M * N * sizeof(int), c);

Функция enqueueNDRangeKernel() принимает вторым параметром смещение (в данном случае ноль, то есть вычисления с начала массива), третьим параметром размерность обрабатываемых данных (данном случае это двумерный массив), а четвертый параметр это размерность группы. В нашем случае это массив 16 на 16, то есть 256 тредов в каждой группе. Более подробно можно почитать здесь.

На стороне CPU произведение матриц выполнено предельно просто:

1. void seqMultiplyMatrices(int* a, int* b, int* c,
2. const int M,
3. const int N,
4. const int K) {
5. for (int i = 0; i < M; i++) {
6. for (int j = 0; j < N; j++) {
7. int sum = 0;
8. for (int k = 0; k < K; k++) {
9. sum += a[i*K + k] * b[j + k * N];
10. }
11. c[i*N + j] = sum;
12. }
13. }
14. }

То есть просто последовательно перемножаем строки и столбцы. Я с делал улучшеный вариант перемножения матриц на CPU и он выглядит так:

1. void parMultiplyMatrices_CPU(int* a, int* b, int* c,
2. const int M,
3. const int N,
4. const int K)
5. {
6. const auto thrCnt = std::thread::hardware_concurrency();
7. std::vector<std::thread> workers;
8. const int rod = N % thrCnt;
9. const unsigned int bColCntPerThread = (N - rod) / thrCnt;
10. for (int thrId = 0; thrId < thrCnt; ++thrId)
11. {
12. workers.push_back(std::thread([thrId, thrCnt, bColCntPerThread, a, b, c, M, N, K]()
13. {
14. std::vector<int> bColBuf(K);
15. auto bColbegIdx = thrId * bColCntPerThread;
16. bool isLastThread = (thrId == thrCnt - 1);
17. auto bColEndIdx = isLastThread ? N : bColbegIdx + bColCntPerThread;
18. //each thread iterate some amount of columns in B matrix
19. for (auto bColIdx = bColbegIdx; bColIdx < bColEndIdx; ++bColIdx)
20. {
21. //fill column buffer for B matr
22. for (int i = 0; i < K; ++i)
23. {
24. bColBuf[i] = b[i * N + bColIdx];
25. }
26. //multiply all rows of A to current column of B
27. for (int aRowIdx = 0; aRowIdx < M; ++aRowIdx)
28. {
29. int tmp = 0;
30. for (int aColIdx = 0; aColIdx < K; ++aColIdx)
31. {
32. tmp += a[aRowIdx * M + aColIdx] * bColBuf[aColIdx];
33. }
34. c[aRowIdx * M + bColIdx] = tmp;
35. }
36. }
37. }));
38. }
39.  
40. for (std::thread& t : workers)
41. {
42. t.join();
43. }
44. workers.clear();
45. }
46.  

Кроме многопоточности, я добавил еще одно важное улучшение: кеширование столбца одной из перемножаемых матриц. Почему именно столбца? - потому что строки вычитываются из памяти очень просто - они последовательно расположены в RAM друг за другом, а вот если идти по столбцам, то для считывания одного элемента надо перескакивать всю строку из сотен или тысяч байт, а это очень невыгодный режим работы памяти. Поэтому я в строках 22-25 сначала считываю весь столбец матрицы в буфер целиком, а потом многократно использую его для вычисления результата в строке 32. Результаты вот такие:

Самый медленный способ это ожидаемо наивная реализация перемножения матриц на процессоре (3 сек), на втором месте кешированное и распараллеленное перемножение на процессоре (0.163 сек), а на первом месте по скорости ожидаемо перемножение матриц на видеокарте (0.019 сек). Но все же, честно говоря, ожидал большей производительности от RTX 3060, которая оказалась быстрее процессора всего в 8.6 раз, а вот "кешированный параллельный" способ оказался в 18.7 раз быстрее, чем "наивная" реализация. Интересно, а что будет, если оставить только кеширование на CPU и полностью убрать многопоточность? (Что я сделал thrCnt = 1 в строке 6). А вот что:

И тут у нас вместо 0.163 сек получается 0.367 сек, то есть распараллеливание давало прирост всего лишь в 2.3 раза! А основной "ускоритель" вычисления - это именно использования кеша для столбцов. Что еще раз говорит нам об важности алгоритмов и понимания работы подсистемы памяти. Сам по себе мощный многоядерный процессор ничего не ускорит.

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

Код этого пример можно найти здесь.