Тестирование скорости перемножения матриц 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 

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