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

Назначение

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

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

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

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

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

Разложение Холецкого используется для факторизации матрицы в виде множителей 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  PDPBSV ( UPLO, N, BW, NRHS, A, JA, DESCA,
                                        B, IB, DESCB, WORK, LWORK, INFO) 

   Параметры
UPLO - если UPLO = ' U ', сохраняется верхний треугольник матрицы
A (1 : N, JA : JA + N - 1);
если UPLO = ' L ', сохраняется нижний треугольник матрицы
A (1 : N, JA : JA + N - 1)
(глобальный входной параметр, тип символьный);
N - порядок распределенной подматрицы A (1 : N, JA : JA + N - 1), N і 0 (глобальный входной параметр, тип целый);
BW - число поддиагоналей в L или U, 0 Ј BW Ј N - 1 (глобальный входной параметр, тип целый);
NRHS - число правых частей, т.е. число столбцов распределенной подматрицы
B ( IB : IB + N - 1, 1 : NRHS), NRHS і 0 (глобальный входной параметр, тип целый);
A - указатель на локальную память, занимаемую массивом двойной точности с первым размером LLD_A і (bw + 1), хранящимся в DESCA;
на входе - этот массив содержит локальные куски факторизуемой распределенной матрицы A (1 : N, JA : JA + N - 1). Эта локальная порция хранится в запакованном ленточном формате, используемом в пакете LAPACK. Более подробно см. в разделе документации "Схема размещения в локальной памяти и блочно - столбцовое разбиение ленточных матриц";
на выходе - массив A содержит информацию, касающуюся подробностей факторизации; заметим, что перестановки выполняются над матрицей A так, что возвращаемые множители отличаются от тех, которые возвращаются в пакете LAPACK;
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 і ( NB + 2 * bw) * bw + max ((bw * NRHS), bw * bw)
(локальный или глобальный входной параметр, тип целый);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если i - ый фактический параметр является массивом и его j - ый элемент имеет недопустимое значение,
то INFO = - ( i * 100 + j ),
если i - ый фактический параметр является скаляром и имеет недопустимое значение,
то INFO = - i;
> 0 - если INFO = K Ј NPROCS, то подматрица, хранящаяся на процессоре INFO и локально факторизуемая, не является положительно определенной, и факторизация не может быть завершена;
если INFO = K > NPROCS, то подматрица, хранящаяся на процессоре
INFO - NPROCS, и представляющая взаимодействие с другими процессорами, не является положительно определенной, и факторизация не может быть завершена 
   Версии
PSPBSV -   PCPBSV     PZPBSV     решение системы с симметричной (эрмитовой) положительно определенной ленточной матрицей методом Холецкого для случаев вещественной одинарной точности, комплексной одинарной точности и комплексной двойной точности матриц соответственно 
   Вызываемые подпрограммы

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

PDPBTRF - PZPBTRF    Треугольное разложение симметричной (эрмитовой) положительно определенной ленточной матрицы методом Холецкого
PDPBTRS - PZPBTRS    Решение системы линейных уравнений с симметричной (эрмитовой) положительно определенной ленточной матрицей, на основе разложения Холецкого, полученного подпрограммой PDPBTRF(PZPBTRF) 
   Замечания по использованию
  В подпрограммах PSPBSV, PCPBSV и PZPBSV параметры A, B и WORK имеют REAL, COMPLEX и DOUBLE COMPLEX соответственно
  Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), PXERBLA ( из пакета PBLAS)
  Для того чтобы избежать дублирования, наибольшая часть проверок ошибок в этой подпрограмме не выполняется, эти проверки выполняются в подпрограммах PDPBTRF(PZPBTRF) и PDPBTRS(PZPBTRS)
  Имеются следующие ограничения на входные параметры.
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.

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

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

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

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

      CHARACTER   UPLO
      PARAMETER   ( UPLO = 'L' )

      INTEGER   ICTXT, 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, PDPBSV, 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) = 2
      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  PDPBSV( UPLO, N, BW, 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, MXLLDA, 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 ) = T
         A( 2 ) = MO
         A( 3 ) = T
         A( 4 ) = MO
         A( 5 ) = T
         A( 6 ) = 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 ) = T
         A( 2 ) = MO
         A( 3 ) = T
         A( 4 ) = MO
         A( 5 ) = T
         A( 6 ) = 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 ) = T
         A( 2 ) = MO
         A( 3 ) = T
         A( 4 ) = MO
         A( 5 ) = 1.0D0
         A( 6 ) = 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