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

Назначение

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

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

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

         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.

Используется разложение Холецкого для представления матрицы sub (A) в виде 

         sub(A) = UT * U    или    sub(A) = L * LT ,

  где  U - верхняя треугольная матрица,
          L - нижняя треугольная матрица.

Полученное разложение используется затем для решения указанной выше системы уравнений.

После этого производится итерационное уточнение полученного решения и делается оценка границ ошибок.

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

Литература:

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

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

         CALL  PDPOSV2 ( UPLO, N, NRHS, A, IA, JA, DESCA,
                                          AF, IAF, JAF, DESCAF, B, IB, JB, DESCB,
                                          X, IX, JX, DESCX, FERR, BERR, WORK, LWORK,
                                          IWORK, LIWORK, INFO)

   Параметры
UPLO - 
 если UPLO = ' U ', сохраняется верхняя треугольная
                                подматрица sub(A),
 если UPLO = ' L ', сохраняется нижняя треугольная
                                подматрица sub(A)
(глобальный входной параметр, тип символьный);
N - порядок распределенной подматрицы sub (A), N і 0 (глобальный входной параметр, тип целый);
NRHS - число правых частей, т.е. столбцов распределенной подматрицы sub (B), NRHS і 0 (глобальный входной параметр, тип целый);
A - указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_A, LOCc ( JA + N - 1));
на входе - этот массив содержит локальные части исходной квадратной ( N на N ) симметричной распределенной матрицы sub (A);
если UPLO = ' U ', главная (ведущая) верхняя треугольная часть матрицы sub (A) размера N на N содержит верхнюю треугольную часть распределенной матрицы, а ее строго нижняя треугольная часть (под главной диагональю) не используется;
если UPLO = ' L ', ведущая нижняя треугольная часть матрицы sub (A) размера N на N содержит нижнюю треугольную часть распределенной матрицы, а ее строго верхняя треугольная часть (над главной диагональю) не используется;
на выходе A не изменяется (локальный входной параметр);
IA - номер строки в глобальном массиве A, указывающий на первую строку подматрицы sub (A) (глобальный входной параметр, тип целый);
JA - номер столбца в глобальном массиве A, указывающий первый столбец подматрицы sub (A) (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN_);
DESCA содержит информацию о размещении A в памяти; полное описание DESCA см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр);
AF - указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_AF, LOCc ( JA + N - 1));
на выходе, если INFO = 0, массив AF содержит локальные части распределенных множителей U или L, полученных при разложении Холецкого sub (A) = UT * U или sub (A) = L * LT (локальный выходной параметр);
IAF - номер строки в глобальном массиве AF, указывающий на первую строку подматрицы sub (AF) (глобальный входной параметр, тип целый);
JAF - номер столбца в глобальном массиве AF, указывающий на первый столбец подматрицы sub (AF) (глобальный входной параметр, тип целый);
DESCAF - дескриптор распределенной матрицы AF (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
B - указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_B, LOCc ( JB + NRHS - 1));
на входе - это локальные куски правой части распределенной матрицы sub (B);
на выходе B не изменяется (локальный входной параметр);
IB - номер строки в глобальном массиве B, указывающий на первую строку подматрицы sub (B) (глобальный входной параметр);
JB - номер столбца в глобальном массиве B, указывающий на первый столбец подматрицы sub (B) (глобальный входной параметр);
DESCB - дескриптор распределенной матрицы B (одномерный массив длины DLEN_);
DESCB содержит информацию о размещении B в памяти; полное описание DESCB см. в разделе документации "Дескрипторы глобальных массивов" (глобальный входной параметр);
X - указатель на локальную память, занимаемую массивом двойной точности размерности ( LLD_X, LOCc ( JX + NRHS - 1));
на выходе (если INFO = 0) это - распределенная матрица X, столбцы которой образуют NRHS векторов решений исходной системы, полученных в результате итерационного уточнения (локальный выходной параметр);
IX - номер строки в глобальном массиве X, указывающий на первую строку подматрицы sub ( X ) (глобальный входной параметр, тип целый);
JX - номер столбца в глобальном массиве X, указывающий на первый столбец подматрицы sub ( X ) (глобальный входной параметр, тип целый);
DESCX - дескриптор распределенной матрицы X (одномерный массив длины DLEN_); подробное описание см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
FERR - одномерный массив длины LOCc ( N_B), каждый элемент FERR ( j ) которого содержит оценку границы относительной ошибки вычисленного вектора решений X ( j ) (т.е.  j - го столбца матрицы X ( IX : IX + N - 1, JX : JX + NRHS - 1) ), полученную в соответствии с описанием в http://www.netlib.org/scalapack/slug/, раздел 4.6 (см. также [12, 25, 27, 28]) (локальный выходной параметр, тип DOUBLE PRECISION );
BERR - одномерный массив длины LOCc ( N_B); содержащий покомпонентные оценки границ при обратном анализе относительных ошибок для каждого вектора решений X ( j ), т.е. наименьшее относительное возмущение в любом элементе матрицы A ( IA : IA + N - 1, JA : JA + N - 1) или B ( IB : IB + N - 1, JB : JB + NRHS - 1), которое делает X ( j ) точным решением (см. описание в http://www.netlib.org/scalapack/slug/, раздел 4.6 а также [12, 25, 27, 28]) (локальный выходной параметр, тип DOUBLE PRECISION );
WORK - одномерный рабочий массив двойной точности длины LWORK;
на выходе - WORK (1) содержит минимальную требуемую величину LWORK (локальное рабочее пространство, локальный выходной параметр);
LWORK - размер рабочего пространства (массива) WORK;  LWORK - локальный входной параметр и должен быть не меньше, чем
LWORK = PDPORFS ( LWORK ) + LOCr ( N_A);
где PDPORFS ( LWORK ) означает величину рабочего массива WORK, которая требуется при работе вызываемой подпрограммы PDPORFS (см.ниже);
смысл обозначений 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 - то ведущий (угловой) минор порядка k матрицы
A ( IA : IA + k - 1, JA : JA + k - 1) не является положительно определенным, и разложение матрицы не может быть завершено, а решение системы не может быть найдено; 
   Версии
PSPOSV2 -   PCPOSV2     PZPOSV2     решение системы с симметричной (эрмитовой) положительно определенной матрицей методом Холецкого с итерационным уточнением решения и оценкой границ ошибок для вещественных данных одинарной точности, комплексных данных одинарной точности и комплексных данных двойной точности соответственно 
   Вызываемые подпрограммы

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

PDPOTRF - PZPOTRF    Треугольное разложение симметричной (эрмитовой) положительно определенной матрицы методом Холецкого
PDPOTRS - PZPOTRS    Решение системы с симметричной (эрмитовой) положительно определенной матрицей с использованием разложения Холецкого, полученного подпрограммой PDPOTRF(PZPOTRF)
PDPORFS - PZPORFS    Выполнение итерационного уточнения решения системы A X = B, полученного подпрограммой PDPOTRS (PZPOTRS), и оценка границ ошибок полученного решения
PDLANSY - вычисление значения 1 - ой нормы, нормы Фробениуса, бесконечной нормы или наибольшего абсолютного значения любого элемента вещественной симметричной матрицы 
   Замечания по использованию
  1.  В подпрограммах PSPOSV2, PCPOSV2, PZPOSV2 параметры A, B, AF, X и WORK имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно
  2.  В версии подпрограммы PZPOSV2 вместо параметров IWORK и LIWORK целого типа используются параметры с именами RWORK и LRWORK. При этом параметр RWORK имеет тип DOUBLE PRECISION, а параметр LRWORK і 2 * LOCc (N_A).
  3.  Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), PXERBA ( из пакета PBLAS), CHK1MAT, PCHK2MAT, INDXG2P, INFOG2L, LSAME, NUMROC ( из библиотеки ScaLAPACK_TOOLS)

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

Решается система уравнений 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   |

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

0 1 2
0  2  -1
-1   2 		 0   0
 0   0 		 0   0
-1   0 		 0
 0 		 0   0
 0   0
 0   0
 0   0 		 0   0
-1   0 		 0  -1
 0   0 		 0
 0 		 2  -1
-1   2
 0   0  0  -1  0   0  1  0   0
1  0  -1
 0   0 		 0   0
 0   0 		 2  -1
-1   2 		 0
 0 		 0   0
-1   0
 0   0
 0   0 		 2  -1
-1   2 		 0   0
 0   0 		 0
-1 		 0  -1
 0   0

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

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

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

      PROGRAM TDPOSV2
      include 'mpif.h'
      INTEGER     DLEN_, IA, JA, IB, JB, IAF, JAF, IX, JX,
     $                     MB, N, NB, RSRC, CSRC,
     $                     MXLLDA, MXLLDB, MXLLDAF, MXLLDX, NRHS, NBRHS,
     $                     NOUT, MXLOCR, MXLOCC, MXRHSC, LWORK, LIWORK
      PARAMETER   ( DLEN_ = 9, IA = 1, JA = 1, IB = 1, JB = 1,
     $                           IAF = 1, JAF = 1, IX = 1, JX = 1,
     $                           N = 9, MB = 2, NB = 2, RSRC = 0, CSRC = 0,
     $                           MXLLDA = 5, MXLLDB = 5, MXLLDAF = 5, MXLLDX = 5,
     $                           NRHS = 1, NBRHS = 1, NOUT = 6, MXLOCR = 5,
     $                           MXLOCC = 4, MXRHSC = 1, LIWORK = 5, LWORK = 15 )
*     $                           MXLOCC = 4,MXRHSC = 1, LIWORK = -1,LWORK = -1 )

      CHARACTER   UPLO
      PARAMETER   ( UPLO = 'U' )

      INTEGER   ICTXT, INFO, MYCOL, MYROW, NPCOL, NPROW
*
      INTEGER            DESCA( DLEN_ ), DESCB( DLEN_ ),
     $                           DESCAF( DLEN_ ), DESCX( DLEN_ ),
     $                           IWORK( LIWORK )
*     $                           IWORK( 5 )
      DOUBLE PRECISION   A( MXLLDA, MXLOCC ), B( MXLLDB, MXRHSC ),
     $                                      AF( MXLLDAF, MXLOCC ), X( MXLLDX, MXRHSC ),
     $                                      FERR(1), BERR(1),
     $                                      WORK1( NB ), WORK( LWORK )
*     $                                      WORK1( NB ), WORK( 15 )
*
      EXTERNAL   BLACS_EXIT, BLACS_GRIDEXIT, BLACS_GRIDINFO,
     $                       DESCINIT, PDMATINIT, PDPOSV2, 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, MB, NB, RSRC, CSRC, ICTXT, MXLLDA,
     $                               INFO )
      CALL  DESCINIT( DESCB, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
     $                               MXLLDB, INFO )
      CALL DESCINIT( DESCAF, N, N, MB, NB, RSRC, CSRC, ICTXT, MXLLDAF,
     $                               INFO )
      CALL DESCINIT( DESCX, N, NRHS, NB, NBRHS, RSRC, CSRC, ICTXT,
     $                               MXLLDX, INFO )
*
      CALL  PDMATINIT( N, A, IA, JA, DESCA, B, IB, JB, DESCB, INFO )
*
      CALL  PDPOSV2( UPLO, N, NRHS, A, IA, JA, DESCA, AF, IAF, JAF,
     $                               DESCAF, B, IB, JB, DESCB, X, IX, JX, DESCX,
     $                               FERR, BERR, 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            DESCA( * ), DESCB( * ), N, IA, JA, IB, JB,
     $                            INFO
      DOUBLE PRECISION   A( * ), B( * )
*
      INTEGER            I, J
      DOUBLE PRECISION   T, MO, Z
*
      EXTERNAL           PDELSET
*
      T = 2.0D0
      MO = -1.0D0
      Z = 0.0D0
*
      DO 20  J = 1, N
      DO 10  I = 1, N
         IF( I .EQ. J ) THEN
             CALL  PDELSET( A, I, J, DESCA, T )
         ELSE
             CALL  PDELSET( A, I, J ,DESCA, Z )
         END IF
 10 CONTINUE
 20 CONTINUE

      CALL  PDELSET( A, 9, 9, DESCA, 1.0D0 )
*
      DO 40  I = 1, N-1
         J = I + 1
         CALL  PDELSET( A, I, J, DESCA, MO )
 40 CONTINUE

      DO 50  J = 1, N-1
         I = J + 1
         CALL  PDELSET( A, I, J, DESCA, MO )
 50 CONTINUE

      J = 1
      DO 30  I = 1, N
        CALL  PDELSET( B, I, J, DESCB, Z )
 30 CONTINUE

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

Результаты:

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

  X = ( 1.D0, 2.D0, 3.D0, 4.D0, 5.D0, 6.D0, 7.D0, 8.D0, 9.D0 )

 Значение   INFO  =  0,

    FERR  =  0.12745360E-12,   BERR  =  0.11102230E-15