Нет ускорения для векторных сумм с потоками

У меня есть программа на С++, которая в основном выполняет некоторые вычисления матрицы. Для них я использую LAPACK/BLAS и обычно ссылаюсь на MKL или ACML в зависимости от платформы. Многие из этих матричных вычислений работают на разных независимых матрицах, и поэтому я использую std:: thread, чтобы эти операции выполнялись параллельно. Тем не менее, я заметил, что у меня нет ускорения при использовании большего количества потоков. Я проследил эту проблему до простой процедуры Blas. Кажется, что если два потока используют эту процедуру параллельно, каждый поток занимает в два раза больше времени, хотя оба потока работают на разных массивах.

Следующее, что я попробовал, - это написать новый простой метод для выполнения векторных дополнений для замены подпрограммы daxpy. С помощью одного потока этот новый метод выполняется так же быстро, как и процедура BLAS, но при компиляции с gcc он испытывает те же проблемы, что и процедура BLAS: удвоение числа потоков, выполняемых параллельно, также удваивает количество времени, которое требуется каждому потоку, поэтому ускорение не достигается. Однако с использованием компилятора Intel С++ эти проблемы исчезают: с увеличением числа потоков время, которое требуется одному потоку, постоянно.

Однако мне нужно также скомпилировать системы, где нет компилятора Intel. Поэтому мои вопросы: почему нет ускорения работы с gcc и есть ли возможность улучшить производительность gcc?

Я написал небольшую программу для демонстрации эффекта:

// $(CC) -std=c++11 -O2 threadmatrixsum.cpp -o threadmatrixsum -pthread

#include <iostream>
#include <thread>
#include <vector>

#include "boost/date_time/posix_time/posix_time.hpp"
#include "boost/timer.hpp"

void simplesum(double* a, double* b, std::size_t dim);

int main() {

    for (std::size_t num_threads {1}; num_threads <= 4; num_threads++) {
        const std::size_t N { 936 };

        std::vector <std::size_t> times(num_threads, 0);    

        auto threadfunction = [&](std::size_t tid)
        {
            const std::size_t dim { N * N };
            double* pA = new double[dim];
            double* pB = new double[dim];

            for (std::size_t i {0}; i < N; ++i){
                pA[i] = i;
                pB[i] = 2*i;
            }   

            boost::posix_time::ptime now1 = 
                boost::posix_time::microsec_clock::universal_time();    

            for (std::size_t n{0}; n < 1000; ++n){
                simplesum(pA, pB, dim);
            }

            boost::posix_time::ptime now2 = 
                boost::posix_time::microsec_clock::universal_time(); 
            boost::posix_time::time_duration dur = now2 - now1; 
            times[tid] += dur.total_milliseconds(); 
            delete[] pA;
            delete[] pB;
        };

        std::vector <std::thread> mythreads;

        // start threads
        for (std::size_t n {0} ; n < num_threads; ++n)
        {
            mythreads.emplace_back(threadfunction, n);
        }

        // wait for threads to finish
        for (std::size_t n {0} ; n < num_threads; ++n)
        {
            mythreads[n].join();
            std::cout << " Thread " << n+1 << " of " << num_threads 
                      << "  took " << times[n]<< "msec" << std::endl;
        }
    }
}

void simplesum(double* a, double* b, std::size_t dim){

    for(std::size_t i{0}; i < dim; ++i)
    {*(++a) += *(++b);}
}

Outout с gcc:

Thread 1 of 1  took 532msec
Thread 1 of 2  took 1104msec
Thread 2 of 2  took 1103msec
Thread 1 of 3  took 1680msec
Thread 2 of 3  took 1821msec
Thread 3 of 3  took 1808msec
Thread 1 of 4  took 2542msec
Thread 2 of 4  took 2536msec
Thread 3 of 4  took 2509msec
Thread 4 of 4  took 2515msec

Вывод с icc:

Thread 1 of 1  took 663msec
Thread 1 of 2  took 674msec
Thread 2 of 2  took 674msec
Thread 1 of 3  took 681msec
Thread 2 of 3  took 681msec
Thread 3 of 3  took 681msec
Thread 1 of 4  took 688msec
Thread 2 of 4  took 689msec
Thread 3 of 4  took 687msec
Thread 4 of 4  took 688msec

Итак, с icc время, необходимое для одного потока, выполняется, вычисления постоянны (как я и ожидал, мой процессор имеет 4 физических ядра), а с gcc время для одного потока увеличивается. Замена процедуры simplesum BLAS:: daxpy дает те же результаты для icc и gcc (неудивительно, так как большинство времени проводится в библиотеке), которые почти такие же, как и указанные выше результаты gcc.

Ответы

Ответ 1

Ответ довольно прост: ваши потоки борются за пропускную способность памяти!

Считайте, что вы выполняете одно добавление с плавающей запятой на 2 магазина (одна инициализация, одна после добавления) и 2 чтения (в дополнение). Большинство современных систем, обеспечивающих множественный процессор, фактически должны совместно использовать контроллер памяти из нескольких ядер.

В системе с двумя физическими гнездами процессора и 12 ядрами (24 с HT) выполнялось следующее. В вашем исходном коде показана ваша проблема:

Thread 1 of 1  took 657msec
Thread 1 of 2  took 1447msec
Thread 2 of 2  took 1463msec
[...]
Thread 1 of 8  took 5516msec
Thread 2 of 8  took 5587msec
Thread 3 of 8  took 5205msec
Thread 4 of 8  took 5311msec
Thread 5 of 8  took 2731msec
Thread 6 of 8  took 5545msec
Thread 7 of 8  took 5551msec
Thread 8 of 8  took 4903msec

Однако, просто увеличивая арифметическую плотность, мы можем наблюдать значительное увеличение масштабируемости. Чтобы продемонстрировать, я изменил вашу процедуру добавления, чтобы также выполнить возведение в степень: *(++a) += std::exp(*(++b));. Результат показывает почти идеальное масштабирование:

Thread 1 of 1  took 7671msec
Thread 1 of 2  took 7759msec
Thread 2 of 2  took 7759msec
[...]
Thread 1 of 8  took 9997msec
Thread 2 of 8  took 8135msec
Thread 3 of 8  took 10625msec
Thread 4 of 8  took 8169msec
Thread 5 of 8  took 10054msec
Thread 6 of 8  took 8242msec
Thread 7 of 8  took 9876msec
Thread 8 of 8  took 8819msec

Но как насчет ICC?

Во-первых, ICC inlines simplesum. Доказательство того, что вложение происходит просто: с помощью icc я отключил многопрофильную межпроцедурную оптимизацию и переместил simplesum в свою собственную единицу перевода. Разница поразительна. Производительность от

Thread 1 of 1  took 687msec
Thread 1 of 2  took 688msec
Thread 2 of 2  took 689msec
[...]
Thread 1 of 8  took 690msec
Thread 2 of 8  took 697msec
Thread 3 of 8  took 700msec
Thread 4 of 8  took 874msec
Thread 5 of 8  took 878msec
Thread 6 of 8  took 874msec
Thread 7 of 8  took 742msec
Thread 8 of 8  took 868msec

Для

Thread 1 of 1  took 1278msec
Thread 1 of 2  took 2457msec
Thread 2 of 2  took 2445msec
[...]
Thread 1 of 8  took 8868msec
Thread 2 of 8  took 8434msec
Thread 3 of 8  took 7964msec
Thread 4 of 8  took 7951msec
Thread 5 of 8  took 8872msec
Thread 6 of 8  took 8286msec
Thread 7 of 8  took 5714msec
Thread 8 of 8  took 8241msec

Это уже объясняет, почему библиотека плохо работает: ICC не может встроить ее, и, следовательно, независимо от того, что еще может заставить ICC работать лучше, чем g++, этого не произойдет.

Он также дает подсказку относительно того, что ICC может делать прямо здесь... Что, если вместо выполнения simplesum 1000 раз, меняет местами петли, чтобы он

  • Загружает два дубликата
  • Добавляет их 1000 раз (или даже выполняет a = 1000 * b)
  • Сохраняет два двухлокальных номера

Это увеличит арифметическую плотность без добавления каких-либо экспонент к функции... Как это доказать? Ну, для начала давайте просто реализовать эту оптимизацию и посмотреть, что произойдет! Для анализа мы рассмотрим производительность g++. Напомним наши результаты:

Thread 1 of 1  took 640msec
Thread 1 of 2  took 1308msec
Thread 2 of 2  took 1304msec
[...]
Thread 1 of 8  took 5294msec
Thread 2 of 8  took 5370msec
Thread 3 of 8  took 5451msec
Thread 4 of 8  took 5527msec
Thread 5 of 8  took 5174msec
Thread 6 of 8  took 5464msec
Thread 7 of 8  took 4640msec
Thread 8 of 8  took 4055msec

И теперь, давайте обменяться

for (std::size_t n{0}; n < 1000; ++n){
    simplesum(pA, pB, dim);
}

с версией, в которой внутренний цикл был сделан внешним контуром:

double* a = pA; double* b = pB;
for(std::size_t i{0}; i < dim; ++i, ++a, ++b)
{
    double x = *a, y = *b;
    for (std::size_t n{0}; n < 1000; ++n)
    {
        x += y;
    }
    *a = x;
}

Результаты показывают, что мы находимся на правильном пути:

Thread 1 of 1  took 693msec
Thread 1 of 2  took 703msec
Thread 2 of 2  took 700msec
[...]
Thread 1 of 8  took 920msec
Thread 2 of 8  took 804msec
Thread 3 of 8  took 750msec
Thread 4 of 8  took 943msec
Thread 5 of 8  took 909msec
Thread 6 of 8  took 744msec
Thread 7 of 8  took 759msec
Thread 8 of 8  took 904msec

Это доказывает, что оптимизация обмена петлями действительно является основным источником превосходных характеристик ICC здесь.

Обратите внимание, что ни один из проверенных компиляторов (MSVC, ICC, g++ и clang) не заменит цикл умножением, что улучшает производительность на 200x в однопоточном и 15x в 8-поточных случаях. Это связано с тем, что численная неустойчивость повторных дополнений может приводить к дико отличающимся результатам при замене одним умножением. При тестировании с целыми типами данных вместо типов данных с плавающей запятой эта оптимизация происходит.

Как мы можем заставить g++ выполнить эту оптимизацию?

Интересно, что истинный убийца для g++ не является невозможностью выполнить обмен петлями. Когда вызывается с -floop-interchange, g++ может также выполнить эту оптимизацию. Но только когда шансы значительно сложены в свою пользу.

  • Вместо std::size_t все оценки были выражены как int s. Не long, а не unsigned int, но int. Мне все еще трудно поверить, но, похоже, это сложное требование.

  • Вместо инкрементных указателей индексируйте их: a[i] += b[i];

  • g++ нужно сообщить -floop-interchange. Простой -O3 недостаточно.

Когда все три критерия выполнены, производительность g++ похожа на то, что обеспечивает ICC:

Thread 1 of 1  took 714msec
Thread 1 of 2  took 724msec
Thread 2 of 2  took 721msec
[...]
Thread 1 of 8  took 782msec
Thread 2 of 8  took 1221msec
Thread 3 of 8  took 1225msec
Thread 4 of 8  took 781msec
Thread 5 of 8  took 788msec
Thread 6 of 8  took 1262msec
Thread 7 of 8  took 1226msec
Thread 8 of 8  took 820msec

Примечание. Версия g++, используемая в этом эксперименте, - 4.9.0 для x64 Arch linux.

Ответ 2

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

summutex.lock();
simplesum(pA, pB, dim);
summutex.unlock();

решает проблему промахов кэша, но, конечно, не дает оптимального ускорения. Во всяком случае, поскольку теперь другие потоки заблокированы, метод simplesum может также использовать все доступные потоки для суммы

void simplesum(double* a, double* b, std::size_t dim, std::size_t numberofthreads){

    omp_set_num_threads(numberofthreads);
    #pragma omp parallel 

    {
    #pragma omp for
    for(std::size_t i = 0; i < dim; ++i)
    {
      a[i]+=b[i];

    }
    }
}

В этом случае все потоки работают на одном блоке в памяти: он должен находиться в кэше процессора, и если процессору необходимо загрузить некоторые другие части памяти в свой кеш, другие потоки извлекают выгоду из этого всего колодца (в зависимости от будь то кеш L1 или L2, но я считаю, что детали не имеют большого значения для этого обсуждения).

Я не утверждаю, что это решение идеально или почти оптимально, но, похоже, оно работает намного лучше, чем исходный код. И он не полагается на некоторые трюки переключения циклов, которые я не могу сделать в своем фактическом коде.