БЧА НИВЦ МГУ. AGG3C_P. Обобщенная проблема собственных значений для матриц общего вида

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

Тексты тестовых примеров
tagg3c_p.zip

Подпрограмма: AGG3C (модуль AGG3C_p)

Назначение

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

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

Подпрограмма AGG3C вычисляет все собственные значения обобщенной проблемы Ax = λBx двух комплексных матриц А и B размера N на N с помощью LZ - алгоритма.

Информация о вычисленных собственных значениях выдается в вектоpах ALFA и BETA длины N, причем K - ое собственное значение λ_k определяется по формуле:

                        λ_k = ALFA(k) / BETA(k), при BETA(k) ≠ 0,
   (1)                 λ_k = ∞                            , при BETA(k) = 0,          
                        λ_k = любое число          , при ALFA(k) = BETA(k) = 0.

R.S.Martin, J.H.Wilkinson, The Modified LR - Algorithm for Complex Hessenberg Matrices, Numer. Math., 12, 1968.

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

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

procedure AGG3C(var A :Array of Complex; var B :Array of Complex;
                var ALFA :Array of Complex;
                var BETA :Array of Complex; N :Integer;
                var IERR :Integer);

   Параметры

A, B -	комплексные двумерные массивы размера N на N, содержащие исходные матрицы;
ALFA - BETA	комплексные векторы длины N, в которых содержится информация о вычисленных собственных значениях; при этом собственные значения определяются формулой (1);
N -	порядок исходных матриц (тип: целый);
IERR -	целая переменная, служащая для сообщения об ошибках, обнаруженых в ходе работы подпрограммы; значение IERR полагается равным номеpу собственного значения, для вычисления которого потребовалось более 30 итераций; при этом собственные значения с индексами IERR + 1, ..., N вычислены правильно, а с индексами 1, ..., IERR не вычисляются.

   Версии :   нет

   Вызываемые подпрограммы

AFG3C -	приведение пары комплексных матриц к веpхнему почти треугольному и к веpхнему треугольному виду с помощью LZ - алгоритма.
AGT3C -	вычисление всех собственных значений в обобщенной проблеме Ax = λBx для комплексных верхней почти треугольной матрицы A и веpхней треугольной матрицы B с помощью LR - алгоритма.
UTAG10 -	подпрограмма выдачи диагностических сообщений при работе подпрограммы AGG3C.

   Замечания по использованию

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

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

Unit TAGG3C_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AGG3C_p;

function TAGG3C: String;

implementation

function TAGG3C: String;
var
J,I,N,K,IERR :Integer;
ALFA :Array [0..2] of Complex;
ВЕТА :Array [0..2] of Complex;
LAMBDA :Array [0..2] of Complex;
const
A :Array [0..8] of Complex = ( ( re:1.0; im:0.0 ),( re:-10.0; im:0.0 ),( 
re:5.0; im:0.0 ),( re:0.5; im:0.0 ),( re:2.0; 
im:0.0 ),( re:1.0; im:0.0 ),( re:0.0; im:0.0 ),( 
re:0.0; im:0.0 ),( re:0.5; im: 0.0 ) );
B :Array [0..8] of Complex = ( ( re:0.5; im:0.0 ),( re:3.0; im:0.0 ),( re:4.0; 
im:0.0 ),( re:0.0; im:0.0 ),( re:3.0; im:0.0 ),( 
re:0.5; im:0.0 ),( re:0.0; im:0.0 ),( re:0.0; 
im:0.0 ),( re:1.0; im: 0.0 ) );
label
_11,_10;
begin
Result := '';  { результат функции }
Result := Result + Format('%s',
 [' ОПРЕДЕЛЕНИЕ СОБСТВЕННЫХ ЗНАЧЕНИЙ ' + #$0D#$0A +
 ' YРАВНЕНИЯ BИДA:AX=LAMBDA*BX, ГДЕ A - КОМПЛЕКСНАЯ  MATPИЦA/ И B - КОМПЛЕКСНАЯ MATPИЦA']) + #$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 %20.16f ',
 [A[(I-1)+(J-1)*3].re,A[(I-1)+(J-1)*3].im]) + #$0D#$0A;
   end;
 end;
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 %20.16f ',
 [B[(I-1)+(J-1)*3].re,B[(I-1)+(J-1)*3].im]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
N := 3;
AGG3C(A,B,ALFA,BETA,N,IERR);
Result := Result + Format('%s',[' PEЗYЛЬTAT' + #$0D#$0A]) + #$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 %20.16f ',
 [A[(I-1)+(J-1)*3].re,A[(I-1)+(J-1)*3].im]) + #$0D#$0A;
   end;
 end;
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 %20.16f ',
 [B[(I-1)+(J-1)*3].re,B[(I-1)+(J-1)*3].im]) + #$0D#$0A;
   end;
 end;
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;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' BETA' + #$0D#$0A]);
for I:=1 to 3 do
 begin
  Result := Result + Format('%20.16f %20.16f ',
 [BETA[I-1].re,BETA[I-1].im]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' IERR' + #$0D#$0A]);
Result := Result + Format('%3d ',[IERR]) + #$0D#$0A;
for K:=1 to 3 do
 begin
  if ( CReal(BETA[K-1]) = 0.0 ) 
   then goto _11;
  LAMBDA[K-1] := DivC(ALFA[K-1],BETA[K-1]);
  goto _10;
_11:
  LAMBDA[K-1] := Cmplx(3.4E38,3.4E38);
_10:
 end;
Result := Result + Format('%s',[' LAMBDA=ALFA/BETA' + #$0D#$0A]);
for I:=1 to 3 do
 begin
  Result := Result + Format('%20.16f %20.16f ',
 [LAMBDA[I-1].re,LAMBDA[I-1].im]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TAGG3C',Result);  { вывод результатов в файл TAGG3C.res }
exit;
end;

end.

Результаты:

      IERR  =  0

                      |  2.79963-3.65018*i |
      ALFA  =  |  0.42142+0.54945*i |
                      | -1.09877+0.0*i        |

                      | -1.05899-1.84748*i |
      BETA  =  | -0.15941+0.27809*i |
                      | -2.19753+0.0*i        |

Собственные значения   λ_k = ALFA(k) / BETA(k)

      λ₁  =  0.8333 + 1.993i
      λ₂  =  0.8333 - 1.993*i
      λ₃  =  0.5000 + 0.0*i