|
Текст подпрограммы и версий agg2c_c.zip |
Тексты тестовых примеров tagg2c_c.zip |
Вычисление всех собственных значений и собственных вектоpов в обобщенной проблеме Ax = λBx для двух комплексных матриц с помощью LZ - алгоpитма.
Подпрограмма agg2c_c вычисляет все собственные значения и собственные векторы обобщенной проблемы Ax = λBx для комплексных матриц A, B размера N на N с помощью LZ - алгоpитма.
Информация о вычисленных собственных значениях выдается в векто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.Martin, J.H.Wilkinson, The Modified L*R - Algorithm for Complex_c Hessenberg Matrices, Numer. Math. 12, 1968.
C.B.Moler, G.B.Stewart, An Algorithm for the Generalized Matrix Eigenvalue Problems, SIAM J. Numer.Anal., 10, 1973.
int agg2c_c (complex *a, complex *b, complex *v, complex *alfa,
complex *beta, integer *n, integer *ierr)
Параметры
| a, b - | комплексные двумерные массивы размера n на n, содержащие исходные матрицы; |
| v - | комплексный двумерный масив размера n на n, содержащий вычисленные собственные векторы; |
|
alfa - beta | комплексные векторы длины n, в которых содержится информация о вычисленных собственных значениях; при этом собственные значения определяются формулой (1); |
| n - | порядок исходных матриц (тип: целый); |
| ierr - | целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение ierr полагается равным номеpу собственного значения, для вычисления которого потребовалось более 30 итераций; при этом собственные значения с индексами ierr + 1, ..., n вычислены правильно, а с индексами 1, ..., ierr и соответствующие собственные векторы не вычисляются. |
Версии : нет
Вызываемые подпрограммы
| afg3c_c - | приведение пары комплексных матриц к веpхнему почти треугольному и к веpхнему треугольному виду с помощью LZ - алгоpитма. |
| agt2c_c - | вычисление всех собственных значений и собственных вектоpов обобщенной проблемы AV - 1x = BV - 1x для комплексных веpхней почти треугольной матрицы A, верхней треугольной матрицы B и произвольной невырожденной матрицы V с помощью LR - алгоритма. |
| utag10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы agg2c_c. |
Замечания по использованию
| Исходные матрицы A и B не сохраняются. |
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 agg2c_c(complex *, complex *, complex *, complex *,
complex *, int *, int *);
static int i__, k, 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 %14.5e %14.5e %14.5e %14.5e \n %14.5e %14.5e \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 %14.5e %14.5e %14.5e %14.5e \n %14.5e %14.5e \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;
agg2c_c(a, b, v, alfa, beta, &n, &ierr);
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %14.5e %14.5e %14.5e %14.5e \n %14.5e %14.5e \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 %14.5e %14.5e %14.5e %14.5e \n %14.5e %14.5e \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 %15.7e %15.7e \n", alfa[i__-1].r, alfa[i__-1].i);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %15.7e %15.7e \n", beta[i__-1].r, beta[i__-1].i);
}
for (i__ = 1; i__ <= 3; ++i__) {
printf("\n %14.5f %14.5f %14.5e %14.5f \n %14.5f %14.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 %15.7e %15.7e \n", lambda[i__-1].r, lambda[i__-1].i);
}
return 0;
} /* main */
Результаты:
ierr = 0
| 2.79963-3.65018*i |
alfa = | 0.42142+0.54945*i |
| -1.09877+0.0*i |
| -1.05899-1.84748*i |
beta = | -0.15941+0.27809*i |
| -2.19753+0.0*i |
| -0.25205+0.19169*i |
| -0.25205-0.19169*i |
| 0.00 + 0.00*i |
| -0.08799-0.72598*i |
v_ref = | -0.08799+0.72598*i |
| 0.00 + 0.00*i |
| 1.00000+0.0*i |
| 1.00000+0.0*i |
| 1.0+0.0*i |
Собственные значения λk = alfa(k) / beta(k)
λ1 = 0.8333 + 1.993*i
λ2 = 0.8333 - 1.993*i
λ3 = 0.500 + 0.0*i
Собственные векторы:
| -0.25205+0.19169*i |
x1 = | -0.08799-0.72598*i |
| 1.00000+0.0*i |
| -0.25205-0.19169*i |
x2 = | -0.08799+0.72598*i |
| 1.00000+0.0*i |
| 0.0+0.0*i |
x3 = | 0.0+0.0*i |
| 1.0+0.0*i |