3. Дескриптор вектора (матрицы) правой части системы уравнений.

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

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

Для случая ленточных и трехдиагональных глобальных матриц  A, блочно - столбцовое распределение которых по процессам описывается с помощью дескрипторов типа 501 (см.), распределение глобального вектора (матрицы) правых частей  B должно быть описано с помощью дескриптора типа 502. Он, также как и дескриптор типа 501, состоит их 7 элементов, имеющих аналогичные закрепленные символические имена. Ниже приведена таблица для элементов дескриптора типа 502.

В отличие от глобальной ленточной матрицы  A, которая подвергается блочно - столбцовому распределению (описываемому дескриптором типа 501), глобальный вектор (матрица)  B подвергается так называемому блочно - строчному распределению. Это означает, что вектор (матрица)  B разбивается на блоки по строкам (каждая строка целиком попадает в локальную память одного из процессов). В том случае, если B - вектор, строка вырождается в один элемент, т.е. блокируются несколько соседних элементов вектора  B. Подробнее см. в п.5.4. Конкретный пример распределения вектора  B по процессам можно найти в разделе "Пример использования" описания одной из подпрограмм для решения систем с ленточной матрицей  A (см.).

Инициализацию дескрипторов глобальных матриц  A и  B необходимо сделать до обращения к подпрограммам Комплекса. Для плотных матриц, распределение которых описывается дескрипторами типа 1, инициализацию проще всего выполнить с помощью специальной служебной подпрограммы Комплекса с именем DESCINIT, которой в качестве параметров необходимо передать значения всех элементов дескриптора. В первой части пособия в п.4.2. приводится пример обращения к этой подпрограмме при инициализации дескриптора DESCA для плотной матрицы коэффициентов системы  A.

Приведём здесь обращение к этой подпрограмме при инициализации дескриптора DESCB для матрицы (вектора) правых частей системы уравнений  B с плотной матрицей системы A. Подобные обращения можно найти в текстах тестовых примеров к программам решения системы уравнений с плотными матрицами ( TPDGESV.f, TPDPOSV.f ).

     CALL DESCINIT (DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC,
                                    ICTXT, MXLLDA, INFO)

Этот вызов подрограммы DESCINIT эквивалентен следующим операторам присваивания:

        DESCB(1) = 1
        DESCB(2) = ICTXT
        DESCB(3) = N
        DESCB(4) = NRHS
        DESCB(5) = NB
        DESCB(6) = NBRHS
        DESCB(7) = RSRC
        DESCB(8) = CSRC
        DESCB(9) = MXLLDB

Здесь первый элемент обозначает тип дескриптора, второй - служебный параметр подпрограммы BLACS'а (контекст).

Примеры инициализации дескрипторов DESCA и DESCB, имеющих типы 501 и 502, можно найти в разделе "Пример использования" описания одной из подпрограмм для решения систем с ленточной матрицей  A (см.).

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

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

Пример демонстрирует решение системы уравнений  A * X = B с матрицей  A размера 9 * 9. Ниже показана матричная форма записи этой системы.

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

Сначала мы должны задать решетку процессов и распределить по ней исходные глобальные матрицы  A и  B.

Предположим, что матрица  A разделена на блоки размером MB * NB, где  MB = NB = 2. Выбрана решетка процессов размером 2 * 3, т.е. NPROW = 2 и NPCOL = 3. Таким образом имеет место случай распределения матрицы A, приведенный в п.6.2 в первой части методического пособия. После подстановки конкретных значений элементов матрицы  A получаем следующую картинку распределения ее элементов по процессам. (Вверху указаны номера столбцов решетки процессов, слева - номера строк решетки процессов.)

Из последнего рисунка видно, что процесс с координатами (0, 0) содержит локальный массив (часть матрицы  A) размера (5, 4).

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

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

После вызова подпрограммы PDGESV и выполнения всех необходимых вычислений процесс (0, 0) будет содержать на месте локальной части вектора  B локальную часть глобального вектора решений  X. Ниже изображено, на месте каких элементов локального вектора  B (обозначенных как ^~b_i ) с локальными индексами, какие элементы глобального вектора  X с глобальными индексами оказываются расположенными по окончании вычислений. Это расположение соответствует расположению элементов исходного глобального вектора правых частей  B, показанному на предыдущем рисунке. Ниже показаны также вычисленные значения элементов вектора  X.

Аналогичное соответствие и полученные результаты изображены далее для процесса с координатами (1, 0).

Таким образом глобальный выходной вектор решений  X имеет следующее значение:

Кроме того, выполняется проверка точности полученных результатов вычислением нормализованной невязки по формуле:

                            || A * x - b ||
                     ______________________   .
                     ( || x || * || A || * eps * N )

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

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

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

      AIJ = 19.0D0

      DO 20 J = 1, N
         DO 10 I = 1, N
            IF( I.LE.J ) THEN
                CALL PDELSET( A, I, J, DESCA, AIJ )
            ELSE
                CALL PDELSET( A, I, J ,DESCA, -AIJ )
            END IF
 10 CONTINUE
         IF( J.EQ.1 .OR. J.EQ.7 ) AIJ = 3.0D0
         IF( J.EQ.2 .OR. J.EQ.4 .OR. J.EQ.6 ) AIJ = 1.0D0
         IF( J.EQ.3 ) AIJ = 12.0D0
         IF( J.EQ.5 ) AIJ = 16.0D0
         IF( J.EQ.8 ) AIJ = 11.0D0
 20 CONTINUE

      CALL PDELSET( A, 8, 2, DESCA, 3.0D0 )

      BIJ = 0.0D0

      J = 1
      DO 30 I = 1, N
        CALL PDELSET( B, I, J, DESCB, BIJ )
 30 CONTINUE

      CALL PDELSET( B, 3, 1, DESCB, 1.0D0 )

Полностью весь фортранный текст рассматриваемого примера приведен в Приложении 2_2.

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