Десять простых шагов к стройной прикладной программе для гибридного вычислительного комплекса «К-100».
- Введение.
- Шаг 1. Выбираем вариант инструментальной среды. Осваиваем команды трансляции и запуска.
- Шаг 2. Учимся писать на С «настоящие Фортрановские программы».
- Шаг 3. Выбираем численный метод для реализации в качестве модельного гибридного приложения.
- Шаг 4. Пишем параллельную программу без использования ускорителей CUDA.
- Шаг 5. Выделяем в программе вычислительное ядро.
- Шаг 6. Реализуем коммуникации между процессором и ускорителем в пределах вычислительного узла.
- Шаг 7. Переписываем текст вычислительного ядра для процессора с использованием алгоритма, эффективно реализуемого на ускорителе.
- Шаг 8. Проектируем вычислительное ядро для мультипроцессора ускорителя CUDA.
- Шаг 9. Реализуем вычислительное ядро для ускорителя.
- Шаг 10. Переводим дух и думаем, что дальше.
- Приложение 1. Обзор базового программного обеспечения разработки гибридных параллельных приложений, доступного на K-100.
- Приложение 2. Что читать, чтобы не только успешно разбирать демонстрационные примеры, но и программировать самостоятельно.примеры, но и программировать самостоятельно.
- Исходные тексты примеров
Введение.
Настоящий документ сознательно написан в жанре «HOWTO», то есть объединяет в себе все основные сведения о разработке прикладных программ для К-100, без отсылок к отдельным руководствам по разным аспектам рассматриваемого вопроса. Тем не менее, изучение этого документа все же предполагает некоторый объем предварительных знаний у читателя.
В настоящем документе НЕ РАССМАТРИВАЮТСЯ вопросы организации удаленного доступа к ГВК, вроде порядка редактирования текстовых файлов, их передачи с локального компьютера на ГВК и обратно. Не рассматриваются также вопросы взаимодействия пользователя с системой очередей (СУППЗ), установленной на К-100. Предполагается, что этими знаниями и навыками читатель уверенно владеет, почерпнув их, при необходимости, из других документов.
Также предполагается, что читатель уверенно программирует на С и, скорее всего, имеет некоторое представление о разработке параллельных программ для традиционного вычислительного кластера.
Шаг 1. Выбираем вариант инструментальной среды. Осваиваем команды трансляции и запуска.
К-100 оснащен большим набором трансляторов стандартно используемых языков (С, С++, Фортран), а также большим набором коммуникационных библиотек, отличающихся как по своим возможностям, так и по конкретному исполнению. Имеется также несколько вариантов выбора средств программирования ускорителей CUDA, которыми оснащены вычислительные узлы. Разные коммуникационные библиотеки (и их варианты исполнения) используют разные коммуникационные сети для взаимодействия вычислительных модулей ГВК в процессе счета. В настоящем документе предполагается использование коммуникационных библиотек MPI и OSHCC. OSHCC или OSH++ - это Shmem с открытым исходным кодом. Ранее используемая библиотека Shmem (продукт фирмы Qlogic) имела закрытый исходный код, поэтому когда в ней что-то сломалось, исправить не удалось, и несколько лет назад она была признана неработающей. Для oshcc все исправлено.
Проведенные нами многочисленные и подробные серии экспериментов позволяют рекомендовать в качестве основного варианта следующую конфигурацию инструментальной среды:
- трансляторы Intel,
- - коммуникационные библиотеки – MPI (реализация OpenMPI для своего варианта Infiniband), и oshcc/osh++ (реализация openSHMEM с открытым исходным кодом),
- - средства программирования ускорителей CUDA – стандартный CUDA Toolkit от nVidia.
Конкретные соображения о том, почему мы рекомендуем именно такую инструментальную среду, можно прочитать в Приложении 1.
Для того, чтобы настроиться на конкретный вариант инструментальной среды, необходимо скопировать в свою домашнюю директорию специально разработанные для этого варианта стандартные конфигурационные файлы .bash_profile и .bashrc. Впоследствии их, конечно, можно расширять и модифицировать, но с сохранением исходных настроек. Чтобы скопированные в домашнюю директорию конфигурационные файлы начали действовать, необходимо завершить сеанс – «выйти», и «войти» снова.
Стандартные комплекты конфигурационных файлов для разных вариантов инструментальной среды хранятся в директории /common/profile.versions. Каждая поддиректория этой директории соответствует некоторому варианту. Директория для варианта, который мы выбрали, называется openshmemmpi.
В разных вариантах инструментальной среды команды трансляции и запуска прикладных программ могут несколько отличаться.
В выбранном нами варианте команда запуска на счет всегда одна. Это команда mpirun, стандартная для СУППЗ, применяемой на К-100, как и на всех прочих машинах серии МВС. В отличие от более ранних вариантов этой команды, обязательным является ключ «-ppn <число процессов на узле>».
Например, команда:
mpirun -ppn 5 -np 25 myapplicationзапускает приложение myapplication в виде параллельной программы из 25 процессов, располагая их по 5 процессов на одном вычислительном узле. Максимальное задаваемое число процессов на узле равно 12. При подсчете числа свободных процессоров СУППЗ считает, что на каждом вычислительном узле процессоров именно 12.
Правильное задание ключа «-ppn» бывает особенно важно для гибридных приложений. Например, если требуется, как в примере, разбираемом ниже, чтобы у каждого процесса параллельной программы был в распоряжении отдельный ускоритель CUDA, программу придется запускать с ключом «-ppn 3», поскольку в состав каждого вычислительного узла К-100 входит 3 ускорителя.
В отличие от команды запуска на счет, команды трансляции параллельных приложений могут быть разными, в зависимости от используемых в приложении языков и/или коммуникационных библиотек.
Если в программе используется shmem (не важно, только shmem или в смеси с MPI), то транслировать программу следует командой oshcc (для программ на С), или же командой osh++ (для программ на С++). Использование shmem в программах на Фортране планируется, но пока не реализовано.
Если в программе используется только MPI (без shmem), то команды трансляции – стандартные для свободно распространяемых версий MPI: mpicc, mpiCC, mpif77, mpif90 и т. п.
Наконец, если программа гибридная, то ее следует поделить (разместив в разных файлах с исходным текстом) на 2 части: коммуникационную и ускорительную.
Файлы коммуникационной части должны иметь тип «.cpp» и не должны содержать языковых конструкций CUDA.
Файлы ускорительной части должны иметь тип «.cu». В них должны быть сосредоточены все необходимые конструкции CUDA – расширений, спрятанные в функции со стандартным (без CUDA) внешним интерфейсом. Функции ускорительной части должны вызываться по мере надобности из коммуникационной части программы.
Файлы исходных текстов ускорительной части (.cu) транслируются в объектные модули (.o) командой:
nvcc -arch=compute_20 -c <имя файла cuda>.cu
Файлы исходных текстов коммуникационной части транслируются, как указано выше, и компонуются с объектными модулями ускорительной части, например:
nvcc -arch=compute_20 –c mykernel.cu mpiCC –o myapp myapp.cpp mykernel.o
Даже если Вы планируете писать коммуникационную часть программы на обычном базовом С, без «++», следует, тем не менее, файлу с исходным текстом дать расширение «.cpp», и использовать команды трансляции/компоновки mpiCC или shmemCC, то есть формально превратить Вашу программу в программу на С++. Это необходимо для правильной компоновки с объектными модулями, которые делает nvcc.
Шаг 2. Учимся писать на С «настоящие Фортрановские программы».
Как читатель уже, наверное, понял из сказанного выше, мы рекомендуем начинать осваивать разработку гибридных приложений на С, а не на Фортране. Такой выбор, на наш взгляд, имеет целый ряд преимуществ, как общего характера, так и технических, применительно к текущему состоянию программного обеспечения для гибридных машин. Однако, имеет он также и весьма существенный недостаток. В отличие от Фортрана, канонический С (и С++) лишен средств работы с многомерными массивами, размеры которых задаются динамически, то есть не константами, а переменными, значения которых заранее не известны. Во многих программах вычислительного характера эта возможность – главное и чуть ли не единственное, что требуется программисту от языка. Вероятно, именно из-за наличия в языке этой возможности многие программисты вычислительных задач выбирают Фортран, утверждая, что язык С, этой возможности лишенный, «плохо приспособлен для вычислительных задач». Другие программисты вычислительных приложений все же не отказываются от С, но вынужденно пользуются теми или иными языковыми ухищрениями (массивы, представленные как иерархические структуры указателей, специальные типы «вектор» и «матрица», организованные средствами С++, системы макросов для «многомерного» доступа к одномерным массивам и т. п.). Все это представляется не очень правильным по трем следующим причинам:
- Простое должно делаться просто, сложностей нам и так хватит, когда будем изучать использование CUDA,
- Многие из упомянутых выше приемов могут заметно снижать быстродействие универсального процессора, особенно в «тугих» критических циклах работы с такими «искусственно многомерными» массивами,
- «Искусственные» способы организации многомерных массивов зачастую плохо соотносятся с использованием CUDA. Ускоритель не очень приспособлен для обработки сложных ссылочных структур. Данные, копируемые из памяти универсального процессора в память ускорителя, обычно требуется располагать в ней в виде плотных массивов, а необходимость предварительно «добывать» куски таких массивов из экземпляров организованных средствами С++ типов «вектор» или «матрица» отнюдь не упростит нашу и без того довольно сложную задачу освоения гибридного программирования.
Словом, перед тем, как перейти к поэтапному построению примера гибридной программы, нам необходимо выбрать раз и навсегда простой и эффективный способ организации многомерных массивов, и попрактиковаться в его использовании на тривиальных примерах.
Наша задача несколько облегчается тем, что в языке С99, являющимся расширением базового С, массивы с динамически задаваемыми размерами поддерживаются, хотя и с некоторыми ограничениями. Расширенные возможности С99 официально не включены в стандарт С++, но фактически традиционно присутствуют в трансляторах основных производителей, в том числе – в трансляторах Intel и PGI, используемых на К-100. Из минимального подмножества этих расширений и будем создавать необходимые нам «настоящие Фортрановские» многомерные массивы.
Напомним кратко основные сведения о многомерных массивах в С. Многомерный массив в С – это массив массивов на единицу меньшей размерности. Так, двумерный массив, объявленный как
double da[100][30]- это, в действительности, одномерный массив длиной 100, элементами которого являются одномерные массивы вещественных чисел длиной 30. Из этого следует, что хранятся многомерные массивы в памяти «по строкам», то есть последний индекс меняется быстрее всех, а первый – медленнее всех, то есть не как в Фортране, а в точности наоборот. В отличие от Фортрана, нельзя задать нижнюю границу изменения индекса: индексация по каждому измерению всегда начинается с нуля.
Определимся теперь конкретнее с набором возможностей, которых нам не хватает. Нам необходимо уметь:
- Объявить массив, размеры которого по всем, вообще говоря, измерениям задаются не константами, а текущими (на некоторый фиксированный момент) значениями переменных,
- Отвести память под такой массив,
- Передать такой массив в функцию в качестве аргумента, причем значения размеров переданного массива по переменным измерениям передать в качестве аргументов этой же функции, как это обычно делается в программах на Фортране.
Если мы раз и навсегда научимся делать это, причем просто и единообразно, наш выбор языка С для создания гибридных приложений вычислительного характера ни в чем не будет уступать выбору Фортрана.
Предлагаемые возможности проиллюстрируем тривиальным примером (текст – в файле t2.cpp). Тестовая программа создает двумерный массив динамически задаваемого размера, заполняет его некоторым шаблоном, печатает как одномерный, чтобы можно было убедиться, что данные в нем расположены по строкам, а затем передает в функцию, в которой он, в свою очередь, частично печатается. Разберем этот пример строка за строкой:
#include <stdio.h>
#include <stdlib.h>
#include "dimension.h"
int main( void )
{
int n = 4;
{
DIM2( float, d, n );
Заголовочный файл «dimension.h» содержит описания макросов объявления многомерных массивов, из которых в данной программе используется макрос «DIM2» (объяснение см. ниже).
В инициализирующей части программы объявляется и получает значение переменная n, в которой хранится задаваемый динамически размер строки создаваемого двумерного массива d. В нашем тестовом случае переменной n просто присваивается константа, хотя в реальной программе значение n, скорее всего, будет вводиться из файла исходных данных или вычисляться на основе значений, вводимых из файла исходных данных. Как только значение n стало известным, можно объявлять массив d, строки которого имеют размер n.
В каноническом С принято писать в теле функции сначала все объявления переменных, а затем – все выполняемые операторы. В данном случае мы физически не можем так поступить – переменная n может получить требуемое значение только в результате работы выполняемых операторов, да и само объявление массива с динамическим размером становится, в известном смысле, выполняемым – оно «срабатывает» в некоторый момент выполнения программы, как какой-нибудь оператор присваивания. Момент же этот наступает после того, как получила свое значение переменная n. Не вдаваясь в многословные рассуждения о том, какие диалекты С и/или С++ позволяют нарушать канонический стиль записи программы, и как именно, условимся поступать описанным ниже простым и заведомо работоспособным способом.
В том месте программы, где все необходимые динамические размеры уже известны, просто откроем новый блок (в нашем примере – открывающая фигурная скобка после строки «int n = 4;»). В пределах этого блока будем поступать так, как мы привыкли поступать в теле функции – сначала все объявления, потом все выполняемые операторы. В некоторых случаях, для некоторых вариантов языка (С++ или же базовый С), для некоторых трансляторов дополнительный блок можно не открывать, но это – не очень хороший стиль программирования, и мы условимся впредь, для ясности, его избегать.
Руководствуясь правилом «сначала объявления», в самом начале вновь открытого блока объявим указатель на двумерный массив с размером строки, равным текущему значению n. Тип элементов массива – float, имя – d. Именно так читается строка: «DIM2( float, d, n );», завершающая приведенный выше фрагмент тестовой программы. Цифра «2» в названии макроса означает размерность массива. Поскольку, как мы уже отмечали выше, двумерный массив хранится в памяти по строкам, при объявлении указателя на него обязательно надо знать размер строки, но не требуется знать число строк. Впоследствии, при реальном отведении памяти под массив, который будет адресоваться этим указателем, мы решим, сколько строк нам потребуется.
Теперь разберем следующий фрагмент:
float *pd;
int i, j;
extern void tf( void *f, int k, int l );
printf( "%ld %ld\n", sizeof(d), sizeof(*d) );
d = (typeof(d))malloc( 3*sizeof(*d) );
for ( i = 0; i < 3; i++ )
for ( j = 0; j < 4; j++ )
d[i][j] = i*10+j;
pd = (float*)d;
for ( i = 0; i < 12; i++ ) printf( " %f", pd[i] );
printf( "\n" );
tf( d, n, 2 );
}
return 0;
}
Первые две строки этого фрагмента – объявления переменных, локальных в блоке. Их вполне можно было бы объявить и выше по тексту, непосредственно в теле функции main(), расположение этих объявлений именно здесь никакого «сокровенного смысла» не имеет. Далее следует определение прототипа функции, аргументами которой являются указатель на двумерный массив с динамическими размерами, и сами размеры. Тип массива с динамическими размерами нельзя (по правилам С) использовать в прототипе функции, поэтому аргумент – указатель на массив объявлен как «void *». Внутри функции мы снова приведем его к надлежащему виду (см. ниже).
Печать двух целочисленных значений в следующей строке призвана продемонстрировать, что мы не ошиблись в конструировании типа переменной d. Значение sizeof(d) должно быть равно 8, поскольку d – всего лишь указатель. Значение же sizeof(*d) – размер типа, адресуемого указателем d, должно быть равно размеру строки двумерного массива в байтах, то есть массива элементов типа float длиной n, то есть, в данном случае, 16.
В следующей строке (с обращением к функции «malloc()») под двумерный массив отводится память, в количестве 3 строк по n вещественных чисел каждая. Указатель d устанавливается на начало этой области памяти. С этого момента мы располагаем полностью готовым к использованию двумерным массивом d, размером 3 на 4. Обе величины (число строк, равное трем, и число столбцов, равное четырем) мы задавали динамически. Приведение типа (typeof(d)) перед обращением к функции malloc() обязательно – ведь мы, рассчитывая в перспективе включить наш код в гибридное приложение, собираемся транслировать его транслятором С++, который в этом месте очень строг.
Остаток текста приведенного выше фрагмента программы тривиален. Массив d заполняется простым регулярным шаблоном. Затем значение d присваивается указателю на одномерный вещественный массив, который печатается. По порядку расположения чисел из шаблона заполнения двумерного массива должно быть видно, что в памяти он, действительно, расположен строками длиной 4, и строк таких 3. Далее массив передается в качестве параметра в функцию, два других параметра которой – размер строки массива и число строк. Поскольку массив хранится в памяти по строкам, число строк можно передать меньшее, чем фактически есть в массиве, что и делается (передается значение 2, а не 3).
Теперь рассмотрим завершающий фрагмент программы – текст функции, получающей массив с динамическими размерами и сами размеры в качестве аргументов.
void tf( void *vf, int k , int l )
{
DIM2( float, f, k ) = (typeof(f)) vf;
int i, j;
printf( "k=%d\n", k );
for ( i = 0; i < l; i++ )
{
for ( j = 0; j < k; j++ )
{
printf( " %f", f[i][j] );
}
printf( "\n" );
}
}
Тип аргумента, в действительности являющегося массивом с динамическими размерами, указан и в заголовке функции как «void*», поскольку, как мы уже говорили выше, С не позволяет использовать типы, определенные лишь внутри некоторого блока, в описаниях и объявлениях функций. Однако, ничто не мешает нам уже опробованным в основной программе способом снова породить, внутри тела функции, нужный нам «внутриблочный» тип, и привести к этому типу значение, переданное в качестве аргумента. Именно это делается в первой после открывающей тело функции фигурной скобке строке. Далее переданный в функцию массив печатается как двумерный, чтобы продемонстрировать, что передан он правильно.
Мы намеренно потратили так много места на очень подробное, возможно, даже излишне подробное, разъяснение работы с многомерными массивами «в Фортрановском стиле». Призываем читателя не пожалеть, если потребуется, даже одного – двух часов на то, чтобы ДЕЙСТВИТЕЛЬНО ПОНЯТЬ все изложенное в этом разделе. Уверенное владение многомерными массивами в описанном стиле гигантски упростит нам как разъяснение примера модельного гибридного приложения, так и последующую работу над реальными программами.
В заключение отметим, что, по аналогии с макросом «DIM2()» для объявления указателей на строки двумерных массивов, в заголовочном файле «dimension.h» определен макрос «DIM3()» для объявления указателей на слои трехмерных массивов. Тестовая программа, демонстрирующая использование этого макроса, очень похожая на только что разобранный пример, находится в файле «t3.cpp». По аналогии с DIM2() и DIM3(), предусмотрены также макросы для массивов больших размерностей, вплоть до DIM7() включительно.
Шаг 3. Выбираем численный метод для реализации в качестве модельного гибридного приложения.
Наша конечная цель состоит в том, чтобы разработать модельную программу, использующую все основные виды параллелизма, предусмотренного в К-100. Программа эта должна быть параллельной в традиционном «кластерном» смысле, то есть состоять из ветвей – процессов, взаимодействующих посредством обращения к MPI и/или shmem. Однако, этого мало. Примерно 80% вычислительной мощности К-100 сосредоточено в его ускорителях CUDA. Разработав «кластерную» параллельную программу, мы затем модифицируем код отдельного процесса таким образом, чтобы он использовал не только центральный процессор (процессоры) вычислительного узла, но и ускоритель.
Выполняя всю эту непростую цепочку постепенного усложнения и совершенствования программы, мы должны постоянно иметь перед глазами образец, «спецификацию» того численного метода, к гибридно – параллельной реализации которого мы стремимся.
Выберем в качестве модельного приложения решение двумерной задачи Дирихле для уравнения теплопроводности на прямоугольной равномерной индексной сетке итерационным методом Якоби, без контроля невязки, то есть с заданием фиксированного числа итераций. Ниже приводится текст программы, реализующей этот расчет в последовательном режиме на универсальном процессоре.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dimension.h"
// Change here the number of steps, the cell geometry, etc
#define NITER 5000
#define STEPITER 1000
#define delx 0.5
#define dely 0.25
// end change here.
/***/
int main( int argc, char **argv )
{
int i, j, n, mx1, my1, mx, my;
FILE *fp;
double rdx2, rdy2, beta;
/***/
if ( argc != 3 )
{
fprintf( stderr, "Usage: %s <nrows> <ncolumns>\n", argv[0] );
return (-1);
}
mx = (int)atol( argv[1] );
my = (int)atol( argv[2] );
mx += 2;
my += 2;
/* Here we know the array sizes, so make the arrays themselves: */
{
DIM2( double, f, my );
DIM2( double, newf, my );
DIM2( double, r, my );
typeof( f ) pf[2];
int curf;
/***/
f = (typeof(f))malloc( 2*mx*sizeof(*f) );
r = (typeof(r))malloc( mx*sizeof(*r) );
if ( (!f) || (!r) )
{
fprintf( stderr, "Cannot allocate, exiting\n" );
return( -1 );
}
curf = 0;
pf[0] = f;
pf[1] = f + mx;
newf = pf[1];
rdx2 = 1./delx/delx;
rdy2 = 1./dely/dely;
beta = 1.0/(2.0*(rdx2+rdy2));
printf( "Solving task on %d by %d grid\n", mx-2, my-2 );
fflush( stdout );
for ( i = 0; i < mx; i++ )
{
for ( j = 0; j < my; j++ )
{
if ( (i==0) || (j==0) || (i==(mx-1)) || (j==(my-1)) )
newf[i][j] = f[i][j] = 1.0;
else
newf[i][j] = f[i][j] = 0.0;
r[i][j] = 0.0;
}
}
mx1 = mx - 1;
my1 = my - 1;
/* Iteration loop: */
for ( n = 0; n < NITER; n++ )
{
if ( !(n%STEPITER) ) printf( "Iteration %d\n", n );
/* Step of calculation starts here: */
f = pf[curf];
newf = pf[1-curf];
for ( i = 1; i < mx1; i++ )
{
for ( j = 1; j < my1; j++ )
{
newf[i][j] = ((f[i-1][j]+f[i+1][j])*rdx2
+(f[i][j-1]+f[i][j+1])*rdy2-r[i][j])*beta;
}
}
/* swap the halves: */
curf = 1 - curf;
}
fp = fopen("foutput.dat", "w" );
for ( i = 1; i < (mx-1); i++ )
{
fwrite( &(newf[i][1]), my-2, sizeof(newf[0][0]), fp );
}
fclose( fp );
}
return 0;
}
Несколько кратких пояснений.
Температура нормализована к диапазону от 0.0 до 1.0. Уравнение решается с правой частью, представленной массивом r, который в данном расчете заполняется нулями, но в пересчете значений температуры, тем не менее, участвует.
Массив значений температуры отводится двойного размера. Имеется вспомогательный массив указателей pf, размером 2, причем pf[0] постоянно указывает на первую половину массива значений температуры, а pf[1] – на вторую. Также имеются два указателя, f и newf, которые указывают на разные половины массива значений температуры, меняясь ролями на каждой итерации (если на текущей итерации f указывает на первую половину массива значений температуры, а newf – на вторую, то на следующей итерации будет наоборот). При этом половина массива значений температуры, на которую на данной итерации указывает f, играет роль «старых» значений, а половина, на которую указывает newf, роль «новых».
В остальном программа вполне тривиальна. Отныне она является для нас руководством к действию. Естественно, мы будем стремиться к параллельной реализации в точности тех вычислений, которые представлены приведенной выше последовательной программой.
Шаг 4. Пишем параллельную программу без использования ускорителей CUDA.
Общая схема параллельной реализации итерационной процедуры Якоби на индексных прямоугольных сетках тривиальна и хорошо известна, остановимся здесь на ней лишь в общих чертах.
Массив значений температуры поделим между процессами параллельной программы на блоки по строкам. Например, если массив содержит 100 строк, не считая граничных условий, а расчет надо проводить на 10 процессорах, то каждый процесс будет рассчитывать новые значения для блока из 10 строк. Чтобы этого добиться, массив старых значений в каждом процессе должен содержать уже не 10, а 12 строк. Две добавленные строки – одна в начале блока, одна – в конце – представляют собой копии краев блоков, расположенных в «соседних» процессах. Эти дополнительные строки, окаймляющие блок, называются теневыми гранями, а сам блок, расположенный в данном процессе, включая теневые грани, будем называть локальной порцией обрабатываемого массива.
Итерация в рамках каждого процесса, таким образом, складывается из двух фаз: вычислительной и коммуникационной. В вычислительной фазе по старым значениям локальной порции формируются новые значения для всех строк локальной порции, кроме теневых граней. В коммуникационной фазе, которая следует за вычислительной, обновленные значения из краев локальных порций пересылаются в теневые грани «соседних» процессов.
Самое трудное в реализации такой схемы – как ни странно, не реализация коммуникаций, а корректное распределение работы между процессами параллельной программы. Число обновляемых строк массива значений температуры, вообще говоря, не делится нацело на число процессоров, на котором запущена задача. В связи с этим велик соблазн «немного подправить» сетку, то есть провести не в точности тот же, а «из физических соображений, почти тот же» расчет, что был бы выполнен, если бы такие исходные данные были заданы эталонной программе для одного процессора. То есть, в конечном итоге, «слегка» изменить суммарное число строк.
Так поступать нельзя ни в коем случае. В реальных задачах, основанных на той же технике параллельной реализации, которую мы здесь для себя наметили, такое «физически несущественное изменение» может привести, например, к тому, что ранее сходившийся метод начнет просто делить на ноль. В любом случае, кошмар при отладке и последующей эксплуатации при таком подходе гарантирован. Как мы уже отмечали выше, мы будем обеспечивать строго корректную параллельную реализацию.
Критерий корректности параллельной реализации в нашем случае очень простой. Если параллельная реализация выполнена правильно, то файлы с результатами расчета для сетки одного и того же размера, с одним и тем же числом итераций, но посчитанные на разном числе процессоров, должны строго, с точностью до бита, совпадать. В противном случае программа, вольно или невольно, скорее всего, «слегка подправляет» размер расчетной сетки.
Кроме корректности вычисления размера локальной порции, введем еще одно, дополнительное требование. Очевидно, при распределении по процессам числа строк сетки, не делящегося нацело на число процессов, последнему по счету процессу достанется не столько строк, сколько всем остальным. Будем требовать, чтобы число строк сетки, обрабатываемых последним процессом, было не больше, а, если надо, меньше, чем у всех остальных. Это важно для того, чтобы единственный перегруженный работой процесс не замедлял весь расчет. Следует иметь в виду, что при делении работы на очень мелкие части это дополнительное требование может оказаться невыполнимым.
С учетом рассказанного выше, приступим к разбору, фрагмент за фрагментом, «кластерного» варианта параллельной программы.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <shmem.h>
#include "shmemtime.h"
#include "dimension.h"
#define FILENAME "a_cpu.dat"
// Change here the number of steps, the cell geometry, etc
#define NITER 5000
#define STEPITER 1000
#define delx 0.5
#define dely 0.25
// end change here.
/***/
int main( int argc, char **argv )
{
int i, j, n, mx1, mx2, my1, my_number, n_of_nodes, totalmx, partmx, leftmx, mx, my;
FILE *fp;
double howlong;
double rdx2, rdy2, beta;
/***/
shmem_init();
my_number = shmem_my_pe();
n_of_nodes = shmem_n_pes();
if ( argc != 3 )
{
if ( !my_number )
{
fprintf( stderr, "Usage: %s \n", argv[0] );
}
return (-1);
}
totalmx = mx = (int)atol( argv[1] );
my = (int)atol( argv[2] );
if ( my < 1 )
{
if ( !my_number )
{
fprintf( stderr, "Number of columns (%d) should be positive\n", my );
}
return (-1);
}
/* Compute the number of rows per node: */
mx = (totalmx+n_of_nodes-1)/n_of_nodes;
/* This is the number of rows for all but the last: */
partmx = mx;
/* This is the number of rows for the last: */
/* It cannot be greater than partmx, but it can be non-positive: */
leftmx = totalmx - partmx*(n_of_nodes-1);
if ( leftmx < 1 )
{
if ( !my_number )
{
fprintf( stderr, "Cannot distribute, too many processors\n" );
}
return (-1);
}
if ( my_number == (n_of_nodes-1) ) mx = leftmx;
/* End rows distribution. */
partmx += 2;
mx += 2;
my += 2;
/* Here we know the array sizes, so make the arrays themselves: */
{
DIM2( double, f, my );
DIM2( double, newf, my );
DIM2( double, r, my );
typeof( f ) pf[2];
int curf;
Этот фрагмент программы является инициализирующим. Он включает в себя вычисление размера локальной порции по упомянутым выше правилам, и объявление указателей на массивы с динамически задаваемыми размерами.
Обращения к функциям shmem_my_pe() и shmem_n_pes() означают выяснение, средствами oshcc, собственного номера процесса и общего числа процессов, соответственно.
Рассмотрим следующий фрагмент:
f = (typeof(f))shmalloc( 2*partmx*sizeof(*f) );
r = (typeof(r))malloc( mx*sizeof(*r) );
if ( (!f) || (!r) )
{
fprintf( stderr, "Cannot allocate, exiting\n" );
exit( -1 );
}
curf = 0;
pf[0] = f;
pf[1] = f + partmx;
newf = pf[1];
rdx2 = 1./delx/delx;
rdy2 = 1./dely/dely;
beta = 1.0/(2.0*(rdx2+rdy2));
if ( !my_number )
{
printf( "Solving heat conduction task on %d by %d grid by %d processors\n",
totalmx, my-2, n_of_nodes);
fflush( stdout );
}
for ( i = 0; i < mx; i++ )
{
for ( j = 0; j < my; j++ )
{
if ( ((i==0)&&(my_number==0)) || (j==0) || ((i==(mx-1))&&(my_number==(n_of_nodes-1))) || (j==(my-1)) )
newf[i][j] = f[i][j] = 1.0;
else
newf[i][j] = f[i][j] = 0.0;
r[i][j] = 0.0;
}
}
mx1 = mx - 1;
my1 = my - 1;
shmem_barrier_all();
/* Iteration loop: */
howlong = shmem_ttt();
for ( n = 0; n < NITER; n++ )
{
if ( !my_number )
{
if ( !(n%STEPITER) )
printf( "Iteration %d\n", n );
}
/* Step of calculation starts here: */
f = pf[curf];
newf = pf[1-curf];
for ( i = 1; i < mx1; i++ )
{
for ( j = 1; j < my1; j++ )
{
newf[i][j] = ((f[i-1][j]+f[i+1][j])*rdx2
+(f[i][j-1]+f[i][j+1])*rdy2-r[i][j])*beta;
}
}
Этот фрагмент включает в себя вычислительную фазу итерационного процесса. Техника ее реализации при помощи двух массивов, которые, меняясь ролями на каждой итерации, попеременно выполняют роль «нового» и «старого», такая же, как в рассмотренной выше эталонной программе для одного процессора.
Поскольку массив значений температуры будет впоследствии участвовать в коммуникациях при помощи shmem, и не является при этом объявленным статически, память для него надо выделять специальной функцией shmalloc() (а не «обычной» malloc()). Обращение к shmalloc() – коллективная операция, требующая указания одинаковой длины массива во всех процессах. Поэтому при вычислении требуемой длины используется partmx, а не mx, то есть в последнем процессе длина берется, вообще говоря, с запасом. Функция shmem_barrier_all() выполняет барьерную синхронизацию всех процессов. В остальном приведенный выше фрагмент программы в комментариях не нуждается.
Перейдем к завершающему фрагменту.
/* Do all the transfers: */
shmem_barrier_all();
if ( my_number > 0 ) shmem_double_put(
&(newf[partmx-1][1]),
&(newf[1][1]),
my-2, my_number-1 );
if ( my_number < (n_of_nodes-1) ) shmem_double_put(
&(newf[0][1]),
&(newf[mx-2][1]),
my-2, my_number+1 );
shmem_barrier_all();
/* swap the halves: */
curf = 1 - curf;
}
if ( !my_number )
{
printf( "Elapsed time: %f sec\n", (float)(shmem_ttt()-howlong) );
fp = fopen( FILENAME, "w" );
fclose( fp );
}
for ( j = 0; j < n_of_nodes; j++ )
{
shmem_barrier_all();
if ( j == my_number )
{
fp = fopen( FILENAME, "a" );
for ( i = 1; i < (mx-1); i++ )
{
fwrite( &(newf[i][1]), my-2, sizeof(newf[0][0]), fp );
}
fclose( fp );
}
}
}
shmem_finalize();
return 0;
}
Фрагмент начинается с коммуникационной фазы, которая состоит из двух обращений к shmem_double_put(), до и после которых выполняется барьерная синхронизация.
Обращение к функции shmem_double_put( target, source, leng, np ) выполняет запись локального («своего») массива source типа double, длиной leng, в массив target, того же типа и той же длины, расположенный в процессе np.
В данном случае каждый процесс переписывает, обращаясь к этой функции, обновленные на вычислительной фазе края локальной порции в «соседские» теневые грани. Для синхронизации коммуникационной и вычислительной фаз используются обращения к функциям shmem_barrier_all(). Оставшаяся часть программы реализует последовательную, процесс за процессом, запись локальных порций в единый файл с результатами расчета, также синхронизированную барьерами.
Полученная программа целиком размещается в одном файле с исходным текстом, должна транслироваться командой oshcc, а запускаться на счет – командой mpirun. Никаких специальных ограничений на значение, задаваемое при этом с ключом «-ppn», не накладывается.
Шаг 5. Выделяем в программе вычислительное ядро.
Чтобы использовать ускорители CUDA, в которых, как мы знаем, сосредоточена львиная доля быстродействия вычислительного узла К-100, нам необходимо, для начала, определиться с тем, какую именно часть имеющегося кода мы хотели бы «переселить» на ускоритель. Требуется выделить в программе все вычислительно – емкие, простые, не содержащие в себе коммуникаций между вычислительными узлами фрагменты, а затем оформить каждый из них в виде отдельно транслируемой функции с простым внешним интерфейсом. Именно такие функции, которые принято называть вычислительными ядрами, мы и будем переписывать для реализации на ускорителе.
В нашем приложении естественным образом выделяется единственное и вполне очевидное вычислительное ядро – код, реализующий собственно пересчет значений сетки по разностному шаблону:
for (i = 1; i < mx1; i++)
{
for ( j=1; j<my1; j++ )
{
newf[i][j]=((f[i-1][j]+f[i+1][j])*rdx2
+(f[i][j-1]+f[i][j+1])*rdy2-r[i][j])*beta;
}
}
Этот код имеет смысл выделить в отдельную функцию, а в текст программы вставить обращение к ней. Модифицированная таким образом программа находится в файлах progrev_shmem.cpp (основная часть) и itstep.cpp (вычислительное ядро, реализованное пока что на универсальном процессоре).
Шаг 6. Реализуем коммуникации между процессором и ускорителем в пределах вычислительного узла.
На первый взгляд может показаться, что для превращения нашего приложения в гибридное нам осталось только переписать текст вычислительного ядра для выполнения на ускорителе. К сожалению, это неверно. Необходимо решить еще несколько вопросов разной важности и сложности.
В составе вычислительного узла – 3 ускорителя. Означает ли это, что в текст нашей программы требуется явно вводить разделение работы, выполняемой вычислительным ядром, на 3 части?
Нет, не означает. Для решения этой проблемы есть более простой способ.
Потребуем (и будем за этим тщательно следить), чтобы наше будущее гибридное приложение всегда запускалось на счет только с ключом «-ppn 3». Таким образом, на каждом из вычислительных узлов будет запущено по 3 процесса с номерами, идущими подряд. Настроив процесс на работу с ускорителем, номер которого (в пределах узла) равен остатку от деления номера процесса на 3, мы решим проблему, снабдив каждый из ускорителей своим, отдельным управляющим процессом. Перегрузки процессорного ресурса при этом не возникнет, поскольку процессорных ядер на вычислительном узле К-100 гораздо больше, чем 3. Отметим, что настройку процесса на конкретный номер ускорителя в пределах вычислительного узла придется предусмотреть в программе явно – ведь ни MPI, ни shmem, ни СУППЗ «не знают» о взаимоотношениях между процессами, процессорными ядрами и ускорителями в пределах вычислительного узла.
Реализация вычислительного ядра на единственном для каждого процесса ускорителе, действительно, свелась бы к переписыванию текста функции itstep(), если бы ускоритель мог непосредственно оперировать данными в памяти универсального процессора. К сожалению, он этого не умеет. У каждого ускорителя имеется своя память, доступная для копирования данных в/из нее по инициативе универсального процессора. Поэтому прежде, чем мы приступим к переписыванию собственно текста вычислительного ядра для работы на ускорителе, мы можем (и должны) выписать в основной, «коммуникационной» части программы, помимо обменов данными между вычислительными узлами, еще и обмены данными между процессором и ускорителем в пределах вычислительного узла.
Логика взаимодействия программы, выполняющейся на процессоре, с программой, выполняющейся на ускорителе, в общих чертах выглядит так.
Программа на процессоре запрашивает у ускорителя области памяти, необходимые для размещения в ускорителе тех массивов, которые будет обрабатывать программа на ускорителе. Это делается обращением к функции cudaMalloc(). Функция эта делает в отношении памяти ускорителя то же, что функция malloc() делает в отношении памяти процессора, то есть находит (в памяти ускорителя, а не процессора!) свободный кусок требуемого размера и сообщает вызывающей программе адрес его начала. Адрес этот запоминается в переменной типа void*, но с ним нельзя обращаться, как со значением переменной того же типа, полученным из обращения к «обычной» malloc(), поскольку это значение есть адрес в памяти ускорителя, а не процессора. Единственное, что мы можем сделать с этим адресом – это запомнить его, с тем, чтобы когда-нибудь потом передать в качестве аргумента другим функциям работы с ускорителем (см. ниже).
Получив адрес областей памяти, зарезервированных в ускорителе, программа на процессоре может копировать данные как из своих массивов в эти области, так и в противоположном направлении. Это делается обращением к функции cudaMemcpy(). В качестве аргументов в эту функцию передаются, в частности, те самые «не настоящие» значения типа void*, которые мы получили из обращения к функции cudaMalloc(), или некоторые смещения относительно этих значений. Как при резервировании памяти, так и при копировании массивов следует учитывать, что размеры (sizeof) базовых типов вещественных чисел (float, double) в процессоре и в ускорителе совпадают, что упрощает вычисление необходимых смещений при выборочном копировании частей массивов.
Организовав необходимые массивы данных в памяти ускорителя, и зная их адреса (те самые «не настоящие» void*), программа на процессоре может теперь вызвать некоторую функцию на ускорителе, передав ей адреса массивов в памяти ускорителя в качестве аргументов. Весьма вероятно, что среди этих массивов будут выходные, то есть такие, значения в которых появятся или изменятся в результате работы этой самой функции на ускорителе. В этом случае программа на процессоре, очевидно, захочет, после вызова функции, скопировать (при помощи уже известной нам cudaMemcpy()) данные, наработанные в памяти ускорителя, в какие-то свои, находящиеся на процессоре, массивы.
Описанная здесь дисциплина может, на первый взгляд, показаться неоправданно усложненной. В самом деле, почему бы не копировать автоматически в ускоритель все входные массивы перед вызовом «ускорительной» функции, а выходные, соответственно, в обратном направлении после вызова? Так, конечно, было бы гораздо проще для программиста, но совершенно недопустимо с точки зрения эффективности. Копирование данных между памятями процессора и ускорителя – важнейшее «узкое место» в гибридных программах, его следует довольно тщательно оптимизировать. Посмотрим с этой точки зрения на нашу модельную программу.
Для того, чтобы могло работать выделенное нами вычислительное ядро, необходимо поддерживать в памяти ускорителя копии размещенных в процессоре массивов значений температуры, а также массива r. Вызывать выполняющуюся на ускорителе функцию, реализующую это вычислительное ядро, мы будем на каждой итерации. При этом массив r достаточно скопировать в память ускорителя всего один раз, поскольку он вообще не меняется в процессе расчета.
С массивами значений температуры ситуация несколько сложнее. Их надо скопировать в память ускорителя перед началом итерационного процесса, а также скопировать обратно после его завершения. После каждой итерации нам необходимо скопировать из памяти ускорителя в память процессора обновленные программой на ускорителе края локальной порции, переслать их в «соседние» процессы, а затем скопировать, уже из памяти процессора в память ускорителя, присланные «соседними» процессами обновленные значения теневых граней. Таким образом, на каждой итерации копируются между памятями процессора и ускорителя только края локальных порций, а не локальные порции целиком. Изложенная выше дисциплина взаимодействия между программой на процессоре и программой на ускорителе идеально подходит для организации именно такого сценария передачи данных.
Теперь все рассмотренное выше необходимо реализовать. Как легко видеть, необходимо внести изменения в основную, коммуникационную часть программы (файл progrev_shmem.cpp). Чтобы не портить текст «кластерного» варианта программы, который является вполне законченным изделием, прибегнем к технике условной трансляции. Потребуем, чтобы при компиляции файла с текстом коммуникационной части приложения в гибридном варианте был определен макрос FORFERMI, а те фрагменты кода, которые придется добавить к тексту основной части программы, будем брать в скобки вида:
#ifdef FORFERMI cudaMalloc(......) ...... #endif
Таким образом, мы сохраним единый текст основной части программы, компилируемый по-разному для «кластерного» и гибридного вариантов.
Модифицированный указанным способом текст основной части программы находится в файле progrev_shmem.cpp. Текст варианта вычислительного ядра без использования ускорителя CUDA по-прежнему находится в файле itstep.cpp, а скрипт для сборки из этих исходных текстов «кластерного» варианта приложения называется compile_for_cpu. Текст варианта вычислительного ядра для ускорителя CUDA, который мы пока не написали, но собираемся, будет находиться в файле itstep.cu, а скрипт для сборки гибридного варианта приложения называется compile_for_fermi.
После всего сказанного в настоящем разделе текст файла progrev_shmem_fermi.cpp, видимо, в дополнительных комментариях не нуждается. Отметим лишь, что использование в коммуникационной части программы функций cudaSetDevice(), cudaMalloc() и cudaMemcpy() никак не противоречит тому, что было сказано выше, в разделе, посвященном командам компиляции и запуска, о разделении программы на коммуникационную и ускорительную части. Коммуникационная часть программы, действительно, не должна содержать языковых конструкций CUDA, но обращения к упомянутым функциям таковыми и не являются, и никаких изменений в процедуре трансляции не требуют. Специальные языковые конструкции CUDA встретятся нам в тексте вычислительного ядра из файла itstep.cu, к написанию которого мы перейдем теперь уже совсем скоро.
Шаг 7. Переписываем текст вычислительного ядра для процессора с использованием алгоритма, эффективно реализуемого на ускорителе.
Вся подготовительная, организационная работа по проектированию и реализации гибридного приложения, наконец, выполнена. Нельзя не порадоваться тому, что мы написали довольно много строк не вполне тривиального кода, руководствуясь обще-программистскими соображениями, не привлекая в явном виде понятий нетрадиционных вычислительных архитектур, то есть, в нашем случае, практически не читая программистских руководств nVidia. Этот факт как нельзя лучше иллюстрирует главный принцип разработки приложений для таких машин, как К-100: «Гибридной машине – гибридное программирование». Пусть мы – разработчики новых суперкомпьютеров – не можем (по вполне объективным, кстати, причинам) сделать вновь создаваемые машины совсем привычными и традиционными. Тогда следует, по крайней мере, позаботиться о том, чтобы так раздражающая пользователей «нетрадиционность» была как можно меньше по объему, и занимала строго определенное место в некоторой структуре, построенной совершенно традиционными, привычными средствами. Вопросы правильной организации приложений не должны сплетаться в тугой клубок с проблемами изучения тонкостей шейдерной архитектуры.
По крайней мере, для нашего «игрушечного» приложения нам пока удавалось следовать этому подходу. Продолжим в том же духе. Хотя технически мы, казалось бы, вполне готовы просто взять и переписать, наконец, текст функции itstep() с использованием технологии CUDA, наберемся терпения и сделаем еще один промежуточный шаг. Попытаемся понять из максимально общих соображений, в терминах традиционного языка С, как именно, за счет чего можно было бы ускорить наше вычислительное ядро при его переносе на ускоритель. Дело в том, что ускоритель – это вовсе не процессор, такой же, как Pentium или Opteron, только более быстрый, а довольно своеобразная параллельная вычислительная система. Надеяться, что формально перенесенный на нее без изменений текст для универсального процессора ускорится сам по себе, было бы так же наивно, как ожидать ускорения последовательной программы для настольного компьютера от механического переноса ее на процессор системы Cray Jaguar.
В основе быстродействия ускорителя CUDA лежит векторный параллелизм. Ускорение достигается за счет одновременной работы большого числа арифметических устройств. Но многие арифметические устройства в состоянии работать одновременно только при наличии памяти, способной их одновременно «прокормить» исходными данными. Та память ускорителя, обмен данными с которой обсуждался в предыдущем разделе, таковой не является. При «лобовом» использовании (внутри программы для ускорителя) данных, расположенных в этой памяти, многочисленные арифметические устройства ускорителя будут не столько работать одновременно, сколько простаивать на конфликтующих друг с другом обращениях к памяти. Значит, мы можем рассчитывать на существенное ускорение нашего вычислительного ядра только в случае, если максимально возможный объем промежуточных данных в процессе вычислений будет размещен в специальной, сверхбыстрой памяти небольшого размера, которая в составе ускорителя CUDA имеется, наряду с «обычной» памятью. В отличие от универсальных процессоров, где роль сверхбыстрой памяти, до некоторой степени, играет кэш, работающий автоматически, в ускорителе CUDA программисту приходится явно указывать в программе, какие именно данные следует размещать в сверхбыстрой памяти. Но в нашем варианте функции itstep() из файла itstep.cpp такие «сверхбыстрые» данные малого объема просто не выделены, при переносе этого текста на ускоритель в сверхбыструю память класть просто нечего.
Для начала, перепишем текст функции itstep()для процессора с явным выделением массивов – кандидатов на размещение в сверхбыстрой памяти. Критерий такого выделения весьма прост. Эти массивы должны быть много меньше по объему, чем все обрабатываемые функцией данные, а попавшие в них значения на протяжении одного вызова функции должны максимально возможное число раз использоваться повторно, без промежуточных записей в «обычную» память.
Посмотрим еще раз внимательно на критический цикл нашего вычислительного ядра:
for( i = 1; i < mx1; i++ )
{
for ( j = 1; j < my1; j++ )
{
newf[i][j]=((f[i-1][j]+f[i+1][j])*rdx2
+(f[i][j-1]+f[i][j+1])*rdy2-r[i][j])*beta;
}
}
В каждом витке внешнего (по строкам) цикла мы имеем дело всего с тремя строками сетки: «верхней» (i-1-й), «средней» (i-й) и «нижней» (i+1-й). Использованные здесь названия соответствуют тому, как двумерные массивы обычно рисуют на картинках. На следующем витке нынешняя «средняя» строка станет «верхней», а «нижняя» – «средней». Таким образом, каждая из строк сетки последовательно играет все три роли – бывает сначала «нижней», потом – «средней» и, наконец, «верхней». Если мы явно предусмотрим в программе использование именно таких трех массивов, то при переносе программы на ускоритель их как раз и можно будет разместить в сверхбыстрой памяти (ведь три строки – это гораздо меньше по объему, чем все поле величин). Коэффициент повторного использования значений двумерного массива при этом будет равен 4 (элементы «средней» строки используются дважды). Модифицированный таким образом текст функции itstep() выглядит так:
#include "dimension.h"
void itstep( int mx, int my, void *pf, void *pnewf, void *pr,
double rdx2, double rdy2, double beta )
{
DIM2( double, f, my ) = (typeof(f))pf;
DIM2( double, newf, my ) = (typeof(newf))pnewf;
DIM2( double, r, my ) = (typeof(r))pr;
double upper[my];
double middle[my];
double lower[my];
int i, j, mx1, my1;
/***/
mx1 = mx - 1;
for ( j = 0; j < my; j++ )
{
upper[j] = f[0][j];
middle[j] = f[1][j];
}
for ( i = 2; i < mx; i++ )
{
for ( j = 0; j < my; j++ ) lower[j] = p[i][j];
for ( j = 1; j < (mx1); j++ )
{
newf[i-1][j] = beta*(
(upper[j]+lower[j])*rdx2
+ (middle[j+1]+middle[j-1])*rdy2
- r[i-1][j]);
}
for ( j = 0; j < my; j++ )
{
upper[j] = middle[j];
middle[j] = lower[j];
}
}
}
Этот вариант текста функции itstep() для универсального процессора интересен даже сам по себе. Например, он делает ненужным использование второго экземпляра сетки (newf). Имея в виду именно этот вариант, можно было бы обойтись единственным массивом f, а итерационная процедура Якоби была бы реализована при этом вполне корректно. Мы, однако, воздержимся от упрощения задним числом коммуникационной части нашей программы – «старый» и «новый» экземпляры сетки нам снова понадобятся впоследствии.
Впрочем, приведенный только что вариант текста вычислительного ядра нужен нам вовсе не сам по себе, а, прежде всего, как образец для подражания при написании варианта itstep() с использованием технологии CUDA. К чему мы (наконец-то!) и приступаем.
Шаг 8. Проектируем вычислительное ядро для мультипроцессора ускорителя CUDA.
Начнем с некоторых важных сведений о CUDA. Мы уже упоминали о том, что ускоритель CUDA – векторный процессор. Это верно, но не вполне точно. Правильнее было бы назвать его векторно – мультитредовым. Мультипроцессор в составе ускорителя выполняет в псевдо-параллельном режиме на сравнительно небольшом (всего 8 или 16) числе одновременно работающих арифметических устройств большое (тысячи) количество одинаковых процессов (тредов), каждый из которых обрабатывает свою маленькую порцию данных. Например, в случае нашего вычислительного ядра в программе для треда, скорее всего, не было бы внутренних (по j) циклов прохода по строке двумерного массива. Вместо этого программа обрабатывала бы всего один – j-й – столбец этого массива, но тредов, выполняющих эту программу, было бы запущено столько, сколько элементов в строке массива (my). Здесь просматривается прямая аналогия с традиционной, «кластерной» параллельной программой, использующей MPI, shmem или подобные им коммуникационные библиотеки. В самом деле, в обоих случаях:
- пишется одна программа,
- ее запускают в некотором числе экземпляров,
- каждый экземпляр запущенной программы, узнав свой номер и количество себе подобных, настраивается на предназначенную ему часть работы.
Одно из несущественных различий – в том, что аналог команды «mpirun», запускающий блок тредов в CUDA, не набирается на клавиатуре, а инициируется программно, оператором вызова «ускорительной» функции из программы на универсальном процессоре. В полной аналогии с командой «mpirun», в операторе вызова такой функции обязательно указывается, в каком количестве тредов эту функцию запустить.
Для чего надо запускать не 8 или 16 (по количеству реально имеющихся в мультипроцессоре арифметических устройств), а сотни или тысячи тредов, работающих псевдо-параллельно? Неужели только для упрощения программы за счет ликвидации в ней внутренних циклов? Конечно, нет. Мультитредовость выполняет совершенно конкретную функцию, очень важную для быстродействия, а не только для упрощения текста программы. Как бы мы ни старались использовать в программе треда сверхбыструю память, совсем без обращений к «обычной», медленной памяти нам не обойтись. Каждое такое обращение очень дорого, равно времени выполнения нескольких сотен арифметических операций с вещественными числами. В состав мультипроцессора CUDA входит специальное устройство – мультитредовый планировщик, который позволяет треду, выполнившему такое обращение, «заснуть» и уступить арифметическое устройство другому треду, обратившемуся к медленной памяти ранее, и как раз готовому «проснуться». В результате арифметические устройства мультипроцессора гораздо меньше простаивают по причине задержек на обращении к памяти (хотя проблема суммарной производительности памяти, вызванная многочисленностью арифметических устройств, остается все равно). Чтобы эффект мультитредовой нейтрализации задержек при обращении к памяти проявился в полной мере, тредов должно быть хотя бы раза в 2 больше, чем величина задержки в тактах (а она, напомним, имеет порядок тысячи). Так что исключение внутренних циклов путем «поручения» каждому треду своего, отдельного столбца сетки – не просто программистский изыск, а прием повышения производительности.
Приведем фрагмент текста вычислительного ядра, написанного с учетом сказанного выше. Будем иметь в виду, что представленный ниже текст – это текст программы для одного треда, обрабатывающего один столбец сетки. Вопрос о том, как правильно вызвать эту функцию из основной, «процессорной» части нашей программы, пока отложим:
#include <cuda.h>
extern __shared__ double pmiddle[];
__global__ void iteration( double *pf, double *pfn, double *pr,
int mx, int my, double rdx2, double rdy2,
double beta )
{
int i, n;
double *middle = pmiddle;
double upper, lower;
/***/
/* pipeline the grid portion through the shared memory "vector registers": */
upper = pf[threadIdx.x];
pf += my;
pfn += my;
middle[threadIdx.x] = pf[threadIdx.x];
pr += my;
for ( i = 2; i < mx; i++ )
{
lower = pf[threadIdx.x+my];
__syncthreads();
if ( (threadIdx.x > 0) && (threadIdx.x < (blockDim.x-1)) )
{
pfn[threadIdx.x] = beta*(
rdx2*(upper+lower)
+ rdy2*(middle[threadIdx.x-1]+middle[threadIdx.x+1])
- pr[threadIdx.x]);
}
pr += my;
pf += my;
pfn += my;
__syncthreads();
upper = middle[threadIdx.x];
middle[threadIdx.x] = lower;
}
}
Спецификатор __global__ при объявлении функции относится к языковым конструкциям CUDA и говорит о том, что данная функция будет выполняться на ускорителе, но вызываться из «процессорной» части программы. О том, откуда взялся массив middle, кто и когда отвел для него память, и какую именно, поговорим немного позже. Обратим внимание на то, что upper и lower, бывшие в прошлой процессорной версии программы массивами, стали теперь одиночными переменными. Это не удивительно – как легко видеть, программа, обрабатывающая единственный, i-й столбец сетки, всегда имеет дело с upper[i] и lower[i], и ни с какими другими элементами этих массивов. Если учесть, что в функциях CUDA одиночные переменные у каждого треда свои, а массивы, если специально не стараться сделать иначе, у всех тредов общие, то становится ясно, почему upper и lower стали одиночными переменными, а middle – нет. Дело в том, что программа обработки i-го столбца записывает только i-й элемент массива middle, но использует (читает) его соседние элементы, те, что записывают «соседние» треды. Чтобы это стало возможным, middle пришлось оставить массивом, общим для всех тредов, а upper и lower вполне логично стали одиночными переменными, которые у каждого треда свои.
threadIdx.x и blockDim.x – встроенные переменные, относятся к языковым конструкциям CUDA, означают в данном случае собственный порядковый номер треда и количество тредов, соответственно. Массивы pf, pfn и pr, общие для всех тредов, переданы как одномерные. Поскольку массивы в С хранятся по строкам, а наша программа обрабатывает столбец, по массивам этим она идет с шагом, равным размеру строки (my). Встроенная функция CUDA syncthreads() выполняет барьерную синхронизацию всех тредов. Больше никаких пояснений по тексту программы, вроде бы, не требуется.
Но где же обещанная сверхбыстрая память? Она уже давно здесь. Дело в том, что переменные и массивы, которыми оперирует наша программа, в действительности, расположены в памяти трех разных видов, два из которых – сверхбыстрые. Остановимся на этом вопросе подробнее.
Из того, что указатели на массивы pf, pfn и pr переданы из «процессорной» части программы как аргументы функции, мы можем заключить, что массивы эти лежат в медленной памяти, поскольку адресов никакой другой памяти внутри ускорителя «процессорная» часть программы «не знает». Однако, все одиночные переменные, к которым относятся, в частности, upper и lower, по умолчанию размещаются в самой быстрой, какая только есть в ускорителе, так называемой регистровой, памяти. На аппаратном уровне регистровая память своя у каждого арифметического устройства. На программном уровне мы видим это через тот факт, что одиночные переменные у каждого треда свои, доступные только этому треду и никакому другому.
Совершенно особый случай – массив middle. Благодаря квалификатору __shared__ (языковая конструкция CUDA), этот массив размещается в так называемой разделяемой памяти. Память эта очень широкополосная, с почти такой же малой задержкой на обращение, как регистровая память, но притом доступная всем тредам, как медленная память. (Поэтому, кстати, на программном уровне массивы __shared__ - общие для всех тредов, подобно массивам, расположенным в медленной памяти). Объем этой памяти совсем невелик: всего лишь первые десятки килобайт, что может быть даже меньше, чем суммарный объем регистровой памяти в мультипроцессоре, и что заведомо гораздо меньше, чем объем медленной памяти (в К-100 – по 1ГБ на ускоритель). Разделяемая память – очень ценный и дефицитный ресурс. К счастью, в нашем модельном приложении каждому треду потребовалось разместить в ней всего по одному значению типа double.
Важное отличие как регистровой, так и разделяемой памяти, от памяти медленной, в том, что ни регистровая, ни разделяемая память не сохраняют своих значений между вызовами «ускорительной» функции. В медленную память «процессорная» часть программы может положить какие-то данные при помощи cudaMemcpy(), и эти данные будут там лежать, пока какая-нибудь «ускорительная» функция, получив при обращении указатель на эти данные в качестве аргумента, их не изменит. В отношении регистровой и разделяемой памяти это неверно. Эти виды памяти не доступны напрямую из «процессорной» части программы. Расположенные в них одиночные переменные и __shared__ - массивы при каждом новом вызове «ускорительной» функции порождаются заново, и в качестве значений до первого присваивания содержат «мусор».
В нашем модельном приложении мы решили упростить себе жизнь, порождая столько тредов, какова длина строки расчетной сетки. Программа в результате упростилась, из нее исчез цикл обхода сетки по строке. Но как быть в случае реальной программы, если в ней потребуется предусмотреть обработку очень широких сеток? Как мы уже отмечали выше, регистровая и, в особенности, разделяемая память – ресурс дефицитный. К сожалению, между тредами, сколько бы их ни было, он распределяется статически. Что будет, если, например, регистры просто кончатся, если их перестанет хватать для всех одиночных переменных каждого треда? Вполне очевидно, что транслятор просто начнет размещать одиночные переменные в медленной памяти. Для ускорителей модели Fermi, установленных на К-100, это до некоторой, очень небольшой, степени терпимо, поскольку эта модель CUDA снабжена еще и кэш-памятью, подобно универсальному процессору. Однако, по мере роста числа не уместившихся одиночных переменных быстродействие начнет весьма заметно уменьшаться. Вопрос о том, как определить, что регистры кончились, имеет несложное решение, но в настоящем тексте не обсуждается. Проще всего просто посчитать на калькуляторе необходимые объемы ресурсов, предварительно узнав на nvidia.com размерные параметры интересующего Вас устройства (Для К-100 – это ускоритель модели Fermi C2050). В заключение отметим, что в реальных программах необходимость соблюдать правильный баланс между желанием, с одной стороны, увеличить количество тредов (чтобы компенсировать задержку при обращении к медленной памяти), а с другой – уменьшить его (чтобы хватило регистров и разделяемой памяти) – важнейшая проблема, которую приходится решать программисту, и важнейший ключ к значительному повышению быстродействия.
Шаг 9. Реализуем вычислительное ядро для ускорителя.
Казалось бы, уж теперь-то мы сделали практически все. Осталось написать «связку» между основной, «процессорной» частью нашего приложения и фрагментом кода для ускорителя, подробно рассмотренным в предыдущем разделе – и гибридное приложение - наконец – готово. К сожалению, для достижения этой благородной цели нам потребуется сделать еще один не вполне тривиальный шаг (теперь уже – точно всего один).
Разобранный в предыдущем разделе фрагмент кода для ускорителя предназначен для выполнения на мультипроцессоре. Но в составе ускорителя мультипроцессоров много (иногда - десятки)! Если мы ограничимся корректным вызовом из основной части приложения того варианта вычислительного ядра, который был разобран в предыдущем разделе, мы используем всего один мультипроцессор, и многократно потеряем в быстродействии. Чтобы использовать все ресурсы устройства, этот код требуется модифицировать для использования на многих мультипроцессорах. Чтобы сделать это, нам, в свою очередь, придется понять, как именно многие мультипроцессоры объединяются в один ускоритель.
Схема такого объединения на аппаратном уровне, действительно, очень проста. Мультипроцессоры внутри ускорителя имеют общую медленную память, а все остальное у них раздельное.
Треды из разобранного в предыдущем разделе примера могли выполнять быструю барьерную синхронизацию и взаимодействовать путем доступа к общей для всех быстрой разделяемой (__shared__) памяти только потому, что выполнялись они на одном мультипроцессоре. Треды на разных мультипроцессорах взаимодействовать не могут, а синхронизируются только завершением вызова «ускорительной» функции из «процессорной» части программы (если вызов завершился, значит, все без исключения треды отработали, все выполненные ими изменения в медленной памяти реально произошли).
На программном уровне это означает, что при запуске «ускорительной» функции «процессорная» часть программы не просто задает требуемое количество тредов, а формирует из них двухуровневую иерархию. Треды объединяются в блоки, при запуске задается как требуемое количество блоков, так и требуемое количество тредов в блоке. Треды внутри блока взаимодействуют так, как мы подробно разбирали выше, в предыдущем разделе. Треды разных блоков не взаимодействуют вообще. Каждый блок выполняется на одном мультипроцессоре, так что при задании числа блоков, меньшего, чем реально имеется мультипроцессоров в ускорителе, ускоритель будет с гарантией загружен не полностью. Реально блоков может быть (и часто бывает) гораздо больше, чем мультипроцессоров, что гарантирует полную загрузку оборудования.
Отсутствие взаимодействия блоков позволяет ускорителю выполнять их либо на разных мультипроцессорах одновременно, либо на одном мультипроцессоре по очереди, либо как-то еще – как получится. О том, в каком ритме и порядке треды разных блоков будут видеть изменения в общей для всех медленной памяти, вносимые друг другом, ничего определенного сказать нельзя. Могут увидеть еще в процессе выполнения вызова функции, если ускоритель решил выполнить соответствующие блоки на разных мультипроцессорах физически одновременно, а кэши случайно вытолкнулись, а могут и не увидеть, если блоки, например, запланированы для выполнения на одном и том же мультипроцессоре по очереди. Гарантируется только итоговая синхронизация по завершению вызова функции: к моменту следующего вызова этой функции выполненные всеми тредами всех блоков изменения в медленной памяти реально вступят в силу.
Как работать в таких невыносимых условиях? Очевидно, программу следует строить так, чтобы:
- разные треды на протяжении одного вызова «ускорительной» функции писали по разным адресам медленной памяти,
- результаты записи новых значений в медленную память, выполненной при текущем вызове функции, читались (использовались) тредами других блоков не раньше, чем при следующем вызове.
Попробуем поделить работу между блоками тредов с учетом этих ограничений.
Поступим примерно так же, как при разделении работы между процессами в основной, «коммуникационной» части программы (см. выше, в Разделе 4). Поделим сетку на блоки строк примерно одинакового размера, обработку каждого такого блока строк «поручим» отдельному блоку тредов. Специально выписывать пересылку данных из краев блока строк в «соседские» теневые грани нам не придется, поскольку все строки лежат в общедоступной медленной памяти.
Убедимся в том, что взаимодействие тредов из разных блоков корректно, поскольку вписывается в сформулированные выше ограничения.
Корректность обеспечивает наша схема работы с двумя версиями сетки, «старой» и «новой», которые при каждой итерации меняются местами.
В самом деле, на каждой итерации треды каждого блока обновляют содержимое только своего блока строк (первое условие соблюдено), и используют значения «соседских» краев, записанные на прошлой итерации, при прошлом вызове функции (второе условие соблюдено).
Разработанный нами в позапрошлом разделе вариант вычислительного ядра с использованием буферных массивов upper, middle и lower позволял нам ограничиться единственной версией сетки, но мы предусмотрительно решили не менять ранее написанную программу. Теперь мы видим, что решили правильно: сейчас пришлось бы менять обратно.
Вот теперь вычислительное ядро, действительно, написано. Проанализируем его текст из предварительной версии файла itstep.cu.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cuda.h>
#define GRIDSIZE 32
extern __shared__ double pmiddle[];
__global__ void iteration( double *f, double *fn, double *rhs,
int nrows, int ncols,
double rdx2, double rdy2, double beta2 )
{
int i, n;
double *pf, *pfn, *prhs;
double *middle = pmiddle;
double upper, lower;
/***/
/* calculate this grid start & size: */
n = nrows/gridDim.x;
pf = f + n*ncols*blockIdx.x;
pfn = fn + n*ncols*blockIdx.x;
prhs = rhs + n*ncols*blockIdx.x;
if ( blockIdx.x == (gridDim.x-1) ) n = nrows - n*(gridDim.x-1);
/* make the shadow edges: */
if ( blockIdx.x == 0 )
{
n++;
}
else
{
pf -= ncols;
pfn -= ncols;
prhs -= ncols;
if ( blockIdx.x == (gridDim.x-1) )
n++;
else
n += 2;
}
/* pipeline the grid portion through the shared memory "vector registers": */
upper = pf[threadIdx.x];
pf += ncols;
pfn += ncols;
middle[threadIdx.x] = pf[threadIdx.x];
prhs += ncols;
for ( i = 2; i < n; i++ )
{
lower = pf[threadIdx.x+ncols];
__syncthreads();
if ( (threadIdx.x > 0) && (threadIdx.x < (blockDim.x-1)) )
{
pfn[threadIdx.x] = beta2*(
rdx2*(upper+lower)
+ rdy2*(middle[threadIdx.x-1]+middle[threadIdx.x+1])
- prhs[threadIdx.x] );
}
prhs += ncols;
pf += ncols;
pfn += ncols;
__syncthreads();
upper = middle[threadIdx.x];
middle[threadIdx.x] = lower;
}
}
void itstep( int nrows, int ncols, void *f, void *fn, void *r,
double rdx, double rdy, double beta )
{
int n, nl;
/***/
n = nrows/GRIDSIZE;
nl = nrows - n*(GRIDSIZE-1);
if ( (n < 1) || (nl < 1) )
{
fprintf( stderr,
"Too many grid blocks (%d) for the number of rows per device (%d), exiting\n",
GRIDSIZE, nrows );
exit( -1 );
}
iteration<<<GRIDSIZE,ncols,ncols*sizeof(double)>>>(
(double*)f, (double*)fn, (double*)r, nrows, ncols, rdx, rdy, beta
);
}
Нам удобно начать анализ этого текста с конца, с функции «itstep()», которая как раз и является «связкой» между «процессорной» и «ускорительной» частью кода. Сама эта функция выполняется на универсальном процессоре, но ее текст содержит языковую конструкцию CUDA – вызов «ускорительной» функции «iteration()». Этот вызов – последний оператор в функции «itstep()». От обычного вызова функции он отличается тем, что после имени функции в тройных угловых скобках указаны параметры запуска тредов. В нашем случае параметры эти имеют вид целых чисел. Первое число задает число блоков, второе – число тредов в блоке, третье – размер разделяемой (__shared__) памяти, отводимое в каждом блоке тредов в дополнение к __shared__ - массивам, объявленным со статическими размерами.
Число блоков не должно быть меньше числа мультипроцессоров, и очень желательно из соображений быстродействия, чтобы оно было кратно 32. Число тредов в блоке у нас, как мы уже отмечали, равно строго длине строки обрабатываемой сетки. Разделяемая память нам требуется в количестве по одному элементу типа double на тред. Перед вызовом ускорительной функции проверяется осуществимость процедуры разделения работы между блоками тредов – каждому блоку должна достаться хотя бы одна строка сетки. От требования, чтобы последнему блоку досталось не больше строк, чем прочим, откажемся – здесь оно часто оказывается трудно выполнимым.
Перейдем к анализу текста функции iteration(). Эта функция выполняется на ускорителе и содержит код для одного треда. Объявление
extern __shared__ double pmiddle[];
не содержит указания размера массива. Это говорит о том, что массив pmiddle, расположенный в сверхбыстрой разделяемой памяти, имеет размер, равный (в байтах) величине, указанной в третьем параметре запуска тредов при вызове функции. Допускается только одно такое объявление массива с неизвестным при компиляции размером. Если требуется организовать несколько __shared__ - массивов неизвестного заранее размера, память под них надо распределять вручную, отсчитывая необходимые смещения от единственного указателя extern __shared__. Если в третьем параметре запуска тредов при вызове функции запрошено больше сверхбыстрой разделяемой памяти, чем реально есть в мультипроцессоре, эффект вызова непредсказуем.
gridDim.x и blockIdx.x – встроенные переменные CUDA, означают общее число блоков и собственный порядковый номер блока, соответственно. Все остальное в тексте должно быть ясно из разобранного выше.
Теперь попробуем все же избавиться от ситуации, когда заданный при запуске программы размер строки сетки оказывается больше максимального допустимого числа тредов в блоке или, что даже более вероятно, оказывается просто настолько большим, что исчерпывается регистровая и/или разделяемая память. К счастью, для решения этой проблемы достаточно внести совсем незначительные изменения в текст функции itstep(). Смысл изменений сводится к тому, что задается предельно допустимое число тредов (константа BLOCKSIZE), которое является одновременно предельно допустимым размером расположенного в сверхбыстрой разделяемой памяти массива middle, выраженным в элементах типа double. Если заданный при запуске программы размер строки сетки не превосходит BLOCKSIZE, то работа программы происходит так, как описано выше. В противном случае «ускорительная» функция iteration() вызывается несколько раз, при каждом вызове обрабатывая не более BLOCKSIZE очередных столбцов сетки. Никаких изменений в тексте самой функции iteration() при этом не требуется. Доработанный таким образом текст функции itstep(), взятый из окончательной версии файла itstep.cu, выглядит так:
void
itstep(intn rows, int ncols, void *f, void *fn, void *r,
double rdx, double rdy, double beta )
{
int n, nl;
/***/
n=nrows/GRIDSIZE;
nl=nrows - n*(GRIDSIZE-1);
if((n<1) || (nl<1))
{
fprintf(stderr,"Too many grid blocks (%d) for the number of rows per device (%d),
exiting\n",
GRIDSIZE,nrows );
exit(-1);
}
n=ncols;
while( n > 2 )
{
nl=n;
if( nl > BLOCKSIZE ) nl=BLOCKSIZE;
iteration<<<GRIDSIZE,nl,nl*sizeof(double)>>>(
(double*)f, (double*)fn,(double*)r,
nrows, ncols, rdx, rdy, beta
);
f = ((double*)f)+(nl-2);
fn = ((double*)fn)+(nl-2);
r = ((double*)r)+(nl-2);
n -= (nl-2);
}
}
Гибридное приложение, таким образом, разработано до конца.
Шаг 10. Переводим дух и думаем, что дальше.
Начнем с краткого перечисления принципиальных возможностей общей организации гибридных приложений, которые мы в данном, модельном случае для простоты и ясности опустили.
Наше гибридное приложение фактически не использует быстродействия универсальных процессоров в составе вычислительного узла ГВК. Между тем, оно составляет хотя и меньшую, но все же вполне заметную (порядка четверти, а то и трети) долю в суммарной производительности узла.
Действительно качественное гибридное приложение должно гораздо лучше, чем это сделали мы, использовать универсальные процессоры. Это можно сделать двумя способами. Первый – вероятно, в реальных приложениях наиболее простой и естественный – состоит в том, чтобы некоторые виды не совсем тривиальной в вычислительном отношении работы выполнять на универсальных процессорах одновременно с работой ускорителя. Для этого надо уметь запускать «ускорительную» функцию асинхронно, с немедленным завершением запуска, а потом отдельным оператором CUDA ждать завершения ранее запущенной на ускорителе работы. Современные версии программного обеспечения CUDA такими возможностями обладают.
Второе, что может потребоваться при организации совместной работы ускорителя и процессорных ядер – это умение писать для процессора параллельные программы в модели параллелизма с общей памятью в пределах одного MPI- (или shmem-) процесса. Поясним, как это связано с гибридным характером приложения.
Суммарное число процессорных ядер, доступных приложению, на узле К-100 равно 12. Реализуя «чисто кластерное», без использования ускорителя, модельное приложение, многоядерность узла можно, в первом приближении, не учитывать вообще. Современные версии MPI и shmem достаточно эффективно реализуют коммуникации между процессами, вне зависимости от того, как они сгруппированы по узлам. Управлять же этой группировкой можно при помощи ключа «-ppn», задаваемого в команде «mpirun» при запуске задачи на счет.
Иное дело – приложение гибридное. На примере нашего приложения мы имели случай убедиться, что правильный ключ при запуске таких приложений – «-ppn 3». Если мы будем разрабатывать приложение в расчете на запуск, скажем, с «-ppn 12», то нам придется сильно усложнить программу, учитывая, что из каждых 12 процессов с номерами, идущими подряд, только три являются гибридными, а остальные 8 – «чисто кластерными». Перспектива делать это вручную вряд ли способна кого-либо вдохновить. При использовании же «удобного» значения «-ppn 3» мы теряем быстродействие 8 процессорных ядер, что обидно. Это быстродействие можно было бы использовать, целиком или частично, если бы каждый из MPI- или shmem- процессов мог, выполняя интенсивную вычислительную работу в «процессорной» части программы, разветвиться, например, на 3 или 4 ветви. Самый правильный способ этого добиться – использовать OpenMP. Впрочем, современные компиляторы позволяют сделать это еще проще, указывая с помощью директив #pragma те циклы, которые следует выполнять в параллельном режиме, с использованием дополнительных процессорных ядер, так же или почти так же, как это делается при использовании OpenMP. Конкретный вид таких частично заменяющих OpenMP директив в разных компиляторах разный, но во всех качественных коммерческих компиляторах они, в том или ином виде, присутствуют.
Подлинным «высшим пилотажем», конечно же, является разработка гибридных приложений, в которых универсальные процессоры и ускорители используются для выполнения одной и той же, по характеру, вычислительной работы, деля ее пропорционально фактически достигаемому на данной работе быстродействию. Впрочем, трудно сказать, оправдывает ли дополнительный скоростной выигрыш от такой организации приложений то весьма значительное усложнение разработки, которое в данном случае неизбежно.
Еще один нетронутый пласт возможностей связан с умением некоторых ускорителей использовать в качестве медленной памяти непосредственно специально выделенную область в памяти универсального процессора. Эта возможность могла бы значительно упростить разработку гибридных приложений, одновременно повысив их быстродействие, если бы не недостаточная производительность интерфейса связи процессора с ускорителем, особенно – в случаях, когда ускорителей в одном вычислительном узле несколько. Недостаточная для того, чтобы «прокормить» вычислитель, пропускная способность канала доступа в память – вещь фундаментальная. Отказываться от трех дополнительных, работающих физически независимо, к тому же – довольно широкополосных устройств памяти, видимо, неразумно.
Словом, эту новую возможность при случае, возможно, стоит испытать, но вряд ли она в обозримом будущем заменит тот способ работы с памятью, который был продемонстрирован выше.
Наиболее общий вывод из всего вышеизложенного состоит в том, что самое сложное и трудоемкое в гибридном приложении – это вовсе не программирование вычислительных ядер для вычислителей нетрадиционной архитектуры как таковое, а аккуратное разделение работы, и главное – организация перемещения данных по иерархии сложной, многоуровневой и гетерогенной системе памяти. Эту работу вполне можно было бы выполнить раз и навсегда, реализовав для нескольких основных шаблонов параллелизма библиотеку коллективных операций работы со сложно, гетерогенно распределенными массивами. Что-то вроде функций работы с распределенными векторами и матрицами известной «кластерной» библиотеки PETSc, только, видимо, более простое. Имея такую библиотеку, скрывающую внутри себя все коммуникации, как между процессами, так и между «процессорной» и «ускорительной» частями программы, можно было бы свести работу прикладного программиста, фактически, ко вводу исходных данных, выводу результатов и написанию вычислительных ядер. Понимая всю важность и неотложность проблемы, мы, тем не менее, не прогнозируем появления таких библиотек в самом ближайшем будущем, поскольку не только реализация, но даже и само проектирование такой библиотеки – дело очень сложное и крайне трудоемкое.
Приложение
1.
Обзор базового программного обеспечения разработки гибридных
параллельных приложений, доступного на К-100.
Из систем программирования, то есть наборов компиляторов С и Фортрана, поставляемых комплектно и совместимых между собой, на К-100 представлены 3: GNU, Intel и Portland Group (PGI).
Компиляторы GNU распространяются свободно. Отсутствует, по крайней мере, в гарантированно пригодном и надежном виде, поддержка OpenMP и Фортрана 90. Качество кода хуже, чем у коммерческих компиляторов. Использование для разработки вычислительно-емких программ на К-100 не рекомендуется.
Компиляторы Intel – символически платные (стоят недорого и легко приобретаются). Характеризуются высоким качеством кода, поддержкой OpenMP, наличием Фортрана 95. Имеются хорошо документированные средства (специальные директивы, добавляемые в текст программы) для автоматизированного распараллеливания циклов на несколько процессорных ядер («Упрощенное подобие OpenMP»). С компиляторами поставляется коммерческая реализация MPI. Имеется реализация Co-Array Фортрана, принудительно привязанная именно к этой реализации MPI (реализованная на ее базе, и не способная работать ни на какой другой реализации MPI). Компиляторы Intel совместимы с компиляторами GNU (исходные тексты на С и С++, часть которых транслировалась компилятором GNU, а часть – компилятором Intel, без проблем собираются друг с другом). Имеется глубоко оптимизированная математическая библиотека, с возможностью вызова вычислительно – емких функций из последовательной программы в параллельном, многоядерном варианте.
Средства программирования CUDA представлены, прежде всего, штатным свободно распространяемым ПО фирмы nVidia – CUDA Toolkit. В комплект входят готовые «ускорительные» библиотеки cuBLAS и cuFFT, практическая применимость которых вызывает большие сомнения из общих соображений, а также собственно система программирования CUDA. Это низкоуровневое ПО для довольно трудоемкого ручного кодирования гибридных приложений на языке С. Разбору примера именно его использования посвящена основная часть настоящего документа. Поддержка Фортрана в CUDA Toolkit ограниченная – «процессорную» часть приложения можно написать на Фортране, но «ускорительную» придется писать на С. Транслятор CUDA от nVidia вызывает в качестве транслятора «процессорной» части кода компилятор GNU С. Поэтому CUDA Toolkit от nVidia можно использовать с компиляторами Intel (они совместимы с GNU), но нельзя – с компиляторами PGI.
Компиляторы PGI имеют собственную встроенную поддержку работы с ускорителями, причем двух разных типов.
С одной стороны, реализован CUDA Фортран – полностью Фортрановский вариант CUDA Toolkit, позволяющий делать все то же самое, с точностью до обозначений, что система программирования CUDA от nVidia, но только на Фортране. CUDA C отсутствует.
С другой стороны, как в Фортран, так и в С встроены средства автоматизированного построения гибридных приложений по тексту последовательной программы, к которому добавлены ориентированные на технологию CUDA директивы, призванные помочь компилятору самостоятельно разделить код на «процессорную» и «ускорительную» части.
Эта возможность из общих соображений крайне привлекательна, однако на практике попытки воспользоваться ей из программ, содержащих вызовы MPI или shmem, к успеху пока не привели.
Как было отмечено выше, в первом разделе основной части настоящего документа, в качестве основного варианта инструментальной среды сегодня рекомендуется сочетание следующих компонентов:
- компиляторы Intel,
- коммуникационные библиотеки OpenMPI и Qlogic shmem,
- CUDA Toolkit от nVidia.
Несколько слов о том, почему мы рекомендуем именно это.
Компиляторы Intel обеспечивают заведомо очень хорошее качество кода, совместимы с компиляторами GNU (и, тем самым, с CUDA Toolkit и любым другим «заточенным» под GNU свободно распространяемым ПО). Выбирая их, мы лишаемся средств автоматизированного построения гибридных программ от PGI, которые пока все равно использовать не удается, а также CUDA Фортрана от PGI.
Тем, кто твердо намерен использовать CUDA Фортран, можно рекомендовать другой вариант инструментальной среды, состоящий из тех же коммуникационных библиотек и компиляторов PGI. Такой вариант подготовлен и испытан, устанавливается Системным Администратором по Вашему запросу.
Отказ от компиляторов PGI означает также отказ от использования HPF. Примеров практического использования HPF вообще известно не очень много. Продвижение этой технологии на К-100 к первоочередным задачам группы поддержки не относится. К тому же неизвестно, как он сочетается с CUDA Фортраном.
Коммуникационная библиотека OpenMPI рекомендована фирмой Qlogic для использования на ее коммуникационном оборудовании совместно с библиотекой shmem. Это активно развивающаяся, известная, качественная свободно распространяемая реализация. Включает в себя поддержку прозрачной организации контрольных точек для параллельных программ, что планируется освоить в среднесрочной перспективе.
Отказ от Intel MPI означает отказ от Co-Array Фортрана, принятие же Intel MPI на вооружение означает отказ от shmem. Учитывая, что данная реализация Co-Array Фортрана нами не тестировалась, в отличие реализации Qlogic shmem, а также то, что shmem предоставляет программисту, с точностью до формы записи, те же возможности, что и Co-Array Фортран, и еще много других, но при этом допускает использование из С/С++, а не только из Фортрана, мы выбрали shmem.
Отметим, что попытки заставить работать средства автоматизированного построения гибридных программ в составе трансляторов PGI будут продолжены. Если на этом направлении будет достигнут успех, компиляторы PGI станут основными. На организации межузловых коммуникаций в прикладных программах это не отразится, поскольку реализации коммуникационных библиотек планируется использовать те же самые. Не изменятся, скорее всего, и команды запуска программ на счет.
Приложение
2.
Что читать, чтобы не только успешно разбирать демонстрационные
примеры, но и программировать самостоятельно.
В основной части настоящего документа мы сознательно избегали отсылок к другим источникам информации, вроде «подробнее об этом можно прочитать в …». Задачей этого текста было не повышение образовательного уровня читателя, а скорейшее прохождение полного пути разработки параллельного гибридного приложения. Если этот марафон впечатлил читателя настолько, что возникло желание разрабатывать свои гибридные приложения в предложенном стиле, рекомендуем систематически освоить следующие вопросы:
- Возможности многоядерного распараллеливания кода (и вообще оптимизации) в компиляторах Intel (или в компиляторах PGI, если Вы выбрали их). Возможности поставляемых с выбранными Вами компиляторами оптимизированных математических библиотек.
- Возможности коммуникационной библиотеки shmem. На наш взгляд, это более простая, а во многих случаях – и более эффективная альтернатива MPI. Ее изучение, по крайней мере, в общих чертах, вряд ли займет у Вас больше одного дня.
- Возможности программирования вычислительных ядер с использованием технологии CUDA. Те знания об этой технологии, которые читатель мог почерпнуть из разобранного выше примера, крайне поверхностны и фрагментарны.
В нашем обзоре мы постарались показать две вещи:
- Написанием кода для ускорителя нетрадиционной архитектуры разработка гибридного приложения не столько начинается, сколько заканчивается,
- Само написание этого кода подчиняется некоторой простой, обозримой и довольно легко объяснимой логике.
Несмотря на безусловную, как нам кажется, верность этих двух утверждений, необходимости все же, рано или поздно, писать текст вычислительного ядра никто пока не отменял, а понимание общей логики и системы не заменяет конкретных знаний, согласно этой системе организованных. Поэтому мы призываем тех, кто все же собирается разрабатывать вычислительные ядра для ускорителей CUDA, потратить некоторое количество времени и сил на более серьезное освоение этой технологии, в первую очередь, по материалам nVidia.
В заключение хотелось бы обратить внимание будущих пользователей на менее очевидное, но, видимо, самое важное направление накопления знаний, необходимых для успешного использования гибридных суперкомпьютеров.
Чтобы дело вообще дошло до написания текста вычислительного ядра, его выделение в используемом численном методе должно быть возможно. Мы постарались, на примере разобранного модельного приложения, показать, что такое вычислительное ядро, как его использование вписывается в систему гетерогенных коммуникаций, какова роль локальности обращений к памяти и ее переиспользования в критическом цикле. Все это необходимо учитывать еще на этапе выбора численного метода. Не будучи сами специалистами в этой области, никаких конкретных реализаций будущим пользователям мы здесь дать не можем. Можем только поблагодарить за внимание, и пожелать успеха.
Исходные тексты примеров
compile_for_cpucompile_for_fermi
dimension.h
shmemtime.h
itstep.cpp
itstep.cu
progrev_shmem.cpp
progrev_shmem_fermi.cpp
run_test
run_test_fermi
Все примеры одним файлом