Библиотека ScaLAPACK требует, чтобы все объекты (векторы и матрицы), являющиеся параметрами подпрограмм, были предварительно распределены по процессорам. Исходные объекты классифицируются как глобальные объекты, и параметры, описывающие их, хранятся в специальном описателе объекта - дескрипторе. Дескриптор некоторого распределенного по процессорам глобального объекта представляет собой массив целого типа, в котором хранится вся необходимая информация об исходном объекте. Части этого объекта, находящиеся в памяти какого-либо процессора, и их параметры являются локальными данными.

Для того, чтобы воспользоваться драйверной или вычислительной подпрограммой из библиотеки ScaLAPACK, необходимо выполнить 4 шага:

1. инициализировать сетку процессоров;
2. распределить матрицы на сетку процессоров;
3. вызвать вычислительную подпрограмму;
4. освободить сетку процессоров.

Библиотека подпрограмм ScaLAPACK поддерживает работу с матрицами трех типов:

1. заполненными (не разреженными) прямоугольными или квадратными матрицами;
2. ленточными матрицами (узкими);
3. заполненными матрицами, хранящимися на внешнем носителе.

Для каждого из этих типов используется различная топология сетки процессоров и различная схема распределения данных по процессорам. Эти различия обуславливают использование четырех типов дескрипторов. Дескриптор - массив целого типа, состоящий из 9 или 7 элементов в зависимости от типа дискриптора. Тип дескриптора, который используется для описания объекта, хранится в первой ячейке самого дескриптора. В таблице 1.1 приведены возможные значения типов дескрипторов.

Инициализация сетки процессоров

Для объектов, которые описываются дескриптором типа 1, распределение по процессорам производится на двумерную сетку процессоров. На рис. 1.2 представлен пример двумерной сетки размера 2 х 3 из 6 процессоров.

Рис. 1.2. Пример двумерной сетки из 6 процессоров.

При таком распределении процессоры имеют двойную нумерацию:

Связь между сквозной и координатной нумерациями определяется параметром процедуры при инициализации сетки (нумерация вдоль строк - row-major или вдоль столбцов - column-major).

Инициализация сетки процессоров выполняется с помощью подпрограмм из библиотеки BLACS. Ниже приводится формат вызовов этих подпрограмм:

CALL BLACS_PINFO (IAM, NPROCS) -

инициализирует библиотеку BLACS, устанавливает некоторый стандартный контекст для ансамбля процессоров (аналог коммуникатора в MPI), сообщает процессору его номер в ансамбле (IAM) и количество доступных задаче процессоров (NPROCS).

CALL BLACS_GET (-1, 0, ICTXT) -

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

CALL BLACS_GRIDINIT (ICTXT, 'Row-major', NPROW, NPCOL) -

выполняет инициализацию сетки процессоров размера NPROW x NPCOL с нумерацией вдоль строк.

CALL BLACS_GRIDINFO(ICTXT, NPROW, NPCOL, MYROW, MYCOL) -

сообщает процессору его позицию (MYROW, MYCOL) в сетке процессоров.

Для ленточных и трехдиагональных матриц (дескриптор 501) используется одномерная сетка процессоров (1 x NPROCS), т.е. параметр NPROW = 1, а NPCOL = NPROCS.

Для объектов, описываемых дескриптором 502 (правые части уравнений с ленточными и трехдиагональными матрицами), используется транспонированная одномерная сетка (NPROCS x 1). Работа с матрицами на внешних носителях (дескриптор 601) в данном пособии не рассматривается.

Распределение матрицы на сетку процессоров

Способ распределения матрицы на сетку процессоров также зависит от типа распределяемого объекта. Для заполненных матриц (тип дескриптора 1) принят блочно-циклический способ распределения данных по процессорам. При таком распределении матрица разбивается на блоки размера MB x NB, где MB - размер блока по строкам, а NB - по столбцам, и эти блоки циклически раскладываются по процессорам. Проиллюстрируем это на конкретном примере распределения матрицы общего вида A(9,9), т.е. M=9, N=9, на сетке процессоров NPROW x NPCOL = 2 х 3 при условии, что размер блока MB x NB = 2 х 2.

Разбиение исходной матрицы на блоки требуемого размера будет выглядеть следующим образом (рис. 1,3):

Рис. 1.3. Разбиение исходной матрицы на блоки размера 2x2.

После циклического распределения блоков вдоль строк и столбцов сетки процессоров получим следующее распределение матричных элементов по процессорам (рис. 1,4):

Рис. 1.4. Распределение элементов матрицы на сетке процессоров 2x3.

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

NP = NUMROC (M, MB, MYROW, RSRC, NPROW)
NQ = NUMROC (N, NB, MYCOL, CSRC, NPCOL)

Здесь
NP - число строк в локальных подматрицах в строке процессоров MYROW;
NQ - число столбцов в локальных подматрицах в столбце процессоров MYCOL.

Входные параметры:

M,	N	-	число строк и столбцов исходной матрицы;
MB,	NB	-	размеры блоков по строкам и по столбцам;
MYROW,	MYCOL	-	координаты процессора в сетке процессоров;
IRSRC,	ICSRC	-	координаты процессора, начиная с которого выполнено распределение матрицы (подразумевается возможность распределения не по всем процессорам);
NPROW,	NPCOL	-	число строк и столбцов в сетке процессоров.

Важно также уметь оценить верхние значения величин NP и NQ для правильного описания размерностей локальных массивов. Можно использовать, например, такие формулы:

NROW = (M-1)/NPROW+MB
NCOL = (N-1)/NPCOL+NB

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

Вызов подпрограммы из PTOOLS, выполняющей аналогичное действие по формированию дескрипторов, имеет вид:

CALL DESCINIT(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, NROW, INFO)

Параметр INFO - статус возврата. Если INFO = 0, то подпрограмма выполнилась успешно; INFO = -I означает, что I-й параметр некорректен.

В расчетных формулах при решении задач линейной алгебры, как правило, фигурируют выражения, содержащие матричные элементы, идентифицируемые их глобальными индексами (I, J). После распределения матрицы по процессорам, а точнее, при заполнении локальных матриц на каждом процессоре, мы имеем дело с локальными индексами (i,j), поэтому очень важно знать связь между локальными и глобальными индексами. Так, если в сетке процессоров NPROW строк, а размер блока вдоль строк равен MB, то i-я строка локальных подматриц в MYROW строке сетки процессоров должна быть заполнена I-ми элементами строки исходной матрицы.

Формула для вычисления индекса I имеет вид:
I = MYROW*MB + ((i -1)/MB)*NPROW*MB + mod(i-1, MB) + 1,
Аналогично можно записать формулу для индесов столбцов:
J = MYCOL*NB + ((j -1)/NB)*NPCOL*NB + mod(j - 1,NB) + 1.
Здесь деление подразумевает деление нацело, а функция mod вычисляет остаток от деления.

Симметричные и эрмитовы матрицы на сетке процессоров хранятся в виде полных матриц, однако используется при этом либо верхний, либо нижний треугольники. Поэтому в процедурах, работающих с такими матрицами, указывается параметр UPLO = 'U' для верхнего треугольника или UPLO = 'L' для нижнего треугольника.

Совершенно другой принцип распределения по процессорам принят для ленточных и трехдиагональных матриц. В этом случае используется блочный принцип распределения. В каждом процессоре хранится только один блок, размер которого выбирается из соображений равномерности распределения, т.е. NB N/NPROCS. При этом объекты левой части уравнения распределяются по столбцам, а правой - по строкам.

На рис. 1.5 представлен пример ленточной матрицы A(7,7).

A =

Рис. 1.5. Пример ленточной матрицы 7х7.

Для характеристики ленточных матриц вводится параметр ширины ленты BW для симметричных матриц и два параметра BWL (ширина поддиагонали) и BWU (ширина наддиагонали) для несимметричных матриц. Так для рассмотренного на рис. 1.5 примера:

• если матрица A симметричная положительно определенная, то BW = 2;
• если матрица A несимметричная, то BWL = 2 и BWU = 2.

Кроме того, если матрица несимметричная, то возможны два варианта факторизации этой матрицы: без выбора главного элемента и с выбором главного элемента столбца. В обоих случаях матрицы хранятся в виде прямоугольных матриц и распределяются по столбцам, однако в первом случае количество строк матрицы равно BWL+BWU + 1, а во втором 2*(BWL+BWU) +1.

Схема хранения матричных элементов несимметричной матрицы при использовании процедур "без выбора главных элементов" в 3-х процессорах при условии NB = 3 представлена на рис. 1.6. Звездочками отмечены неиспользуемые позиции.

Рис. 1.6. Распределение по процессорам несимметричной ленточной матрицы "без выбора главного элемента".

При использовании процедур с "выбором главного элемента" необходимо предусмотреть выделение дополнительной памяти для рабочих ячеек. Схема распределения матричных элементов в этом случае представлена на рис. 1.7, где через F обозначена дополнительная память.

Рис. 1.7. Распределение по процессорам несимметричной ленточной матрицы "с выбором главного элемента".

В случае, когда матрица A является симметричной положительно определенной, достаточно хранить либо нижние поддиагонали (рис. 1.8), либо верхние (рис. 1.9).

Рис. 1.8. Распределение по процессорам симметричной положительно определенной ленточной матрицы (UPLO = 'L').

Рис. 1.9. Распределение по процессорам симметричной положительно определенной ленточной матрицы (UPLO = 'U').

Трехдиагональные матрицы хранятся в виде трех векторов (DU,D,DL) в случае несимметричных матриц (рис. 1.10) и двух векторов (D,E) для симметричных матриц (рис. 1.11, 1.12).

Рис. 1.10. Пример распределения по процессорам несимметричной трехдиагональной матрицы

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

Рис. 1.11. Пример распределения по процессорам симметричной трехдиагональной матрицы (UPLO = 'L').

Рис. 1.12. Пример распределения по процессорам симметричной трехдиагональной матрицы (UPLO = 'U').

Для всех рассмотренных выше матриц в качестве дескриптора используется массив целого типа из 7 элементов. Его заполнение выполняется напрямую с помощью операторов присваивания в соответствии со следующей таблицей:

Дескриптор для правых частей (табл. 1.4) имеет такую же структуру, как и для левых. Отличие состоит в том, что разложение по процессорам производится не по столбцам, а по строкам. Кроме того, во избежание лишних пересылок данных его параметры должны соответствовать параметрам дескриптора для левых частей:

Обращения к подпрограммам библиотеки ScaLAPACK

Обращение к подпрограммам ScaLAPACK рассмотрим на примере вызова подпрограммы решения систем линейных алгебраических уравнений с матрицами общего вида. Имя подпрограммыа PDGESV указывает, что:

• тип матриц - double precision (второй символ D)
• матрица общего вида, т.е. не имеет ни симметрии, ни других специальных свойств (3-й и 4-й символы GE);
• подпрограмма выполняет решение системы линейных алгебраических уравненийаааа A * X = B (последние символы SV).

Обращение к подпрограмме имеет вид:

CALL PDGESV(N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB, INFO)

Здесь

N	-	размерность исходной матрицы A (полной);
NRHS	-	количество правых частей в системе (сколько столбцов в матрице B);
А	-	на входе локальная часть распределенной матрицы A, на выходе локальная часть LU разложения;
IA, JA	-	индексы левого верхнего элемента подматрицы матрицы А, для которой находится решение (т.е. подразумевается возможность решать систему не для полной матрицы, а для ее части);
DESCA	-	дескриптор матрицы А (подробно рассмотрен выше);
IPIV	-	рабочий массив целого типа, который на выходе содержит информацию о перестановках в процессе факторизации. Длина массива должна быть не меньше количества строк в локальной подматрице;
B	-	на входе локальная часть распределенной матрицы B, на выходе локальная часть полученного решения (если решение найдено);
IB, JB	-	то же самое, что IA, JA для матрицы А;
DESCB	-	дескриптор матрицы B;
INFO	-	целая переменная, которая на выходе содержит информацию о том, успешно или нет завершилась подпрограмма, и причину аварийного завершения.

Примерно такой же набор параметров используется для всех драйверных и вычислительных подпрограмм. Назначение их достаточно понятно, следует только внимательно следить за тем, какие параметры относятся к исходным (глобальным) объектам, а какие имеют локальный характер (в первую очередь это касается размерностей массивов и векторов). Подробную информацию обо всех подпрограммах ScaLAPACK можно найти на WWW серверах [2].

Освобождение сетки процессоров

Программы, использующие пакет ScaLAPACK, должны заканчиваться закрытием всех BLACS-процессов. Это осуществляется с помощью вызова подпрограмм:

CALL BLACS_GRIDEXIT (ICTXT) - закрытие контекста;
CALL BLACS_EXIT (0) - закрытие всех BLACS процессов.

Примечание: Кроме перечисленных подпрограмм в библиотеку BLACS входят подпрограммы пересылки данных, синхронизации процессов, выполнения глобальных операций по всему ансамблю процессоров или его части (строке или столбцу). Это позволяет обойтись без использования каких-либо других коммуникационных библиотек, тем не менее, допускается и совместное использование BLACS с другими коммуникационными библиотеками. Полное описание подпрограмм библиотеки BLACS включено в документацию по пакету ScaLAPACK [2].