О пирамидальном сложении на параллельной архитектуре

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

Для второго шага reduce обычно используется пирамидальная схема: сначала в N/2 потоков сложим N результатов попарно, затем сложим N/2 в N/4 и так далее. Число итераций равно, очевидно, log2N.

Возникает вопрос, «сколько данных складывать на каждой итерации?» Ведь можно складывать в N/4-N/16-N/256 кучек, можно по 1/8-64-512 и так далее. Из общих соображений, складывать по несколько лучше чем по два. Конечно, потоков получается меньше, но меньше и оверхед на создание-завершение потока.

Для NVidia CUDA идея "делать не по 2", выбирая все динамически, оказалась очень плохой. Да, с одной стороны мы действительно имеем оптимум при сложении по 8 или по 16. С другой стороны, код для вычисления содержит больше условного исполнения, отчего все ухудшается:

  • С одной стороны, "умный" код для сложения по 8 примерно втрое быстрее, чем при сложении им же "по 2".
  • С другой стороны, код рассчитанный только на сложение по 2 - в полтора раза быстрее чем умный складывает по восемь
Приведу на всякий случай быстрый код, взятый из примера scalarprod (CUDA SDK)
for(int stride = N / 2; stride > 0; stride >>= 1){
     __syncthreads();
     for(int iAccum = threadIdx.x; iAccum < stride; iAccum += blockDim.x)
               accumResult[iAccum] += accumResult[stride + iAccum];
}
Здесь вложенный for - это на самом деле такой IF, который для рассматриваемого случая, когда N == blockDim.x, удовлетворяется только для части потоков, а выполняется тело цикла для этих потоков только один раз. Свой "умный" код не привожу, слишком уж умный.

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

for(int stride = blockDim.x / 4; stride > 0; stride >>=2){
  __syncthreads();
   for(int iAccum = threadIdx.x; iAccum < stride;
         iAccum += blockDim.x)
        {
                        data[iAccum] += data[stride + iAccum];
                        data[iAccum] += data[stride*2 + iAccum];
                        data[iAccum] += data[stride*3 + iAccum];
        }
}
Да, дальнейший unroll дает падение производительности.

Comments

Вообще это широкоизвестный алгоритм. Но есть и более дешевые в плане сложности - поищите различные реализации Reduce для MPI.

По поводу возможности использование GPU для проведения параллельных вычислений - реальное приложение GPGPU.

Перед мной стоит задача задействовать ресурсы графического процессора для проведения параллельных вычислений.

У меня решается задача, требующая очень большого количества массивных вычислений.

На данный момент существует программа, позволяющее распараллеливать проведение расчетов на ПРОИЗВОЛЬНОЕ количество процессов.

Программа позволяет также осуществлять запуск этих нескольких процессов на разных ядрах в многоядерных CPU с целью увеличения производительности вычислений.

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

Для решения этой задачи, по всей видимости, мне нужно будет создание программы-конвертера, которая

- работала бы в фоновом режиме, во время работы расчетной программы

- позволяла бы осуществлять перехватывание команд между расчетной программой и CPU

- перенаправляла бы вычисления из CPU в GPU.

Сложность заключается в том, что все это нужно мне сделать БЕЗ МОДИФИКАЦИИ ИСХОДНОГО ТЕКСТА ПРОГРАММЫ из-за того, что разработчик уехал, а работоспособного кода не осталось, т.к. полетел винт на его бывшей машине.

Поэтому, мне хотелось бы узнать по поводу способов решения этой задачи!

Мне кажется, что наиболее эффективным способом решения этой проблемы является задействование ресурсов GPU с помощью использования ШЕЙДЕРОВ на языке их описания OpenGL Shading Language (GLSL).

В синтаксисе GLSL существуют возможность выполнять в GPU такие же никоуровневые команды как и в ассемблере для CPU (ADD, SUB и т.д.).

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

Таким образом, необходимо решить 2 следующие задачи:

1. Найти способ перехватывания инструкций между расчетной программой и CPU для каждого из запущенных процессов.

2. Найти способ перенаправлять вычисления из CPU в GPU и в обратную сторону.

Было бы интересно услышать, какие способы можно использовать для решения описанной выше задачи и какие методы решения существуют.

Понимаете какое дело, если вы будете перехватывать команды CPU по одной, то никакого выигрыша не получите, только проигрыш.

Более того, бинарная трансляция (предположим что мы ее сделали) на уровне отдельных команд тоже не спасет, архитектура принципиально разная.

Придется переписывать.

Дело в том, что программа состоит из нескольких исполняемых модулей, которые выполняют следующие функции:

1. PreProcess модуль для запуска интерфейса подготовки исходных данных (компонентный состав, давление, температура и т.д.)

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

3. Solver программа-решатель, которая запускается из под модуля Setup в заданном количестве экземпляров-процессов.

4. Monitor программа работающая в фоновом режиме и работающая арбитром при раздаче очередных данных для каждого потока (экземпляры программы решателей Solver)

5.PostProccess модуль для вывода расcчитанных данных.

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

В принципе, при помощи инструментов из MS DDK можно эмулировать заданное количество процессоров, а затем использовать их для отладки многопоточных приложений.

В то же время, набор низкоуровневых команд GPU не аналогичен набору команд CPU.

Значит, нужен системный драйвер для эмуляции процессора это раз.

И механизм, который мог бы обеспечить исполнение основных команд на центральном процессоре, а арифмитических на GPU.

В GLSL есть возможность использования низкоуровневых шейдеров, которые могут оперировать с неким подобием ассемблерного кода, выполняя операции ADD, SUB и некотрые др. (см. гл.20, стр. 877 из кники OpenGL- Суперкнига, 3-е изд. под ред. Ричарда Райта-мл. и Бенджамина Липчача)

Я дисассемблировал исполняемый файл расчетной программы, о которой идет речь.

Можно ознакомиться с первыми несколькими страницами asm-овского кода, который я выложил на нижеприведенном ресурсе:

http://downfile.narod.ru/solver.txt

Программа действительно распараллеливает вычисления на ПРОИЗВОЛЬНОЕ количество процессов.

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

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

Вопрос заключается, собственно говоря, в том как это можно сделать?

Ну вы не можете кормить GPU "по одной команде".

И архитектурно не можете, ну да это ладно, но и по факту это не имеет смысла. Ну допустим вы выполнили одну команду на CPU, а вторую хотите на GPU, данные так и будем по байтику пересылать ?

Садитесь и пишите свой Solver для GPU, деваться вам некуда.

Добрый день Алексей,
Затронутая в статье проблема относится к классу проблем редукции. Интересное решение дал Mark Harris в SC07_CUDA_5_Optimization_Harris.pdf.
Развивая эти идеи я пришел к возможности реализации некоторых из методов метапрограммирования в CUDA, в частности удалось пирамиду представить в виде квази-рекурсии (здесь слово квази употреблено поскольку настоящей рекурсии нет, а используются новые instance той же функции). Детали смотрите в http://www.gamedev.ru/code/forum/?id=80317

Да, текст от Харриса реально вставляет.

Но вот что я жутко не люблю, так это современные приемы оптимизации, когда руками анроллим, руками проверяем крайние случаи, руками реализуем логарифмическую пирамиду. Дурдом. Программа должна быть красивой.

В тот-то и фокус что не руками. Это все можно поручить компилятору. Вот так

// f = 2^log2(N) ÷
template
struct Exp2Log2{ enum { value = 2*Exp2Log2::value };};
template<>
struct Exp2Log2<1> { enum { value = 1 };};

// RecursivSum, " "
// sdata blockSize
template
__device__ void RecursivSum(__shared__ myType *sdata)
{
if (threadIdx.x < blockSize-blockSizeAligned)
sdata[threadIdx.x] += sdata[threadIdx.x + blockSizeAligned];
__syncthreads();
// instance , blockSize, blockSizeAligned/2:
if (blockSizeAligned >= 2) //
RecursivSum(sdata);
};

// g_idata n, g_odata
template
__global__ void VectorSum(const myType *g_idata, myType *g_odata, const unsigned int n)
{
__shared__ int sdata[blockSize];
sdata[threadIdx.x] = 0;
for (unsigned int i = threadIdx.x; i < n; i += blockSize)
sdata[threadIdx.x] += g_idata[i];
__syncthreads();
RecursivSum::value>(sdata);
if (threadIdx.x == 0) g_odata[0] = sdata[0];
}

Что-то здесь двиг форума исказил код программы, в template нет параметров. Правильный код смотрите в http://www.gamedev.ru/code/forum/?id=80317.

Ну да, такие вещи обычно на препроцессоре делаются, но темплейтер ничем не хуже (считать ли параметризацию с типом и авто-fallback удобнее или нет - дело вкуса)

>Ну да, такие вещи обычно на препроцессоре делаются, но темплейтер ничем не хуже

"обычно"? в CUDA?

Очень интересно. Не могли бы Вы дать ссылку на CUDA программу с таким примером?

Не в CUDA, а вообще вычислителями.

Но насколько я помню исходники CuBLAS, там все в макросах (я туда сунул нос и сбежал)

Исходники CuBLAS я тоже смотрел. Тоже не понравилось. Те же макросы.
Может быть каждая из CuBLAS функций и эффективна, но сама концепция порочна. Программу для GPU нельзя собирать как кубики из CuBLAS. Она по определению будет не эффективной.
Оказывается только вызов пустой глобал функции стоит 0.16 мсек. Имхо, нвидиа разработала CuBLAS на скорую руку, причем в Си стиле и предназначает в основном для использования в старых фортран-пакетах.

Empty call на 160 микросекунд - это что-то очень странное. Люди получают цифры в районе 10-30 usec. Подозреваю, что вы неправильно меряете.

С другой стороны, гонять данные на видеокарту, чтобы посчитать 2+2 или даже просто сумму элементов вектора - категорически невыгодно. Вы же этот вектор туда льете со скоростью, ну три гигабайта в секунду, ну может 4 (надо, кстати, на новой материнке померять). А на CPU это будет молотить со скоростью памяти т.е. в разы выше.

Конечно, всякий анроллинг и "рекурсию" делать через макроподстановки (неважно, препроцессором или темплейтером) - удобно.

А что касается CuBLAS - так это паттерн такой. Берем 200 мегабайтные матрицы и умножаем, умножаем, умножаем...

Мерил Cuda функциями
CUT_SAFE_CALL(cutCreateTimer(&timer));
CUT_SAFE_CALL(cutStartTimer(timer));
//......................
CUT_SAFE_CALL(cutStopTimer(timer));
CUT_SAFE_CALL(cutDeleteTimer(timer));
Такой результат можно получить на любой глобал функции (т.е. вызываемой из хоста, но работающей в девайсе).

Ну вот я взял пример template из CUDA SDK 2.0, перенес там startTimer прямо перед вызовом kernel, stopTimer - сразу после.

16 микросекунд - если пустой kernel. В 10 раз меньше чем у вас.
kernel с работой - 30 микросекунд.

Может быть вы вместе с заливкой-выливкой данных меряете время ?

Вот такой тест

__global__ void Empty(void)
{
}

вызывается из

void runTimeTest(void)
{
CUT_DEVICE_INIT();
unsigned int n =2;

// create and start timer
unsigned int timer = 0;
CUT_SAFE_CALL(cutCreateTimer(&timer));
CUT_SAFE_CALL(cutStartTimer(timer));

// setup execution parameters
dim3 threads(2);
dim3 grid(n / threads.x);

Empty <<< grid, threads >>>();

// check if kernel execution generated and error
CUT_CHECK_ERROR("Kernel execution failed");

// stop and destroy timer
CUT_SAFE_CALL(cutStopTimer(timer));
printf("Processing time: %f (ms) Bandwidth %f gb/sek \n", cutGetTimerValue(timer),
1000*(4*n)/(1024*1024*1024*cutGetTimerValue(timer)));
};

У меня в релиз делается за 18 ms.

Милли или микро ?

Хороший вопрос.
Метод измерения времени был без изменений взят из NVIDIA CUDA SDK\projects\matrixMul
На сколько я могу судить (100% гарантии дать не могу), там речь идет все же о миллисекундах.

Ну вот если у matrixMul мерять время только по kernel, без копирования данных туда-сюда, то для матрицы 16x16 получается 0.035 тамошних
единиц т.е. 35 микросекунд.

Сейчас в matrixMul.h поставил соответствующие параметры
//matrixMul.h
#ifndef _MATRIXMUL_H_
#define _MATRIXMUL_H_

// Thread block size
#define BLOCK_W_SIZE 16
#define BLOCK_H_SIZE 16

// Matrix dimensions
// (chosen as multiples of the thread block size for simplicity)
#define WA (1 * BLOCK_W_SIZE) // Matrix A width
#define HA (1 * BLOCK_H_SIZE) // Matrix A height
#define WB (1* BLOCK_H_SIZE ) // Matrix B width
#define HB WA // Matrix B height
#define WC WB // Matrix C width
#define HC HA // Matrix C height

#define NBX (WC / BLOCK_W_SIZE)

#endif // _MATRIXMUL_H_

Т.е. на один блок 16х16, показал Processing time: 0.063115 (ms).
Карта gf8800gt.
Значит все-таки
CUT_SAFE_CALL(cutStartTimer(timer));
//......................
CUT_SAFE_CALL(cutStopTimer(timer));
дают время в миллисекундах.

В топике 10.05.08 02:27 дан конечно совершенно абсурдный результат 18 мс, я конечно это понимал, поэтому и привел весь код.
Почему так хотел бы знать.

Черт, опять забыл про двиг сайта. Нужно конечно вызывать
Empty {{{ grid, threads }}}();
Здесь я угловые скобки заменил фигурными.

Впрочем, некоторые из CuBLAS функций могут быть полезны (ну там инициализация, аллокаторы, диагностика ошибок и т.п.), но "тяжелые" матричные и векторные вычисления наверное все же лучше делать самому, причем как девайс функции.

Ну да, и если вы больше ничего не меняли, то эти 60 микросекунд - это вместе с I/O

Там было только с выводом, без вывода 40.

Там было только с выводом, без вывода 40.

Ну вот, 40, а не 160.

помоги пожалуйста с одним вопросом:
есть 2 потока. Как вручную распределить потоки по разным ядрам?
Среда программирования Delphi
процессор 2-х или 4-х ядерный

В дельфи - понятия не имею. Вообще в Win32 - через SetThreadAffinityMask