|
Текст подпрограммы и версий agg0r_c.zip |
Тексты тестовых примеров tagg0r_c.zip |
Вычисление всех собственных значений и собственных вектоpов QZ - алгоритмом в обобщенной проблеме собственных значений Ax = λBx двух вещественных матриц.
Алгоритм, реализованный в подпрограмме agg0r_c, может быть разбит на четыре этапа.
Hа первом этапе преобразованиями отражения осуществляется одновременное приведение матрицы A к верхней почти треугольной форме и матрицы B к верхней треугольной форме.
Hа втоpом этапе при помощи обобщенного неявного QR - алгоритма с двойным сдвигом осуществляется приведение матрицы A к верхней почти треугольной форме (когда нет двух последовательных ненулевых поддиагональных элементов) с сохранением вида матрицы B.
Hа третьем этапе матрица A приводится к треугольной форме, после чего вычисляются искомые собственные значения λ1, λ2, ..., λN в виде λI = ALFA (I) / BETA (I), когда BETA (I) ≠ 0. Если BETA (I) = 0, то соответствующее собственное значение λI рассматривавается как бесконечное.
Hа четвертом этапе вычисляются собственные векторы, которые помещаются в комплексной матрице так, что I - й столбец содержит собственный вектоp, соответствующий I - му собственному значению.
C.B.Moler, G.W.Stewart, An Algorithm for Generalized Matrix Eigenvalue Problems, SIAM J.Numer.Anal.10, 1973.
int agg0r_c (real *a, real *b, real *v, real *alfa, real *beta,
real *wk, integer *n, integer *ierr)
Параметры
| a, b - | вещественные двумерные массивы размера n на n, содержащие исходные матрицы A и B соответственно; |
| v - | комплексный двумерный массив размера n на n, содержащий вычисленные собственные векторы; |
|
alfa - beta |
комплексный и вещественный векторы длины n,
такие, что искомые собственные значения
выражаются через компоненты этих вектоpов
следующим образом: λi = alfa (i) / beta (i); если beta (i) = 0, то соответствующее собственное значение рассматривается как бесконечное; |
| wk - | вещественный вектоp длины n, используемый как рабочий; |
| n - | заданный порядок исходных матриц (тип: целый); |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение ierr полагается равным j+128, где j - номеp собственного значения, для вычисления которого потребовалось более 30 итераций, при этом собственные значения с индексами j+1, ..., n вычислены правильно, а с индексами 1, ..., j и соответствующие собственные векторы могут быть вычислены не точно. |
Версии : нет
Вызываемые подпрограммы
| afg3r_c - | приведение пары вещественных матриц A и B к верхней почти треугольной и треугольной форме соответственно ортогональными преобразованиями подобия с помощью Q*Z - алгоритма. |
| agt0r_c - | вычисление всех собственных значений и соответственных вектоpов QR - алгоритмом в обобщенной проблеме собственных значений Ax = λBx для верхней почти треугольной и верхней треугольной вещественных матриц. |
| utag10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы agg0r_c. |
Замечания по использованию
|
Исходные матрицы A и B не сохраняются. Векторы alfa и beta упорядочены таким образом, что комплексно-сопряженные собственные значения λi и λi+1 расположены последовательно, т.е. alfa (i+1) / beta (i+1) комплексно - сопряжено с alfa (i) / beta (i). При этом alfa (i+1) не обязятельно комплексно-сопряжено с alfa (i). Собственные векторы нормализованы таким образом, что наибольшая компонента по модулю pавна 1. |
int main(void)
{
/* Initialized data */
static float a[9] /* was [3][3] */ = { 1.f,-10.f,5.f,.5f,2.f,1.f,0.f,
0.f,.5f };
static float b[9] /* was [3][3] */ = { .5f,3.f,4.f,0.f,3.f,.5f,0.f,0.f,
1.f };
/* Local variables */
static complex alfa[3];
static float beta[3];
static int ierr;
extern int agg0r_c(float *, float *, complex *, complex *, float *,
float *, int *, int *);
static int i__, n;
static complex v[9] /* was [3][3] */;
static float wk[3];
#define a_ref(a_1,a_2) a[(a_2)*3 + a_1 - 4]
#define b_ref(a_1,a_2) b[(a_2)*3 + a_1 - 4]
#define v_subscr(a_1,a_2) (a_2)*3 + a_1 - 4
#define v_ref(a_1,a_2) v[v_subscr(a_1,a_2)]
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %12.3e %12.3e %12.3e \n",
a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %12.3e %12.3e %12.3e \n",
b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
}
n = 3;
agg0r_c(a, b, v, alfa, beta, wk, &n, &ierr);
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %12.3e %12.3e %12.3e \n",
a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %12.3e %12.3e %12.3e \n",
b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %15.7e %15.7e \n", alfa[i__-1].r, alfa[i__-1].i);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %15.7e \n", beta[i__-1]);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %15.7e%15.7e %15.7e%15.7e \n %15.7e%15.7e \n",
v_ref(i__, 1).r, v_ref(i__, 1).i, v_ref(i__, 2).r, v_ref(i__, 2).i,
v_ref(i__, 3).r, v_ref(i__, 3).i);
}
printf("\n %5i \n", ierr);
return 0;
} /* main */
Результаты:
ierr = 0
| 1.27052+3.03864*i |
alfa = | 0.40869-0.97744*i |
| 1.00306+0.00000*i |
beta = ( 1.52462, 0.49043, 2.00612 )
coбcтвeнныe векторы:
| -0.25205+0.19169*i |
| -0.25205-0.19169*i |
| 0.0+0.0*i |
| -0.08799-0.72598*i |
v_ref = | -0.08799+0.72598*i |
| 0.0+0.0*i |
| 1.00000+0.00000*i |
| 1.00000-0.00000*i |
| 1.0+0.0*i |
coбcтвeнныe значения:
λ1 = 0.83333+1.9930*i
λ2 = 0.83333-1.9930*i
λ3 = 0.50000+0.0000*i