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

Назначение

Решение системы  A * X = B или  AT * X = B с матрицей общего вида методом Гаусса с выбором ведущего элемента по столбцу и оценка обратного числа обусловленности

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

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

(1)         sub(A) * X   =  sub(B)     или

(2)         sub(A)T * X  =  sub(B) ,

  где  sub(A) = A ( IA : IA+N-1, JA : JA+N-1) - квадратная распределенная матрица
                                                                           порядка N,
                  X -  распределенная матрица размеров N на NRHS,
         sub(B) = B ( IB : IB+N-1, JB : JB+NRHS-1) -  распределенная матрица
                                                                                 размеров N на NRHS. 

Используется LU - разложение методом Гаусса с выбором ведущего элемента по столбцу для факторизации матрицы sub(A) в виде

         sub(A) = P * L * U ,

  где  P - матрица перестановок,
         L - нижняя треугональная матрица с единичными диагональными элементами,
         U - верхняя треугональная матрица. 

Матрицы L и U помещаются в памяти на место матрицы sub (A). Полученное LU - разложение используется для решения системы  sub (A) * X = sub (B)  или  sub (A) T * X = sub (B).

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

Литература:

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

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

         CALL  PDGESV2 ( TRANS, N, NRHS, A, IA, JA, DESCA, IPIV,
                                          B, IB, JB, DESCB, RCOND, WORK, LWORK, IWORK, LIWORK, INFO)

Параметры

TRANS - если TRANS = ' N ', решается система (1),
если TRANS = ' T ', решается система (2),
(глобальный входной параметр, тип символьный);
N - порядок распределенной подматрицы sub (A), N ≥ 0 (глобальный входной параметр, тип целый);
NRHS - число правых частей, т.е. число столбцов распределенной подматрицы sub (B), NRHS ≥ 0 (глобальный входной параметр, тип целый);
A - указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_A, LOCc ( JA + N - 1));
на входе - это локальные части факторизуемой распределенной матрицы sub(A) порядка N;
на выходе - массив A содержит локальные части множителей L и U, полученные при разложении (факторизации) подматрицы в виде sub (A) = P * L * U; единичные диагональные элементы матрицы L не хранятся в памяти (локальный входной и локальный выходной параметр);
IA - номер строки в глобальном массиве A, указывающий на первую строку подматрицы sub (A) (глобальный входной параметр, тип целый);
JA - номер столбца в глобальном массиве A, указывающий на первый столбец подматрицы sub (A) (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
IPIV - одномерный массив длины ( LOCr (M_A) + MB_A), содержащий информацию о выборе ведущего элемента на каждом шаге исключения Гаусса;
элемент массива IPIV ( i ) содержит глобальный номер строки, которая переставляется со строкой, имеющей локальный номер  i;
этот массив содержит информацию о перестановках строк распределенной матрицы A (локальный выходной параметр, тип целый);
B - указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_B, LOCc ( JB + NRHS - 1));
на входе - это распределенная матрица sub (B) правой части;
на выходе (если INFO = 0) на месте sub (B) находится распределенная матрица X, столбцы которой образуют NRHS векторов решений исходной системы (локальный входной и локальный выходной параметр);
IB - номер строки в глобальном массиве B, указывающий на первую строку подматрицы sub (B) (глобальный входной параметр, тип целый);
JB - номер столбца в глобальном массиве B, указывающий на первый столбец подматрицы sub (B) (глобальный входной параметр, тип целый);
DESCB - дескриптор распределенной матрицы B (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
RCOND - оценка обратного числа обусловленности матрицы
A ( IA : IA + N - 1, JA : JA + N - 1);
если RCOND меньше машинной точности (машинного epsilon) (в частности, если RCOND = 0), матрица A является вырожденной с машинной точностью; в этом случае код возврата  INFO > 0 (глобальный выходной параметр, тип DOUBLE PRECISION);
WORK - одномерный рабочий массив двойной точности длины LWORK; на выходе - WORK (1) содержит минимальную требуемую величину значение LWORK (локальное рабочее пространство, локальный выходной параметр);
LWORK - размер рабочего пространства (массива) WORK;  LWORK - локальный входной параметр и должен быть не меньше, чем
LWORK = PDGECON ( LWORK ) + LOCr ( N_A);
где PDGECON ( LWORK ) означает величину рабочего массива WORK, которая требуется при работе вызываемой подпрограммы PDGECON (см.ниже);
смысл обозначений LOCr (...) и LOCc (...) можно найти в разделах документации "Схема размещения в локальной памяти и блочно - циклическое отображение плотных матриц" и "Пример блочно - циклического распределения плотной матрицы по решетке процессов";
минимальная требуемая величина LWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LWORK = - 1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LWORK возвращается в элементе массива WORK (1); при этом сообщений об ошибках не выдается (локальный или глобальный входной параметр);
IWORK - одномерный рабочий массив целого типа длины LIWORK; на выходе IWORK (1) содержит минимальную требуемую величину LIWORK (локальное рабочее пространство, локальный выходной параметр);
LIWORK - размер рабочего пространства (массива) IWORK; LIWORK - локальный входной параметр и должен быть не меньше, чем LIWORK = LOCr ( N_A);
см. раздел документации "Схема размещения в локальной памяти и блочно - циклическое отображение плотных матриц";
минимальная требуемая величина LIWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LIWORK = - 1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LIWORK возвращается в элементе массива IWORK (1); при этом сообщений об ошибках не выдается (локальный или глобальный входной параметр, тип INTEGER);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если i - ый фактический параметр является массивом и его j - ый элемент имеет недопустимое значение,
то INFO = - ( i * 100 + j ),
если i - ый фактический параметр является скаляром и имеет недопустимое значение,
то INFO = - i;
> 0 - если INFO = K (и K ≤ N ), то элемент U(IA+K-1, JA+K-1) является точным нулем; разложение завершено, но матрица U вырождена, и поэтому решение системы не вычисляется;
если INFO = N+1, RCOND меньше машинного eps; разложение завершено, но матрица sub(A) вырождена с машинной точностью, и поэтому решение системы не вычисляется

Версии

PSGESV2 -   PCGESV2     PZGESV2     решение системы с матрицей общего вида методом Гаусса с выбором ведущего элемента по столбцу для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно

Вызываемые подпрограммы

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

PDGETRF - PZGETRF    LU - разложение матрицы общего вида методом Гаусса с выбором ведущего элемента по столбцу
PDGETRS - PZGETRS    Решение системы A X = B, AT X = B или AH X = B на основе LU - разложения, полученного подпрограммой PDGETRF(PZGETRF)
PDGECON - PZGECON    Оценка обратного числа обусловленности матрицы общего вида

Замечания по использованию

  1.  В подпрограммах PSGESV2, PCGESV2, PZGESV2 параметры A и B имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно
  2.  В версии подпрограммы PZGESV2 вместо параметров IWORK и LIWORK целого типа используются параметры с именами RWORK и LRWORK. При этом параметр RWORK имеет тип DOUBLE PRECISION, а параметр LRWORK ≥ 2 * LOCc (N_A).
  3.  Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), PXERBA ( из пакета PBLAS), CHK1MAT, PCHK2MAT, INDXG2P ( из библиотеки ScaLAPACK_TOOLS)

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

Решение системы A X = B , где A - матрица порядка 9
Пусть матрица А системы имеет вид:

 19  3  1  12  1  16  1  3 11
-19  3  1  12  1  16  1  3 11
-19 -3  1  12  1  16  1  3 11
-19 -3 -1  12  1  16  1  3 11
-19 -3 -1 -12  1  16  1  3 11
-19 -3 -1 -12 -1  16  1  3 11
-19 -3 -1 -12 -1 -16  1  3 11
-19  3 -1 -12 -1 -16 -1  3 11
-19 -3 -1 -12 -1 -16 -1 -3 11

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

 0
 0
 1
 0
 0
 0
 0
 0
 0

Матрица разбивается на блоки размеров 2 на 2.
Осуществляется пропуск программы на 6 процессах и используется двумерная решетка процессов 2 на 3.
Ниже показано, как разбитая на блоки матрица A распределяется (блочно-циклически) по решетке процессов 2 * 3.

0 1 2
0  19   3
-19   3
 1   3
 1   3
 1   12
 1   12
11
11
 1   16
 1   16
-19  -3
-19  -3
 1   3
 1   3
-1  -12
-1  -12
11
11
 1   16
-1   16
-19  -3 -1  -3 -1  -12 11 -1  -16
1 -19  -3
-19  -3
 1   3
 1   3
 1   12
-1   12
11
11
 1   16
 1   16
-19  -3
-19   3
 1   3
-1   3
-1  -12
-1  -12
11
11
-1  -16
-1  -16

Распределение матрицы B (блочно-циклическое) по процессам

0 1 2
0 0
0
0
0
0
1 1
0
0
0

Фрагмент фортранного текста вызывающей программы
Полный текст теста можно получить в tdgesv2.zip, а принятые обозначения - посмотреть в разделе документации "Обозначения и упрощения в примерах по использованию подпрограмм комплекса"; в тесте используются: служебная подпрограмма PDELSET (для распределения элементов глобальной матрицы/вектора по паралельным процессам) и служебная подпрограмма PDLAPRNT (для выборки и распечатки элементов распределенной матрицы/вектора единым массивом)

      PROGRAM TDGESV2
      include 'mpif.h'

      INTEGER     DLEN_, IA, JA, IB, JB, N, NB, RSRC,
     $                     CSRC, MXLLDA, MXLLDB, NRHS, NBRHS, NOUT,
     $                     MXLOCR, MXLOCC, MXRHSC, LWORK, LIWORK
      PARAMETER   ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
     $                           N = 9, NB = 2, RSRC = 0,
     $                           CSRC = 0, MXLLDA = 5, MXLLDB = 5, NRHS = 1,
     $                           NBRHS = 1, NOUT = 6, MXLOCR = 5, MXLOCC = 4,
     $                           MXRHSC = 1, LWORK = 24, LIWORK = 5 )

      CHARACTER   TRANS
      PARAMETER   ( TRANS = 'N' )

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

      DOUBLE PRECISION   RCOND

      INTEGER   DESCA( DLEN_ ), DESCB( DLEN_ ), IPIV( MXLOCR+NB ), IWORK( LIWORK )

      DOUBLE PRECISION   A( MXLLDA, MXLOCC ), B( MXLLDB, MXRHSC ), WORK( LWORK )

      EXTERNAL   BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
     $                       DESCINIT, PDMATINIT, PDGESV2, SL_INIT

      DATA   NPROW / 2 / , NPCOL / 3 /

      CALL  SL_INIT( ICTXT, NPROW, NPCOL )
      CALL  BLACS_GRIDINFO( ICTXT, NPROW, NPCOL, MYROW, MYCOL )

      IF( MYROW.EQ.-1 ) GO TO 10

      CALL  DESCINIT( DESCA, N, N, NB, NB, RSRC, CSRC, ICTXT, MXLLDA,
     $                               INFO )
      CALL  DESCINIT( DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
     $                               MXLLDB, INFO )

      CALL  PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )

      CALL  PDGESV2( TRANS, N, NRHS, A, IA, JA, DESCA, IPIV,
     $                               B, IB, JB, DESCB, RCOND, WORK, LWORK, IWORK, LIWORK, INFO )

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

      SUBROUTINE PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
* PDMATINIT генерирует и распределяет матрицы  A и B
*           по решетке процессов 2 x 3
* Параметры
      INTEGER            IA, IB, INFO, JA, JB, N
      INTEGER            DESCA( * ), DESCB( * )
      DOUBLE PRECISION   A( * ), B( * )
*
* Локальные переменные
      INTEGER            I, J
      DOUBLE PRECISION   AIJ, BIJ
*
* Внешние подпрограммы
      EXTERNAL           PDELSET
*
* Выполняемые операторы
      INFO = 0
*
      IF( IA.NE.1 ) THEN
         INFO = -3
      ELSE IF(JA .NE.1) THEN
         INFO = -4
      END IF
*
      AIJ = 19.0D0
*
      DO 20 J = 1, N
         DO 10 I = 1, N 
            IF( I.LE.J ) THEN
                CALL PDELSET( A, I, J, DESCA, AIJ )
            ELSE
                CALL PDELSET( A, I, J ,DESCA, -AIJ )
            END IF
   10 CONTINUE
         IF( J.EQ.1 .OR. J.EQ.7 ) AIJ = 3.0D0
         IF( J.EQ.2 .OR. J.EQ.4 .OR. J.EQ.6 ) AIJ = 1.0D0
         IF( J.EQ.3 ) AIJ = 12.0D0
         IF( J.EQ.5 ) AIJ = 16.0D0
         IF( J.EQ.8 ) AIJ = 11.0D0
   20 CONTINUE

      CALL PDELSET( A, 8, 2, DESCA, 3.0D0 )
*
      BIJ = 0.0D0
*
      J = 1
      DO 30 I = 1, N
        CALL PDELSET( B, I, J, DESCB, BIJ )
   30 CONTINUE

      CALL PDELSET( B, 3, 1, DESCB, 1.0D0 )
*
      RETURN
      END

Результаты:

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

 X = ( 0.D0, - 0.16666666666667D0, 0.5D0, 0.D0, 0.D0, 0.D0, - 0.5D0, 0.16666666666667D0, 0.D0 )

 Значение   INFO  =  0 ,      RCOND = 0.68027211 E-2