5. Блочно-столбцовое разбиение и схема размещения в локальной памяти
ленточных матриц с примером реализации размещения на языке ФОРТРАН.

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

Размер блока выбирается исходя из числа столбцов N глобальной матрицы  A и числа процессоров P, используемых для решения задачи. При этом используется одномерная решетка процессов размером 1 * P. Каждый процесс должен получить один из блоков матрицы. Таким образом, размер блока NB определяется из формулы: NB і [ N / P ]. Если N делится на P нацело, все блоки, распределяемые на все процессы, имеют одинаковый размер. В противном случае последний блок, распределяемый на последний процесс, имеет размер < NB. При этом K - й столбец матрицы A попадает в память процесса с номером [ K / NB ]. Как уже говорилось, процессы нумеруются от 0 до (P - 1). Координаты первого процесса в такой решетке есть (0,0), а последнего - (0, P - 1). На рисунке условно показано разбиение некоторой глобальной матрицы A размера M * N на 3 блока по столбцам для распределения ее по трем процессам (координаты процессов в решетке указаны под соответствующими блоками).

Такой способ распределения матрицы  A позволяет использовать в подпрограммах высоко эффективные алгоритмы, называемые "разделяй и властвуй" [см.].

Рассмотрим схему размещения элементов глобальной ленточной матрицы  A в локальной памяти параллельных процессов на конкретных примерах.

Пусть матрица  A имеет размеры 7 * 7 и ширину ленты, равную 2 :

Пусть матрица распределяется по решетке процессов 1 * 3 с размером блоков NB_A = 3.

Рассмотрим несколько случаев распределения по решетке процессов несимметричных и симметричных положительно определенных ленточных матриц.

5.1. Случай несимметричной ленточной матрицы, при факторизации которой используются алгоритмы без выбора ведущего элемента.

Пусть число поддиагоналей BWL матрицы  A равно 2, а число наддиагоналей BWU равно 2.

Ниже на рисунке представлено распределение такой матрицы по процессам. Символом "_*" отмечены элементы в двумерном локальном массиве, которые в дальнейших вычислениях не используются.

Из этого рисунка видно, что ведущая размерность этих локальных массивов (т.е. число строк в блоках глобальной матрицы  A, которые помещаются в каждый процессор) должна быть равна BWL + 1 + BWU (т.к. элементы изображенных массивов заносятся в локальную память по столбцам). О ведущей размерности локальных массивов см. в первой части пособия в п.5.2.2.

5.2. Случай несимметричной ленточной матрицы, при факторизации которой используются алгоритмы с выбором ведущего элемента по столбцам.

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

В рассматриваемом нами примере имеем BWL = 2 и BWU = 2. Ниже на рисунке представлено распределение глобальной ленточной матрицы по процессам с выделением необходимого дополнительного пространства. Дополнительно резервируемые элементы в каждом локальном массиве обозначены буквами "F".

Из этого рисунка видно, что ведущая размерность такого локального массива должна быть равна 2 * (BWL + BWU) + 1.

5.3. Случай симметричной положительно определенной ленточной матрицы, у которой в памяти сохраняется только нижний треугольник (параметр UPLO в подпрограммах обработки таких матриц полагается равным "L").

Пусть ширина ленты такой матрицы BW = 2. Ниже на рисунке представлено распределение такой матрицы ( 7 * 7 ) по той же решетке процессов ( 1 * 3 ).

5.4. Случай симметричной положительно определенной ленточной матрицы, у которой в памяти сохраняется только верхний треугольник (параметр UPLO в подпрограммах обработки таких матриц полагается равным "U").

Ниже на рисунке представлено распределение такой матрицы по той же решетке процессов ( 1 * 3 ).

Из двух последних рисунков видно, что ведущая размерность такого локального массива равна BW + 1.

Кроме этого, при решении систем линейных уравнений предполагается, что вектор (матрица) правых частей B (размером  N * NRHS) распределяется по той же решетке процессов, что и описанные выше матрицы  A, но только не блочно - столбцовым, а блочно - строчным образом. Это означает, что в блоки объединяются соседние строки матрицы B (или просто соседние элементы, если B - вектор). Ниже на рисунке условно представлено блочно - строчное разбиение матрицы (вектора)  B для распределения ее по той же решетке процессов ( 1 * 3 ), что и соответствующая ленточная (или трехдиагональная) матрица  A.

               1                     NRHS              Процессы
               _ _ _ _ _ _ _ _ _ _ _
            1 |                               |
               |                               |                  (0,0)
               |                               |
               _ _ _ _ _ _ _ _ _ _ _
               |                               |
               |                               |                  (0,1)
               |                               |
               _ _ _ _ _ _ _ _ _ _ _
               |                               |
               |                               |                  (0,2)
           N |                               |
               _ _ _ _ _ _ _ _ _ _ _

Распределение элементов вектора B по этим процессам будет следующим.

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

Приведем здесь фортранные операторы реализующие 2 - ой способ распределения матрицы  A вида

с помощью служебной подпрограммы PDELSET.

6. Схема размещения в локальной памяти трехдиагональных матриц
с примером реализации размещения на языке ФОРТРАН.

Рассмотрим пример трехдиагональной глобальной матрицы  A размером 7 * 7 ( N_A = 7)

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

6.1. Если матрица  A несимметричная, она представляется в памяти в виде трех отдельных векторов (D, DL, DU), каждый из которых содержит элементы одной из диагоналей.

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

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

                                          Процессы

Символом "_*" отмечены элементы, которые не используются при вычислениях.

6.2. Если трехдиагональная матрица  A является симметричной положительно определенной, она должна быть представлена в памяти в виде двух векторов, один из которых (D), содержит элементы главной диагонали, а другой (E) - элементы одной из кодиагоналей.

Если предполагается, что в памяти сохраняется нижняя диагональ матрицы  A (параметр UPLO в соответствующей подпрограмме полагается равным "L"), то вектор E содержит элементы нижней диагонали.

Если предполагается, что в памяти сохраняется верхняя диагональ матрицы  A (параметр UPLO полагается равным "U"), то вектор E содержит элементы нижней диагонали. Оба вектора D и E должны иметь одинаковую длину, равную длине вектора D, т.е. быть выровнены. При этом элементы верхней или нижней диагонали сдвигаются к началу вектора E, т.е. последний элемент кодиагонали становится предпоследним элементом вектора E.

Сказанное проиллюстрировано на двух рисунках ниже (первый - для UPLO = "L", второй - для UPLO = "U"). На них также показано распределение векторов, разбитых на блоки длиной NB_A = 3, в локальную память каждого из трех процессов.

                                          Процессы

                                          Процессы

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

Приведем здесь фортранные операторы реализующие 2 - ой способ распределения матрицы  A вида:

с помощью служебной подпрограммы PDELSET.

      T = 2.0D0
      MO = -1.0D0
      Z = 0.0D0

      DO 10 J = 1, N-1
       CALL PDELSET( D, 1, J, DESCD, T )
 10 CONTINUE

      DO 30 J = 2, N
      CALL PDELSET( DL, 1, J, DESCD, MO )
 30 CONTINUE

      DO 20 J = 1, N-1
      CALL PDELSET( DU, 1, J, DESCD, MO )
 20 CONTINUE

      CALL PDELSET( D, 1, N, DESCD, 1.D0 )