Текст подпрограммы и версий
afg7r_c.zip  afg7d_c.zip  afg7c_c.zip  afg7p_c.zip 		Тексты тестовых примеров
tafg7r_c.zip  tafg7d_c.zip  tafg7c_c.zip  tafg7p_c.zip

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

Назначение

Приведение вещественной матрицы к верхней форме Хессенберга методом отражений.

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

Подпрограмма afg7r_c осуществляет приведение вещественной матрицы  А порядка  N к верхней форме Хессенберга  Н = Р А P, где  Р - произведение матриц отражения.

Подпрограмма afg7r_c требует задания чисел low и IGН, удовлетворяющих условию 

 ai j  =  0 ,   ecли  i > j   и  ecли  1 ≤ j < low  или  igh < i ≤ N ,

которое означает, что первые (LОW - 1) столбцов и последние (N - IGН) строк матрицы  А имеют верхнюю треугольную форму.

Тогда достаточно будет привести к верхней форме Хессенберга только подматрицу матрицы  А, расположенную в строках и столбцах с номерами от low до igh.

Приведение к форме Хессенберга осуществляется с помощью следующей последовательности преобразований отражения 

     Ai  =  Pi Ai -1 Pi  ,      i = LOW, low +1, ..., igh - 2 ,
 где
   ALOW - 1  =  A  ,
   Pi  =  I - ui uiH / hi  ,      hi  =  uiH ui / 2  ,

а векторы  ui выбираются таким образом, чтобы первые  i столбцов матрицы  Аi имели бы форму Хессенберга.
При этом у векторов  ui отличными от нуля являются только компоненты с номерами  i + 1, i + 2, ..., IGН.

Информация о выполненных преобразованиях, а именно векторы  ui, порождающие матрицы отражения  Рi, запоминается и впоследствии может быть использована для восстановления собственных векторов исходной матрицы.

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

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

    int afg7r_c (integer *nm, integer *n, integer *low,
                 integer *igh, real *a, real *ort)

   Параметры
nm - число строк двумерного массива  a, указанное при описании этого массива в вызывающей подпрограмме (тип: целый);
n - порядок исходной матрицы,  n ≤ nm (тип: целый);
      low -
      igh   		заданные граничные индексы строк и столбцов подматрицы исходной матрицы, которую требуется привести к форме Хессенберга (тип: целый); если матрица масштабировалась, то low, igh - выходные параметры подпрограммы amb1r_c , в общем случае можно взять low = 1,  igh = n;
a - вещественный двумерный массив размерности nm * n, содержащий в своих первых  n строках исходную матрицу; в результате работы подпрограммы  a содержит вычисленную матрицу Хессенберга, а в остальной части массива  a в столбцах с номерами  i = low, low + 1, ..., igh - 2 запоминаются ненулевые компоненты векторов  ui (см. математическое описание), начиная со 2 - ой ненулевой, т.е. с ( i + 2) - ой компоненты;
ort - вещественный вектор длины igh, содержащий в результате работы подпрограммы оставшуюся часть информации о векторах  ui, при этом компонента ort ( i ) содержит  i - ую компоненту вектора  ui - 1,  i = low + 1, ..., igh - 1, т.е. первую ненулевую компоненту вектора  ui - 1;
подпрограмма использует только компоненты вектора ort с номерами low + 1, low + 2, ..., igh. 
   Версии
afg7d_c - приведение к верхней форме Хессенберга методом отражений вещественной матрицы, заданной с удвоенной точностью.
afg7c_c - приведение к верхней форме Хессенберга методом отражений комплексной матрицы (см. замечания по использованию).
afg7p_c - приведение к верхней форме Хессенберга методом отражений комплексной матрицы, заданной с удвоенной точностью (см. замечания по использованию). 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1. 

В подпрограмме afg7d_c параметры a, ort имеют тип double.

  2. 

В подпрограмме afg7c_c исходная комплексная матрица порядка  n задается в двух вещественных массивах ar и ai размеров nm * n, содержащих в своих первых  n строках ее вещественные и мнимые части соответственно.
Параметр nm определяет число строк двумерных массивов ar и ai, указанное при описании этих массивов в вызывающей подпрограмме.
В результате работы подпрограммы afg7c_c массивы ar и ai содержат соответственно вещественную и мнимую части вычисленной матрицы Хессенберга, а в столбцах с номерами  i = low, low + 1, ..., igh - 2 под поддиагональными элементами - соответственно вещественные и мнимые части векторов  ui, начиная с ( i + 2) - ой компоненты.
Для запоминания первых ненулевых компонент векторов  ui используются два вещественных вектора ortr и orti длины igh, при этом  i - ые компоненты этих векторов содержат соответственно вещественную и мнимую части  i - ой компоненты вектора  ui - 1.

Первый оператор подпрограммы afg7c_c имеет вид: 

int afg7c_c(integer *nm, integer *n, integer *low,
            integer *igh, real *ar, real *ai, real *ortr, real *orti)
  3.  Подпрограмма afg7p_c имеет такие же параметры, как и подпрограмма afg7c_c, только при этом ar, ai, ortr и orti имеют тип double.

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

int main(void)
{
    /* Initialized data */
    static float a[16] /* was [4][4] */ = { 2.f,0.f,3.f,4.f,1.f,1.f,-.6f,
                             -.8f,1.f,-.6f,1.64f,-.48f,1.f,-.8f,-.48f,1.36f };
    /* System generated locals */
    int i__1, i__2;

    /* Local variables */
    extern int afg7r_c(int *, int *, int *, int *, float *, float *);
    static int i__, n;
    static float ort[4];

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

    n = 4;
    ort[0] = 0.f;
    afg7r_c(&c__4, &c__4, &c__1, &c__4, a, ort);

    for (i__ = 1; i__ <= 4; ++i__) {
        printf("\n  %15.7e %15.7e %15.7e %15.7e \n",
             a_ref(i__, 1), a_ref(i__, 2), a_ref(i__, 3), a_ref(i__, 4));
    }
     printf("\n  %15.7e %15.7e %15.7e %15.7e \n",
             ort[0], ort[1], ort[2], ort[3]);
    return 0;
} /* main */


Результаты:

                  |    2     - 1.4      1    - 0.22
                  |  - 5       1.        1    - 9.1e - 13
    a_ref  =  |    3       1.        1    + 9.1e - 13
                  |    4     - 0.8      0       2.

      (ort( i ),  i  =  1, 3)  =  0., 5., - 1.6

 Это означает, что

             |    2    - 1.4      1    - 0.2
             |  - 5      1         1    - 9.1e - 13
    h   =  |    0      1         1       9.1e - 13
             |    0      0         0       2.

    pi  =  i - ui uit / hi  ,    hi  =  uiT ui / 2  ,   i  =  1, 2  , 

   где

    u1T  =  ( 0., 5., 3., 4 )  ,
    u2T  =  ( 0., 0., - 1.6, - 0.8 )  .