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

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

Назначение

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

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

Обобщенная проблема Ax = λBx с якобиевой матрицей А и диагональной положительной матрицей B последовательностью преобразований сводится к стандартной задаче Qz = λz. Здесь Q - симметричная трехдиагональная матрица, coбственныe значeния которой совпадают с собственными значениями исходной обобщенной проблемы, а собственные векторы  z  связаны с собственными векторами  x  обобщенной проблемы соотношением x = Cz, где матрица C определяется последовательностью выполненных преобразований.

Для решения стандартной проблемы применяется неявный QL - алгоритм со сдвигом. Собственные векторы  z  нормированы таким образом, что  zTz = E,  где E - единичная матрица.

Дж.Х. Уилкинсон, Алгебраическая проблема собственных значений, "Hаука", M., 1970.

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

    int agj0r_c (real *a, real *b, real *v, real *ev, real *rab,
            integer *n, integer *ierr)

   Параметры
a - вещественный двумерный массив размера n на 3, содержащий якобиеву матрицу A в компактной форме записи для ленточных матриц;
b - вещественный вектоp длины n, содержащий элементы диагональной матрицы B;
v - вещественный двумерный массив размера n на n, в столбцах которого содержатся вычисленные оpтонормированные собственные векторы;
ev - вещественный вектоp длины n, содержащий вычисленные в возрастающем порядке собственные значения;
rab - вещественный вектоp длины 2 на n, используемый как рабочий;
n - порядок исходных матриц (тип: целый);
ierr - целая переменная служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение ierr полагается равным:
 

номеpу собственного значения, для вычисления которого потребовалось более 30 итераций, при этом собственные значения с индексами 1, 2, ..., ierr - 1 и соответствующие собственные векторы вычислены правильно, но собственные значения неупорядочены;

ierr = 65 - когда исходная диагональная матрица B не является положительно определенной;

ierr = 66 - когда в исходной матрице Якоби не все попарные произведения соответствующих произведений побочных диагоналей неотрицательны;

ierr = 67 - когда в матрице Якоби имеется pавное нулю произведение соответствующих элементов побочных диагоналей, причем сомножители равны нулю неодновpеменно; в этом случае нет симметризации, необходимой для правильного вычисления собственных вектоpов;

ierr = 68 - когда матрица Q получена, но матрица преобразования C не вычислена, так как она алгоритмически вырождена; в этом случае собственные векторы не вычисляются. 
   Версии :   нет 
   Вызываемые подпрограммы
afj0r_c - приведение матрицы Якоби к симметрической трехдиагональной матрице преобразованием подобия.
afe1r_c - одновременное приведение вещественной симметрической трехдиагональной матрицы и диагональной положительно определенной матрицы двусторонним диагональным преобразованием соответственно к вещественной симметрической трехдиагональной и единичной матрицам.
utag10_c -
utag11_c   		подпрограммы выдачи диагностических сообщений при работе подпрограммы agj0r_c. 
   Замечания по использованию
  Подпрограмма agj0r_c сохраняет матрицу A.

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

int main(void)
{
    /* Initialized data */
    static float a[15]   /* was [5][3] */ = { 0.f,1.f,1.f,1.f,1.f,0.f,0.f,0.f,
                                              0.f,0.f,4.f,4.f,4.f,4.f,0.f };
    static float b[5] = { 2.f,2.f,2.f,2.f,2.f };

    /* Local variables */
    static int ierr;
    extern int agj0r_c(float *, float *, float *, float *, float *,
                       int *, int *);
    static int i__, n;
    static float v[25] /* was [5][5] */, ev[5], rab[10];

#define a_ref(a_1,a_2) a[(a_2)*5 + a_1 - 6]
#define v_ref(a_1,a_2) v[(a_2)*5 + a_1 - 6]

    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e %15.7e %15.7e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e \n", b[i__-1]);
    }
    n = 5;
    agj0r_c(a, b, v, ev, rab, &n, &ierr);

    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e %15.7e %15.7e \n",
                a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3));
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %15.7e \n", b[i__-1]);
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %14.5e \n", ev[i__-1]);
    }
    for (i__ = 1; i__ <= 5; ++i__) {
         printf("\n %14.5e %14.5e %14.5e %14.5e %14.5e \n",
                v_ref(i__, 1), v_ref(i__, 2), v_ref(i__, 3),
                v_ref(i__, 4), v_ref(i__, 5));
    }
    printf("\n %5i \n", ierr);
    return 0;
} /* main */


Результаты:

                   | -1.7320508 |                       |  0.2041241 |
                   | -1.0000000 |                       | -0.1767767 |
      ev    =  |  0.0000000 |    ,        v1  =  |  0.1020621 |    , 
                   |  1.0000000 |                        | -0.0441942 |
                   |  1.7320508 |                        |  0.0127578 |

                 | -0.3535534 |                       | -0.4082483 |
                 |  0.1767767 |                        | -0.0000000 |
      v2  =  |  0.0000000  |    ,        v3  =  |  0.1020621 |    , 
                 | -0.0441942 |                        |  0.0000000 |
                 |  0.0220971 |                        | -0.0255155 |

                 | -0.3535534 |                       | -0.2041241 |
                 | -0.1767767 |                       | -0.1767767 |
      v4  =  | -0.0000000 |    ,        v5  =  | -0.1020621 |    , 
                 |  0.0441942 |                       | -0.0441942 |
                 |  0.0220971 |                       | -0.0127578 |

      ierr  =  0.