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

Назначение

Решение системы с симметричной положительно определенной трехдиагональной матрицей методом Холецкого

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

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

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

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

Используется разложение Холецкого для факторизации матрицы в виде множителей L и L '.

Матрицы 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  PDPTSV ( N, NRHS, D, E, 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 (глобальный входной параметр, тип целый);
D - указатель на локальную часть глобального вектора двойной точности, в котором хранится главная диагональ матрицы A;
на выходе, этот массив содержит информацию, касающуюся множителей, на которые разложена матрица A; необходимо, чтобы размер (size) был ≥ DESCA (NB_) (локальный входной и локальный выходной параметр);
E - указатель на локальную часть глобального вектора двойной точности, в котором хранится верхняя диагональ матрицы A; глобально, на DU (n) нет ссылок, и DU должно быть выровнено с D;
на выходе, этот массив содержит информацию о множителях, на которые разложена матрица A; необходимо, чтобы размер (size) был ≥ DESCA (NB_) (локальный входной и локальный выходной параметр);
JA - номер столбца в глобальном массиве A, указывающий на начало матрицы (которая может быть либо всей матрицей A, либо подматрицей матрицы A) (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN); содержит информацию о размещении A в памяти; полное описание DESCA см. в разделе документации "Дескрипторы глобальных массивов";
если тип дескриптора одномерный ( DTYPE_A = 501),  DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_A = 1),      DLEN ≥ 9;
(глобальный и локальный входной параметр, тип целый);
B - указатель на локальную память, занимаемую массивом двойной точности с локальной ведущей размерностью lld_b ≥ NB;
на входе - этот массив содержит локальные куски правых частей
B ( IB : IB + N - 1, 1 : NRHS);
на выходе - он содержит локальный кусок распределенной матрицы решений X;
IB - номер строки в глобальном массиве B, указывающий на первую строку матрицы (которая может быть либо всей матрицей B, либо подматрицей матрицы B) (глобальный входной параметр, тип целый);
DESCB - дескриптор распределенной матрицы B (одномерный массив длины DLEN); содержащий информацию о размещении матрицы B в памяти; полное описание DESCB см. в разделе документации "Дескрипторы глобальных массивов";
если тип дескриптора одномерный ( DTYPE_B = 502),  DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_B = 1),      DLEN ≥ 9;
(глобальный и локальный входной параметр, тип целый);
WORK - временное рабочее пространство двойной точности; это пространство может изменять свое содержимое между вызовами подпрограмм; WORK должно иметь длину, задаваемую параметром LWORK;
на выходе - WORK (1) содержит минимальное значение LWORK (локальное рабочее пространство, локальный рабочий параметр);
LWORK - размер отводимого пользователем рабочего пространства с именем WORK; если LWORK слишком мало, минимальный приемлемый (допустимый) размер будет возвращен в элементе WORK (1) и будет выдан код ошибки;
LWORK ≥ (12 * NPCOL + 3 * NB) + max ((10 + 2 * min (100, NRHS)) *
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, и представляющая взаимодействие с другими процессорами, не является положительно определенной, и факторизация не может быть завершена. 
   Версии
PSPTSV -   PCPTSV     PZPTSV     решение системы с симметричной положительно определенной трехдиагональной матрицей методом Холецкого для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно 
   Вызываемые подпрограммы

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

PDPTTRF - Треугольное разложение симметричной положительно определенной трехдиагональной матрицы методом Холецкого
PDPTTRS - Решение системы с симметричной положительно определенной трехдиагональной матрицей на основе разложения Холецкого, полученного подпрограммой PDPTTRF 
   Замечания по использованию
  В подпрограммах PSPTSV, PCPTSV и PZPTSV параметры B, D, E и WORK имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно
  Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), PXERBLA ( из пакета PBLAS)
  Для того чтобы избежать дублирования, большинство проверок ошибок в этой подпрограмме не выполняется, эти проверки выполняются в подпрограммах PDPTTRF и PDPTTRS
 
Список параметров подпрограммы PZPTSV имеет следующий вид:

PZPTSV ( UPLO, N, NRHS, D, E, JA, DESCA, B, IB, DESCB, WORK, LWORK, INFO),
где дополнительный параметр UPLO означает: 
   если UPLO = 'U' - задается верхняя диагональ матрицы  A;
   если UPLO = 'L' - задается нижняя диагональ матрицы  A
   (глобальный входной параметр, тип символьный);
  Имеются следующие ограничения на входные параметры.
1. Не циклическое ограничение
P * NB ≥ mod ( JA - 1, NB) + N
2. Ни один из процессоров не может содержать больше одного блока глобальной матрицы.
3. Размер блока (NB) не может быть слишком маленьким:
если матрица обрабатывается более чем на одном процессоре, то должно быть выполнено следующее неравенство:
NB ≥ 2
Это объясняется тем, что на каждом процессоре действия выполняются над блоками матрицы размера 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.

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

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

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

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

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

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

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

      INTRINSIC   DBLE
      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) = 3
      DESCB(7) = 0

      CALL  MATINIT( D, E, DESCA, B, DESCB )

      CALL  PDPTSV( N, NRHS, D, E, JA, DESCA, B, IB,
     $                             DESCB, WORK, LWORK, INFO)

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

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

      INTEGER   DESCA( * ), DESCB( * )
      DOUBLE PRECISION   D( * ), E( * ), 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
         D(1) = T
         D(2) = T
         D(3) = T
         E(1) = MO
         E(2) = MO
         E(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
         D(1) = T
         D(2) = T
         D(3) = T
         E(1) = MO
         E(2) = MO
         E(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
         D(1) = T
         D(2) = T
         D(3) = 1.D0
         E(1) = MO
         E(2) = MO
         E(3) = 0.D0

         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