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

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

Назначение

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

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

Алгоритм, реализованый в подпрограмме agg1r_c, может быть разбит на три этапа.

Hа первом этапе преобразованиями отражения осуществляется одновременное приведение матрицы A к верхней почти треугольной форме и матрицы B к верхней треугольной форме.

Hа втоpом этапе при помощи обобщенного неявного QR - алгоритма с двойным сдвигом осуществляется приведение матрицы A к верхней почти треугольной форме (когда нет двух последовательных ненулевых поддиагональных элементов) с сохранением вида матрицы B.

Hа третьем этапе матрица A приводится к треугольной форме, после чего вычисляются искомые собственные значения λ1, λ2, ..., λN в виде λI = ALFA (I) / BETA (I), когда BETA (I) ≠ 0. Если BETA (I) = 0, то соответствующее собственное значение λI рассматривается как бесконечное.

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

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

    int agg1r_c (real *a, real *b, real *alfa, real *beta,
            real *wk, integer *n, integer *ierr)

   Параметры
a, b - вещественные двумерные массивы размера n на n, содержащие исходные матрицы A и B соответственно;
  alfa -
  beta   		комплексный и вещественный векторы длины n, такие, что искомые собственные значения выражаются через компоненты этих вектоpов следующим образом: λi = alfa (i) / beta (i); если beta (i) = 0, то соответствующее собственное значение рассматривается как бесконечное;
wk - вещественный вектоp длины n, используемый как рабочий;
n - заданный порядок исходных матриц (тип: целый);
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение ierr полагается равным j+128, где j - номеp собственного значения, для вычисления которого потребовалось более 30 итераций, при этом собственные значения с индексами j+1, ..., n вычислены правильно, а с индексами 1, ..., j могут быть вычислены не точно. 
   Версии :   нет 
   Вызываемые подпрограммы
afg3r_c - приведение пары вещественных матриц A и B к верхней почти треугольной и треугольной форме соответственно ортогональными преобразованиями подобия с помощью QZ - алгоритма.
agt1r_c - вычисление всех собственных значений QR - алгоритмом в обобщенной проблеме собственных значений Ax = λBx для верхней почти тpеугольной и верхней треугольной вещественных матриц.
utag10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы agg1r_c. 
   Замечания по использованию
 

Исходные матрицы A и B не сохраняются.

Векторы alfa и beta упорядочены таким образом, что комплексно-сопряженные собственные значения λi и λi+1 расположены последовательно, т.е. alfa (i+1) / beta (i+1) комплексно - сопряжено с alfa (i) / beta (i). При этом alfa (i+1) не обязательно комплексно - сопряжено с alfa (i).

Пример использования 

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 complex alfa[3];
    static float beta[3];
    static int ierr;
    extern int agg1r_c(float *, float *, complex *, float *, float *,
                       int *, int *);
    static int i__, n;
    static float wk[3];

#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 %12.3e %12.3e %12.3e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
    }
    n = 3;
    agg1r_c(a, b, alfa, beta, wk, &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 %12.3e %12.3e %12.3e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
    }
    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 \n", beta[i__-1]);
    }
    printf("\n %5i \n", ierr);
    return 0;
} /* main */


Результаты:

      ierr  =  0

                      | 1.27052+3.03864*i |
      alfa  =  | 0.40869-0.97744*i |
                      | 1.00306+0.00000*i | 

      beta  =  ( 1.52462,  0.49043,  2.00612 )

coбcтвeнныe значения: 

      λ1  =  0.83333+1.9930*i
      λ2  =  0.83333-1.9930*i
      λ3  =  0.50000+0.0000*i