Многопоточная сказка о потерянном времени

    В публикации «Сказка о потерянном времени» пользователь crea7or рассказал, как он опровергал Гипотезу Эйлера на современном CPU.

    Мне же было интересно узнать как покажет себя GPU, и я сравнил однопоточный код с многопоточным для CPU и совсем многопоточным для GPU, с помощью архитектуры параллельных вычислений CUDA.

    Железо и компиляторы
    CPU Core i5 4440
    GPU GeForce GTX 770
    MSVS 2015 update 3
    tookit CUDA 8.
    Конечно, сборка была релизная и 64 битная, с оптимизацией /02 и /LTCG.
    Отключение сброса видеоадаптера
    Так как вычисления длятся более двух секунд, видеодрайвер завершает CUDA-программу. Чтобы этого не происходило, нужно указать достаточное время до сброса через ключ реестра
    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
    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 29

      0
      До какого N возможно проверить решения, скажем, за час работы?
        0
        Сложность алгоритма O(n^4). И эта теория отлично согласуется с практикой, увеличение N в 2 раза увеличивает время в 16 раз.
        Если взять лучшее время 12.8 секунд для CPU #2 на размере N=1500, то увеличение N в 4 раза увеличит время в 256 раз, до 54,5 минуты.
        Т.е. примерно за час на CPU можно посчитать для N=6000.
        0
        Приготовился читать эпопею борьбы с CPU, GPU и компиляторами, попытки победить тот-самый-старый-баг (ну, у всего ПО есть старый-добрый-баг, который никак не добьют годами), и был приятно удивлен длиной текста. Действительно — сказка!

        Самая короткая сказка о потерянном времени.
          +1

          одно из чисел будет кратно 11

            0
            Есть доказательство?
              0

              (0..10) ^5 mod 11== 0, 1, 10


              0 0 0 1 10 =0
              0 1 1 10 10 = 0
              значит минимум 1 кратен 11

                0

                т.е. первый параметр берешь с шагом 11 0 11 22 33 итд
                все остальное так же


                переписываешь в виде e^5-a^5-b^5-c^5=d^5
                e берешь с шагом в 11


                ускорение в 5.5 раза 11/2

                  0
                  27^5 + 84^5 + 110^5 + 133^5 = 144^5
                  здесь кратно 11 число 110.

                  55^5 + 3183^5 + 28969^5 + 85282^5 = 85359^5
                  а здесь — 55

                  e нельзя перебирать с шагом 11.
                  И заранее не известно, какое можно.
                    0

                    первое ты перебираешь с шагом в 11, а порядок остальных делаешь независимым от него.

                      0
                      Как минимум надо два варианта — паребирать с шагом 11 значение или e, или одно из слагаемых.
                      Потом — оптимальнее всего перебирать от больших к меньшим (чтобы сузить диапазон самых вариабельных значений).
                      Перебор одного из слагаемых в полном диапазоне может легко съесть всю экономию.
                        0

                        можно разбить на 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";
                        }
                        }


                        }

                          0
                          Во первых, у Вас перепутано условие — если is11 == true, то перебирать d надо полностью, а не с шагом 11, т.к. одно из предыдущих слагаемых уже кратно 11, и на d такое ограничение уже не распространяется.
                          Во вторых если is11 == false — то перебирать d надо все равно полностью, т.к. нет гарантии что кратным 11 числом является не e, а d.
              +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% случаев можно вообще не выполнять…
              0

              Ускорим во много раз
              (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

                0
                Ускорим во много раз

                Операция определения остатка весьма тормозная.
                Настолько тормозная, что практически равна по скорости с поиском i^5 в таблице 5-х степеней. Если поиск сделан быстрым, то проверка остатка только добавит задержку, если она отсекает менее 90-95% вариантов.
                +1
                 if (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 и сравнить значение по итератору с искомым.
                  0
                  Попадание внутрь if случается так редко, что это не оказывает влияния на производительность
                  0
                  А как быстро будет работать 128-битная версия? Такие вычисления в GPU вообще предусмотрены?
                  А то 64 бита заканчиваются очень быстро.
                    0
                    128bit можно: http://stackoverflow.com/questions/6162140/128-bit-integer-on-cuda
                      0

                      зачем? достаточно float
                      ДробнаяЧасть( r1/p1+r2/p2+r3/p3+.....)<1/sqrt(p1p2p3*… ,5)


                      f_p — таблица преобразования по p


                      r_p= F_p^-1[F_p[a]+F_p[b]+F_p[c]+F_p[d]]

                        0
                        Что такое rX и pX?
                          0

                          p — простое число для которого 5 степень биективна
                          r_p остаток деления на p

                            0
                            А можете объяснить Вашу идею более детально?
                            Сорри, но из Ваших формул мне вообще не понятно что с чем сравнивать и почему это будет работать (
                    0

                    Получить по сравнению с одним ядром 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 вроде бы это сейчас поддерживает.
                      0
                      Все правильно говорите. Можно сделать лучше.
                      Я же просто продемонстрировал что даже скопировав код, не углубляясь в тонкости GPU, можно получить ускорение по времени в 10 раз.
                        0

                        А попробуйте внешний цикл на CPU вынести?


                        не углубляясь в тонкости GPU

                        Это зря. Не рекомендую такой подход. У GPU совершенно другая архитектура, CPU-код очень плохо ложится на неё. Если задача не разбита на несколько тысяч одновременно работающих нитей (в несколько раз больше, чем ядер) — видеокарта будет работать неэффективно. Туда же относятся особенности ветвлений.


                        В общем случае, архитектуру "железа" следует иметь в виду.


                        (Конечно, есть случаи, когда 10х будет достаточно, тогда можно просто скопировать код. Но чаще получить 200х вместо 10х будет лучше)

                        0
                        Там внутренний цикл 64 бита. Производительность 770 в 64-битных вычислениях около 145 ГФлопс.
                        У Core i5-4440 99 ГФлопс с AVX. Так что чудес от видеокарты не ждите.
                        0

                        А для GPU не пытались проверить полосу (bandwidth) и сравнить с максимальной?


                        Я боюсь, что задачка для CUDA не очень показательная: объем данных небольшой, проблем с обращениями к памяти не будет.

                          0
                          Так просто, замечание. Мне стоит самому попробовать с вашим кодом, конечно. Использованный пул потоков (https://github.com/progschj/ThreadPool/blob/master/ThreadPool.h) довольно неэффективен, безотносительно зависимости потоков по данным. Немного больше труда составит изготовить пул потоков с task stealing. Как и почему — очень хорошо есть здесь:


                          Можно прямо код брать.

                          Only users with full accounts can post comments. Log in, please.