Текст подпрограммы и версий
afh5c_c.zip , afh5p_c.zip 		Тексты тестовых примеров
tafh5c_c.zip , tafh5p_c.zip

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

Назначение

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

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

Подпрограмма afh5c_c приводит исходную эрмитову матрицу А порядка N к вещественной симметричной матрице Т = VРАР*V*, где V - диагональная унитарная матрица, а Р - унитарная матрица.
Матрица Р имеет вид: 

        P  =  PN-1 * PN-2* ...* P2 * P1 ,
где 
       Pi  =  I - ui ui* / hi  ,       hi = ui* ui / 2  ,       i = 1, 2, ..., N - 1 ,
ui -  N - мерный вектор, выбираемый таким образом, что матрица PiPi - 1 ... P1AP1 ... Pi - 1Pi имеет трехдиагональную структуру для последних  i  строк и столбцов, при этом последние  i  компонент вектора  ui  равны нулю ;
I -  единичная матрица порядка N .

Информация о матрице V и векторах  ui  ,  i = 1, ..., N - 1, запоминается и может быть впоследствии использована для восстановления собственных векторов исходной матрицы А.

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

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

    int afh5c_c (integer *nm, integer *n, real *a, real *d,
            real *e, real *e2, real *tau)

   Параметры
nm - число строк двумерного массива c, указанное при описании этого массива в вызывающей подпрограмме (тип: целый);
n - порядок исходной матрицы, n ≤ nm (тип: целый);
c - вещественный двумерный массив размерности nm*n, содержащий на входе в подпрограмму в своих первых n строках информацию об исходной комплексной эрмитовой матрице A, при этом:

ci j = re( ai j ) ,   для i ≥ j ,

ci j = im( aj i ) ,   для i < j ;

на выходе из подпрограммы в первых n строках массива c содержится информация об унитарных матрицах Pi, i = 1,2,...,n - 1 (см. математическое описание), при этом вещественные и мнимые части первых n - 1 элементов вектора ui* (последние i компонент равны нулю) хранятся соответственно в (n - i + 1) - ой строке слева от диагонального элемента и в (n - i + 1) - ом столбце над диагональным элементом, а на месте диагонального элемента (n - i + 1) - ой строки запоминается величина √ h i ;
d - вещественный вектор длины n, содержащий на выходе из подпрограммы вычисленные диагональные элементы трехдиагональной матрицы T;
e - вещественный вектор длины n, содержащий на выходе из подпрограммы в своих последних n - 1 компонентах поддиагональные элементы матрицы  T ,  e ( 1 ) = 0;
e2 - вещественный вектор длины n, на выходе из подпрограммы e2 ( i ) = (e ( i ))2 ;
tau - вещественный двумерный массив размерности 2*n, содержащий на выходе из подпрограммы в своей первой строке вещественные части диагональных элементов матрицы V, а во второй строке - мнимые части диагональных элементов матрицы V. 
   Версии
afh5p_c - приведение комплексной эрмитовой матрицы, заданной в компактной форме, к симметричной трехдиагональной матрице унитарными преобразованиями подобия с удвоенной точностью. 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1. 

Фактические параметры, соответствующие формальным параметрам e и e2, могут совпадать. В этом случае будут вычислены только поддиагональные элементы.

  2.  В подпрограмме afh5p_c параметры c, d, e, e2, tau имеют тип double.

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

int main(void)
{
    /* Initialized data */
    static float a[9] /* was [3][3] */ = { 1.f,3.f,0.f,-4.f,1.f,0.f,1.f,0.f,
                                           1.f };
    /* Local variables */
    extern int afh5c_c(int *, int *, float *, float *, float *, float *,
                       float *);
    static float d__[3], e[3];
    static int i__, n;
    static float e2[3], tau[6]   /* was [2][3] */;

#define a_ref(a_1,a_2) a[(a_2)*3 + a_1 - 4]
#define tau_ref(a_1,a_2) tau[(a_2)*2 + a_1 - 3]

    afh5c_c(&c__3, &c__3, a, d__, e, e2, tau);

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


Результаты:

                |  0.   -8.                         1. |
      c  =   | -6.    7.0710678119     0.  |
                |  0.    1.                         1. |

      d__   =   (1., 1., 1.)

      e     =   (0., 5., 1.)

      e2    =   (0., 25., 1.)

                       | -0.6     -1.     1. |
      tau_ref  =   | -0.8      0.      0. |

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

      t  =  v * p2 * p1 * a * p1 * p2 * v* ,
 
               | 1.   5.   0. |                |    1.          3 + 4i   -1. |
      t  =  | 5.   1.   1. |  ,    a  =  |    3-4i      1.           0. |
               | 0.   1.   1. |                |    1.          0.          1. |

      pi = I - ui*ui* / hi ,     i = 1, 2  , 
 где
      u1* =   ( i, 1, 0 ) ,     √ h1  =  1

      u2*  =  ( -6-8i, 0, 0 ) ,     √ h2  =  7.0710678119

      v  =  diag ( -0.6-0.8i, -1., 1. )