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

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

Назначение

Вычисление всех собственных значений и собственных вектоpов в обобщенной проблеме AV - 1x = λBV - 1 для вещественных верхней почти треугольной матрицы A произвольной невырожденной матрицы V и верхней треугольной матрицы B с помощью QR - алгоритма.

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

Подпрограмма agt0r_c вычисляет все собственные значения и собственные векторы обобщенной проблемы

     AV - 1x = λBV - 1x 

для вещественных верхней почти треугольной матрицы A, верхней треугольной матрицы B и произвольной невырожденной матрицы V с помощью QR - алгоритма. Матрицы A, B и V имеют размеры N на N.

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

      λk = ALFR(K) / BETA(K) + i ALFI(K) / BETA(K),  при BETA(K) ≠ 0,
(1)  λk = ∞                                                          , при BETA(K) = 0,
      λk = любое число                                        , при ALFR(K) = ALFI(K) = 
                                                                                    = BETA(K) = 0. 

Вычисленные собственные векторы Xk обобщенной проблемы помещаются на место матрицы V таким образом, что вещественному значению λk соответствует k - ый столбец Vk матрицы V, т.е. Xk = Vk, а паре комплексно - сопряженных собственных значений λk и λk+1 соответствуют столбцы Vk и Vk+1 матрицы V. При этом

     Xk = Vk + i Vk+1,    Xk+1 = Vk - i Vk+1.  

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

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

    int agt0r_c (real *a, real *b, real *v, real *alfr,
             real *alfi, real *beta, integer *n, integer *ierr)

Параметры

a, b - вещественные двумерные массивы размера n на n, содержащие соответственно верхнюю почти треугольную и верхнюю треугольную матрицы;
v - вещественный двумерный массив размера n на n, содержащий заданную невырожденную матрицу V; в результате работы подпрограммы в столбцах массива v содержатся вычисленные собственные векторы Xk таким образом, что вещественному собственному значению λk соответствует k - ый столбец vk, т.е Xk = vk, а паре комплексно - сопряженных собственных значений λk и λk+1 соответствуют два столбца vk и vk+1; при этом Xk = vk + i vk+1,  а  Xk+1 = vk - i vk+1;
  alfr -
 alfi  
  beta  
вещественные векторы длины n, содержащие информацию о собственных значениях λk обобщенной проблемы; при этом собственные значения определяются формулой (1);
n - заданный порядок исходных матриц A, B и V (тип: целый);
ierr - целая переменная, служащая для сообщения об ошибках, обнаруженых в ходе работы подпрограммы; значение ierr полагается равным j+128, где j - номеp собственного значения, для вычисления которого потребовалось более 30 итераций, при этом собственные значения с индексами j+1, ..., n вычислены правильно, а с индексами 1, ..., j и собственные векторы могут быть вычислены не точно.

Версии : нет

Вызываемые подпрограммы

utag10_c - подпрограмма выдачи диагностических сообщений при работе подпрограммы agt0r_c.

Замечания по использованию

  1. 

Подпрограмма agt0r_c не сохраняет исходную информацию.

  2. 

B массиве v должна задаваться единичная матрица, если решается обобщенная проблема Ax = λBx для вещественных верхней почти треугольной матрицы A и верхней треугольной матрицы B.

  3.  Вычисленные собственные векторы нормализованы так, что наибольшая компонента по модулю pавна 1.

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

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 float alfi[3], beta[3], alfr[3];
    static int ierr;
    extern int afg3r_c(float *, float *, float *, int *, int *),
               agt0r_c(float *, float *, float *, float *, float *, float *,
                       int *, int *);
    static int i__, m, n;
    static float v[9] /* was [3][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]
#define v_ref(a_1,a_2) v[(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 %13.3e %13.3e %13.3e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
    }
    n = 3;
    m = 0;
    afg3r_c(a, b, v, &n, &m);

    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 %13.3e %13.3e %13.3e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
    }
    for (i__ = 1; i__ <= 3; ++i__) {
         printf("\n %14.5e %14.5e %14.5e \n",
                v_ref(i__, 1), v_ref(i__, 2), v_ref(i__, 3));
    }
    agt0r_c(a, b, v, alfr, alfi, beta, &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 %13.3e %13.3e %13.3e \n",
                b_ref(i__, 1), b_ref(i__, 2), b_ref(i__, 3));
    }
    for (i__ = 1; i__ <= 3; ++i__) {
         printf("\n %14.5e \n", alfr[i__-1]);
    }
    for (i__ = 1; i__ <= 3; ++i__) {
         printf("\n %14.5e \n", alfi[i__-1]);
    }
    for (i__ = 1; i__ <= 3; ++i__) {
         printf("\n %14.5e \n", beta[i__-1]);
    }
    for (i__ = 1; i__ <= 3; ++i__) {
         printf("\n %14.5e %14.5e %14.5e \n",
                v_ref(i__, 1), v_ref(i__, 2), v_ref(i__, 3));
    }
    printf("\n %5i \n", ierr);
    return 0;
} /* main */


Результаты:

       ierr  =  0
       alfr  =   ( 1.27052,   0.40869,   1.00306)
       alfi   =   ( 3.03864,  -0.97744,   0.00000)
       beta  =   ( 1.52462,   0.49043,   2.00612)

                | -0.25205   0.19169   0.00000 |
       v  =  | -0.08799  -0.72598   0.00000 |
                |  1.00000    0.00000   1.00000 |

Собственные значения
   λk  =   (alfr(k) + i alfi(k)) / beta(k),  k  =  1, 2, 3:

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

Собственные векторы:

                | -0.25205 + 0.19169i |
      x1  =  | -0.08799 - 0.72598i |
                |  1.00000 + 0.00000i |

                 | -0.25205 - 0.19169i |
      x2  =  | -0.08799 + 0.72598i |
                 |  1.00000 - 0.00000i |

                 |  0.00000 + 0.00000i |
      x3  =  |  0.00000 + 0.00000i |
                 |  1.00000 + 0.00000i |