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

Назначение

Вычисление собственных значений вещественной симметричной матрицы, принадлежащих заданному интервалу индексов (IL, IU), и собственных векторов, соответствующих этим собственным значениям

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

На первом этапе подпрограмма приводит симметричную матрицу A к трехдиагональной форме методом отражений:
A = Q*D*QT, где Q - ортогональная матрица, а D - симметричная трехдиагональная матрица.
На втором этапе неявным QL или QR алгоритмом находятся собственные значения матрицы D, принадлежащие заданному интервалу индексов (IL, IU), которые совпадают с собственными значениями матрицы A, поскольку матрицы A и D ортогонально подобны.
На третьем этапе вычисляются собственные векторы трехдиагональной матрицы D при помощи метода обратных итераций. Программа не гарантирует ортогональность тех из собственных векторов, которые соответствуют близко расположенным собственным значениям (см. параметры ORFAC, ICLUSTR, GAP). Однако по желанию пользователя (см. параметр ORFAC) может осуществляться процесс переортогонализации собственных векторов с помощью модифицированного метода Грама - Шмидта (см., например, [1] ).
Для получения собственных векторов матрицы A, матрица Z, состоящая из собственных векторов - столбцов трехдиагональной матрицы D, умножается слева на ортогональную матрицу Q.

1. Воеводин В.В. Вычислительные основы линейной алгебры, М., "Наука", 1977, стр. 128 - 135.

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

         CALL  PDSYEV6 (UPLO, N, A, IA, JA, DESCA, IL, IU,
                                        ABSTOL, M, NZ, W, ORFAC, Z, IZ, JZ,
                                        DESCZ, WORK, LWORK, IWORK, LIWORK, IFAIL,
                                        ICLUSTR, GAP, 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 см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
IL - заданная нижняя граница интервала индексов, в котором ищутся собственные значения;
предполагается, что собственные значения располагаются по возрастанию; IL ≥ 1; (глобальный входной параметр, тип целый);
IU - заданная верхняя граница интервала индексов, в котором ищутся собственные значения;
предполагается, что собственные значения располагаются по возрастанию; min ( IL, N) ≤ IU ≤ N; (глобальный входной параметр, тип целый);
ABSTOL - заданная абсолютная точность, с которой должны быть вычислены собственные значения;
заданная точность считается достигнутой, если собственное значение лежит в интервале [a, b], ширина которого не превосходит ABSTOL + EPS * MAX ( | a |, | b | ), где EPS - машинная точность;
если ABSTOL ≤ 0, то ABSTOL полагается равным EPS * norm ( D ), где norm ( D ) - первая норма трехдиагональной матрицы  D, полученной посредством преобразования матрицы A к трехдиагональной форме;
собственные значения будут вычислены с наибольшей точностью, если ABSTOL = 2 * PDLAMCH ('S'), где PDLAMCH ('S') = sfmin - минимальное вещественное число, при котором 1 / sfmin не вызывает переполнения (глобальный входной параметр, тип DOUBLE PRECISION);
M - общее количество собственных значений, найденных в заданном интервале индексов; 0 ≤ M ≤ N; (глобальный выходной параметр, тип целый);
NZ - общее количество найденных собственных векторов; 0 ≤ NZ ≤ M; число столбцов матрицы Z, содержащих вычисленные собственные векторы;
NZ = M, если пользователь задал рабочее пространство достаточных размеров и подпрограмма была в состоянии определить это до начала вычислений; чтобы получить все требующиеся собственные векторы, пользователь должен задать как достаточных размеров матрицу Z, так и достаточное рабочее пространство для процесса вычисления векторов (cм. LWORK);
(глобальный выходной параметр, тип целый);
W - одномерный массив двойной точности длины N;
если INFO = 0, первые M элементов массива содержат вычисленные собственные значения в возрастающем порядке (глобальный выходной параметр);
ORFAC - определяет, какие собственные векторы должны быть подвергнуты процессу переортогонализации [1]; собственные векторы, которые соответствуют собственным значениям, находящимся на расстоянии  tol = ORFAC * norm (A) друг от друга, должны быть переортогонализованы; однако если рабочего пространства недостаточно (см. LWORK),  tol может быть уменьшено, чтобы все собственные векторы, которые должны быть переортогонализованы, могли быть заданы в одном процессе;
если ORFAC = 0, переортогонализация не выполняется; если ORFAC задан отрицательным (< 0), то он полагается равным 10 - 3
(глобальный входной параметр, тип DOUBLE PRECISION);
Z - массив двойной точности с глобальной размерностью (N, N) и с локальной размерностью (LLD_Z, LOCc ( JZ + N - 1));
если INFO = 0, первые M столбцов матрицы Z содержат вычисленные ортонормированные собственные векторы симметричной матрицы A, соответствующие вычисленным собственным значениям; если какой - либо собственный вектор не может быть вычислен с заданной точностью, то соответствующий столбец матрицы Z содержит последнее приближение к этому собственному вектору, а его индекс сохраняется в IFAIL;
(локальный выходной параметр);
IZ - глобальный номер строки матрицы Z, который указывает на начало строк подматрицы sub (Z) (глобальный входной параметр, тип целый);
JZ - глобальный номер столбца матрицы Z, который указывает на начало столбцов подматрицы sub (Z) (глобальный входной параметр, тип целый);
DESCZ - дескриптор распределенной матрицы Z (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_Z = 501),  то DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_Z = 1),      то DLEN ≥ 9;
DESCZ содержит информацию о размещении Z в памяти; полное описание DESCZ см. в разделе документации "Дескрипторы глобальных массивов"; DESCZ (CTXT_) должен быть равен DESCA (CTXT_) (глобальный и локальный входной параметр, тип целый);
WORK - одномерный рабочий массив двойной точности длины LWORK;
на выходе в элементе WORK (1) возвращается необходимая длина рабочего пространства, при которой, однако, не гарантируется ортогональность вычисленных собственных векторов. Процесс переортогонализации выполняется только в том случае, если указанная величина рабочего пространства будет увеличина на (L - 1) * N, где L - наибольшее число собственных значений, входящих в группу подряд идущих близко расположенных собственных значений, а N - порядок матрицы A (см. параметры ORFAC, ICLUSTR);
(локальное рабочее пространство, локальный выходной параметр);
LWORK - задаваемая длина рабочего пространства WORK;
минимальная требуемая величина LWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LWORK = - 1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LWORK возвращается в элементе массива WORK (1); (локальный или глобальный входной параметр, тип целый);
IWORK - одномерный рабочий массив целого типа длины LIWORK;
на выходе в элементе IWORK (1) содержится необходимая длина рабочего пространства IWORK (локальное рабочее пространство, локальный выходной параметр);
LIWORK - задаваемая длина рабочего пространства IWORK;
LIWORK ≥ 6 * NNP, где NNP = MAX ( N, NPROW * NPCOL + 1, 4 );
минимальная требуемая величина LIWORK может быть вычислена самой подпрограммой, если к ней обратиться со значением LIWORK = - 1; при этом подпрограмма не производит никаких других вычислений, а требуемое значение LIWORK возвращается в элементе массива IWORK (1);
(локальный или глобальный входной параметр, тип целый);
IFAIL - одномерный массив целого типа длины N;
если INFO = 0, первые M элементов массива равны 0;
если на выходе (MOD ( INFO, 2) .NE. 0), то IFAIL содержит индексы собственных векторов, которые не были вычислены с заданной точностью;
(глобальный выходной параметр);
ICLUSTR - одномерный массив целого типа длины (2 * NPROW * NPCOL);
в соответствии с параметром ORFAC (> 0) среди вычисленных собственных значений выделяются одна или несколько групп близко расположенных подряд идущих собственных значений; пусть  I обозначает номер такой группы среди других выделенных групп;
целый массив ICLUSTR содержит пары чисел, обозначающих индексы первого и последнего из собственных значений, входящих в каждую из выделенных групп собственных значений, которым соответствуют собственные векторы, не подвергавшиеся процессу переортогонализации из - за недостаточных размеров выделенного рабочего пространства; следовательно, собственные векторы с индексами от ICLUSTR (2 * I - 1) до ICLUSTR (2 * I ) могут быть не ортогональными;
( ICLUSTR (2 * K) .NE. 0 .AND. ICLUSTR (2 * K + 1) .EQ. 0), если и только если K равно числу групп;
(глобальный выходной параметр);
GAP - одномерный массив двойной точности длины NPROW * NPCOL;
на выходе этот массив содержит промежуток (разность) между собственными значениями, собственные векторы которых не могут быть переортогонализованы;
значения в этом массиве соответствуют группам собственных значений, индексы которых содержатся в массиве ICLUSTR; в результате скалярное произведение собственных векторов, соответствующих I - ой группе, может быть определено величиной (C * n) / GAP ( I ), где  C - небольшая константа;
(глобальный выходной параметр);
INFO - целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр)
= 0 - успешное завершение работы;
< 0 - если i - ый фактический параметр подпрограммы является массивом и его j - ый элемент имеет недопустимое значение, тогда INFO = - ( i * 100 + j );
если i - ый фактический параметр является скаляром и имеет недопустимое значение, тогда INFO = - i;
> 0 - если (MOD ( INFO,2) .NE. 0), то один или более собственных векторов не могут быть вычислены с заданной точностью; их индексы запоминаются в массиве IFAIL; проверьте, что ABSTOL = 2.0 * PDLAMCH('U');
если (MOD ( INFO/2, 2) .NE. 0), то собственные векторы, соответствующие одной или более группам собственных значений, не могут быть переортогонализованы из-за недостаточной величины рабочего пространства; индексы этих групп запоминаются в массиве ICLUSTR;
если (MOD ( INFO/4, 2) .NE. 0), то ограниченность (размеры) рабочего пространства не дает подпрограмме вычислить все собственные векторы в интервале индексов (IL, IU); число вычисленных собственных векторов задается в параметре NZ. 
   Версии
PSSYEV6 -   PCHEEV6     PZHEEV6     вычисление собственных значений и собственных векторов симметричной матрицы, принадлежащих заданному интервалу индексов, для случаев вещественных данных одинарной точности, комплексных данных одинарной точности (эрмитова матрица), комплексных данных двойной точности (эрмитова матрица) соответственно 
   Вызываемые подпрограммы

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

PDSYNTRD - PZHENTRD    приведение симметричной (эрмитовой) матрицы к трехдиагональной форме методом отражений
PDSTEBZ - вычисление собственных значений (всех, в заданном интервале, или в интервале индексов) симметричной трехдиагональной матрицы
PDSTEIN - PZSTEIN    вычисление собственных векторов симметричной (эрмитовой) трехдиагональной матрицы посредством обратной итерации
PDORMTR - PZUNMTR    умножение матрицы общего вида на матрицу отражения, вычисленную подпрограммой PDSYTRD (PZHETRD)
PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой
PDLASET - внедиагональные элементы матрицы полагаются равными ALFA, а диагональные элементы - равными BETA
PDLANSY - PZLANHE    вычисление значений 1 - ой нормы, нормы Фробениуса, бесконечной нормы или наибольшего абсолютного значения любого элемента вещественной симметричной, комплексной симметричной или эрмитовой матрицы
PDLASCL - PZLASCL    умножение вещественной (эрмитовой) матрицы на вещественный скаляр
PJLAENV - выбирает проблемно - зависимые параметры для локального окружения 
   Замечания по использованию
  1.  В подпрограммах PSSYEV6, PCHEEV6, PZHEEV6 параметры A, WORK и Z имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно, а параметры ABSTOL, ORFAC, W - тип REAL, REAL и DOUBLE PRECISION соответственно
  2.  Список параметров подпрограммы PZHEEV6 имеет следующий вид:

PZHEEV6 (UPLO, N, A, IA, JA, DESCA, IL, IU, ABSTOL, M, NZ, W, ORFAC, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, IWORK, LIWORK, IFAIL, ICLUSTR, GAP, INFO),
где дополнительные параметры RWORK и LRWORK означают: 
   RWORK - одномерный рабочий массив двойной точности длины LRWORK;
                      на выходе в элементе RWORK (1) возвращается необходимая
                      длина рабочего массива RWORK (локальное рабочее пространство,
                      локальный выходной параметр);
 LRWORK - задаваемая длина рабочего пространства RWORK;
                       (локальный входной параметр, тип целый);
  3.  Используются подпрограммы BLACS_GRIDINFO, DGEBS2D, DGEBR2D (из пакета BLACS), DCOPY, DSCAL (из пакета BLAS), LSAME, INDXG2P, NUMROC, CHK1MAT, PCHK2MAT, PXERBLA, SL_GRIDRESHAPE, PDELGET(PZELGET) (из библиотеки ScaLAPACK_TOOLS)

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

Вычисление собственных значений, принадлежащих заданному интервалу индексов (IL, IU), и соответствующих им собственных векторов вещественной симметричной матрицы 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, IL = 1, IU = 3. Тогда матрица 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 = 2.

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

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

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

      PROGRAM TDSYEV6
      include 'mpif.h'
      INTEGER                         LWORK, LIWORK, MAXN, LDA, MAXPROCS, NOUT
      DOUBLE PRECISION    ZERO, MONE
      PARAMETER                  ( LDA = 100, LWORK = 5592, LIWORK = 30,
     $                                        MAXN = 100,
     $                                        MAXPROCS = 64, 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, IZ, JZ, IU, IL
      DOUBLE PRECISION     ABSTOL

      INTEGER                          DESCA( 9 ), DESCZ( 9 ), IWORK(LIWORK),
     $                                         IFAIL(4), ICLUSTR(8)
      DOUBLE PRECISION     A( LDA, LDA ), W( MAXN ), GAP(4),
     $                                         WORK( LWORK ), Z( LDA, LDA ), WORK1( 4 )

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

      N = 4
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
      IZ = 1
      JZ = 1
      IL = 1
      IU = 3
      ORFAC = 0.0D0

      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
      ABSTOL = 2.0D0* PDLAMCH(CTXT,'U')

      CALL DESCINIT( DESCA, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
      CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
*
* Построение матрицы A
*
      CALL PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* Вычисление собственных значений и собственных векторов
*
      CALL  PDSYEV6( UPLO, N, A, IA, JA, DESCA, IL, IU, ABSTOL, M, NZ, W,
     $                              ORFAC, Z, IZ, JZ, DESCZ, WORK, LWORK, IWORK, LIWORK,
     $                              IFAIL, ICLUSTR, GAP, 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
                         M  =  3
                       NZ  =  3

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

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

 Собственные векторы:

 -.485377931880-01    .912072227031-01    .495193043046+00

 -.122222718667+00   .426620092092+00  -.798642041262+00

 -.282485135303+00  -.877038303257+00  -.298844131980+00

  .950214627338+00  -.201197301593+00  -.166273644422+00