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

Назначение

Вычисление всех собственных значений вещественной симметричной матрицы

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

На первом этапе подпрограмма приводит симметричную матрицу A к трехдиагональной форме методом отражений:
A =  Q*D*QT, где Q - ортогональная матрица, а D - симметричная трехдиагональная матрица. На втором этапе неявным QL или QR алгоритмом находятся собственные значения матрицы D, которые совпадают с собственными значениями матрицы A, поскольку матрицы A и D ортогонально подобные. Собственные значения вычисляются с машинной точностью.

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

         CALL  PDSYEV1 ( UPLO, N, A, IA, JA, DESCA,
                                         W, WORK, LWORK, INFO)

   Параметры
UPLO - если UPLO = ' U ' - задается верхний треугольник матрицы A;
если UPLO = ' L ' - задается нижний треугольник матрицы A
(глобальный входной параметр, тип символьный);
N - порядок распределенной подматрицы sub (A); N ≥ 0 (глобальный входной параметр, тип целый);
A - массив двойной точности, распределенный по процессам блочно - циклическим образом (см.), глобальная размерность которого (N, N), а локальная размерность ( LLD_A, LOCc ( JA + N - 1)),
на входе - это симметричная матрица A;
если UPLO = ' U ', то используется только верхняя треугольная часть матрицы A;
если UPLO = ' L ', то используется только нижняя треугольная часть матрицы A;
на выходе - нижний треугольник (при UPLO = ' L ' ) или верхний треугольник (при UPLO = ' U ' ) матрицы A, включая диагональ, не сохраняется (локальный входной параметр, рабочее пространство);
IA - глобальный номер строки матрицы A, который указывает на начало подматрицы, которая должна быть обработана (глобальный входной параметр, тип целый);
JA - глобальный номер столбца матрицы A, который указывает на начало подматрицы, которая должна быть обработана (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_A = 501),  DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_A = 1),      DLEN ≥ 9;
DESCA содержит информацию о размещении A в памяти; полное описание DESCA см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
W - одномерный массив двойной точности длины N; если  INFO = 0, собственные значения располагаются в нем по возрастанию (глобальный выходной параметр);
WORK - одномерный рабочий массив двойной точности длины LWORK;
на выходе, в элементе WORK (1) возвращается необходимая длина рабочего пространства
(локальное рабочее пространство, локальный выходной параметр);
LWORK - задаваемая длина рабочего пространства WORK;
LWORK ≥ 5 * N + MAX ( NB * (NP + 1), 3 * NB) + 1, где
NP = NUMROC ( N, NB, MYROW, IAROW, NPROW) - число строк локальной части матрицы sub (A);
NB = DESCA (NB_);
Минимальная требуемая величина LWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LWORK = -1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LWORK возвращается в элементе массива WORK (1);
(локальный или глобальный входной параметр, тип целый);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если i - ый фактический параметр подпрограммы является массивом и его j - ый элемент имеет недопустимое значение, тогда INFO = - ( i * 100 + j ), если i - ый фактический параметр является скаляром и имеет недопустимое значение, тогда INFO = - i;
> 0 - если INFO = i, где 1 ≤ i ≤ N, то i - ое собственное значение не сходится после выполнения 30 * N итераций;
если INFO = N+1, то подпрограмма не гарантирует точности вычисленных собственных значений. 
   Версии
PSSYEV1 -   PCHEEV1     PZHEEV1     вычисление всех собственных значений симметричной (эрмитовой) матрицы для случаев вещественных данных одинарной точности, комплексных данных одинарной точности, комплексных данных двойной точности соответственно 
   Вызываемые подпрограммы

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

PDSYTRD - PZHETRD    приведение симметричной (эрмитовой) матрицы к трехдиагональной форме методом отражений
DSTEQR2 - ZSTEQR2    вычисление всех собственных значений и, возможно, собственных векторов симметричной трехдиагональной матрицы
PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой
PDLASET - внедиагональные элементы матрицы полагаются равными ALFA, а диагональные элементы - равными BETA
PDLANSY - вычисление значений 1 - ой нормы, нормы Фробениуса, бесконечной нормы или наибольшего абсолютного значения любого элемента вещественной симметричной, комплексной симметричной или эрмитовой матрицы
PDLASCL - умножение вещественной матрицы на вещественный скаляр 
   Замечания по использованию
  1.  В подпрограммах PSSYEV1, PCHEEV1, PZHEEV1 параметры A и WORK имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно, а параметр W - REAL, REAL и DOUBLE PRECISION соответственно
  2.  Список параметров подпрограммы PZHEEV1 имеет следующий вид:

PZHEEV1 (UPLO, N, A, IA, JA, DESCA, W, WORK, LWORK, RWORK, LRWORK, INFO),
где дополнительные параметры RWORK и LRWORK означают: 
   RWORK - одномерный рабочий массив двойной точности длины LRWORK;
                      на выходе в элементе RWORK (1) возвращается необходимая
                      длина рабочего массива RWORK (локальное рабочее пространство,
                      локальный выходной параметр);
 LRWORK - задаваемая длина рабочего пространства RWORK;
                       LRWORK ≥ 2 * N (локальный входной параметр, тип целый);
  3.  Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS),
DCOPY, DSCAL ( из пакета BLAS),
LSAME, INDXG2P, NUMROC, CHK1MAT, PCHK1MAT, PXERBLA ( из библиотеки ScaLAPACK_TOOLS)

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

1. Вычисление всех собственных значений симметричной матрицы A, которая имеет вид: 

           |   1   0   0   0   0   0   1   |
           |   0   1   0   0   0   0   2   |
           |   0   0   1   0   0   0   3   |
           |   0   0   0   1   0   0   4   |
           |   0   0   0   0   1   0   5   |
           |   0   0   0   0   0   1   6   |
           |   1   2   3   4   5   6   7   |

Порядок матрицы N = 7.
Пусть матрица разбивается на блоки с размером блоков NB = 2.

Осуществляется пропуск программы на 4 процессах, которые образуют решетку 2 на 2.

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

Матрица A инициализируется посредством обращения к подпрограмме PDELSET из библиотеки ScaLAPACK_TOOLS. 

      PROGRAM TDSYEV11
      include 'mpif.h'
      INTEGER                         LWORK, MAXN, LDA, MAXPROCS, NOUT
      DOUBLE PRECISION    ZERO, MONE
      PARAMETER                  ( LDA = 100, LWORK = 264, MAXN = 100,
    $                                         MAXPROCS = 512, NOUT =6, ZERO = 0.0D+0,
    $                                         MONE = -1.0D+0 )
      CHARACTER                 UPLO
      PARAMETER                 ( UPLO = 'U' )
      INTEGER                        CTXT, I, IAM, INFO, MYCOL, MYROW, N, NB,
    $                                        NPCOL, NPROCS, NPROW, IA, JA
*     ..
      INTEGER                        DESCA( 9 )
      DOUBLE PRECISION   A( LDA, LDA ), W( MAXN ),
    $                                        WORK( LWORK ), WORK1( 7 )
*     ..
      EXTERNAL                    BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
    $                                        BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
    $                                        BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT,
    $                                        PDSYEV1
*     ..
      N = 7
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
*
      CALL  BLACS_PINFO( IAM, NPROCS )
      IF( ( NPROCS .LT. 1 ) ) THEN
         CALL  BLACS_SETUP( IAM, NPROW*NPCOL )
      END IF
*
      CALL  BLACS_GET( -1, 0, CTXT )
      CALL  BLACS_GRIDINIT( CTXT, 'R', NPROW, NPCOL )
      CALL  BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )

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

      CALL  DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )

* Построение матрицы A

      CALL  PDMATINIT( N, A, IA, JA, DESCA, INFO )

* Вычисление собственных значений

      CALL  PDSYEV1( UPLO, N, A, IA, JA, DESCA, W, WORK, LWORK, INFO )

      CALL  BLACS_GRIDEXIT( CTXT )
 20 CONTINUE

      CALL  BLACS_EXIT( 0 )
      STOP
      END
*
      SUBROUTINE  PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* PDMATINIT генерирует и распределяет матрицу A по решетке процессов
*
      DOUBLE PRECISION   ONE
      PARAMETER                 ( ONE = 1.0D+0 )
*
      INTEGER                        IA, INFO, JA, N
      INTEGER                        DESCA( * )
      DOUBLE PRECISION   A( * )
      INTEGER                        I, J, MYCOL, MYROW, NPCOL, NPROW
      EXTERNAL                    PDELSET
      INTRINSIC                     DBLE
*
      INFO = 0
*
      IF( IA .NE. 1 ) THEN
         INFO = -3
      ELSE IF( JA .NE. 1 ) THEN
         INFO = -4
      END IF
*
      DO 2 J = 1, N
      DO 1 I = 1, N
      CALL  PDELSET( A, I, J, DESCA, 0.0D0)
   1 CONTINUE
   2 CONTINUE
      DO 4 J = 1, N
      DO 3 I = 1, N
         IF( I .EQ. J ) THEN
            CALL  PDELSET( A, I, J, DESCA, 1.0D0)
         END IF
   3 CONTINUE
   4 CONTINUE
      J = N
      DO 5 I = 1, N
      CALL  PDELSET( A, I, J, DESCA, DBLE(I) )
   5 CONTINUE
      I = N         
      DO 6 J = 1, N
      CALL  PDELSET( A, I, J, DESCA, DBLE(J) )
   6 CONTINUE
*
      RETURN
      END

Результаты:

 Значение  INFO  =   0

 Собственные значения:

  W(1) = -6.D+00
  W(2) =   1.D+00
  W(3) =   1.D+00
  W(4) =   1.D+00
  W(5) =   1.D+00
  W(6) =   1.D+00
  W(7) = 14.D+00

2. Вычисление всех собственных значений симметричной матрицы A, элементы которой задаются формулами: 


                                      1
          A( I, J )   =   --------------                                       при  I  ≠  J
                                 I + J - 1


                                      1                      N - I + 1
          A( I, J )   =   --------------    +    ---------------       при  I  =  J
                                 I + J - 1                     N

Пусть порядок матрицы N = 4. Тогда матрица A имеет вид 

            |    2            1/2            1/3                 1/4        |

            |   1/2    (1/3 + 3/4)       1/4                1/5        |
   
            |   1/3          1/4       (1/5 + 1/2)          1/6        |

            |   1/4          1/5             1/6         (1/7 + 1/4)   |

Пусть матрица разбивается на блоки с размером блоков NB = 1 (т.е. каждый элемент матрицы является отдельным блоком).

Осуществляется пропуск программы на 4 процессах, которые образуют решетку 2 на 2.

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

      PROGRAM TDSYEV1
      include 'mpif.h'
      INTEGER                        LWORK, MAXN, LDA, MAXPROCS, NOUT
      DOUBLE PRECISION   ZERO, MONE
      PARAMETER                 ( MAXN = 100, LDA = 100, LWORK = 264,
    $                                        MAXPROCS = 512, NOUT =6, ZERO = 0.0D+0,
    $                                        MONE = -1.0D+0 )
      CHARACTER                 UPLO
      PARAMETER                 ( UPLO = 'U' )

      INTEGER                       CTXT, I, IAM, INFO, MYCOL, MYROW, N, NB,
    $                                        NPCOL, NPROCS, NPROW, IA, JA

      INTEGER                        DESCA( 9 )
      DOUBLE PRECISION   A( LDA, LDA ), W( MAXN ),
    $                                        WORK( LWORK ), WORK1( 4 )

      EXTERNAL                    BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
    $                                        BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
    $                                        BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT,
    $                                        PDSYEV1

      N = 4
      NB = 1
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1

      CALL  BLACS_PINFO( IAM, NPROCS )
      IF( ( NPROCS .LT. 1 ) ) THEN
         CALL  BLACS_SETUP( IAM, NPROW*NPCOL )
      END IF

      CALL  BLACS_GET( -1, 0, CTXT )
      CALL  BLACS_GRIDINIT( CTXT, 'R', NPROW, NPCOL )
      CALL  BLACS_GRIDINFO( CTXT, NPROW, NPCOL, MYROW, MYCOL )

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

      CALL  DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
*
* Построение матрицы A
*
      CALL  PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* Вычисление собственных значений
*
      CALL  PDSYEV1( UPLO, N, A, IA, JA, DESCA, W, WORK, LWORK, INFO )
*
      CALL  BLACS_GRIDEXIT( CTXT )
 20 CONTINUE

      CALL  BLACS_EXIT( 0 )
      STOP
      END
*
      SUBROUTINE PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* PDMATINIT генерирует и распределяет матрицу A по решетке процессов
*
      DOUBLE PRECISION   ONE
      PARAMETER                 ( ONE = 1.0D+0 )
      INTEGER                        IA, INFO, JA, N
      INTEGER                        DESCA( * )
      DOUBLE PRECISION   A( * )
      INTEGER                        I, J, MYCOL, MYROW, NPCOL, NPROW
      EXTERNAL                    BLACS_GRIDINFO, PDELSET
      INTRINSIC                     DBLE
*
      INFO = 0
*
      IF( IA .NE. 1 ) THEN
         INFO = -3
      ELSE IF( JA .NE. 1 ) THEN
         INFO = -4
      END IF
*
      DO 20 J = 1, N
         DO 10 I = 1, N
           IF( I .EQ. J ) THEN
             CALL  PDELSET( A, I, J, DESCA,
     $                                 ( DBLE( N-I+1 ) ) / DBLE( N )+ONE /
     $                                 ( DBLE( I+J )-ONE ) )
           ELSE
             CALL  PDELSET( A, I, J, DESCA, ONE / ( DBLE( I+J )-ONE ) )
           END IF
 10    CONTINUE
 20 CONTINUE
      RETURN
      END

Результаты:

 Значение  INFO  =   0

 Собственные значения:

  W(1) = 0.30481403606055D+00
  W(2) = 0.58196120214202D+00
  W(3) = 0.90849811000513D+00
  W(4) = 0.23809171279828D+01