Как уже говорилось выше, перед обращением к одной из программ Комплекса пользователю необходимо разделить исходную(ые) матрицу(ы) на блоки и распределить её(их) по параллельным процессам. Для матриц разных видов приняты свои способы разбиения и распределения матриц (см.,кроме этого п., п.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.
I | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 |
p | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
m | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
x | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
m*NB + x | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 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.
I | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 |
p | 1 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 |
m | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 |
x | 1 | 2 | 3 | 1 | 2 | 3 | 1 | 2 | 3 | 1 | 2 | 3 | 1 | 2 | 3 | 1 |
m*NB + x | 1 | 2 | 3 | 1 | 2 | 3 | 4 | 5 | 6 | 4 | 5 | 6 | 7 | 8 | 9 | 7 |
В двумерном случае предположим, что матрица разделена на 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). Локальные элементы каждого столбца матрицы хранятся рядом в памяти каждого процесса.
a11 a12 a21 a22 |
a13 a14 a23 a24 |
a15 a25 | |
a31 a32 a41 a42 |
a33 a34 a43 a44 |
a35 a45 | |
a51 a52 | a53 a54 | a55 |
0 | 1 | ||
0 |
a11 a12 a21 a22 |
a15 a25 |
a13 a14 a23 a24 |
a51 a52 | a55 | a53 a54 | |
1 |
a31 a32 a41 a42 |
a35 a45 |
a33 a34 a43 a44 |
Цифры 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
Заметим, что эти вычисления могут привести к очень большой переоценке величины действительно требующегося пространства.
Приводимый в настоящем разделе пример является простой иллюстрацией блочно - циклического распределения плотной глобальной матрицы A по двумерной решетке процессов.
В соответствии с принятой схемой, сначала плотная матрица A размера M * N разделяется на блоки размера MB * NB начиная с левого верхнего угла этой матрицы. Эти блоки затем равномерно распределяются по каждому измерению решетки процессов. Точные математические соотношения, соответствующие такой схеме распределения, приводятся в п.6.1.
Таким образом каждый процесс владеет набором блоков, которые расположены в его локальной памяти рядом, в двумерном массиве, хранящемся по столбцам (см. ниже).
Ниже на рисунке показано разделение матрицы A размером 9 * 9 на блоки размером 2 * 2.
a11 a12 a21 a22 |
a13 a14 a23 a24 |
a15 a16 a25 a26 |
a17 a18 a27 a28 |
a19 a29 | |
a31 a32 a41 a42 |
a33 a34 a43 a44 |
a35 a36 a45 a46 |
a37 a38 a47 a48 |
a39 a49 | |
a51 a52 a61 a62 |
a53 a54 a63 a64 |
a55 a56 a65 a66 |
a57 a58 a67 a68 |
a59 a69 | |
a71 a72 a81 a82 |
a73 a74 a83 a84 |
a75 a76 a85 a86 |
a77 a78 a87 a88 |
a79 a89 | |
a91 a92 | a93 a94 | a95 a96 | a97 a98 | a99 |
Далее показано отображение полученных блоков на решетку процессов размером 2 * 3 (двумерное блочно - циклическое распределение данных).
Чтобы понять, как распределятся изображенные выше блоки по процессам решетки, сначала напишем на каждой из клеток (блоков) рисунка координаты процессов, куда они должны быть распределены в соответствии с установленным блочно - циклическим распределением, описанным в п.6.1.
(0,0) | (0,1) | (0,2) | (0,0) | (0,1) | |
(1,0) | (1,1) | (1,2) | (1,0) | (1,1) | |
(0,0) | (0,1) | (0,2) | (0,0) | (0,1) | |
(1,0) | (1,1) | (1,2) | (1,0) | (1,1) | |
(0,0) | (0,1) | (0,2) | (0,0) | (0,1) |
Пусть первый блок распределяется в процесс (0, 0). Тогда все блоки, расположенные с ним в той же строке будут распределяться в ту же строку решетки процессов, т.е. первая координата в этих клетках будет равна 0.
Вторая же координата (столбца решетки) будет циклически изменяться: 0, 1, 2, 0, 1, т.к. столбцов в решетке только 3 (с координатами 0, 1, 2).
Все блоки, расположенные с первым блоком в том же самом столбце, будут распределяться в тот же самый столбец решетки процессов. Т.е. вторая координата в этих клетках будет равна 0. Первая же координата (строки решетки) будет циклически изменяться: 0, 1, 0, 1, 0, т.к. число строк в решетке только 2 (с координатами 0 и 1).
Другими словами, каждая координата независимо от другой (в строках - слева направо, а в столбцах - сверху вниз), циклически изменяется.
Если теперь мы соединим вместе (пристыкуем) все клетки (блоки), имеющие одинаковые координаты процесса, то получим массив элементов, который и должен расположиться в локальной памяти процесса с такими координатами.
Сборка клеток (блоков) с одинаковыми координатами делается так:
Тем самым получаем единый двумерный массив, расположенный в локальной памяти процесса с указанными в клетках координатами.
Ниже представлено расположение элементов исходной матрицы во всех процессах решетки после завершения процесса блочно - циклического распределения. Вверху указаны номера столбцов решетки процессов, слева - номера строк решетки процессов.
0 | 1 | 2 | |||
0 |
a11 a12 a21 a22 |
a17 a18 a27 a28 |
a13 a14 a23 a24 |
a19 a29 |
a15 a16 a25 a26 |
a51 a52 a61 a62 |
a57 a58 a67 a68 |
a53 a54 a63 a64 |
a59 a69 |
a55 a56 a65 a66 | |
a91 a92 | a97 a98 | a93 a94 | a99 | a95 a96 | |
1 |
a31 a32 a41 a42 |
a37 a38 a47 a48 |
a33 a34 a43 a44 |
a39 a49 |
a35 a36 a45 a46 |
a71 a72 a81 a82 |
a77 a78 a87 a88 |
a73 a74 a83 a84 |
a79 a89 |
a75 a76 a85 a86 |
Ниже в таблице указаны характеристики локальных массивов для каждого из процессов решетки.
Координаты процесса 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.
В настоящем пособии рассматриваются три способа реализации на языке ФОРТРАН рассылки элементов исходных глобальных матриц по параллельным процессам, поскольку каждый из них удобнее использовать в трех разных ситуациях.
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.