Мне же было интересно узнать как покажет себя GPU, и я сравнил однопоточный код с многопоточным для CPU и совсем многопоточным для GPU, с помощью архитектуры параллельных вычислений CUDA.
GPU GeForce GTX 770
MSVS 2015 update 3
tookit CUDA 8.
Конечно, сборка была релизная и 64 битная, с оптимизацией /02 и /LTCG.
HKEY_LOCAL_MACHINE\SYSTEM\CurrentControlSet\Control\GraphicsDrivers\TdrDelay и перезагрузить компьютер.
CPU
Для начала я распараллелил алгоритм для CPU. Для этого в функцию передается диапазон для перебора значений внешнего цикла a. Затем весь диапазон 1..N разбивается на кусочки и скармливается ядрам процессора.
void Algorithm_1(const uint32_t N, const std::vector<uint64_t>& powers, const uint32_t from, const uint32_t to) {
for (uint32_t a = from; a < to; ++a) {
for (uint32_t b = 1; b < a; ++b) {
for (uint32_t c = 1; c < b; ++c) {
for (uint32_t d = 1; d < c; ++d) {
const uint64_t sum = powers[a] + powers[b] + powers[c] + powers[d];
if (std::binary_search(std::begin(powers), std::end(powers), sum)) {
const auto it = std::lower_bound(std::begin(powers), std::end(powers), sum);
uint32_t e = static_cast<uint32_t>(std::distance(std::begin(powers), it));
std::cout << a << " " << b << " " << c << " " << d << " " << e << "\n";
}
}
}
}
}
}
В моем случае было 4 ядра, и вычисления я выполнял через чуть-чуть доработанный пул потоков, который взял за основу тут. Так как никаких зависимостей между данными нет, скорость возрастает практически пропорционально количеству ядер.
GPU
Для GPU получается немного сложнее. Внешние два цикла перебора (a и b) будут развернуты, а в функции останется только перебор двух внутренних циклов (c и d). Можно себе представить что программа будет выполняться параллельно для всех значений а = 1..N и b = 1..N. На самом деле конечно это не так, все исполняться параллельно они все-таки не смогут, но это уже забота CUDA максимально распараллелить выполнение, этому можно помочь, если правильно настроить конфигурацию блоков и потоков.
int blocks_x = N;
int threads = std::min(1024, static_cast<int>(N));
int blocks_y = (N + threads - 1) / threads;
dim3 grid(blocks_x, blocks_y);
NaiveGPUKernel<<<grid, threads>>>(N);
blocks_x — это диапазон перебора первого цикла a.
А вот перебор второго цикла b пришлось сделать сложнее из-за ограничения видеокарты на количество потоков (1024), и поместить счетчик одновременно в threads и blocks_y:
a и b вычисляется так:
const int a = blockIdx.x + 1;
const int b = threadIdx.x + blockIdx.y * blockDim.x + 1;
Это точно не оптимальная конфигурация, так как получается много холостых циклов. Но дальше подстраивать значения я не стал.
Код для GPU получается вполне прямолинейный и очень похожий на CPU, только бинарный поиск приходится сделать руками.
__constant__ uint64_t gpu_powers[8192];
inline __device__ int gpu_binary_search(const uint32_t N, const uint64_t search) {
uint32_t l = 0;
uint32_t r = elements_count - 1;
uint32_t m;
while (l <= r) {
m = (l + r) / 2;
if (search < gpu_powers[m])
r = m - 1;
else if (search > gpu_powers[m])
l = m + 1;
else
return l;
}
return -1;
}
__global__ void NaiveGPUKernel(const uint32_t N) {
const int a = blockIdx.x + 1;
const int b = threadIdx.x + blockIdx.y * blockDim.x + 1;
if (b >= a)
return;
for (uint32_t c = 1; c < b; ++c) {
for (uint32_t d = 1; d < c; ++d) {
const uint64_t sum = gpu_powers[a] + gpu_powers[b] + gpu_powers[c] + gpu_powers[d];
const auto s = gpu_binary_search(N, sum);
if (s > 0) {
printf("%d %d %d %d %d\n", a, b, c, d, s);
}
}
}
}
Так же я реализовал оптимизации, подсказанные тут.
Замеры скорости
Замеры для CPU делались по два раза, так как они оказались намного стабильнее GPU, которых я делал уже по семь и выбирал лучшее время. Разброс для GPU мог быть двукратным, это я объясняю тем что для CUDA-расчетов использовался единственный в системе видеоадаптер.
| N | CPU время, мс | CPU (4 потока) время, мс | GPU время, мс | |
|---|---|---|---|---|
| @antonkrechetov | 100 | 58.6 | 16.7 | 14.8 |
| Базовый | 100 | 45.3 | 13.0 | 10.7 |
| Базовый оптимизация #1 | 100 | 6.3 | 2.1 | 12.7 |
| Базовый оптимизация #2 | 100 | 1.4 | 0.7 | 0.8 |
| @antonkrechetov | 250 | 2999.7 | 782.9 | 119.0 |
| Базовый | 250 | 2055.6 | 550.1 | 90.9 |
| Базовый оптимизация #1 | 250 | 227.2 | 60.3 | 109.2 |
| Базовый оптимизация #2 | 250 | 42.9 | 11.9 | 6.0 |
| @antonkrechetov | 500 | 72034.2 | 19344.1 | 1725.83 |
| Базовый | 500 | 38219.7 | 10200.8 | 976.7 |
| Базовый оптимизация #1 | 500 | 3725.1 | 926.5 | 1140.36 |
| Базовый оптимизация #2 | 500 | 630.7 | 170.2 | 48.7 |
| @antonkrechetov | 750 | 396566.0 | 105003.0 | 11521.2 |
| Базовый | 750 | 218615.0 | 57639.2 | 5742.5 |
| Базовый оптимизация #1 | 750 | 19082.7 | 4736.8 | 6402.1 |
| Базовый оптимизация #2 | 750 | 3272.0 | 846.9 | 222.9 |
| Базовый оптимизация #2 | 1000 | 10204.4 | 2730.3 | 1041.9 |
| Базовый оптимизация #2 | 1250 | 25133.1 | 6515.3 | 2445.5 |
| Базовый оптимизация #2 | 1500 | 51940.1 | 14005.0 | 4895.2 |
А теперь вколем немного адреналина в CPU!
И этим адреналином будет Profile-guided optimization.
Для PGO я привожу только время CPU, потому что GPU мало изменилось, и это ожидаемо.
| N | CPU время, мс |
CPU (4 потока) время, мс |
CPU+PGO время, мс |
CPU+PGO (4 потока) время, мс |
|
|---|---|---|---|---|---|
| @antonkrechetov | 100 | 58.6 | 16.7 | 55.3 | 15.0 |
| Базовый | 100 | 45.3 | 13.0 | 42.2 | 12.1 |
| Базовый оптимизация #1 | 100 | 6.3 | 2.1 | 5.2 | 1.9 |
| Базовый оптимизация #2 | 100 | 1.4 | 0.7 | 1.3 | 0.8 |
| @antonkrechetov | 250 | 2999.7 | 782.9 | 2966.1 | 774.1 |
| Базовый | 250 | 2055.6 | 550.1 | 2050.2 | 544.6 |
| Базовый оптимизация #1 | 250 | 227.2 | 60.3 | 200.0 | 53.2 |
| Базовый оптимизация #2 | 250 | 42.9 | 11.9 | 40.4 | 11.4 |
| @antonkrechetov | 500 | 72034.2 | 19344.1 | 68662.8 | 17959.0 |
| Базовый | 500 | 38219.7 | 10200.8 | 38077.7 | 10034.0 |
| Базовый оптимизация #1 | 500 | 3725.1 | 926.5 | 3132.9 | 822.2 |
| Базовый оптимизация #2 | 500 | 630.7 | 170.2 | 618.1 | 160.6 |
| @antonkrechetov | 750 | 396566.0 | 105003.0 | 404692.0 | 103602.0 |
| Базовый | 750 | 218615.0 | 57639.2 | 207975.0 | 54868.2 |
| Базовый оптимизация #1 | 750 | 19082.7 | 4736.8 | 15496.4 | 4082.3 |
| Базовый оптимизация #2 | 750 | 3272.0 | 846.9 | 3093.8 | 812.7 |
| Базовый оптимизация #2 | 1000 | 10204.4 | 2730.3 | 9781.6 | 2565.9 |
| Базовый оптимизация #2 | 1250 | 25133.1 | 6515.3 | 23704.3 | 6244.1 |
| Базовый оптимизация #2 | 1500 | 51940.1 | 14005.0 | 48717.5 | 12793.5 |
Видно что PGO смогло ускорить оптимизацию #1 на целых 18%, для остальных прирост оказался скромнее.
Тут водятся драконы
Забыл упомянуть забавную особенность. Сначала я искал первое решение, и поставил return сразу после printf. Так вот это снижало производительность на порядок. Когда стал искать все решения, то производительность резко скакнула вверх.
__global__ void NaiveGPUKernel(const uint32_t N) {
const int a = blockIdx.x + 1;
const int b = threadIdx.x + blockIdx.y * blockDim.x + 1;
if (b >= a)
return;
for (uint32_t c = 1; c < b; ++c) {
for (uint32_t d = 1; d < c; ++d) {
const uint64_t sum = gpu_powers[a] + gpu_powers[b] + gpu_powers[c] + gpu_powers[d];
const auto s = gpu_binary_search(N, sum);
if (s > 0) {
printf("%d %d %d %d %d\n", a, b, c, d, s);
return; <- это дало замедление на порядок
}
}
}
}
Что можно сделать еще
Для CPU попробовать развернуть циклы, считать по два числа за раз, использовать векторные команды процессора.
Для GPU поиграть с порядком вычислений, попробовать по-экономить на регистрах, лучше подобрать параметры блоков и потоков.
Выводы
1. Писать на CUDA просто.
2. Без ухищрений, прямолинейный код на CUDA оказался быстрее CPU реализации в десятки раз.
3. В числодробильных задачах GPU сильно быстрее CPU за те же деньги.
4. GPU требуется много времени на «раскачку» и для совсем коротких задач ускорения может не быть.
5. Более сложный алгоритм может оказаться быстрее для CPU, но медленнее для GPU.
6. Профильная оптимизация хорошо ускоряет код и уменьшает размер файла одновременно.
Исходники
Традиционно можно потрогать код на GitHub
Комментарии (28)
dmitry_ch
22.12.2016 08:25Приготовился читать эпопею борьбы с CPU, GPU и компиляторами, попытки победить тот-самый-старый-баг (ну, у всего ПО есть старый-добрый-баг, который никак не добьют годами), и был приятно удивлен длиной текста. Действительно — сказка!
Самая короткая сказка о потерянном времени.
erwins22
22.12.2016 09:04+1одно из чисел будет кратно 11
drBasic
22.12.2016 09:04Есть доказательство?
erwins22
22.12.2016 09:34(0..10) ^5 mod 11== 0, 1, 10
0 0 0 1 10 =0
0 1 1 10 10 = 0
значит минимум 1 кратен 11erwins22
22.12.2016 13:30т.е. первый параметр берешь с шагом 11 0 11 22 33 итд
все остальное так же
переписываешь в виде e^5-a^5-b^5-c^5=d^5
e берешь с шагом в 11
ускорение в 5.5 раза 11/2
zuborg
22.12.2016 20:0427^5 + 84^5 + 110^5 + 133^5 = 144^5
здесь кратно 11 число 110.
55^5 + 3183^5 + 28969^5 + 85282^5 = 85359^5
а здесь — 55
e нельзя перебирать с шагом 11.
И заранее не известно, какое можно.erwins22
22.12.2016 20:15первое ты перебираешь с шагом в 11, а порядок остальных делаешь независимым от него.
zuborg
22.12.2016 20:57Как минимум надо два варианта — паребирать с шагом 11 значение или e, или одно из слагаемых.
Потом — оптимальнее всего перебирать от больших к меньшим (чтобы сузить диапазон самых вариабельных значений).
Перебор одного из слагаемых в полном диапазоне может легко съесть всю экономию.erwins22
22.12.2016 21:29можно разбить на 2 варианта
for (uint32_t a = from; a < to; ++a) {
if (a%11=0) is11=true else is11=false
for (uint32_t b = 1; b < a; ++b) {
if (b%11=0) is11=true
for (uint32_t c = 1; c < b; ++c) {
if (c%11=0) is11=true
if(is11)
{
for (uint32_t d = 11; d < c; d+=11) {
const uint64_t sum = powers[a] + powers[b] + powers[c] + powers[d];
if (std::binary_search(std::begin(powers), std::end(powers), sum)) {
const auto it = std::lower_bound(std::begin(powers), std::end(powers), sum);
uint32_t e = static_cast<uint32_t>(std::distance(std::begin(powers), it));
std::cout << a << " " << b << " " << c << " " << d << " " << e << "\n";
}
}
else
{
for (uint32_t d = 1; d < c; ++d) {
const uint64_t sum = powers[a] + powers[b] + powers[c] + powers[d];
if (std::binary_search(std::begin(powers), std::end(powers), sum)) {
const auto it = std::lower_bound(std::begin(powers), std::end(powers), sum);
uint32_t e = static_cast<uint32_t>(std::distance(std::begin(powers), it));
std::cout << a << " " << b << " " << c << " " << d << " " << e << "\n";
}
}
}
zuborg
22.12.2016 23:49Во первых, у Вас перепутано условие — если is11 == true, то перебирать d надо полностью, а не с шагом 11, т.к. одно из предыдущих слагаемых уже кратно 11, и на d такое ограничение уже не распространяется.
Во вторых если is11 == false — то перебирать d надо все равно полностью, т.к. нет гарантии что кратным 11 числом является не e, а d.
paluke
22.12.2016 11:56+1Если посмотреть вот так: a^5+b^5+c^5 = e^5-d^5
— то разность пятых степеней может принимать не любые значения по модулям 11, 31, 41 и 61… Так что проверив сумму a^5+b^5+c^5 по модулю 11*31*41*61, четвертый цикл примерно в 80% случаев можно вообще не выполнять…
erwins22
22.12.2016 09:48Ускорим во много раз
(a^5) mod 7 =
таблица соответствия
t f[t]
0 0
1 1
2 4
3 5
4 2
5 3
6 6
зная первые 4 последнее находим как (f[a]+f[b]+f[c]+f[d])mod 7 =f[e]
аналогично для 13
t R[r]
0 0
1 1
2 6
3 9
4 10
5 5
6 2
7 11
8 8
9 3
10 4
11 7
12 12
зная первые 4 последнее находим как (R[a]+R[b]+R[c]+R[d])mod 13 =R[e]
аналогично для 17
zuborg
22.12.2016 21:04Ускорим во много раз
Операция определения остатка весьма тормозная.
Настолько тормозная, что практически равна по скорости с поиском i^5 в таблице 5-х степеней. Если поиск сделан быстрым, то проверка остатка только добавит задержку, если она отсекает менее 90-95% вариантов.
mrMOD
22.12.2016 12:07+1if (std::binary_search(std::begin(powers), std::end(powers), sum)) { const auto it = std::lower_bound(std::begin(powers), std::end(powers), sum);
вы тут при попадании внутрь if второй раз бинарный поиск делаете (aka lower_bound), зачем?
нужно сделать 1 раз lower_bound и сравнить значение по итератору с искомым.paluke
22.12.2016 12:25Попадание внутрь if случается так редко, что это не оказывает влияния на производительность
Akon32
22.12.2016 14:54Получить по сравнению с одним ядром CPU 10-кратное ускорение на видеокарте с 1500 ядрами как-то слишком неэффективно. На других задачах у меня, как правило, получалось 20х на 200 ядрах, без особых оптимизаций. Согласитесь, 150-кратное замедление по сравнению с "теорией" — это много.
Нужно попробовать учесть специфику CUDA (выполнение программы варпами и некоторую нелюбовь к ветвлениям):
- Для определения, является ли число sum пятой степенью целого числа, вместо бинарного поиска можно бы попробовать
double r = pow(sum,0.2) ; bool isPower5 = abs((uint32)(r+0.5)-r) < EPS;Не уверен, что возведение в степень действительно будет быстрее, но вдруг. - Самое важное: нити в вашей программе получают блоки работы разного объёма, поэтому в конце будут работать лишь несколько (а не несколько тысяч, как поддерживает GPU) варпов, на которые свалилась бОльшая часть работы. Остальная часть GPU будет простаивать. Нужно разбить задачу на блоки примерно равного объёма, например, как-то пронумеровавать их, или обсчитывать меньшими кусками (вызывать ядро много раз), чтобы равномерно загрузить GPU. На OpenMP было бы достаточно написать
#pragma omp sсhedule(dynamic,n), на CUDA, наверно, можно попробовать рекурсивно вызывать ядро из того же ядра, CUDA вроде бы это сейчас поддерживает.
drBasic
22.12.2016 16:40Все правильно говорите. Можно сделать лучше.
Я же просто продемонстрировал что даже скопировав код, не углубляясь в тонкости GPU, можно получить ускорение по времени в 10 раз.Akon32
22.12.2016 17:16А попробуйте внешний цикл на CPU вынести?
не углубляясь в тонкости GPU
Это зря. Не рекомендую такой подход. У GPU совершенно другая архитектура, CPU-код очень плохо ложится на неё. Если задача не разбита на несколько тысяч одновременно работающих нитей (в несколько раз больше, чем ядер) — видеокарта будет работать неэффективно. Туда же относятся особенности ветвлений.
В общем случае, архитектуру "железа" следует иметь в виду.
(Конечно, есть случаи, когда 10х будет достаточно, тогда можно просто скопировать код. Но чаще получить 200х вместо 10х будет лучше)
vlanko
27.12.2016 14:51Там внутренний цикл 64 бита. Производительность 770 в 64-битных вычислениях около 145 ГФлопс.
У Core i5-4440 99 ГФлопс с AVX. Так что чудес от видеокарты не ждите.
- Для определения, является ли число sum пятой степенью целого числа, вместо бинарного поиска можно бы попробовать
nikolay_karelin
23.12.2016 16:46А для GPU не пытались проверить полосу (bandwidth) и сравнить с максимальной?
Я боюсь, что задачка для CUDA не очень показательная: объем данных небольшой, проблем с обращениями к памяти не будет.
zuborg
До какого N возможно проверить решения, скажем, за час работы?
drBasic
Сложность алгоритма O(n^4). И эта теория отлично согласуется с практикой, увеличение N в 2 раза увеличивает время в 16 раз.
Если взять лучшее время 12.8 секунд для CPU #2 на размере N=1500, то увеличение N в 4 раза увеличит время в 256 раз, до 54,5 минуты.
Т.е. примерно за час на CPU можно посчитать для N=6000.