Вычисление всех собственных значений и собственных векторов в обобщенной проблеме собственных значений для вещественных симметричных матриц A и B
Данная подпрограмма реализует алгоритм вычисления всех собственных значений и собственных векторов уравнений вида:
(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 PDSYGV2 ( IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB,
W, Z, IZ, JZ, DESCZ, WORK, LWORK, IFAIL, 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 см. в разделе документации "Дескрипторы глобальных массивов" (глобальный и локальный входной параметр, тип целый); |
| 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 ≥ 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); (локальный входной параметр, тип целый); |
| IFAIL - |
обеспечивает дополнительную информацию, когда
INFO ≠ 0; если MOD (INFO/16,2) ≠ 0, то IFAIL содержит порядок наименьшего минора матрицы B, который не является положительно определенным (глобальный выходной параметр, тип целый); |
| INFO - | целая переменная, диагностирующая результат работы подпрограммы (глобальный выходной параметр) |
| = 0 - | успешное завершение работы; |
| < 0 - | если i - ый фактический параметр подпрограммы является массивом и его j - ый элемент имеет недопустимое значение, тогда INFO = - ( i * 100 + j ), если i - ый фактический параметр является скаляром и имеет недопустимое значение, тогда INFO = - i; |
| > 0 - |
если (MOD(INFO/8,2).NE.0), то значит PDSTEBZ не смогла вычислить собственные
значения; если (MOD(INFO/16,2).NE.0), то значит матрица B не является положительно определенной, и IFAIL содержит порядок наименьшего минора, который не является положительно определеннным. |
Версии
| PSSYGV2 - PCHEGV2 PZHEGV2 | вычисление всех собственных значений в обобщенной проблеме собственных значений для вещественных симметричных (эрмитовых) матриц A и B для случаев вещественных данных одинарной точности, комплексных данных одинарной точности, комплексных данных двойной точности соответственно |
Вызываемые подпрограммы
Здесь указаны только базовые подпрограммы (2 - ого уровня), которые вызываются из целевой подпрограммы (1 - ого уровня).
| PDSYNGST - PZHENGST | приведение обобщенной проблемы собственных значений к стандартной форме |
| PDSYEV2 - PZHEEV2 | вычисление всех собственных значений и собственных векторов симметричной (эрмитовой) матрицы в линейной проблеме собственных значений |
| PDPOTRF - PZPOTRF | приведение симметричной (эрмитовой) матрицы к трехдиагональной форме методом отражений |
| PDLAMCH - | вычисление машинных параметров для арифметики с плавающей запятой |
Замечания по использованию
| 1. | В подпрограммах PSSYGV2, PCHEGV2, PZHEGV2 параметры A, B, Z и WORK имеют тип REAL, COMPLEX и DOUBLE COMPLEX соответственно, а параметр W - REAL, REAL и DOUBLE PRECISION соответственно | |
| 2. |
Список параметров подпрограммы PZHEGV2 имеет следующий вид: PZHEGV2 (IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB, W, Z, IZ, JZ, DESCZ, WORK, LWORK, RWORK, LRWORK, IFAIL, 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 ( из пакета BLACS), DSCAL ( из пакета BLAS), LSAME, INDXG2P, NUMROC, CHK1MAT, PCHK1MAT, PCHK2MAT, PXERBLA ( из библиотеки ScaLAPACK_TOOLS) |
Вычисление всех собственных значений и собственных векторов в обобщенной проблеме собственных значений (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.
Фрагмент фортранного текста вызывающей программы.
(полный текст теста можно получить в tdsygv21.zip,
а принятые обозначения - посмотреть в разделе документации
"Обозначения и упрощения в примерах
по использованию подпрограмм комплекса").
Матрицы A и B инициализируются посредством обращения к подпрограмме PDELSET из библиотеки ScaLAPACK_TOOLS.
PROGRAM TDSYGV21
include 'mpif.h'
INTEGER LWORK, MAXN, LDA, MAXPROCS, NOUT
PARAMETER ( LDA = 100, LWORK = 146, MAXN = 100,
$ 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, IFAIL
*
INTEGER DESCA( 9 ), DESCB(9), DESCZ(9)
DOUBLE PRECISION A( LDA, LDA ), B( LDA, LDA ), Z( LDA, LDA ), W( MAXN ),
$ WORK( LWORK )
*
EXTERNAL BLACS_EXIT, BLACS_GET, BLACS_GRIDEXIT,
$ BLACS_GRIDINFO, BLACS_GRIDINIT, BLACS_PINFO,
$ BLACS_SETUP, DESCINIT, PDMATINIT, PDLAPRNT,
$ PDSYGV2
*
N = 5
NB = 2
NPROW = 2
NPCOL = 2
IA = 1
JA = 1
IB = 1
JB = 1
IZ = 1
JZ = 1
IBTYPE = 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( 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 PDSYGV2( IBTYPE, UPLO, N, A, IA, JA, DESCA, B, IB, JB, DESCB,
W, Z, IZ, JZ, DESCZ, WORK, LWORK, IFAIL, 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
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
Собственные значения:
W(1) = 0.432787211016963
W(2) = 0.663662748392315
W(3) = 0.943859004668386
W(4) = 1.10928454001752
W(5) = 1.49235323254300
Собственные векторы:
Z(1, 1) = 0.134590573961348964D+00
Z(2, 1) = -0.612947224715858080D-01
Z(3, 1) = -0.157902562211254954D+00
Z(4, 1) = 0.109465787723999147D+00
Z(5, 1) = -0.414730117966489037D-01
Z(1, 2) = 0.829198064861523532D-01
Z(2, 2) = 0.153148395665945264D+00
Z(3, 2) = -0.118603667911389210D+00
Z(4, 2) = -0.182813041785799935D+00
Z(5, 2) = 0.356172036818147315D-02
Z(1, 3) = 0.191710031573882250D+00
Z(2, 3) = -0.158991211513700409D+00
Z(3, 3) = 0.748390709386720365D-01
Z(4, 3) = -0.137468929466840417D+00
Z(5, 3) = 0.889778923495063573D-01
Z(1, 4) = 0.142011959884590616D+00
Z(2, 4) = 0.142419950546665203D+00
Z(3, 4) = 0.120997623004495347D+00
Z(4, 4) = 0.125531015188065975D+00
Z(5, 4) = 0.769220728307276440D-02
Z(1, 5) = -0.763867178777808409D-01
Z(2, 5) = 0.170980018713419862D-01
Z(3, 5) = -0.666645336709070502D-01
Z(4, 5) = 0.860480093056616990D-01
Z(5, 5) = 0.289433414168931813D+00