Вычисление всех собственных значений и собственных векторов вещественной симметричной матрицы
На первом этапе подпрограмма приводит симметричную матрицу A к
трехдиагональной форме методом отражений:
A = Q*D*QT, где Q - ортогональная матрица,
а D - симметричная трехдиагональная матрица. На втором этапе неявным
QL или QR алгоритмом находятся собственные значения матрицы D, которые
совпадают с собственными значениями матрицы A, поскольку матрицы A и D
ортогонально подобны. Собственные значения и собственные векторы вычисляются
с машинной точностью.
CALL PDSYEV2 ( UPLO, N, A, IA, JA, DESCA, W,
Z, IZ, JZ, DESCZ, 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, собственные значения располагаются в нем по возрастанию (глобальный выходной параметр); |
| 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 вычисляется самой подпрограммой, если к ней обратиться со значением 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, то подпрограмма не гарантирует точности вычисленных собственных значений. |
Версии
| PSSYEV2 - PCHEEV2 PZHEEV2 | вычисление всех собственных значений и собственных векторов симметричной (эрмитовой) матрицы для случаев вещественных данных одинарной точности, комплексных данных одинарной точности, комплексных данных двойной точности соответственно |
Вызываемые подпрограммы
Здесь указаны только базовые подпрограммы (2 - ого уровня), которые вызываются из целевой подпрограммы (1 - ого уровня).
| PDSYTRD - PZHETRD | приведение симметричной (эрмитовой) матрицы к трехдиагональной форме методом отражений |
| DSTEQR2 - ZSTEQR2 | вычисление всех собственных значений и, возможно, собственных векторов симметричной трехдиагональной матрицы |
| PDORMTR - PZUNMTR | умножение матрицы общего вида на матрицу отражения, вычисленную подпрограммой PDSYTRD (PZHETRD) |
| PDLAMCH - | вычисление машинных параметров для арифметики с плавающей запятой |
| PDLASET - | внедиагональные элементы матрицы полагаются равными ALFA, а диагональные элементы - равными BETA |
| PDLANSY - | вычисление значений 1 - ой нормы, нормы Фробениуса, бесконечной нормы или наибольшего абсолютного значения любого элемента вещественной симметричной, комплексной симметричной или эрмитовой матрицы |
| PDLASCL - | умножение вещественной матрицы на вещественный скаляр |
Замечания по использованию
| 1. | В подпрограммах PSSYEV2, PCHEEV2, PZHEEV2 параметры A, WORK и Z имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно, а параметр W - REAL, REAL и DOUBLE PRECISION соответственно | |
| 2. |
Список параметров подпрограммы PZHEEV2 имеет следующий вид:
PZHEEV2 (UPLO, N, A, IA, JA, DESCA, W, Z, IZ, JZ, DESCZ,
WORK, LWORK, RWORK, LRWORK, INFO),
где дополнительные параметры RWORK и LRWORK означают:
RWORK - одномерный рабочий массив двойной точности длины LRWORK;
на выходе в элементе RWORK (1) возвращается необходимая
длина рабочего массива RWORK (локальное рабочее пространство,
локальный выходной параметр);
LRWORK - задаваемая длина рабочего пространства RWORK;
LRWORK ≥ 2 * N + 2 * N - 2 (локальный входной параметр, тип целый);
| |
| 3. |
Используются подпрограммы BLACS_GRIDINFO ( из пакета BLACS), DCOPY, DSCAL ( из пакета BLAS), LSAME, INDXG2P, NUMROC, CHK1MAT, PCHK2MAT, PXERBLA, SL_GRIDRESHAPE, PDELGET ( из библиотеки ScaLAPACK_TOOLS) |
Вычисление всех собственных значений и собственных векторов симметричной матрицы 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.
Фрагмент фортранного текста вызывающей программы.
(полный текст теста можно получить в tdsyev2.zip,
а принятые обозначения - посмотреть в разделе документации
"Обозначения и упрощения в примерах
по использованию подпрограмм комплекса").
Матрица A инициализируется посредством обращения к подпрограмме PDELSET из библиотеки ScaLAPACK_TOOLS.
PROGRAM TDSYEV2
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, IZ, JZ
INTEGER DESCA( 9 ), DESCZ( 9 )
DOUBLE PRECISION A( LDA, LDA ), Z( 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,
$ PDSYEV2
N = 4
NB = 1
NPROW = 2
NPCOL = 2
IA = 1
JA = 1
IZ = 1
JZ = 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 )
CALL DESCINIT( DESCZ, N, N, NB, NB, 0, 0, CTXT, LDA, INFO )
*
* Построение матрицы A
*
CALL PDMATINIT( N, A, IA, JA, DESCA, INFO )
*
* Вычисление всех собственных значений и собственных векторов
*
CALL PDSYEV2( UPLO, N, A, IA, JA, DESCA, W,
$ Z, IZ, JZ, DESCZ, 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 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
Собственные векторы:
-.485377931880-01 .912072227031-01 .495193043046+00 .862617629821+00
-.122222718667+00 .426620092092+00 -.798642041262+00 .406482218545+00
-.282485135303+00 -.877038303257+00 -.298844131980+00 .248391118466+00
.950214627338+00 -.201197301593+00 -.166273644422+00 .170190725350+00