Текст подпрограммы и версий
afh5r_p.zip , afh5e_p.zip 		Тексты тестовых примеров
tafh5r_p.zip , tafh5e_p.zip

Подпрограмма:  AFH5R (модуль AFH5R_p)

Назначение

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

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

Подпрограмма АFН5R осуществляет приведение вещественной матрицы А порядка N к симметричной трехдиагональной матрице Т с помощью следующей последовательности ортогональных преобразований подобия 

            Ai = Pi*Ai-1*Pi ,     i = 1, 2, ..., N-2 ,
 где
           А0 = А ,
           Рi =  I - uiuiT / hi ,      hi = uiTui / 2 ,

а векторы ui ∈ ЕN выбираются таким образом, чтобы последние  i строк и столбцов матрицы Аi имели бы трехдиагональную структуру. При этом у векторов ui отличными от нуля являются только компоненты с номерами 1, 2, ..., N - 1 . I - единичная матрица порядка N .

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

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

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

procedure AFH5R(N :Integer; NV :Integer; var A :Array of Real;
                var D :Array of Real; var E :Array of Real;
                var E2 :Array of Real);


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

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

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

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

Unit tafh5r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AFH5R_p;

function tafh5r: String;

implementation

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

end.

Результаты:

      A   =  (0., 0., 0., -1., 1., 1., 4., 0., 4., 4.)T
      D   =  (2., 1., 1., 2.)T
      E   =  (0., 0., -1., -4.)T
      E2  =  (0., 0., 1., 16.)T

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

              |  2     0     0     0    |
              |  0     1    -1     0    |
      T  =  |  0    -1     1    -4   |
              |  0     0    -4     2    |

      Pi = I - uiuiT / hi ,   i  =  1, 2 , 
 где
      u1  =  (4., 0., 4., 0.)T ,      √ h 1   =  4

      u2  =  (-1., 1., 0., 0.)T ,     √ h 2   =  1