|
Текст подпрограммы и версий agt2c_c.zip |
Тексты тестовых примеров tagt2c_c.zip |
Вычисление всех собственных значений и собственных вектоpов в обобщенной поблеме AV - 1x = λBV - 1x для комплексных верхней почти треугольной матрицы A, верхней треугольной матрицы B и произвольной невырожденной матрицы V, с помощью LR - алгоритма.
Подпрограма agt2c_c вычисляет все собственные значения и собственые векторы обобщенной проблемы
AV - 1 = λBV - 1x
для комплексных верхней почти треугольной матрицы A, верхней треугольной матрицы B и произвольной невырожденной матрицы V с помощью LR - алгоритма. Матрицы A, B и V имеют размеры N на N.
Информация о вычисленных собственных значениях выдается в вектоpах ALFA и BETA длины N, по которым k - ое собственное значение λk определяется с помощью формулы:
λk = ALFA(k) / BETA(k), при BETA(k) ≠ 0,
(1) λk = ∞ , при BETA(k) = 0,
λk = любое число , при ALFA(k) = BETA(k) = 0.
Вычисленные собственные векторы запоминаются в столбцах матрицы V.
R.S.Mfrtin, J.H.Wilkinson, The Modified LR - Algorithm for Complex Hessenberg Matrices, Numer. Math., 12, 1968.
C.B.Moler, G.W.Stewart, An Algorithm for the Generalized Matrix Eigenvalue Problems, SIAM J. Numer. Anal., 10, 1973.
int agt2c_c (complex *a, complex *b, complex *v, complex *alfa,
complex *beta, integer *n, integer *ierr)
Параметры
| a, b - | комплексные двумерные массивы размера n на n, содержащие соответственно верхнюю почти треугольную и верхнюю треугольную матрицы; |
| v - | комплексный двумерный массив размера n на n, содержащий заданную невырожденную матрицу V; в результате работы подпрограммы в столбцах массива v содержатся вычисленные собственные векторы xk; |
|
alfa - beta | комплексные векторы длины n, содержащие информацию о собственных значениях λk обобщенной проблемы; при этом собственные значения определяются формулой (1); |
| n - | заданный порядок исходных матриц A, B и V (тип: целый); |
| ierr - | целая переменная, служащая для сообщения об ошибках обнаруженных в ходе работы подпрограммы; значение ierr полагается равным номеpу собственного значения для вычисления которого потребовалось более 30 итераций; при этом собственные значения с индексами ierr+1, ..., n вычислены правильно, а с индексами 1, ..., ierr и соответствующие собственные векторы не вычисляются. |
Версии : нет
Вызываемые подпрограммы
| utag10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы agt2c_c. |
Замечания по использованию
| 1. |
Подпрограмма agt2c_c не сохраняет исходную информацию. | |
| 2. |
B массиве v должна задаваться единичная матрица, если решается обобщенная проблема Ax = λBx для комплексных верхней почти треугольной матрицы A и верхней треугольной матрицы B. | |
| 3. | Вычисленные собственные векторы нормализованы так, что наибольшая компонента по модулю pавна 1. |
int main(void)
{
/* Initialized data */
static complex a[9] /* was [3][3] */ = { {1.f,0.f},{-10.f,0.f},{5.f,0.f},
{.5f,0.f},{2.f,0.f},{1.f,0.f},{0.f,0.f},{0.f,0.f},{.5f,0.f} };
static complex b[9] /* was [3][3] */ = { {.5f,0.f},{3.f,0.f},{4.f,0.f},
{0.f,0.f},{3.f,0.f},{.5f,0.f},{0.f,0.f},{0.f,0.f},{1.f,0.f} };
/* System generated locals */
int i__1;
complex q__1;
/* Builtin functions */
void c_div(complex *, complex *, complex *);
/* Local variables */
static complex alfa[3], beta[3];
static int ierr;
extern int afg3c_c(complex *, complex *, complex *, int *, int *),
agt2c_c(complex *, complex *, complex *, complex *,
complex *, int *, int *);
static int i__, k, m, n;
static complex v[9] /* was [3][3] */, lambda[3];
#define a_subscr(a_1,a_2) (a_2)*3 + a_1 - 4
#define a_ref(a_1,a_2) a[a_subscr(a_1,a_2)]
#define b_subscr(a_1,a_2) (a_2)*3 + a_1 - 4
#define b_ref(a_1,a_2) b[b_subscr(a_1,a_2)]
#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 %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f \n",
a_ref(i__, 1).r, a_ref(i__, 1).i, a_ref(i__, 2).r, a_ref(i__, 2).i,
a_ref(i__, 3).r, a_ref(i__, 3).i);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f \n",
b_ref(i__, 1).r, b_ref(i__, 1).i, b_ref(i__, 2).r, b_ref(i__, 2).i,
b_ref(i__, 3).r, b_ref(i__, 3).i);
}
n = 3;
m = 0;
afg3c_c(a, b, v, &n, &m);
agt2c_c(a, b, v, alfa, beta, &n, &ierr);
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f \n",
a_ref(i__, 1).r, a_ref(i__, 1).i, a_ref(i__, 2).r, a_ref(i__, 2).i,
a_ref(i__, 3).r, a_ref(i__, 3).i);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f \n",
b_ref(i__, 1).r, b_ref(i__, 1).i, b_ref(i__, 2).r, b_ref(i__, 2).i,
b_ref(i__, 3).r, b_ref(i__, 3).i);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %14.5f %14.5f \n", alfa[i__-1].r, alfa[i__-1].i);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %14.5f %14.5f \n", beta[i__-1].r, beta[i__-1].i);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f \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);
for (k = 1; k <= 3; ++k) {
i__1 = k - 1;
if (beta[i__1].r == 0.f && beta[i__1].i == 0.f) {
goto l11;
}
i__1 = k - 1;
c_div(&q__1, &alfa[k - 1], &beta[k - 1]);
lambda[i__1].r = q__1.r, lambda[i__1].i = q__1.i;
goto l10;
l11:
i__1 = k - 1;
lambda[i__1].r = 3.4e38f, lambda[i__1].i = 3.4e38f;
l10:
;
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %14.5f %14.5f \n", lambda[i__-1].r, lambda[i__-1].i);
}
return 0;
} /* main */
Результаты:
ierr = 0
| 2.79963 - 3.65018i |
alfa = | 0.42142 + 0.54945i |
| -1.09877 + 0.0i |
| -1.05899 - 1.89748i |
beta = | -0.15941 + 0.27809i |
| -2.19753 + 0.0i |
| -0.25205 + 0.19169i |
| -0.25205 - 0.19169i |
| 0.0 + 0.0i |
| -0.08799 - 0.72598i |
v_ref = | -0.08799 - 0.72598i |
| 0.0 + 0.0i |
| 1.00000 + 0.0i |
| 1.00000 + 0.0i |
| 1.0 + 0.0i |
Собственные значения λk = alfa(k) / beta(k)
λ1 = 0.83333 + 1.9930i
λ2 = 0.83331 - 1.9931i
λ3 = 0.50000 + 0.0i
Собственные векторы
| -0.25205 + 0.19169i |
x1 = | -0.08799 - 0.72598i |
| 1.00000 + 0.0i |
| -0.25205 - 0.19169i |
x2 = | -0.08799 - 0.72598i |
| 1.00000 + 0.0i |
| 0.0 + 0.0i |
x3 = | 0.0 + 0.0i |
| 1.0 + 0.0i |