6. Алгоритмы и реализация распределения плотных матриц по решетке процессов.

6.1. Блочно-циклическое отображение плотных матриц в локальную память параллельных процессов.

Как уже говорилось выше, перед обращением к одной из программ Комплекса пользователю необходимо разделить исходную(ые) матрицу(ы) на блоки и распределить её(их) по параллельным процессам. Для матриц разных видов приняты свои способы разбиения и распределения матриц (см.,кроме этого п., п.6.2, а также пп. 5, 6 во второй части пособия).

Здесь мы рассмотрим наиболее общий случай прямоугольной плотной матрицы общего вида.

Для равномерной загрузки параллельных процессов и эффективного использования подпрограмм пакета BLAS при выполнении операций над плотными глобальными матрицами были разработаны специальные способы разбиения таких матриц на прямоугольные блоки ( MB * NB ) и распределения этих блоков по двумерной решетке параллельных процессов. Простой пример конкретного разбиения небольшой плотной матрицы на блоки и распределения их по двумерной решетке процессов, называемого блочно - циклическим отображением, приводится в п.6.2.

Здесь же мы рассмотрим, каким образом блоки, отображенные на один и тот же процесс, располагаются и хранятся в локальной памяти этого процесса. Другими словами, мы выпишем точные формулы, которые связывают элемент глобальной матрицы, определяемый глобальными индексами I и J (т.е. A (I, J)), с координатами процесса, владеющего этим элементом, в решетке процессов (Pr, Pc) и с расположением этого элемента матрицы в локальной памяти этого процесса.

Рассмотрим сначала, для простоты, эти связи на примере одномерного случая, т.е. размещение одномерного глобального массива длины N, разбитого на блоки длины NB, по одномерной решетке из P процессов, перенумерованных, начиная с 0 до (P - 1). Элементы самого глобального массива нумеруются от 1 до N. Прежде всего, массив делится на смежные блоки размера NB. Когда N не делится на NB нацело, последний блок массива элементов будет содержать только mod (N, NB) элементов вместо NB. Блоки, на которые разбивается исходный массив, нумеруются (как и процессы), начиная с 0. Они распределяются по процессам подобно тому, как сдается колода карт участникам игры - по кругу (т.е. циклически).

Другими словами, мы предполагаем, что процесс с номером 0 получает первый блок (с номером 0), k - ый блок связывается с процессом, имеющим координату (номер)  mod (k, P) . Блоки, связанные с одним и тем же процессом, хранятся в памяти рядом (в смежных, прилегающих частях). Отображение элемента массива с глобальным индексом  I определяется следующим соотношением:

     I   =   k * NB + x   =   ( m * P + p ) * NB + x  ,

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

Легко установить соотношения между этими переменными:

     p   =   [ ( I - 1) / NB ] mod P ,
     m   =   [ ( I - 1) / ( P * NB ) ] ,
     x   =   mod ( I - 1, NB ) + 1 .

Эти уравнения позволяют определить локальную информацию (т.е. локальный индекс  m * NB + x), так же как и координату процесса  р, соответствующую глобальному элементу, идентифицируемому его глобальным индексом  I, и наоборот.

В  таблице ниже показано отображение в локальную память при разбиении массива на блоки, когда  P = 2 ,  N = 16 и  NB = 8.

Вовсе не обязательно всегда делать распределение блоков, начиная с процесса с номером 0. Иногда, бывает полезно начать распределение данных с процесса имеющего отличную от нуля координату  SRC. В этом случае соотношения становятся такими:

     p   =   ( SRC + [ ( I - 1) / NB ] ) mod P ,
     m   =   [ ( I - 1) / ( P * NB ) ] ,
     x   =   mod ( I - 1, NB ) + 1 .

Ниже в таблице показано размещение блоков при  P = 2,  SRC = 1,  NB = 3 и  N = 16.

В двумерном случае предположим, что матрица разделена на  MB * NB блоки и что первый блок передан процессу с координатами ( RSRC, CSRC ). Формула, приведенная выше для одномерного случая, должна использоваться повторно независимо для каждого измерения решетки процессов  Pr * Pc. Например, элемент матрицы ( I, J ) находится в процессе с координатами ( Pr, Pc) внутри локального блока ( m, n ) в позиции ( x, y ), задаваемой формулами

  ( m, n )   =   ( [ ( I - 1) / ( Pr * MB ) ], [ ( J - 1) / ( Pc * NB ) ] ) ,

  ( Pr, Pc)   =   ( ( RSRC + [ ( I - 1) / MB ] ) mod Pr, ( CSRC + [ ( J - 1) / NB ] ) mod Pc )

  ( x, y )   =   ( mod ( I - 1, MB ) + 1, mod ( J - 1, NB ) + 1 )

Эти формулы показывают, как матрица  А размера M_A * N_A отображается и хранится на решетке процессов. Она сначала разбивается на MB_A * NB_A блоки, начиная с ее верхнего левого угла. Эти блоки затем равномерно распределяются по решетке процессов циклическим образом.

Каждый процесс владеет набором блоков, которые хранятся рядом (смежно) по столбцам в двумерном массиве, расположенном в памяти "по столбцам".

Это соглашение о локальном размещении позволяет эффективно использовать иерархию локальной памяти посредством вызова пакета BLAS для подмассивов, которые могут быть больше, чем один блок размера MB_A * NB_A.

На рисунке ниже представлено отображение матрицы  5 * 5, разделенной на блоки размером  2 * 2 на решетку процессов размером 2 * 2 (т.е. M_A = N_A = 5, Pr = Pc = 2 и MB_A = NB_A = 2). Локальные элементы каждого столбца матрицы хранятся рядом в памяти каждого процесса.

Цифры 0 и 1 в левой и верхней части последней таблицы означают номера строк и столбцов решетки процессов соответственно.

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

Для локального числа строк используется обозначение LOCr ( ), а для локального числа столбцов - LOCc ( ). Значения LOCr ( ) и LOCc ( ), получаемые подпрограммой NUMROC, являются результатом точных вычислений.

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

а) верхняя граница для  LOCr ( ) оценивается по формуле

                             M_A + MB_A - 1
                             ----------------------    +   Pr  -  1
                                  MB_A
    LOCr( )   =   -------------------------------------------    *   MB_A
                                              Pr

или 

    LOCr( )   =   [ [ M_A / MB_A] / Pr ] * MB_A


б) величина LOCc( ) оценивается по формуле

                              N_A + NB_A - 1
                              ----------------------    +   Pc  -  1
                                 NB_A
    LOCc( )   =   --------------------------------------------    *   NB_A
                                                 Pc
 
или

    LOCc( )   =   [ [ N_A / NB_A ] / Pc ] * NB_A

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

6.2. Пример блочно-циклического распределения плотной матрицы по решетке процессов.

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

В соответствии с принятой схемой, сначала плотная матрица  A размера M * N разделяется на блоки размера MB * NB начиная с левого верхнего угла этой матрицы. Эти блоки затем равномерно распределяются по каждому измерению решетки процессов. Точные математические соотношения, соответствующие такой схеме распределения, приводятся в п.6.1.

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

Ниже на рисунке показано разделение матрицы  A размером 9 * 9 на блоки размером 2 * 2.

Далее показано отображение полученных блоков на решетку процессов размером 2 * 3 (двумерное блочно - циклическое распределение данных).

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

Пусть первый блок распределяется в процесс (0, 0). Тогда все блоки, расположенные с ним в той же строке будут распределяться в ту же строку решетки процессов, т.е. первая координата в этих клетках будет равна 0.

Вторая же координата (столбца решетки) будет циклически изменяться: 0, 1, 2, 0, 1, т.к. столбцов в решетке только 3 (с координатами 0, 1, 2).

Все блоки, расположенные с первым блоком в том же самом столбце, будут распределяться в тот же самый столбец решетки процессов. Т.е. вторая координата в этих клетках будет равна 0. Первая же координата (строки решетки) будет циклически изменяться: 0, 1, 0, 1, 0, т.к. число строк в решетке только 2 (с координатами 0 и 1).

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

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

Сборка клеток (блоков) с одинаковыми координатами делается так:

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

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

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

Ниже в таблице указаны характеристики локальных массивов для каждого из процессов решетки.

     Координаты
     процесса             LLD_A       LOCr ( M_A )       LOCc ( N_A )
    _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
         (0,0)                      5                     5                             4
         (0,1)                      5                     5                             3
         (0,2)                      5                     5                             2
 
         (1,0)                      4                     4                             4
         (1,1)                      4                     4                             3
         (1,2)                      4                     4                             2

Число строк LOCr и число столбцов LOCc матрицы  A, которыми владеет конкретный процесс, могут отличаться у разных процессов в решетке. Подобно этому, для каждого процесса в решетке процессов существует ведущая локальная размерность LLD. Ее величина может быть различной для разных процессов в решетке процессов. Например, как мы можем видеть на рисунке выше, локальный массив, хранящийся в строке решетки процессов с номером 0, должен иметь ведущую локальную размерность LLD не меньше 5, а хранящийся в строке с номером 1 - не меньше 4. Подробнее о ведущей локальной размерности LLD см. в п.5.2.

6.3. Программирование блочно-циклического распределения плотных матриц.

В настоящем пособии рассматриваются три способа реализации на языке ФОРТРАН рассылки элементов исходных глобальных матриц по параллельным процессам, поскольку каждый из них удобнее использовать в трех разных ситуациях.

1. Первый способ удобен только для небольших матриц и требует от пользователя самостоятельного размещения элементов исходных матриц в локальной памяти параллельных процессов в соответствии с общими правилами, описанными в пп.6.1, 6.2, а также во второй части пособия в пп. 4, 5, 6. При этом используются обычные операторы присваивания языка ФОРТРАН. Такой способ достаточно трудоемок и может использоваться для целей обучения и усвоения общих правил распределения матриц (см. п.6.1, а также во второй части пособия пп. 5, 6,4). Фортранная подпрограмма такого способа распределения приводится, например, во второй части пособия в конце Приложения 2_1.

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

Этот способ гораздо проще для пользователя, поскольку не требует самостоятельного определения местоположения элементов матриц в локальной памяти каждого из параллельных процессов. Он опирается на использование служебной подпрограммы PDELSET, которая сама (в соответствии с упомянутыми общими правилами, см. п.6.1) заносит элемент с индексами  I, J исходной матрицы  A ( I, J ) в определенное место локальной памяти того или другого из параллельных процессов.

На ФОРТРАНе это выливается в написание некоторого числа циклов, внутри которых стоит обращение к подпрограмме PDELSET.

Фортранные подпрограммы с реализацией такого способа распределения приводятся во второй части методического пособия в конце Приложения 2_2, а также в третьей части в конце Приложения 3_2.

Еще раз подчеркнем, что хотя в конце Приложения 2_1 и Приложения 2_2 приводятся два разных способа распределения одной и той же исходной матрицы, результат распределения будет одним и тем же. Подробный разбор приводимого в этих приложениях примера дается в п.4.

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

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

Пример обращения к такой (служебной) фортранной подпрограмме (PDLAREAD) приводится во второй части пособия в Приложении 2_3. Однако, пользователь может по образу и подобию написать свою аналогичную подпрограмму, если какие - либо детали размещения исходной матрицы в файле отличаются от приводимого примера.

В этом же примере содержится обращение к служебной подпрограмме записи полученных результатов в файл (PDLAWRITE). Тексты этих подпрограмм можно взять здесь PDLAREAD, PDLAWRIT.