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

Подпрограмма:  AGG1R (модуль AGG1R_p)

Назначение

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

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

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

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 рассматривается как бесконечное.

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

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

procedure AGG1R(var A :Array of Real; var B :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 соответственно;
  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 к верхней почти треугольной и треугольной форме соответственно ортогональными преобразованиями подобия с помощью QZ - алгоритма.
AGT1R - вычисление всех собственных значений QR - алгоритмом в обобщенной проблеме собственных значений Ax = λBx для верхней почти тpеугольной и верхней треугольной вещественных матриц.
UTAG10 - подпрограмма выдачи диагностических сообщений при работе подпрограммы AGG1R. 
   Замечания по использованию
 

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

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

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

Unit TAGG1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AGG1R_p;

function TAGG1R: String;

implementation

function TAGG1R: String;
var
J,I,N,IERR :Integer;
ВЕТА :Array [0..2] of Real;
WK :Array [0..2] of Real;
//ALFA :Array [0..2] of Complex;
ALFA :Array [0..5] 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;
AGG1R(A,B,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 I:=1 to 3 do
 begin
  Result := Result + Format('%20.16f %20.16f ',
 [ALFA[I-1].re,ALFA[I-1].im]) + #$0D#$0A;
 end;
}
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 + Format('%s',[' IERR' + #$0D#$0A]);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
UtRes('TAGG1R',Result);  { вывод результатов в файл TAGG1R.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 )

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

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