Подпрограмма:  PDDBSV

Назначение

Решение системы с ленточной матрицей общего вида методом Гаусса без выбора ведущих элементов

Математическое описание

Подпрограмма вычисляет решение вещественной системы линейных уравнений 

        A(1 : N, JA : JA+N-1) * X = B( IB : IB+N-1, 1 : NRHS) ,

  где  A(1 : N, JA : JA+N-1) - квадратная вещественная ленточная распределенная матрица
                                              порядка N с диагональным преобладанием и 
                                              с шириной ленты BWL, BWU.

Используется метод исключения Гаусса без выбора ведущего элемента для разложения матрицы на множители L и U (нижняя и верхняя треугольные матрицы).

Матрицы A и B должны быть заданы в виде массивов двойной точности. Все вычисления проводятся в режиме DOUBLE PRECISION.

Литература:

http://www.netlib.org/scalapack/slug/index.html

A. Cleary and J. Dongarra, Implementation in scalapack of divide - and - conquer algorithms for banded and tridiagonal linear systems, Computer Science Dept. Technical Report CS - 97 - 358, University of Tennessee, Knoxville, TN, April 1997. (Also LAPACK Working Note #125).

Использование 

         CALL  PDDBSV ( N, BWL, BWU, NRHS, A, JA, DESCA, B, IB, DESCB,
                                        WORK, LWORK, INFO)

   Параметры
N - порядок распределенной подматрицы A (1 : N, JA : JA + N - 1), N ≥ 0 (глобальный входной параметр, тип целый);
BWL - число поддиагоналей, 0 ≤ BWL ≤ N - 1 (глобальный входной параметр, тип целый);
BWU - число наддиагоналей, 0 ≤ BWU ≤ N - 1 (глобальный входной параметр, тип целый);
NRHS - число правых частей, т.е. число столбцов распределенной подматрицы B ( IB : IB + N - 1, 1 : NRHS), NRHS ≥ 0 (глобальный входной параметр, тип целый);
A - указатель на локальную память, занимаемую массивом двойной точности с первой размерностью LLD_A ≥ (bwl + bwu + 1) (хранится в DESCA);
на входе - массив А содержит локальные части факторизуемой распределенной матрицы A (1 : N, JA : JA + N - 1). Эта локальная порция хранится в запакованном ленточном формате, используемом в пакете LAPACK. Более подробно см. в разделе документации "Схема размещения в локальной памяти и блочно - столбцовое разбиение ленточных матриц"
на выходе - этот массив содержит информацию о деталях разложения матрицы;
заметим, что перестановки выполняются над матрицей так, что выдаваемые множители отличаются от множителей, выдаваемых LAPACK 'ом (локальный входной и локальный выходной параметр);
JA - номер столбца в глобальном массиве A, указывающий на начало матрицы (которая может быть либо всей матрицей A, либо подматрицей матрицы A) (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_A = 501),  DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_A = 1),      DLEN ≥ 9;
DESCA содержит информацию о размещении A в памяти; полное описание DESCA см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
B - указатель на локальную память, занимаемую массивом двойной точности с локальной ведущей размерностью lld_b ≥ NB;
на входе -этот массив содержит локальные куски правой части
B ( IB : IB + N - 1, 1 : NRHS);
на выходе - массив содержит локальный кусок решений распределенной матрицы X, столбцы которой образуют NRHS векторов решений исходной системы;
IB - номер строки в глобальном массиве B, указывающий на первую строку матрицы (которая может быть либо всей матрицей B, либо подматрицей матрицы B) (глобальный входной параметр, тип целый);
DESCB - дескриптор распределенной матрицы B (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_B = 502),  DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_B = 1),      DLEN ≥ 9;
DESCB содержит информацию о размещении B в памяти; полное описание DESCB см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
WORK - временное рабочее пространство двойной точности; Это пространство может изменять свое содержимое между вызовами подпрограмм; WORK должно иметь длину, задаваемую параметром LWORK;
на выходе - WORK (1) содержит минимальное значение LWORK (локальное рабочее пространство, локальный выходной параметр);
LWORK - размер входного рабочего пространства пользователя с именем WORK; если LWORK слишком мало, минимальный приемлемый (допустимый) размер будет возвращен в элементе WORK (1) и будет выдан код ошибки;
LWORK ≥ NB * (bwl + bwu)  + 6 * max (bwl, bwu) * max (bwl, bwu) +
max ((max (bwl, bwu) * NRHS),  max (bwl, bwu) * max (bwl, bwu))
(локальный или глобальный входной параметр, тип целый);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если i - ый фактический параметр является массивом и его j - ый элемент имеет недопустимое значение,
то INFO = - ( i * 100 + j ),
если i - ый фактический параметр является скаляром и имеет недопустимое значение,
то INFO = - i;
> 0 - если INFO = K ≤ NPROCS, подматрица, хранящаяся на процессоре INFO и подвергающаяся разложению локально, не является подматрицей с диагональным преобладанием, и процесс разложения не был завершен;

если INFO = K > NPROCS, то подматрица, хранящаяся на процессоре
INFO-NPROCS, представляющая взаимодействие с другими процессорами, не была разложена, и процесс разложения не был завершен. 
   Версии
PSDBSV -   PCDBSV     PZDBSV     решение системы с ленточной матрицей общего вида методом Гаусса без выбора ведущих элементов для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно 
   Вызываемые подпрограммы

Здесь указаны только базовые подпрограммы (2 - ого уровня), которые вызываются из целевой подпрограммы (1 - ого уровня).

PDDBTRF - PZDBTRF    LU - разложение ленточной матрицы общего вида методом Гаусса без выбора ведущего элемента
PDDBTRS - PZDBTRS    Решение системы AX = B, ATX = B, или AH = B с ленточной матрицей общего вида на основе LU - разложения, полученного подпрограммой PDDBTRF(PZDBTRF) 
   Замечания по использованию
  В подпрограммах PSDBSV, PCDBSV, PZDBSV параметры A и B имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно
  Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), PXERBLA ( из пакета PBLAS)
  Для того чтобы избежать дублирования, большинство проверок ошибок в этой подпрограмме не выполняется, эти проверки выполняются в подпрограммах PDDBTRF(PZDBTRF) и PDDBTRS(PZDBTRS)
  Имеются следующие ограничения на входные параметры.
1. Не циклическое ограничение
 P * NB ≥  mod ( JA - 1, NB) + N
2. Ни один из процессоров не может содержать больше одного блока глобальной матрицы.
3. Размер блока (NB) не может быть слишком маленьким:
если матрица обрабатывается более чем на одном процессоре, то должно быть выполнено следующее неравенство:
NB >= 2*MAX(BWL, BWU)
Это объясняется тем, что на каждом процессоре действия выполняются над блоками матрицы размера O (NB). Если этот размер слишком мал, то алгоритм "разделяй и властвуй" является неэффективным.
4. Чтобы избежать лишних передач сообщений между процессорами, при работе с подматрицами следует полагать JA = IB.

Пример использования

Решается система уравнений A X = B , где A - квадратная ленточная матрица порядка 9

Пусть матрица А системы имеет вид:

 2 -1  0  0  0  0  0  0  0
-1  2 -1  0  0  0  0  0  0
 0 -1  2 -1  0  0  0  0  0
 0  0 -1  2 -1  0  0  0  0
 0  0  0 -1  2 -1  0  0  0
 0  0  0  0 -1  2 -1  0  0
 0  0  0  0  0 -1  2 -1  0
 0  0  0  0  0  0 -1  2 -1
 0  0  0  0  0  0  0 -1  1

Вектор правых частей В имеет вид:

                        |   0   |
                        |   0   |
                        |   0   |
                        |   0   |
                        |   0   |
                        |   0   |
                        |   0   |
                        |   0   |
                        |   1   |

Пусть матрица разбивается по столбцам на 3 блока с размером блоков NB = 3.
Осуществляется пропуск программы на 3 процессах, и используется одномерная решетка процессов 1 на 3.

Фрагмент фортранного текста вызывающей программы

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

Другой (более удобный) вариант теста к этой подпрограмме с использованием служебной подпрограммы PDELSET для распределения элементов исходной глобальной матрицы по параллельным процессам (как это описано в разделах методического пособия в пп.12, 15, 16) и с использованием служебной подпрограммы PDLAPRNT (для выборки всех элементов распределенной матрицы/вектора со всех процессов и последующей распечатки их единым массивом) можно получить в tddbsv1.zip.

      PROGRAM TPDDBSV
      include 'mpif.h'
      INTEGER     DLEN_, JA, IB, N, NB, RSRC, CSRC,
     $                    MXLLDA, MXLLDB, NRHS, NBRHS, NOUT,
     $                    MXLOCR, MXLOCC, MXRHSC
      PARAMETER   ( DLEN_ = 7, JA = 1, IB = 1, N = 9, NB = 3,
     $                          RSRC = 0, CSRC = 0, MXLLDA = 3, MXLLDB = 3,
     $                          NRHS = 1, NBRHS = 1, MXLOCR = 3, MXLOCC = 3,
     $                          MXRHSC = 1, NOUT = 6, BWL = 1, BWU = 1, LWORK = 15 )

      INTEGER   CTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
      DOUBLE PRECISION  ANORM, BNORM, EPS, RESID, XNORM

      INTEGER  DESCA( DLEN_ ), DESCB( DLEN_ )
      DOUBLE PRECISION  A( MXLLDA, MXLOCC ), B( MXLLDB, MXRHSC ),
     $                                       WORK( LWORK )

      EXTERNAL    BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
     $                       MATINIT, PDDBSV, SL_INIT

      DATA        NPROW / 1 / , NPCOL / 3 /

      CALL  SL_INIT( ICTXT, NPROW, NPCOL )
      CALL  BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )
      IF( MYROW.EQ.-1 ) GO TO 10

      DESCA(1) = 501
      DESCA(2) = ICTXT
      DESCA(3) = N
      DESCA(4) = NB
      DESCA(5) = CSRC
      DESCA(6) = 3
      DESCA(7) = 0

      DESCB(1) = 502
      DESCB(2) = ICTXT
      DESCB(3) = N
      DESCB(4) = NB
      DESCB(5) = RSRC
      DESCB(6) = 3
      DESCB(7) = 0

      CALL  MATINIT( A, DESCA, B, DESCB )

      CALL  PDDBSV( N, BWL, BWU, NRHS, A, JA, DESCA, B,
     $                             IB, DESCB, WORK, LWORK, INFO)

      CALL BLACS_GRIDEXIT( ICTXT )
 10 CONTINUE
      CALL BLACS_EXIT( 0 )
      STOP
      END

      SUBROUTINE  MATINIT( A, DESCA, B, DESCB )
*   Подпрограмма генерации и распределения по процессам матриц  A и B

      INTEGER   DESCA( * ), DESCB( * )
      DOUBLE PRECISION   A( * ), B( * )
      INTEGER   CTXT_
      PARAMETER   ( CTXT_ = 2 )
      INTEGER   ICTXT, MYCOL, MYROW, NPCOL, NPROW
      DOUBLE PRECISION   MO, T, Z
      EXTERNAL  BLACS_GRIDINFO

      ICTXT = DESCA( CTXT_ )
      CALL  BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )

      T = 2.0D0
      Z = 0.0D0
      MO = -1.0D0
*
* Локальные части матриц A и B для процесса (0, 0)
*
      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         A( 1 ) = 0.0D0
         A( 2 ) = T
         A( 3 ) = MO
         A( 4 ) = MO
         A( 5 ) = T
         A( 6 ) = MO
         A( 7 ) = MO
         A( 8 ) = T
         A( 9 ) = MO
*
         B( 1 ) = Z
         B( 2 ) = Z
         B( 3 ) = Z
*
* Локальные части матриц A и B для процесса (0, 1)
*
      ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.1 ) THEN
         A( 1 ) = MO
         A( 2 ) = T
         A( 3 ) = MO
         A( 4 ) = MO
         A( 5 ) = T
         A( 6 ) = MO
         A( 7 ) = MO
         A( 8 ) = T
         A( 9 ) = MO
*
         B( 1 ) = Z
         B( 2 ) = Z
         B( 3 ) = Z
*
* Локальные части матриц A и B для процесса (0, 2)
*
      ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.2 ) THEN
         A( 1 ) = MO
         A( 2 ) = T
         A( 3 ) = MO
         A( 4 ) = MO
         A( 5 ) = T
         A( 6 ) = MO
         A( 7 ) = MO
         A( 8 ) = 1.0D0
         A( 9 ) = Z
*
         B( 1 ) = Z
         B( 2 ) = Z
         B( 3 ) = 1.0D0
      END IF
      RETURN
      END

Результаты:

 Решение системы

  X = ( 1.D0, 2.D0, 3.D0,   4.D0, 5.D0, 6.D0,   7.D0, 8.D0, 9.D0 )

 Значение   INFO  =  0