Приложение 2_1

Пример описания целевой подпрограммы
для решения системы A*X = B с матрицей общего вида.

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

Назначение

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

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

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

         sub(A) * 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).

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

Литература:

http://www.srcc.msu.su/num_anal/
http://www.netlib.org/scalapack/slug/index.html

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

         CALL  PDGESV ( N, NRHS, A, IA, JA, DESCA, IPIV,
                                       B, IB, JB, DESCB, INFO)

   Параметры
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_); подробное описание см. в п.5.2 (глобальный и локальный входной параметр, тип целый);
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_); подробное описание см. в п.3 (глобальный и локальный входной параметр, тип целый);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если i - ый фактический параметр является массивом и его j - ый элемент имеет недопустимое значение,
то INFO = - ( i * 100 + j ),
если i - ый фактический параметр является скаляром и имеет недопустимое значение,
то INFO = - i;
> 0 - если INFO = K, то элемент U( IA+K-1, JA+K-1) является точным нулем; разложение завершено, но матрица U вырождена, и поэтому решение системы не вычисляется 
   Версии
PSGESV -  
PCGESV    
PZGESV     		решение системы с матрицей общего вида методом Гаусса с выбором ведущего элемента по столбцу для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно 
   Вызываемые подпрограммы

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

PDGETRF - PZGETRF    LU - разложение матрицы общего вида методом Гаусса с выбором ведущего элемента по столбцу
PDGETRS - PZGETRS    Решение системы AX = B, AT X = B или AH X = B на основе LU - разложения, полученного подпрограммой PDGETRF(PZGETRF) 
   Замечания по использованию
  1.  В подпрограммах PSGESV, PCGESV, PZGESV параметры A и B имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно
  2.  Используются подпрограммы 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

Фрагмент фортранного текста вызывающей программы
(принятые в тестах обозначения можно посмотреть в п.8.3) 

      PROGRAM TPDGESV
      include 'mpif.h'

      INTEGER     DLEN_, IA, JA, IB, JB, M, N, MB, NB, RSRC,
     $                     CSRC, MXLLDA, MXLLDB, NRHS, NBRHS, NOUT,
     $                     MXLOCR, MXLOCC, MXRHSC
      PARAMETER   ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
     $                           M = 9, N = 9, MB = 2, NB = 2, RSRC = 0,
     $                           CSRC = 0, MXLLDA = 5, MXLLDB = 5, NRHS = 1,
     $                           NBRHS = 1, NOUT = 6, MXLOCR = 5, MXLOCC = 4,
     $                           MXRHSC = 1 )
      DOUBLE PRECISION   ONE
      PARAMETER   ( ONE = 1.0D+0 )

      INTEGER   ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW

      INTEGER   DESCA( DLEN_ ), DESCB( DLEN_ ), IPIV( MXLOCR+NB )
      DOUBLE PRECISION   A( MXLLDA, MXLOCC ), B( MXLLDB, MXRHSC ),
     $                                       WORK( MXLOCR )

      EXTERNAL   BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
     $                       DESCINIT, MATINIT, PDGESV, SL_INIT

      INTRINSIC   DBLE
      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, M, N, MB, NB, RSRC, CSRC, ICTXT, MXLLDA,
     $                               INFO )
      CALL  DESCINIT( DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
     $                               MXLLDB, INFO )

      CALL  MATINIT( A, DESCA, B, DESCB )

      CALL  PDGESV( N, NRHS, A, IA, JA, DESCA, IPIV, B, IB, JB, DESCB,
     $                             INFO )

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

      SUBROUTINE  MATINIT( AA, DESCA, B, DESCB )
* MATINIT генерирует и распределяет матрицы  A и B  по решетке процессов 2 x 3

      INTEGER   DESCA( * ), DESCB( * )
      DOUBLE PRECISION   AA( * ), B( * )

      INTEGER   CTXT_, LLD_
      PARAMETER    ( CTXT_ = 2, LLD_ = 9 )

      INTEGER   ICTXT, MXLLDA, MYCOL, MYROW, NPCOL, NPROW
      DOUBLE PRECISION   A, C, K, L, P, S
      EXTERNAL   BLACS_GRIDINFO

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

      S = 19.0D0
      C = 3.0D0
      A = 1.0D0
      L = 12.0D0
      P = 16.0D0
      K = 11.0D0

      MXLLDA = DESCA( LLD_ )

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

      IF( MYROW.EQ.0 .AND. MYCOL.EQ.0 ) THEN
         AA( 1 ) = S
         AA( 2 ) = -S
         AA( 3 ) = -S
         AA( 4 ) = -S
         AA( 5 ) = -S
         AA( 1+MXLLDA ) = C
         AA( 2+MXLLDA ) = C
         AA( 3+MXLLDA ) = -C
         AA( 4+MXLLDA ) = -C
         AA( 5+MXLLDA ) = -C
         AA( 1+2*MXLLDA ) = A
         AA( 2+2*MXLLDA ) = A
         AA( 3+2*MXLLDA ) = A
         AA( 4+2*MXLLDA ) = A
         AA( 5+2*MXLLDA ) = -A
         AA( 1+3*MXLLDA ) = C
         AA( 2+3*MXLLDA ) = C
         AA( 3+3*MXLLDA ) = C
         AA( 4+3*MXLLDA ) = C
         AA( 5+3*MXLLDA ) = -C

         B( 1 ) = 0.0D0
         B( 2 ) = 0.0D0
         B( 3 ) = 0.0D0
         B( 4 ) = 0.0D0
         B( 5 ) = 0.0D0

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

      ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.1 ) THEN
         AA( 1 ) = A
         AA( 2 ) = A
         AA( 3 ) = -A
         AA( 4 ) = -A
         AA( 5 ) = -A
         AA( 1+MXLLDA ) = L
         AA( 2+MXLLDA ) = L
         AA( 3+MXLLDA ) = -L
         AA( 4+MXLLDA ) = -L
         AA( 5+MXLLDA ) = -L
         AA( 1+2*MXLLDA ) = K
         AA( 2+2*MXLLDA ) = K
         AA( 3+2*MXLLDA ) = K
         AA( 4+2*MXLLDA ) = K
         AA( 5+2*MXLLDA ) = K

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

      ELSE IF( MYROW.EQ.0 .AND. MYCOL.EQ.2 ) THEN
         AA( 1 ) = A
         AA( 2 ) = A
         AA( 3 ) = A
         AA( 4 ) = -A
         AA( 5 ) = -A
         AA( 1+MXLLDA ) = P
         AA( 2+MXLLDA ) = P
         AA( 3+MXLLDA ) = P
         AA( 4+MXLLDA ) = P
         AA( 5+MXLLDA ) = -P

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

      ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.0 ) THEN
         AA( 1 ) = -S
         AA( 2 ) = -S
         AA( 3 ) = -S
         AA( 4 ) = -S
         AA( 1+MXLLDA ) = -C
         AA( 2+MXLLDA ) = -C
         AA( 3+MXLLDA ) = -C
         AA( 4+MXLLDA ) = C
         AA( 1+2*MXLLDA ) = A
         AA( 2+2*MXLLDA ) = A
         AA( 3+2*MXLLDA ) = A
         AA( 4+2*MXLLDA ) = -A
         AA( 1+3*MXLLDA ) = C
         AA( 2+3*MXLLDA ) = C
         AA( 3+3*MXLLDA ) = C
         AA( 4+3*MXLLDA ) = C

         B( 1 ) = 1.0D0
         B( 2 ) = 0.0D0
         B( 3 ) = 0.0D0
         B( 4 ) = 0.0D0

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

      ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.1 ) THEN
         AA( 1 ) = A
         AA( 2 ) = -A
         AA( 3 ) = -A
         AA( 4 ) = -A
         AA( 1+MXLLDA ) = L
         AA( 2+MXLLDA ) = L
         AA( 3+MXLLDA ) = -L
         AA( 4+MXLLDA ) = -L
         AA( 1+2*MXLLDA ) = K
         AA( 2+2*MXLLDA ) = K
         AA( 3+2*MXLLDA ) = K
         AA( 4+2*MXLLDA ) = K

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

      ELSE IF( MYROW.EQ.1 .AND. MYCOL.EQ.2 ) THEN
         AA( 1 ) = A
         AA( 2 ) = A
         AA( 3 ) = -A
         AA( 4 ) = -A
         AA( 1+MXLLDA ) = P
         AA( 2+MXLLDA ) = P
         AA( 3+MXLLDA ) = -P
         AA( 4+MXLLDA ) = -P
      END IF
      RETURN
      END

Результаты:

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

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

 Значение   INFO  =  0