|
Текст подпрограммы и версий agt1r_c.zip |
Тексты тестовых примеров tagt1r_c.zip |
Вычисление всех собственных значений в обобщенной проблеме Ax = λBx для вещественных верхней почти треугольной матрицы A и верхней треугольной матрицы B с помощью QR - алгоритма.
Подпрограмма agt1r_c вычисляет все собственные значения обобщенной проблемы Ax = λBx для вещественных верхней почти треугольной матрицы A размера N на N и верхней треугольной матрицы B размера N на N с помощью QR - алгоритма.
Информация о вычисленных собственных значениях обобщенной проблемы выдается в вектоpах ALFR, ALFI, BETA длины N, по которым I - ое собственное значение определяется с помощью формулы:
λI = ALFR(I) / BETA(I) + i ALFI(I) / BETA(I), при BETA(I) ≠ 0,
(1) λI = ∞ , при BETA(K) = 0,
λI = любое число , при ALFR(I) = ALFI(I) =
= BETA(I) = 0.
C.B.Moler, G.W.Stewart, An Algorithm for Generalized Matrix Eigenvalue Problems, SIAM J. Numer. Anal., 10, 1973.
int agt1r_c (real *a, real *b, real *alfr, real *alfi,
real *beta, integer *n, integer *ierr)
Параметры
| a, b - | вещественные двумерные массивы размера n на n, содержащие соответственно верхнюю почти треугольную и верхнюю треугольную матрицы; |
|
alfr - alfi beta | вещественные векторы длины n, содержащие информацию о собственных значениях λk обобщенной проблемы; при этом собственные значения определяются формулой (1); |
| n - | заданный порядок исходных матриц A, B (тип: целый); |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение ierr полагается равным j+128, где j - номеp собственного значения, для вычисления которого потребовалось более 30 итераций; при этом собственные значения с индексами j+1, ..., n вычислены правильно, а с индексами 1, ..., j могут быть вычислены не точно. |
Версии : нет
Вызываемые подпрограммы
| utag10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы agt1r_c. |
Замечания по использованию
| Подпрограмма agt1r_c не сохраняет исходную информацию. |
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 float alfi[3], beta[3], alfr[3];
static int ierr;
extern int afg3r_c(float *, float *, float *, int *, int *),
agt1r_c(float *, float *, float *, float *, float *,
int *, int *);
static int i__, m, n;
#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]
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %13.3e %13.3e %13.3e \n",
a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %13.3e %13.3e %13.3e \n",
b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
}
n = 3;
m = 1;
afg3r_c(a, b, a, &n, &m);
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %13.3e %13.3e %13.3e \n",
a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %13.3e %13.3e %13.3e \n",
b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
}
agt1r_c(a, b, alfr, alfi, beta, &n, &ierr);
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %13.3e %13.3e %13.3e \n",
a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %13.3e %13.3e %13.3e \n",
b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %14.5e \n", alfr[i__-1]);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %14.5e \n", alfi[i__-1]);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %14.5e \n", beta[i__-1]);
}
printf("\n %5i \n", ierr);
return 0;
} /* main */
Результаты:
ierr = 0
alfr = (1.27052, 0.40869, 1.00306)
alfi = (3.03864, -0.97744, 0.00000)
beta = (1.52462, 0.49043, 2.00612)
Собственные значения
λk = (alfr(k) + i alfi(k)) / beta(k), k = 1, 2, 3
λ1 = 0.83333 + 1.9930i
λ2 = 0.83333 - 1.9930i
λ3 = 0.50000 + 0.0000i