Текст подпрограммы и версий
afh4r_c.zip , afh4d_c.zip 		Тексты тестовых примеров
tafh4r_c.zip , tafh4d_c.zip

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

Назначение

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

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

Подпрограмма afh4r_c реализует алгоритм приведения симметрической матрицы А ортогональными подобными преобразованиями к симметрической трехдиагональной форме F = QTАQ, где Q - пpоизведение оpтогональных симметpических матpиц отpажения.

Уилкинсон, Райнш. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра, М., "Машиностроение", 1976.

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

    int afh4r_c (integer *n, real *a, real *d, real *e)

   Параметры
n - порядок исходной матрицы (тип: целый);
a - вещественный двумерный массив размера n*n, содержащий симметрическую матрицу; в результате работы подпрограммы в поддиагональной части массива a в последовательных строках запоминаются векторы, порождающие матрицы отражения Q1, ..., Qn - 2;
d - вещественный вектор длины n, содержащий диагональные элементы трехдиагональной матрицы;
e - вещественный вектор длины n, содержащий в последних n - 1 компонентах поддиагональные элементы симметрической трехдиагональной матрицы. Первый элемент e произволен. 
   Версии
afh4d_c - приведение симметрической матрицы, заданной с удвоенной точностью, к симметрической трехдиагональной матрице ортогональными преобразованиями подобия. 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
 

Исходную матрицу A можно задавать лишь нижним треугольником.

Подпрограмма afh4r_c сохраняет только полный верхний треугольник массива a, в остальных компонентах этого массива размещается информация об ортогональных преобразованиях.

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

int main(void)
{
    /* Initialized data */
    static float a[16] /* was [4][4] */ = { 1.f,.42f,.54f,.66f,.42f,1.f,.32f,
                                 .44f,.54f,.32f,1.f,.22f,.66f,.44f,.22f,1.f };
    /* Local variables */
    extern int afh4r_c(int *, float *, float *, float *);
    static int i__;
    static float ev[4], rab[4];

#define a_ref(a_1,a_2) a[(a_2)*4 + a_1 - 5]

    for (i__ = 1; i__ <= 4; ++i__) {
        printf("\n  %10.2f %10.2f %10.2f %10.2f \n",
             a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3), a_ref(i__, 4));
    }
    afh4r_c(&c__4, a, ev, rab);

    for (i__ = 1; i__ <= 4; ++i__) {
        printf("\n  %10.2f %10.2f %10.2f %10.2f \n",
             a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3), a_ref(i__, 4));
    }
    printf("\n  %14.7f %14.7f %14.7f %14.7f \n",
             ev[0], ev[1], ev[2], ev[3]);
    printf("\n  %10.3e %10.3e %10.3e %10.3e \n",
             rab[0], rab[1], rab[2], rab[3]);
    return 0;
} /* main */


Результаты:

                      | 1.00   0.42   0.54   0.66 |
      a_ref  =    | 0.18   1.00   0.32   0.44 |
                      | 0.42   0.58   1.00   0.22 |
                      | 0.66   0.44   1.04   1.00 |

                    | 0.6670 |
      d__  =    | 0.6501 |
                    | 1.6829 |
                    | 1.0000 |

                   |  0.0000 |
      ev  =     | -0.0879 |
                   | -0.4614 |
                   | -0.8232 |