Текст подпрограммы и версий
afg7r_p.zip  afg7e_p.zip  afg7c_p.zip  afg7z_p.zip 		Тексты тестовых примеров
tafg7r_p.zip  tafg7e_p.zip  tafg7c_p.zip  tafg7z_p.zip

Подпрограмма:  AFG7R (модуль AFG7R_p)

Назначение

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

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

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

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

 ai j  =  0 ,   ecли  i > j   и  ecли  1 ≤ j < LOW  или  IGH < i ≤ N ,

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

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

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

     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.

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

procedure AFG7R(NM :Integer; N :Integer; LOW :Integer; IGH :Integer;
                var A :Array of Real; var ORT :Array of Real);


   Параметры
NM - число строк двумерного массива  А, указанное при описании этого массива в вызывающей подпрограмме (тип: целый);
N - порядок исходной матрицы,  N ≤ NМ (тип: целый);
      LOW -
      IGH   		заданные граничные индексы строк и столбцов подматрицы исходной матрицы, которую требуется привести к форме Хессенберга (тип: целый); если матрица масштабировалась, то LОW, IGН - выходные параметры подпрограммы AMB1R , в общем случае можно взять LОW = 1,  IGН = N;
A - вещественный двумерный массив размерности NМ * N, содержащий в своих первых  N строках исходную матрицу; в результате работы подпрограммы  А содержит вычисленную матрицу Хессенберга, а в остальной части массива  А в столбцах с номерами  i = LОW, LОW + 1, ..., IGН - 2 запоминаются ненулевые компоненты векторов  ui (см. математическое описание), начиная со 2 - ой ненулевой, т.е. с ( i + 2) - ой компоненты;
ORT - вещественный вектор длины IGН, содержащий в результате работы подпрограммы оставшуюся часть информации о векторах  ui, при этом компонента ОRТ ( i ) содержит  i - ую компоненту вектора  ui - 1,  i = LОW + 1, ..., IGН - 1, т.е. первую ненулевую компоненту вектора  ui - 1;
подпрограмма использует только компоненты вектора ОRТ с номерами LОW + 1, LОW + 2, ..., IGН. 
   Версии
AFG7E - приведение к верхней форме Хессенберга методом отражений вещественной матрицы, заданной с расширенной (Extended) точностью.
AFG7C - приведение к верхней форме Хессенберга методом отражений комплексной матрицы (см. замечания по использованию).
AFG7Z - приведение к верхней форме Хессенберга методом отражений комплексной матрицы, заданной с расширенной (Extended) точностью (см. замечания по использованию). 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
  1. 

В подпрограмме АFG7E параметры А, ОRТ имеют тип Extended.

  2. 

В подпрограмме АFG7С исходная комплексная матрица порядка  N задается в двух вещественных массивах АR и АI размеров NМ * N, содержащих в своих первых  N строках ее вещественные и мнимые части соответственно.
Параметр NМ определяет число строк двумерных массивов АR и АI, указанное при описании этих массивов в вызывающей подпрограмме.
В результате работы подпрограммы АFG7С массивы АR и АI содержат соответственно вещественную и мнимую части вычисленной матрицы Хессенберга, а в столбцах с номерами  i = LОW, LОW + 1, ..., IGН - 2 под поддиагональными элементами - соответственно вещественные и мнимые части векторов  ui, начиная с ( i + 2) - ой компоненты.
Для запоминания первых ненулевых компонент векторов  ui используются два вещественных вектора ОRТR и ОRТI длины IGН, при этом  i - ые компоненты этих векторов содержат соответственно вещественную и мнимую части  i - ой компоненты вектора  ui - 1.

Первый оператор подпрограммы АFG7С имеет вид: 

     procedure AFG7C(NM :Integer; N :Integer; LOW :Integer; IGH :Integer;
                var AR :Array of Real; var AI :Array of Real;
                var ORTR :Array of Real; var ORTI :Array of Real);
  3.  Подпрограмма АFG7Z имеет такие же параметры, как и подпрограмма АFG7С, только при этом АR, АI, ОRТR и ОRТI имеют тип Extended.

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

Unit TAFG7R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AFG7R_p;

function TAFG7R: String;

implementation

function TAFG7R: String;
var
N,J,I,_i :Integer;
ORT :Array [0..3] of Real;
const
A :Array [0..15] of Real = ( 2.0,0.0,3.0,4.0,1.0,1.0,-0.6,-0.8,1.0,-0.6,1.64,
-0.48,1.0,-0.8,-0.48,1.36 );
begin
Result := '';  { результат функции }
N := 4;
ORT[0] := 0.0;
AFG7R(4,4,1,4,A,ORT);
Result := Result + #$0D#$0A;
for I:=1 to N do
 begin
  for J:=1 to N do
   begin
    Result := Result + Format('  A= + #$0D#$0A  %20.16f ',
 [A[(I-1)+(J-1)*4]]) + #$0D#$0A;
   eND;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  ORT=' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 3 do
 begin
  Result := Result + Format('%20.16f ',[ORT[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TAFG7R',Result);  { вывод результатов в файл TAFG7R.res }
exit;
end;

end.


Результаты:

             |    2     - 1.4      1    - 0.22
             |  - 5       1.        1    - 9.1E - 13
    A  =  |    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 )  .