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

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

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

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

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

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

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

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

Рис. 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. Разбиение на блоки будет выглядеть следующим образом:

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

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

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

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

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

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 вычисляет остаток от деления. Подразумевается индексация массивов от 1 (как принято в фортране).

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

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

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

A =

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

Для характеристики этой матрицы вводятся параметры ширины ленты:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1. тип матриц - double precision (второй символ D)
2. матрица общего вида, т.е. не имеет ни симметрии, ни других специальных свойств (3-й и 4-й символы GE);
3. подпрограмма выполняет решение системы линейных алгебраических уравненийаааа 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 сервере [4].

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

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

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

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