Введение в комплекс программ решения линейных систем на распределенной памяти

Одна из задач, решаемых комплексом, - это линейные системы (системы линейных алгебраических уравнений)

(1)                            A X = B ,
             
где A - невырожденная квадратная матрица,
       B - либо вектор правой части, если решается единственная (одна)
             система, либо матрица, состоящая из столбцов правых частей,
             если решается несколько систем с одной и той же матрицей A, 
       X - вектор неизвестных, если B - вектор, или матрица, состоящая
             из векторов неизвестных, если B - матрица. 
Под векторами понимаются векторы - столбцы.

Для решения системы (1) применяются специальные блочные реализации (блочные алгоритмы) для прямого и обратного исключения Гаусса.

 Решение системы разбито на два этапа:
 - факторизация матрицы системы;
 - непосредственно само решение системы с использованием
    результатов, полученных на первом этапе.

Каждый из этапов осуществляется посредством обращения к специальным базовым подпрограммам комплекса.

Виды факторизации, зависящие от свойств матрицы  A, представлены ниже.

1. Для матриц общего вида применяется LU-факторизация с выбором ведущего
    элемента по столбцам:

          A  =  P L U ,

где  P - матрица перестановок;
       L - нижняя треугольная матрица с единичными диагональными
             элементами;
       U - верхняя треугольная матрица.

2. Для симметричных или эрмитовых положительно определенных
    матриц применяется разложение Холецкого:

           A  =  UT U   или   A  =  L LT  (для симметричных) ,
           A  =  UH U   или   A  =  L LH  (для эрмитовых) ,

где  U - верхняя треугольная матрица;
        L - нижняя треугольная матрица.

3. Для ленточных матриц общего вида применяется LU-факторизация с выбором
    ведущего элемента по столбцам:

           A  =  P L U Q ,

где  P  и  Q - матрицы перестановок,
       L  и  U - нижняя треугольная и верхняя треугольная ленточные
                     матрицы соответственно.

4. Для ленточных матриц общего вида с диагональным преобладанием, включая
    трехдиагональные матрицы общего вида, применяется LU-факторизация без выбора
    ведущего элемента:

            A  =  P L U PT ,

где  P        - матрица перестановок,
       L и U - нижняя и верхняя треугольные ленточные матрицы
                    соответственно.

5. Для симметричных и эрмитовых положительно определенных ленточных
    матриц применяется разложение Холецкого:

     A  =  P UT U PT  или  A  =  P L LT PT (для симметричных) ,
     A  =  P UH U PT  или  A  =  P L LH PT (для эрмитовых) ,

где  P        - матрица перестановок,
       U и L - верхняя и нижняя треугольные ленточные матрицы
                    соответственно.

6. Для симметричных и эрмитовых положительно определенных
    трехдиагональных матриц применяется  LDLT-факторизация:

  A  =  P UT D U PT  или  A  =  P L D LT PT (для симметричных) ,
  A  =  P UH D U PT  или  A  =  P L D LH PT (для эрмитовых) ,

где  P        - матрица перестановок,
       U и L - двухдиагональные верхняя и нижняя треугольные матрицы
                    соответственно.

Указанные факторизации выполняются базовыми подпрограммами комплекса:

 - для плотных матриц общего вида и симметричных положительно
    определенных - PDGETRF и PDPOTRF;

 - для ленточных матриц общего вида - PDGBTRF (с выбором ведущего
    элемента по столбцу) и PDDBTRF (без выбора ведущего элемента);

 - для симметричных/эрмитовых ленточных положительно определенных
    матриц - PDPBTRF;

 - для трехдиагональных матриц общего вида - PDDTTRF (без выбора
    ведущего элемента);

 - для симметричных/эрмитовых трехдиагональных положительно
    определенных матриц - PDPTTRF.

Перед решением системы матрица  A  разбивается пользователем на квадратные или прямоугольные блоки. После этого пользователь должен задать число параллельных процессов (процессоров) для решения задачи (1). Затем он должен разместить указанные блоки матрицы  A  по процессам (или в памяти процессоров) с применением блочно - циклического распределения (если  A - плотная матрица) или более простого блочного распределения (если  A - трехдиагональная или ленточная). Похожие действия необходимо выполнить и с вектором (или матрицей)  B.

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

Из сказанного выше следует, что исходную матрицу системы  A  не требуется (и не принято) хранить целиком в оперативной памяти процессора как единый двумерный массив. На практике такие матрицы в виде единого целого могут храниться во внешних файлах (на дисковой памяти). По этой причине такую матрицу системы  A  называют глобальной матрицей системы. Аналогично,  B  называют глобальным вектором правой части (или глобальной матрицей правых частей). Поскольку перед началом работы параллельных подпрограмм комплекса глобальные матрицы (векторы) распределяются по параллельным процессам, то в дальнейшем нередко употребляется термин "распределенная глобальная матрица (вектор)". Часть глобальной матрицы (вектора), представленная в локальной памяти одного из параллельных процессов, называется локальной частью матрицы (вектора). Совокупность всех локальных частей на всех, используемых при решении задачи параллельных процессах, и составляет распределенную глобальную матрицу (вектор).

Подпрограммы комплекса составлены на ФОРТРАНе - 77 в стиле SMPD ( Single Program Multiple Data). Это означает, что одна и та же исполняемая программа загружается во все параллельные процессы, заказанные для решения задачи в команде запуска на счет. В то же время, каждый процесс занимается обработкой своих данных, составляющих некоторую часть общих (глобальных) данных задачи. Эти данные отдельного процесса (их называют локальными данными) непосредственно недоступны другим процессам и хранятся в отдельной (локальной) памяти каждого процесса. Необходимый обмен данными между разными процессами производится только посредством передачи сообщений (технология MPI).

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

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

Сказанное выше относится не только к матрице коэффициентов  A, но также и к вектору (матрице) правых частей  B и вектору (матрице) решений  X.

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

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

По информации в дескрипторах, а также учитывая общие правила распределения блоков матриц (векторов) по параллельным процессам всегда можно определить, в локальной памяти какого процесса и где именно находится элемент глобальной матрицы с глобальными индексами I  и  J. См. подробнее в разделе документации "Разбиение на блоки и распределение исходных матриц в локальной памяти параллельных процессов".

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

Множество параллельных процессов, используемых для решения задачи, образует решетку процессов ( process grid ), которая может быть как двумерной, так и одномерной. Каждый отдельный процесс идентифицируется парой целых чисел (координат), определяющих его положение в решетке процессов: (номер строки, номер столбца). Нумерация начинается с номера 0.

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

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

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

Формальные параметры подпрограмм комплекса могут быть глобальными или локальными. Глобальными называются параметры подпрограмм, которые относятся к глобальным объектам исходной системы уравнений (1). Например, число строк или столбцов исходной матрицы  A  (M, N) или число наддиагоналей ( BWU ) в  A, если она является ленточной.

Локальными называются параметры, которые характеризуют объекты, определяемые в локальной памяти отдельного процесса. Например, указатель на локальную память, занимаемую локальной частью глобальной матрицы  A  или глобального вектора правых частей  B.

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