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

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

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

Подпрограмма: AGT3C (модуль AGT3C_p)

Назначение

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

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

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

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

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

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

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

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

procedure AGT3C(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, содержащие информацию о собственных значениях λ_k обобщенной проблемы; при этом собственные значения определяются формулой (1);
N -	заданный порядок исходных матриц A, B (тип: целый);
IERR -	целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; значение IERR полагается равным номеpу собственного значения, для вычисления которого потребовалось более 30 итераций; при этом собственные значения с индексами IERR+1, ..., N вычислены правильно, а с индексами 1, ..., IERR не вычисляются.

   Версии :   нет

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

UTAG10 -

подпрограмма выдачи диагностических сообщений при работе подпрограммы AGT3C.

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

Подпрограмма AGT3C не сохраняет исходные матрицы A, B.

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

Unit TAGT3C_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AFG3C_p, AGT3C_p;

function TAGT3C: String;

implementation

function TAGT3C: String;
var
J,I,N,M,K,IERR :Integer;
V :Array [0..8] of Complex; //рабочий массив для п/п AFG3C
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 - КОМПЛЕКСНАЯ BEPXHЯЯ' + #$0D#$0A +
 ' ПОЧТИ TPEYГОЛЬНАЯ MATPИЦA,' + #$0D#$0A +
 ' B - КОМПЛЕКСНАЯ BEPXHETPEYГОЛЬНАЯ 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;
M := 1;
AFG3C(A,B,V,N,M);
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;
AGT3C(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('TAGT3C',Result);  { вывод результатов в файл TAGT3C.res }
exit;
end;

end.

Результаты:

      IERR  =  0

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

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

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

      λ₁  =  0.83333 + 1.9930i
      λ₂  =  0.83331 - 1.9931i
      λ₃  =  0.50000 + 0.0i