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

Назначение

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

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

Данная подпрограмма реализует алгоритм вычисления собственных значений, принадлежащих заданному интервалу (VL, VU), и собственных векторов, соответствующих этим собственным значениям, для уравнений вида:

          (1)     Аx  =  λBx   ,

          (2)   АBx  =  λx   ,

          (3)   BАx  =  λx   ,

где  A  и  B - вещественные симметричные матрицы и матрица B положительно определена.

При помощи разложения Холецкого для матрицы B:  В = LLT исходные уравнения (1),(2) или (3) приводятся к стандартному виду Cy = λy, где

    для (1):  C  =  L - 1AL - T,  а   y  =  L-Tx .
    для (2):  C  =  L TAL   ,    а   y  =  L-Tx .
    для (3):  C  =  L TAL   ,    а   y  =  Lx .

Стандартная задача решается путем приведения симметричной матрицы C к трехдиагональному виду и вычисления собственных значений  λ и собственных векторов неявным  QL - или QR - алгоритмом.

Восстановление собственных вектоpов  x  исходных уравнений (1),(2) или (3) из соответствующих вектоpов  у  стандартной задачи осуществляется по формулам:

    x   =   L-T y      для   (1) или (2)  ,

    x   =   L y              для   (3) .

Подробнее см. в //http://srcc.msu.su/num_anal/par_prog/

Собственные значения вычисляются с заданной точностью.

Дж.Х. Уилкинсон, Алгебраическая проблема собственных значений, "Hаука", M., 1970.

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

         CALL  PDSYGV4 ( IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB,
                                         VL, VU, ABSTOL, M, NZ, W, ORFAC, Z, IZ, JZ, DESCZ,
                                         WORK, LWORK, IWORK, LIWORK, IFAIL, ICLUSTR, GAP, INFO)

Параметры

IBTYPE - определяет вид обобщенной проблемы собственных значений:
если IBTYPE = 1 -  вид (1);
если IBTYPE = 2 -  вид (2);
если IBTYPE = 3 -  вид (3);
(глобальный входной параметр, тип целый);
UPLO - если UPLO = ' U ' - задается верхний треугольник матрицы A;
если UPLO = ' L ' - задается нижний треугольник матрицы A
(глобальный входной параметр, тип символьный);
N - порядок распределенной подматрицы sub (A); N ≥ 0 (глобальный входной параметр, тип целый);
A - массив двойной точности, распределенный по процессам блочно - циклическим образом (см.), глобальная размерность которого (N, N), а локальная размерность ( LLD_A, LOCc ( JA + N - 1)),
на входе - это локальная часть распределенной симметричной подматрицы sub(A);
если UPLO = ' U ', то используется только верхняя треугольная часть подматрицы sub(A);
если UPLO = ' L ', то используется только нижняя треугольная часть подматрицы sub(A);
на выходе - нижний треугольник (при UPLO = ' L ' ) или верхний треугольник (при UPLO = ' U ' ) матрицы sub(A), включая диагональ, не сохраняется (локальный входной параметр и локальный выходной параметр);
IA - глобальный номер строки матрицы A, который указывает на начало подматрицы, которая должна быть обработана (глобальный входной параметр, тип целый);
JA - глобальный номер столбца матрицы A, который указывает на начало подматрицы, которая должна быть обработана (глобальный входной параметр, тип целый);
DESCA - дескриптор распределенной матрицы A (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_A = 501),  DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_A = 1),      DLEN ≥ 9;
DESCA содержит информацию о размещении A в памяти; полное описание DESCA см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
B - массив двойной точности, распределенный по процессам блочно - циклическим образом (см.), глобальная размерность которого (N, N), а локальная размерность ( LLD_B, LOCc ( JB + N - 1)),
на входе - это локальная часть распределенной симметричной подматрицы sub(B);
если UPLO = ' U ', то используется только верхняя треугольная часть подматрицы sub(B);
если UPLO = ' L ', то используется только нижняя треугольная часть подматрицы sub(B);
на выходе, если INFO <= N, часть матрицы sub(B) содержит матрицу, замененную на треугольный множитель L из разложения Холецкого  sub(B) = L*L T (локальный входной параметр и локальный выходной параметр);
IB - глобальный номер строки матрицы B, который указывает на начало подматрицы, которая должна быть обработана (глобальный входной параметр, тип целый);
JB - глобальный номер столбца матрицы B, который указывает на начало подматрицы, которая должна быть обработана (глобальный входной параметр, тип целый);
DESCB - дескриптор распределенной матрицы B (одномерный массив длины DLEN);
если тип дескриптора одномерный ( DTYPE_B = 501),  DLEN ≥ 7,
если тип дескриптора двумерный    ( DTYPE_B = 1),      DLEN ≥ 9;
DESCB содержит информацию о размещении B в памяти; полное описание DESCB см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый);
VL - заданная нижняя граница интервала, в котором ищутся собственные значения (глобальный входной параметр, тип DOUBLE PRECISION);
VU - заданная верхняя граница интервала, в котором ищутся собственные значения (глобальный входной параметр, тип DOUBLE PRECISION);
ABSTOL - заданная абсолютная точность, с которой должны быть вычислены собственные значения;
заданная точность считается достигнутой, если собственное значение лежит в интервале [a, b], ширина которого не превосходит ABSTOL + EPS * MAX ( | a |, | b | ), где EPS - машинная точность;
если ABSTOL ≤ 0, то ABSTOL полагается равным EPS * norm ( D ), где norm ( 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, собственные значения располагаются в нем по возрастанию (глобальный выходной параметр);
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)); содержит ортонормированные собственные векторы симметричной матрицы A (локальный выходной параметр);
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) возвращается необходимая длина рабочего пространства
(локальное рабочее пространство, локальный выходной параметр);
LWORK - задаваемая длина рабочего пространства WORK;
LWORK ≥ 5 * N + MAX ( 5 * NN , NP0 * MQ0 + 2 * NB * NB) +
ICEIL( NEIG, NPROW * NPCOL) * NN, где
NP0 = NUMROC ( NN, NB, 0, 0, NPROW) - число строк локальной части матрицы sub (A);
MQ0 = NUMROC ( MAX(NEIG, NB, 2), NB, 0, 0, NPCOL) - число столбцов локальной части матрицы sub (A);
NEIG - число собственных векторов;
NB = DESCA (NB_);
NN = MAX( N, NB, 2 )
Минимальная требуемая величина 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;
если MOD (INFO/16,2) ≠ 0, то IFAIL(1) содержит порядок наименьшего минора матрицы  B, который не является положительно определенным
если (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;
если (MOD ( INFO/2, 2) .NE. 0), то собственные векторы, соответствующие одной или более группам собственных значений, не могут быть переортогонализованы из-за недостаточной величины рабочего пространства; индексы этих групп запоминаются в массиве ICLUSTR;
если (MOD ( INFO/4, 2) .NE. 0), то ограниченность (размеры) рабочего пространства не дает подпрограмме вычислить все собственные векторы в интервале (VL, VU); число вычисленных собственных векторов задается в параметре NZ;
если (MOD(INFO/8,2).NE.0), то значит подпрограмма PDSTEBZ не смогла вычислить собственные значения;
если (MOD(INFO/16,2).NE.0), то значит матрица  B не является положительно определенной, и IFAIL(1) содержит порядок наименьшего минора, который не является положительно определеннным.

Версии

PSSYGV4 -   PCHEGV4     PZHEGV4     вычисление собственных значений, принадлежащих заданному интервалу (VL, VU), и собственных векторов, соответствующих этим собственным значениям, в обобщенной проблеме собственных значений для вещественных симметричных (эрмитовых) матриц A и B для случаев вещественных данных одинарной точности, комплексных данных одинарной точности, комплексных данных двойной точности соответственно

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

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

PDSYNGST - PZHENGST    приведение обобщенной проблемы собственных значений к стандартной форме
PDSYEV4 - PZHEEV4    вычисление собственных значений, принадлежащих заданному интервалу (VL, VU), и собственных векторов, соответствующих этим собственным значениям, симметричной (эрмитовой) матрицы в линейной проблеме собственных значений
PDPOTRF - PZPOTRF    приведение симметричной (эрмитовой) матрицы к трехдиагональной форме методом отражений
PDLAMCH - вычисление машинных параметров для арифметики с плавающей запятой

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

  1.  В подпрограммах PSSYGV4, PCHEGV4, PZHEGV4 параметры A, B, Z и WORK имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно, а параметры ABSTOL, ORFAC, GAP, VL, VU, W - REAL, REAL и DOUBLE PRECISION соответственно
  2.  Список параметров подпрограммы PZHEGV4 имеет следующий вид:

PZHEGV4 (IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB, VL, VU,
    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;
                       LRWORK ≥4*N + MAX( 5*NN, NP0 * MQ0 ) +
                                      ICEIL( NEIG, NPROW*NPCOL)*NN , где

                      NEIG = число требующихся собственных векторов
                      NB = DESCA( NB_ )
                      NN = MAX( N, NB, 2 )
                      NP0 = NUMROC( NN, NB, 0, 0, NPROW )
                      MQ0 = NUMROC( MAX( NEIG, NB, 2 ), NB, 0, 0, NPCOL )
                      (локальный входной параметр, тип целый);
  3.  Используются подпрограммы BLACS_GRIDINFO, DGEBR2D, DGEBS2D ( из пакета BLACS),
DSCAL ( из пакета BLAS),
LSAME, INDXG2P, NUMROC, CHK1MAT, PCHK1MAT, PCHK2MAT, PXERBLA ( из библиотеки ScaLAPACK_TOOLS)

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

Вычисление собственных значений, принадлежащих заданному интервалу (VL, VU), и собственных векторов, соответствующих этим собственным значениям, в обобщенной проблеме собственных значений (1) для вещественных симметричных матриц A и B, которые имеют вид ( заданы только верхние треугольники ):

           |   10    2    3    1    1    |
           |    0   12    1    2    1    |
           |    0    0   11    1   -1    |
           |    0    0    0    9     1    |
           |    0    0    0    0   15    |

           |   12    1   -1    2    1    |
           |    0   14    1   -1    1    |
           |    0    0   16   -1    1    |
           |    0    0    0   12   -1    |
           |    0    0    0    0   11    |

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

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

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

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

      PROGRAM TDSYGV41
      include 'mpif.h'
      INTEGER                         LWORK, LIWORK, MAXN, LDA, MAXPROCS, NOUT
      PARAMETER                  ( LDA = 100, LWORK = 5577, MAXN = 100, LIWORK = 30,
    $                                         MAXPROCS = 512, NOUT = 6 )
      CHARACTER                 UPLO
      PARAMETER                 ( UPLO = 'U' )
      INTEGER                        CTXT, I, IAM, INFO, MYCOL, MYROW, N, NB, IZ, JZ,
    $                                        NPCOL, NPROCS, NPROW, IA, JA, IB, JB, IBTYPE
      DOUBLE PRECISION    ABSTOL, ORFAC, VL, VU
*
      INTEGER                        DESCA( 9 ), DESCB(9), DESCZ(9), IWORK( LIWORK ), IFAIL(5), ICLUSTR(8)
      DOUBLE PRECISION   A( LDA, LDA ),  B( LDA, LDA ), Z( LDA, LDA ), GAP( 4 ), W( MAXN ),
    $                                        WORK( LWORK ), WORK1( 5 )
*
      EXTERNAL                    BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
    $                                        BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
    $                                        BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT,
    $                                        PDSYGV4
*
      N = 5
      NB = 2
      NPROW = 2
      NPCOL = 2
      IA = 1
      JA = 1
      IB = 1
      JB = 1
      IZ = 1
      JZ = 1
      VL = 0.5D0
      VU = 1.5D0
      IBTYPE = 1
      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( DESCB, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
      CALL  DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )

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

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

* Решение обобщенной проблемы собственных значений

      CALL  PDSYGV4( IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB, VL, VU, 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, B, IB, JB, DESCB, INFO )
*
* PDMATINIT генерирует и распределяет матрицы A и B по решетке процессов
*
      INTEGER                        IA, INFO, JA, N, IB, JB
      INTEGER                        DESCA( * ), DESCB( * )
      DOUBLE PRECISION   A( * ), B( * )
      INTEGER                        I, J, MYCOL, MYROW, NPCOL, NPROW
      EXTERNAL                    PDELSET
*
      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)
       CALL  PDELSET( B, I, J, DESCB, 0.0D0)
   1 CONTINUE
   2 CONTINUE
*
      DO 3  I = 1, N-1
          J = I + 1
            CALL  PDELSET( A, I, J, DESCA, 1.0D0)
            CALL  PDELSET( B, I, J, DESCB, 1.0D0)
         END IF
   3 CONTINUE
*
      J = N
      DO 5  I = 1, N-1
       CALL  PDELSET( A, I, J, DESCA, 1.0D0)
       CALL  PDELSET( B, I, J, DESCB, 1.0D0)
   5 CONTINUE
*
       CALL  PDELSET( A, 1, 1, DESCA, 10.0D0)
       CALL  PDELSET( A, 1, 2, DESCA, 2.0D0)
       CALL  PDELSET( A, 1, 3, DESCA, 3.0D0)
       CALL  PDELSET( A, 1, 4, DESCA, 1.0D0)
       CALL  PDELSET( A, 2, 2, DESCA, 12.0D0)
       CALL  PDELSET( A, 2, 4, DESCA, 2.0D0)
       CALL  PDELSET( A, 3, 3, DESCA, 11.0D0)
       CALL  PDELSET( A, 3, 5, DESCA, -1.0D0)
       CALL  PDELSET( A, 4, 4, DESCA, 9.0D0)
       CALL  PDELSET( A, 5, 5, DESCA, 15.0D0)
*
       CALL  PDELSET( B, 1, 1, DESCA, 12.0D0)
       CALL  PDELSET( B, 1, 3, DESCA, -1.0D0)
       CALL  PDELSET( B, 1, 4, DESCA, 2.0D0)
       CALL  PDELSET( B, 2, 2, DESCA, 14.0D0)
       CALL  PDELSET( B, 2, 4, DESCA, -1.0D0)
       CALL  PDELSET( B, 3, 3, DESCA, 16.0D0)
       CALL  PDELSET( B, 3, 4, DESCA, -1.0D0)
       CALL  PDELSET( B, 4, 4, DESCA, 12.0D0)
       CALL  PDELSET( B, 4, 5, DESCA, -1.0D0)
       CALL  PDELSET( B, 5, 5, DESCA, 11.0D0)
*
      RETURN
      END

Результаты:

 Значение  INFO  =   0
                        M = 4

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

  W(1) =  0.663662748392315
  W(2) =  0.943859004668386
  W(3) =  1.10928454001752
  W(4) =  1.49235323254300

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

   Z(1, 1) =  0.829198064861523532D-01
   Z(2, 1) =  0.153148395665945264D+00
   Z(3, 1) = -0.118603667911389210D+00
   Z(4, 1) = -0.182813041785799935D+00
   Z(5, 1) =  0.356172036818147315D-02

   Z(1, 2) =  0.191710031573882250D+00
   Z(2, 2) = -0.158991211513700409D+00
   Z(3, 2) =  0.748390709386720365D-01
   Z(4, 2) = -0.137468929466840417D+00
   Z(5, 2) =  0.889778923495063573D-01

   Z(1, 3) =  0.142011959884590616D+00
   Z(2, 3) =  0.142419950546665203D+00
   Z(3, 3) =  0.120997623004495347D+00
   Z(4, 3) =  0.125531015188065975D+00
   Z(5, 3) =  0.769220728307276440D-02

   Z(1, 4) = -0.763867178777808409D-01
   Z(2, 4) =  0.170980018713419862D-01
   Z(3, 4) = -0.666645336709070502D-01
   Z(4, 4) =  0.860480093056616990D-01
   Z(5, 4) =  0.289433414168931813D+00