Текст подпрограммы и версий
afh5c_p.zip , afh5z_p.zip 		Тексты тестовых примеров
tafh5c_p.zip , tafh5z_p.zip

Подпрограмма:  AFH5C (модуль AFH5C_p)

Назначение

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

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

Подпрограмма AFH5C приводит исходную эрмитову матрицу А порядка 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.

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

procedure AFH5C(NM :Integer; N :Integer; var A :Array of Real;
                var D :Array of Real; var E :Array of Real;
                var E2 :Array of Real; var TAU :Array of Real);


   Параметры
NM - число строк двумерного массива С, указанное при описании этого массива в вызывающей подпрограмме (тип: целый);
N - порядок исходной матрицы, N ≤ NМ (тип: целый);
C - вещественный двумерный массив размерности NМ*N, содержащий на входе в подпрограмму в своих первых N строках информацию об исходной комплексной эрмитовой матрице А, при этом:

Ci j = Re( Аi j ) ,   для i ≥ j ,

Сi j = Im( Аj i ) ,   для i < j ;

на выходе из подпрограммы в первых N строках массива C содержится информация об унитарных матрицах Рi, i = 1,2,...,N - 1 (см. математическое описание), при этом вещественные и мнимые части первых N - 1 элементов вектора ui* (последние i компонент равны нулю) хранятся соответственно в (N - i + 1) - ой строке слева от диагонального элемента и в (N - i + 1) - ом столбце над диагональным элементом, а на месте диагонального элемента (N - i + 1) - ой строки запоминается величина √ h i ;
D - вещественный вектор длины N, содержащий на выходе из подпрограммы вычисленные диагональные элементы трехдиагональной матрицы Т;
E - вещественный вектор длины N, содержащий на выходе из подпрограммы в своих последних N - 1 компонентах поддиагональные элементы матрицы  Т ,  Е ( 1 ) = 0;
E2 - вещественный вектор длины N, на выходе из подпрограммы Е2 ( I ) = (Е ( I ))2 ;
TAU - вещественный двумерный массив размерности 2*N, содержащий на выходе из подпрограммы в своей первой строке вещественные части диагональных элементов матрицы V, а во второй строке - мнимые части диагональных элементов матрицы V. 
   Версии
AFH5Z - приведение комплексной эрмитовой матрицы, заданной в компактной форме, к симметричной трехдиагональной матрице унитарными преобразованиями подобия с расширенной (Extended) точностью. 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1. 

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

  2.  В подпрограмме АFН5Z параметры С, D, Е, Е2, ТАU имеют тип Extended.

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

Unit tafh5c_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AFH5C_p;

function tafh5c: String;

implementation

function tafh5c: String;
var
N,J,I,_i :Integer;
D :Array [0..2] of Real;
E :Array [0..2] of Real;
E2 :Array [0..2] of Real;
TAU :Array [0..5] of Real;
const
A :Array [0..8] of Real = ( 1.0,3.0,0.0,-4.0,1.0,0.0,1.0,0.0,1.0 );
begin
Result := '';  { результат функции }
AFH5C(3,3,A,D,E,E2,TAU);
N := 3;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  A='  ]) + #$0D#$0A;
for J:=1 to N do
 begin
  for I:=1 to N do
   begin
    Result := Result + Format('  %20.16f ',
 [A[(I-1)+(J-1)*3]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   D=']);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[D[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   E=']);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[E[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  E2=']);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[E2[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  TAU='  ]) + #$0D#$0A;
for J:=1 to N do
 begin
  for I:=1 to 2 do
   begin
    Result := Result + Format('  %20.16f ',
 [TAU[(I-1)+(J-1)*2]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
UtRes('tafh5c',Result);  { вывод результатов в файл tafh5c.res }
exit;
end;

end.


Результаты:

                |  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  =   | -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. )