Текст подпрограммы и версий
agg0r_p.zip  		Тексты тестовых примеров
tagg0r_p.zip 

Подпрограмма:  AGG0R (модуль AGG0R_p)

Назначение

Вычисление всех собственных значений и собственных вектоpов QZ - алгоритмом в обобщенной проблеме собственных значений AX = λBX двух вещественных матриц.

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

Алгоритм, реализованный в подпрограмме AGG0R, может быть разбит на четыре этапа.

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

Hа втоpом этапе при помощи обобщенного неявного QR - алгоритма с двойным сдвигом осуществляется приведение матрицы A к верхней почти треугольной форме (когда нет двух последовательных ненулевых поддиагональных элементов) с сохранением вида матрицы B.

Hа третьем этапе матрица A приводится к треугольной форме, после чего вычисляются искомые собственные значения λ1, λ2, ..., λN в виде λI = ALFA (I) / BETA (I), когда BETA (I) ≠ 0. Если BETA (I) = 0, то соответствующее собственное значение λI рассматривавается как бесконечное.

Hа четвертом этапе вычисляются собственные векторы, которые помещаются в комплексной матрице так, что I - й столбец содержит собственный вектоp, соответствующий I - му собственному значению.

C.B.Moler, G.W.Stewart, An Algorithm for Generalized Matrix Eigenvalue Problems, SIAM J.Numer.Anal.10, 1973.

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

procedure AGG0R(var A :Array of Real; var B :Array of Real;
                var V :Array of Real; var ALFA :Array of Real;
                var BETA :Array of Real; var WK :Array of Real;
                N :Integer; var IERR :Integer);


   Параметры
A, B - вещественные двумерные массивы размера N на N, содержащие исходные матрицы A и B соответственно;
V - комплексный двумерный массив размера N на N, содержащий вычисленные собственные векторы;
  ALFA -
  BETA   		комплексный и вещественный векторы длины N, такие, что искомые собственные значения выражаются через компоненты этих вектоpов следующим образом:
λI = ALFA (I) / BETA (I); если BETA (I) = 0, то соответствующее собственное значение рассматривается как бесконечное;
WK - вещественный вектоp длины N, используемый как рабочий;
N - заданный порядок исходных матриц (тип: целый);
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение IERR полагается равным J+128, где J - номеp собственного значения, для вычисления которого потребовалось более 30 итераций, при этом собственные значения с индексами J+1, ..., N вычислены правильно, а с индексами 1, ..., J и соответствующие собственные векторы могут быть вычислены не точно. 
   Версии :   нет 
   Вызываемые подпрограммы
AFG3R - приведение пары вещественных матриц A и B к верхней почти треугольной и треугольной форме соответственно ортогональными преобразованиями подобия с помощью Q*Z - алгоритма.
AGT0R - вычисление всех собственных значений и соответственных вектоpов QR - алгоритмом в обобщенной проблеме собственных значений Ax = λBx для верхней почти треугольной и верхней треугольной вещественных матриц.
UTAG10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы AGG0R. 
   Замечания по использованию
 

Исходные матрицы A и B не сохраняются.

Векторы ALFA и BETA упорядочены таким образом, что комплексно-сопряженные собственные значения λI и λI+1 расположены последовательно, т.е. ALFA (I+1) / BETA (I+1) комплексно - сопряжено с ALFA (I) / BETA (I). При этом ALFA (I+1) не обязятельно комплексно-сопряжено с ALFA (I).

Собственные векторы нормализованы таким образом, что наибольшая компонента по модулю pавна 1.

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

Unit TAGG0R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AGG0R_p;

function TAGG0R: String;

implementation

function TAGG0R: String;
var
J,I,N,IERR :Integer;
ВЕТА :Array [0..2] of Real;
WK :Array [0..2] of Real;
ALFA :Array [0..5] of Real;
V :Array [0..17] of Real;
const
A :Array [0..8] of Real = ( 1.0,-10.0,5.0,0.5,2.0,1.0,0.0,0.0,0.5 );
B :Array [0..8] of Real = ( 0.5,3.0,4.0,0.0,3.0,0.5,0.0,0.0,1.0 );
begin
Result := '';  { результат функции }
Result := Result + Format('%s',
 [' ОПРЕДЕЛЕНИЕ СОБСТВЕННЫХ ЗНАЧЕНИЙ И СОБСТВЕННЫХ BEKTOPOB' + #$0D#$0A +
 ' YРАВНЕНИЯ BИДA:AX=ЛЯMБДA*BX, ГДЕ A И B - ВЕЩЕСТВЕННЫЕ MATPИЦЫ' + #$0D#$0A]) + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' A' + #$0D#$0A]);
for I:=1 to 3 do
 begin
  for J:=1 to 3 do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*3]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' B' + #$0D#$0A]);
for I:=1 to 3 do
 begin
  for J:=1 to 3 do
   begin
    Result := Result + Format('%20.16f ',[B[(I-1)+(J-1)*3]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
N := 3;
AGG0R(A,B,V,ALFA,BETA,WK,N,IERR);
Result := Result + Format('%s',[' PEЗYЛЬTAT' + #$0D#$0A]) + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' A' + #$0D#$0A]);
for I:=1 to 3 do
 begin
  for J:=1 to 3 do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*3]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' B' + #$0D#$0A]);
for I:=1 to 3 do
 begin
  for J:=1 to 3 do
   begin
    Result := Result + Format('%20.16f ',[B[(I-1)+(J-1)*3]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' ALFA' + #$0D#$0A]);
for J:=1 to 3 do
 begin
  for I:=1 to 2 do
   begin
    Result := Result + Format('%20.16f ', [ALFA[(I-1)+(J-1)*2]]);
   end;
  Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' BETA' + #$0D#$0A]);
for I:=1 to 3 do
 begin
  Result := Result + Format('%20.16f ',[BETA[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' V' + #$0D#$0A]);
for J:=1 to 9 do
 begin
  for I:=1 to 2 do
   begin
    Result := Result + Format('%20.16f ', [V[(I-1)+(J-1)*2]]);
   end;
  Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' IERR' + #$0D#$0A]);
Result := Result + Format(' %3d ',[IERR]) + #$0D#$0A;
UtRes('TAGG0R',Result);  { вывод результатов в файл TAGG0R.res }
exit;
end;

end.

Результаты:

      IERR = 0
 
                      | 1.27052+3.03864*i |
      ALFA  =  | 0.40869-0.97744*i |
                      | 1.00306+0.00000*i |
 
      BETA  =  ( 1.52462,  0.49043,  2.00612 )

Собственные векторы: 

               | -0.25205+0.19169*i |
               | -0.25205-0.19169*i |
               |  0.0+0.0*i                 |
               | -0.08799-0.72598*i |
      V  =  | -0.08799+0.72598*i |
               |  0.0+0.0*i                 |
               |  1.00000+0.00000*i |
               |  1.00000-0.00000*i  |
               |  1.0+0.0*i                 |

Собственные значения: 

      λ1  =  0.83333+1.9930*i
      λ2  =  0.83333-1.9930*i
      λ3  =  0.50000+0.0000*i