Текст подпрограммы и версий
aet1r_p.zip , aet1e_p.zip 		Тексты тестовых примеров
taet1r_p.zip , taet1e_p.zip

Подпрограмма:  AET1R (модуль AET1R_p)

Назначение

Вычисление всех собственных значений и собственных векторов вещественной верхней матирцы Хессенберга QR - алгоpитмом с двойным сдвигом.

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

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

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

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

procedure AET1R(NM :Integer; N :Integer; LOW :Integer;
                IGH :Integer; var H :Array of Real;
                var WR :Array of Real; var WI :Array of Real;
                var Z :Array of Real; var IERR :Integer);


   Параметры
NM - число строк двумерных массивов Н и Z, указанное при описании этих массивов в вызывающей подпрограмме (тип: целый);
N - порядок исходной матрицы (тип: целый);
      LOW -
      IGH   		выходные параметры подпрограммы АМВ1R (тип: целый); если матрица не масштабировалась, то можно взять LОW := 1, IGН := N;
H - вещественный двумерный массив размерности NМ на N, содержащий на входе в подпрограмму в своих первых N строках исходную матрицу Хессеберга;
WR, WI - вещественные векторы длины N, содержащие на выходе из подпрограммы соответственно вещественные и мнимые части вычисленных собственных значений, при этом комплексно - сопряженные собственные значения располагаются последовательно, причем первым идет собственное значение с положительной мнимой частью;
Z - вещественный двумерный массив размерности NМ на N, содержащий на выходе из подпрограммы вычисленные собственные векторы, при этом:
если WI (J) = 0 , то соответствующий собственный вектор расположен в J - ом столбце массива Z;
если WI (J) > 0 , то в J - ом и (J + 1) - ом столбцах располагаются соответственно вещественная и мнимая части собственного вектора;
если WI (J) < 0 , то собственный вектор не запоминается, т.к. может быть получен как сопряженный собственному вектору, соответствующему (J - 1) - ому собственному значению;
  собственные векторы не нормируются;
IERR - целочисленная переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; на выходе из подпрограммы IЕRR равно индексу собственного значения для вычисления которого требуется более 30 итераций, при этом собственные значения с индексами IЕRR + 1, IЕRR + 2, ...,N вычислены правильно, а собственные векторы не вычисляются; если вычислены все собственные значения и собственные векторы, то IЕRR := 0. 
   Версии
AET1E - вычисление всех собственных значений и собственных векторов вещественной матрицы Хессенберга, заданной с расширенной (Extended) точностью, с помощью QR - алгоритма с двойным сдвигом. 
   Вызываемые подпрограммы
UTAE12 - подпрограмма выдачи диагностических сообщений при работе подпрограмм АЕТ1R и АЕТ1E. 
   Замечания по использованию
  1. 

В подпрограмме АЕТ1E параметры Н, WR, WI, Z имеют тип Extended.

  2.  Подпрограмма АЕТ1R (АЕТ1E) не сохраняет исходный массив Z.

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

Unit TAET1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AET1R_p;

function TAET1R: String;

implementation

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

end.

Результаты:

       IERR  =  0 ,   WR  =  (2., 4., 2.) ,  WI  =  (0., 0., 0.)

                | 1   4a    0   |
       Z  =  | 0     a   -a   |
                | 0     a    a   |

   где    a = 0.707106781187