Текст подпрограммы и версий
agt3c_c.zip  		Тексты тестовых примеров
tagt3c_c.zip 

Подпрограмма:  agt3c_c

Назначение

Вычисление всех собственных значений в обобщенной проблеме Ax = λBx для комплексных верхней почти треугольной матрицы A и верхней треугольной матрицы B с помощью LR - алгоритма.

Математическое описание

Подпрограмма agt3c_c вычисляет все собственные значения обобщенной проблемы Ax = λBx для комплексных верхней почти треугольной матрицы A размера N на N и верхней треугольной матрицы B размера N на N с помощью LR - алгоритма.

Информация о вычисленных собственных значениях обобщенной проблемы выдается в вектоpах ALFA, BETA длины N, причем k - ое собственное значение вычисляется по формуле: 

                   λk = ALFA(k) / BETA(k),   при BETA(k) ≠ 0,
  (1)            λk = ∞                             ,   при BETA(k) = 0,
                   λk = любое число          ,   при ALFA(k) = BETA(k) = 0.

R.S.Martin and J.H.Wilkinson, The Modified LR - Algorithm for Complex Hessemberg Matrices, Numer. Math., 12, 369-376, 1968.

C.B.Moler, G.W.Stewart, An Algorithm for Generalized Matrix Eigenvalue Problems, SIAM J. Numer. Anal., 10, 1973.

Использование 

    int agt3c_c (complex *a, complex *b, complex *alfa,
        complex *beta, integer *n, integer *ierr)

   Параметры
a, b - комплексные двумерные массивы размера n на n, содержащие соответственно верхнюю почти треугольную и верхнюю треугольную матрицы;
  alfa -
  beta   		комплексные векторы длины n, содержащие информацию о собственных значениях λk обобщенной проблемы; при этом собственные значения определяются формулой (1);
n - заданный порядок исходных матриц A, B (тип: целый);
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение ierr полагается равным номеpу собственного значения, для вычисления которого потребовалось более 30 итераций; при этом собственные значения с индексами ierr+1, ..., n вычислены правильно, а с индексами 1, ..., ierr не вычисляются. 
   Версии :   нет 
   Вызываемые подпрограммы
utag10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы agt3c_c. 
   Замечания по использованию
  Подпрограмма agt3c_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 afg3c_c(complex *, complex *, float *, int *, int *),
               agt3c_c(complex *, complex *, complex *, complex *,
                       int *, int *);
    static int i__, k, m, n;
    static float v;
    static complex 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)]

    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 = 1;
    afg3c_c(a, b, &v, &n, &m);

    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);
    }
    agt3c_c(a, b, 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);
    }
    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.84748i |
      beta  =  | -0.15941 + 0.27809i |
                      | -2.19753 + 0.0i        |

coбcтвeнныe знaчeния   λk = alfa(k) / beta(k)

      λ1  =  0.83333 + 1.9930i
      λ2  =  0.83331 - 1.9931i
      λ3  =  0.50000 + 0.0i