Исследования › Гибридный реконфигурируемый вычислитель › Инженерная методика адаптации приложения к гибридному кластеру с ускорителями на ПЛИС.

Часть 3. Параллельная реализация на процессорах общего назначения. Запись обменов данными для блочного распределения данных с перекрытием.

Первый этап преобразования нашей программы состоит в ее параллельной реализации с использованием MPI. Объяснение того, как при этом делится работа между MPI - процессами, в общих чертах заимствовано из «Десяти простых шагов», из Шага 4.

Массив значений температуры (сетку) поделим между процессами параллельной программы на блоки по строкам. Например, если массив содержит 100 строк, не считая граничных условий, а расчет надо проводить на 10 процессорах, то каждый процесс будет рассчитывать новые значения для блока из 10 строк. Чтобы этого добиться, массив старых значений в каждом процессе должен содержать уже не 10, а 12 строк. Две добавленные строки – одна в начале блока, одна – в конце – представляют собой копии краев блоков, расположенных в «соседних» процессах. Эти дополнительные строки, окаймляющие блок, называются теневыми гранями, а сам блок, расположенный в данном процессе, включая теневые грани, будем называть локальной порцией обрабатываемого массива.

Итерация в рамках каждого процесса, таким образом, складывается из двух фаз: вычислительной и коммуникационной. В вычислительной фазе по старым значениям локальной порции формируются новые значения для всех строк локальной порции, кроме теневых граней. В коммуникационной фазе, которая следует за вычислительной, обновленные значения из краев локальных порций пересылаются в теневые грани «соседних» процессов.

Блочное распределение массивов данных с «нарезкой» и перекрытием блоков по одному измерению — довольно типичный прием параллельной реализации многих численных методов. Простейший вариант таких блочно распределенных массивов реализован в библиотеке VSDA (Very Simple Distribured Arrays, Очень простые распределенные массивы). Смысл этой совсем небольшой библиотеки в том, чтобы, исходя из общей длины распределенного массива и общего числа MPI - процессов, автоматически рассчитать для каждого MPI - процесса размер локальной порции, размер и размещение теневых граней, и, с учетом этого, уметь автоматически выполнять, по запросу вызывающей программы, все необходимые обмены данными. Библиотека VSDA реализована с использованием MPI. Параллельная версия программы, написанная с использованием библиотеки VSDA, приводится ниже. Описание библиотеки VSDA приводится в Приложении 1.

#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "vsda.h"
#include "dimension.h"
// Change here the number of steps, the cell geometry, etc
#define STEPITER 1000
#define delx 0.5
#define dely 0.25
#define NITER 10002
// end change here.
/***/
    extern void itstep( int mx, int my, void *f, void *newf, void *r, 
                                    double rdx2, double rdy2, double beta );
    int main( int argc, char **argv )
{
    int i, j, k, n, mx1, mx2, my1, my_number, n_of_nodes, 
            totalmx, partmx, leftmx, mx, my, niter;
    FILE *fp;
    double howlong;
    double rdx2, rdy2, beta;
    VSDA vsda;
    long *distr_out;
    MPI_Request rq[5];
/***/
    MPI_Init( &argc, &argv );
    MPI_Comm_rank( MPI_COMM_WORLD, &my_number );
    MPI_Comm_size( MPI_COMM_WORLD, &n_of_nodes );
    if ( argc != 3 )
    {
        if ( !my_number )
        {
            fprintf( stderr, "Usage: %s <nrows> <ncolumns>\n", argv[0] );
            fflush( stderr );
        }
        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 );
            fflush( stderr );
        }
        return (-1); 
    }
    distr_out = (long*)malloc( n_of_nodes*sizeof( *distr_out ) );
    if ( !distr_out )
    {
        fprintf( stderr, "No memory\n" );
        fflush( stderr );
        return( -1 );
    }
    totalmx += 2;    
    my += 2;
    niter = NITER;
/* Here we know the array width, so make the arrays themselves: */
    {
        DIM2( double, f, my );
        DIM2( double, newf, my );
        DIM2( double, r, my );
        void *tmpf;
/* Make the distributed array, discover the local portion length: */
        VSDA_Create( MPI_COMM_WORLD, VSDA_BLOCKED, 0, 
                           (long)totalmx, (long)sizeof(*f), (long)1, NULL, distr_out, &vsda );
        mx = distr_out[my_number];
        if ( n_of_nodes > 1 )
        {
            if ( my_number > 0 ) mx ++;
            if ( my_number < (n_of_nodes-1) ) mx ++;
        }
        f = (typeof(f))malloc( mx*sizeof(*f) );
        newf = (typeof(newf))malloc( mx*sizeof(*newf) );
        r = (typeof(r))malloc( mx*sizeof(*r) );
        if ( (!f) || (!newf) || (!r) )
        {
            fprintf( stderr, "Cannot allocate, exiting\n" );
            exit( -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-2, 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;
           }
         }
/* Iteration loop: */
        howlong = MPI_Wtime();
        for ( n = 0; n < niter; n++ )
        {
            if ( !my_number )
            {
                if ( !(n%STEPITER) ) printf( "Iteration %d\n", n );
            }
/* Step of calculation starts here: */
            itstep( mx, my, f, newf, r, rdx2, rdy2, beta );
/* Do all the transfers: */
            VSDA_Update_begin( &vsda, rq, newf );
            VSDA_Update_end( rq );
/* swap the halves: */      
            tmpf = newf;
            newf = f;
            f = (typeof(f))tmpf;
        }
        if ( !my_number )
        {
            printf( "Elapsed time: %f sec\n", (float)(MPI_Wtime()-howlong) );
            fp = fopen( "output_par_1.dat", "w" );
            fclose( fp );
        }
        for ( j = 0; j < n_of_nodes; j++ )
        {
            MPI_Barrier( MPI_COMM_WORLD );
            if ( j == my_number )
            {
                int istart, isize;
                istart = 1;
                isize = mx-1; 
                if ( n_of_nodes > 1 )
                {
                    if ( my_number > 0 ) istart = 1;
                    if ( my_number < (n_of_nodes-1) ) isize = mx-1;
                }
                fp = fopen( "output_par_1.dat", "a" );
                for ( i = istart; i < isize; i++ )
                {
                    fwrite( &(f[i][1]), my-2, sizeof(f[0][0]), fp );
                } 
                fclose( fp );
            }
        }
    }
    MPI_Finalize();
    return 0;
}

Немного пояснений. Программа разбита на две функции - main() и itstep(). Выше приведен текст функции main(). Текст функции itstep(), выполняющей шаг итерационного процесса над локальной порцией массивов f, newf и r, приводится ниже.

Текст функции main() до строки:

       distr_out = (long*)malloc( n_of_nodes*sizeof( *distr_out ) );

Обращением к функции VSDA_Create() создается распределенный массив. Фактически в этой программе используются три распределенных массива, f, newf и r, но, поскольку распределены они одинаково, к VSDA_Create() достаточно обратиться один раз. В обращении:

    VSDA_Create( MPI_COMM_WORLD, VSDA_BLOCKED, 0, 
                           (long)totalmx, (long)sizeof(*f), (long)1, NULL, distr_out, vsda );

1. MPI_COMM_WORLD - коммуникатор, который будет использоваться при реализации обменов краями для этого распределенного массива;
2. VSDA_BLOCKED - тип распределения данных в массиве: блоками с перекрытием. Другие типы пока не поддерживаются;
3. Признак того, что в начале локальной порции нулевого процесса, и в конце локальной порции последнего процесса, имеются не участвующие в обменах искусственные теневые грани. Задан равным нулю, что означает "искусственных теневых граней нет, наличие краевых условий программист учтет самостоятельно";
4. totalmx - размер распределенного массива в распределяемых элементах, в данном случае - в строках по my вещественных чисел в каждой;
5. sizeof(*f) - размер элемента распределенного массива (в данном случае - строки из my вещественных чисел) в байтах;
6. Размер теневой грани в распределяемых элементах (в данном случае - в строках), задан равным единице согласно используемому разностному шаблону;
7. Массив с требуемым распределением элементов распределенного массива по MPI - процессам, NULL означает "вычислить распределение самостоятельно";
8. distr_out - массив с итоговым распределением элементов распределенного массива по MPI - процессам, в i-м элементе этого массива будет находиться число строк каждого из массивов: f, newf или r, в локальной порции i-го MPI - процесса, не считая теневые грани;
9. &vsda - указатель на созданный дескриптор распределенного массива.

Теперь для выполнения обмена краями любого распределенного массива, характер распределения которого таков, как было задано при обращении к VSDA_Create(), достаточно "предъявить" соответствующей функции библиотеки VSDA созданный дескриптор и конкретную локальную порцию.

В следующем фрагменте кода:

    mx = distr_out[my_number];
    if ( n_of_nodes > 1 )
     {
      if ( my_number > 0 ) mx ++;
      if ( my_number < (n_of_nodes-1) ) mx ++;
     }

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

      VSDA_Update_begin( &vsda, rq, newf );
      VSDA_Update_end( rq );

Запрос разбит на две фазы: "начать обмен" и "дождаться конца обмена", но логика нашей программы такова, что мы никак не можем воспользоваться этим для совмещения обменов со счетом, и просто выполняем обе фазы одну за другой. Информация о характере распределения массива, края которого обмениваются, находится в дескрипторе vsda, на локальную порцию указывает newf. Запрос на обмен есть коллективная операция: все MPI - процессы выдают его логически одновременно, по ходу его выполнения происходят обмены при помощи MPI. О том, почему размер массива rq задан равным именно пяти, можно узнать из описания библиотеки VSDA.

После завершения итерационного процесса, результаты счета из всех MPI - процессов записываются в один общий файл на диске в порядке их (процессов) номеров. Делается это следующим образом.

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

Разбор текста функции main() завершен. Ниже приводится текст функции itstep().

#include "dimension.h"
      void itstep( int mx, int my, void *pf, void *pnewf, void *pr, 
                                double rdx2, double rdy2, double beta, int niter )
{
      DIM2( double, f, my ) = (typeof(f))pf;
      DIM2( double, newf, my ) = (typeof(newf))pnewf;
      DIM2( double, r, my ) = (typeof(r))pr;
      int i, j, k, mx1, my1;
      void *tmpf;
/***/
      mx1 = mx - 1;
      my1 = my - 1;
      for ( k = 0; k < niter; k++ )
       {
        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;
           }
         }
       }
}

Пояснять и комментировать в этом тексте, на наш взгляд, нечего.

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

В приведенном выше варианте программы, естественным кандидатом на перенос в сопроцессор (точнее, на аппаратную реализацию в качестве сопроцессора) является функция itstep(). Реализовав внутреннюю логику этой функции аппаратно, мы решили бы задачу создания гибридно - параллельного приложения, если бы у каждого MPI - процесса был в распоряжении собственный сопроцессор. Добиться последнего совсем легко - ведь мы всегда знаем, сколько сопроцессоров подключено к каждому вычислительному узлу нашего кластера, и имеем возможность, при помощи опций команды mpirun, разместить на каждом вычислительном узле именно по столько MPI - процессов.

Считать, что с аппаратной реализацией функции itstep() все проблемы создания гибридно - параллельного приложения решены, нам мешает крайне ограниченный объем внутренней памяти сопроцессора. Сопроцессор может обработать только то, что было предварительно записано в его внутреннюю память, объем которой - от одного до, в лучшем случае, 16 мегабайт. Это на несколько порядков меньше, чем объем памяти типичного вычислительного сервера в расчете на одно процессорное ядро (читателю, знакомому с GPGPU-кластерами, предлагаем представить себе GPGPU, в составе которого вообще нет глобальной памяти, а есть только регистры и shared memory). Чтобы уметь обрабатывать в сопроцессоре объемы данных, сопоставимые с объемом памяти на одно процессорное ядро в типичном сегодняшнем сервере, нам придется разделить массивы f, newf и r на куски, каждый из которых умещается в памяти сопроцессора, и "прогнать" эти куски через сопроцессор по очереди. По сравнению со случаем, когда обрабатываемые массивы умещаются в памяти сопроцессора целиком (так бывает, например, если сопроцессор - GPGPU), при разбиении обрабатываемых данных на куски мы очень сильно увеличиваем объем данных, передаваемых из процессора в сопроцессор и обратно. В самом деле, в первом (гипотетическом, хорошем) случае, нам достаточно было бы один раз скопировать обрабатываемые данные в сопроцессор, на каждой итерации копировать между процессором и сопроцессором только края и теневые грани, а в конце расчета - один раз выгрузить из сопроцессора насчитанные данные. В нашем же (реальном, плохом) случае, на каждой итерации придется, в общей сложности, все обрабатываемые данные один раз записать в сопроцессор, а все насчитанные - один раз выгрузить. Оправдается ли такой огромный накладной расход ускорением обработки внутри сопроцессора? Для того, чтобы грубо оценить это, нам даже не требуется ничего знать о свойствах сопроцессора. Вполне достаточно знаний о свойствах численного метода и процессора общего назначения.

Проанализировав текст itstep(), легко увидеть, что на вычисление каждого нового значения массива температур тратится 7 полезных операций над вещественными числами, и всего 3 обращения к памяти за элементами массивов (число обращений к массиву температур надо делить на 4, поскольку каждое старое значение, будучи один раз скопировано в сопроцессор, используется в вычислении 4-х новых значений). На каждое обращение к памяти, таким образом, приходится чуть больше двух полезных арифметических операций. Но обращение процессора к памяти присутствует и при копировании данных из процессора в сопроцессор, так что его нам ускорить не получится. Имеет ли смысл бороться за ускорение двух арифметических операций, пусть даже в 100 раз, на фоне не ускоренного обращения к памяти? Очевидно, нет. Чтобы копирование данных в сопроцессор ради ускорения их обработки имело смысл, число полезных арифметических операций над каждым значением, единожды скопированным в сопроцессор, придется увеличить многократно. Ради этого нам и потребуется пресловутая "глубокая структурная перестройка программы".

Прежде, чем к ней приступать, сделаем два важных замечания.

Прежде всего, опровергнем популярную гипотезу о том, что ключ наших проблем - в "отсутствии прямого доступа сопроцессора к памяти процессора" и, следовательно, таких проблем не бывает у сопроцессоров, разделяющих с процессором общую память. Оставляя в стороне тот факт, что в действительности такой "прямой доступ" существует, отметим, что для преодоления наших трудностей он совершенно бесполезен. Материал, из которого внутри ПЛИС сделана внутренняя память сопроцессора, обладает целым рядом уникальных свойств, как по скорости доступа к сделанной из него памяти, так и по количеству одновременно работающих каналов такого доступа. Ничего даже близко подобного память вычислительного сервера общего назначения делать не умеет. Верно и более сильное утверждение: если бы умела - никакие гибридно - параллельные вычислители нетрадиционной архитектуры не были бы нужны.

Также легко заметить, что так огорчившая нас слабая насыщенность обращений к данным арифметическими операциями - во многом следствие того, что мы рассматриваем двумерную задачу с очень простым разностным шаблоном, в то время как сегодня трехмерные задачи гораздо актуальнее. Здесь уместно вспомнить, как на заре параллельных вычислений сюрпризом для многих оказалась бессмысленность "распараллеливания" одномерных задач, которые, впрочем, за вполне приемлемое время решались уже тогда на однопроцессорной машине. В самом деле, при использовании трехмерного разностного шаблона с диагональным соседством ячеек сетки и порядком выше первого, мы получили бы совсем иную насыщенность обращений к памяти арифметикой. Однако, в некоторых случаях, даже для трехмерных задач, уровень насыщенности обращений к памяти арифметикой может нас не устроить. Поэтому, не забывая о том, что приложение наше - всего лишь модельное, попробуем на его примере выработать методику повышения уровня локализации вычислений в мелких блоках данных.

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