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

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

Назначение

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

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

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

Информация о вычисленных собственных значениях выдается в векто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.

R.S.Martin, 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 agg3c_c (complex *a, complex *b, complex *alfa,
        complex *beta, integer *n, integer *ierr)

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

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

      λ1  =  0.8333 + 1.993i
      λ2  =  0.8333 - 1.993*i
      λ3  =  0.5000 + 0.0*i