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

Назначение

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

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

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

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

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

Используется метод исключения Гаусса без выбора ведущего элемента для разложения матрицы на множители 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  PDDTSV ( N, NRHS, DL, D, DU, JA, DESCA,
                                        B, IB, DESCB, WORK, LWORK, INFO)

   Параметры
N - порядок распределенной подматрицы A (1 : N, JA : JA + N - 1), N ≥ 0 (глобальный входной параметр, тип целый)
NRHS - число правых частей, т.е. число столбцов распределенной
подматрицы B ( IB : IB + N - 1, 1 : NRHS), NRHS ≥ 0 (глобальный входной параметр, тип целый)
DL - указатель на локальную часть глобального вектора двойной точности, в котором хранится нижняя диагональ матрицы A;
глобально, к DL (1) ссылок нет, и DL должно быть выровнено с D; необходимо, чтобы размер (size) был ≥ DESCA (NB_); на выходе, этот массив содержит информацию, касающуюся множителей, на которые разложена матрица A (локальный входной и локальный выходной параметр)
D - указатель на локальную часть глобального вектора двойной точности, в котором хранится главная диагональ матрицы A;
на выходе, этот массив содержит информацию, касающуюся множителей, на которые разложена матрица A; необходимо, чтобы размер (size) был ≥ DESCA (NB_) (локальный входной и локальный выходной параметр)
DU - указатель на локальную часть глобального вектора двойной точности, в котором хранится верхняя диагональ матрицы A;
глобально, к DU (n) ссылок нет, и DU должен быть выровнен с D; на выходе, этот массив содержит информацию, касающуюся множителей, на которые разложена матрица A; необходимо, чтобы размер был ≥ DESCA (NB_) (локальный входной и локальный выходной параметр)
JA - номер столбца в глобальном массиве A, указывающий на на начало матрицы, которая должна быть обработана (которая может быть либо всей матрицей A, либо подматрицей матрицы A) (глобальный входной параметр, тип целый)
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_A = 501 или 502),  DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_A = 1),                     DLEN ≥ 9.
Содержит информацию о размещении 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.
Содержит информацию о размещении B в памяти. Полное описание DESCB см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый)
WORK - временное рабочее пространство двойной точности. Это пространство может изменять свое содержимое между вызовами подпрограмм. WORK должно иметь длину, задаваемую параметром LWORK; на выходе - WORK (1) содержит минимальное значение LWORK (локальное рабочее пространство, локальный выходной параметр)
LWORK - размер входного рабочего пространства пользователя с именем WORK. Если LWORK слишком мало, минимальный допустимый размер будет возвращен в элементе WORK (1) и будет выдан код ошибки.
LWORK ≥ (12 * NPCOL + 3 * NB) +  max (10 * NPCOL + 4 * NRHS, 8 * NPCOL)
(локальный входной или глобальный выходной параметр, тип целый)
INFO - целая переменная, диагностирующая результат работы подпрограммы (локальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если i - ый фактический параметр является массивом и его j - ый элемент имеет недопустимое значение,
то INFO = - ( i * 100 + j ),
если i - ый фактический параметр является скаляром и имеет недопустимое значение,
то INFO = - i;
> 0 - если INFO = K ≤ NPROCS, подматрица, хранящаяся на процессоре INFO и подвергающаяся разложению локально, не является подматрицей с диагональным преобладанием, и процесс разложения не был завершен;
если INFO = K > NPROCS, то подматрица, хранящаяся на процессоре
INFO-NPROCS, представляющая взаимодействие с другими процессорами, не была разложена, и процесс разложения не был завершен. 
   Версии
PSDTSV -   PCDTSV     PZDTSV     решение системы с трехдиагональной матрицей общего вида методом Гаусса без выбора ведущего элемента для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно 
   Вызываемые подпрограммы

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

PDDTTRF - PZDTTRF    LU - разложение трехдиагональной матрицы общего вида методом Гаусса без выбора ведущего элемента
PDDTTRS - PZDTTRS    Решение системы с трехдиагональной матрицей общего вида на основе LU - разложения, полученного подпрограммой PDDTTRF(PZDTTRF) 
   Замечания по использованию
  В подпрограммах PSDTSV, PCDTSV, PZDTSV параметры B, D, DL, DU, WORK имеют тип REAL, COMPLEX и DOUBLE COMPLEX, соответственно
  Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), PXERBLA ( из пакета PBLAS)
  Для того чтобы избежать дублирования, большинство проверок ошибок в этой подпрограмме не выполняется, эти проверки выполняются в подпрограммах PDDTTRF(PZDTTRF) и PDDTTRS(PZDTTRS)
  Имеются следующие ограничения на входные параметры.
1. Не циклическое ограничение
P * NB ≥  mod ( JA - 1, NB) + N
2. Ни один из процессоров не может содержать больше одного блока глобальной матрицы.
3. Размер блока (NB) не может быть слишком маленьким:
если матрица обрабатывается более чем на одном процессоре, то должно быть выполнено следующее неравенство:
NB ≥ 2*BW
Это объясняется тем, что на каждом процессоре действия выполняются над блоками матрицы размера 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.

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

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

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

      PROGRAM TPDDTSV
      include 'mpif.h'
      INTEGER     DLEN_, JA, IB, N, NB, RSRC, CSRC, MXLLDB,
     $                    NRHS, NBRHS, NOUT, MXRHSC, LWORK
      PARAMETER   ( DLEN_ = 7, JA = 1, IB = 1, N = 9, NB = 3,
     $                           RSRC = 0, CSRC = 0, MXLLDB = 3,
     $                           NRHS = 1, NBRHS = 1,
     $                           MXRHSC = 1, NOUT = 6, LWORK = 79)

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

      INTEGER   DESCA( DLEN_ ), DESCB( DLEN_ )
      DOUBLE PRECISION   B( MXLLDB, MXRHSC ), WORK( LWORK ),
     $                                       DL( NB ), D( NB ), DU( NB )

      EXTERNAL   BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
     $                       MATINIT, PDDTSV, 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) = 0
      DESCA(7) = 0

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

      CALL  MATINIT( DL, D, DU, DESCA, B, DESCB )

      CALL  PDDTSV( N, NRHS, DL, D, DU, JA, DESCA,
     $                              B, IB, DESCB, WORK, LWORK, INFO)

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

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

      INTEGER   DESCA( * ), DESCB( * )
      DOUBLE PRECISION   DL( * ), D( * ), DU( * ), B( * )
      INTEGER   CTXT_
      PARAMETER   ( CTXT_ = 2 )
      INTEGER   ICTXT, MYCOL, MYROW, NPCOL, NPROW

      DOUBLE PRECISION   T, MO, Z

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

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

* Локальные части матриц A и B для процесса (0, 0)

      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         DU(1) = Z
         DU(2) = MO
         DU(3) = MO
         D(1) = T
         D(2) = T
         D(3) = T
         DL(1) = MO
         DL(2) = MO
         DL(3) = MO

         B(1) = Z
         B(2) = Z
         B(3) = Z

* Локальные части матриц A и B для процесса (0, 1)

      ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.1 ) THEN
         DU(1) = MO
         DU(2) = MO
         DU(3) = MO
         D(1) = T
         D(2) = T
         D(3) = T
         DL(1) = MO
         DL(2) = MO
         DL(3) = MO

         B(1) = Z
         B(2) = Z
         B(3) = Z

* Локальные части матриц A и B для процесса (0, 2)

      ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.2 ) THEN
         DU(1) = MO
         DU(2) = MO
         DU(3) = MO
         D(1) = T
         D(2) = T
         D(3) = 1.D0
         DL(1) = MO
         DL(2) = MO
         DL(3) = 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