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-функции.