Подпрограмма:  AGS0R (модуль AGS0R_p)

Назначение

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

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

Метод Ланцоша для обобщенной проблемы собственных значений

(1)       A x  =  λ B x ,  

где  A и  B - заданные вещественные симметричные матрицы порядка  N,  B > 0, и заданного начального вектора  v ∈ E^N,  v ≠ 0, состоит в построении последовательности векторов Ланцоша  v₁, v₂, ..., v_j и последовательности симметричных трехдиагональных матриц  T₁, T₂, ..., T_j по следующим формулам:

         v1  =  v / (vT Bv)1/2 ,    β1  =  0 ,
         ui   =  Avi - βi Bvi - 1 ,
(2)     αi  =  uiT vi ,

         T1  =  α1 ,                           i = 1 ,

                 |               |  0 |
                 |               |  .  |
         Ti =  |    Ti-1     |  0 |  ,          i > 1 ,
                 | _______| βi |
                 | 0 ... 0 βi   αi |

         wi + 1  =  ui - αi Bvi ,

         решение системы    Bui +1  =  wi + 1 ,
         βi + 1  =  (uTi + 1 wi + 1)1/2 ,
         vi + 1  =  ui + 1 / βi + 1             i = 1, ..., j . 

Метод Ланцоша (2) может быть рассмотрен как реализация для задачи (1) метода Ритца, где в качестве подпространства аппроксимации используется подпространство  Крылова   K ^j = span { v₁, B^- 1 Av₁, ..., (B^- 1 A)^j - 1 v₁ }, в котором задано скалярное произведение (x, y)_B = x^T By. При этом, если (θ_i ^j, s_i ^j),  1 ≤ i ≤ j, - собственная пара симметричной трехдиагональной матрицы  T_j, получаемой на  j - ом шаге метода Ланцоша (2), а  V_j - матрица, столбцами которой являются векторы Ланцоша  v₁, v₂, ..., v_j, то (θ_i ^j, V_j s_i ^j) будет соответствующей парой Ритца для пучка (A, B).

Полезной особенностью метода Ланцоша является то, что часто уже при  j ≈ 2 √N крайние пары Ритца оказываются хорошими приближениями к искомым крайним собственным парам пучка (A, B).

В подпрограмме AGS0R используется метод Ланцоша с выборочной ортогонализацией - устойчивая к влиянию ошибок округления модификамция простого метода Ланцоша (2), отличающаяся от него тем, что на некторых шагах процесса (2) перед вычислением  v_i + 1 проводится  B - ортоганализация вектора  u_i + 1 относительно некоторых векторов Ритца.

Как видно из формул (2), метод Ланцоша использует исходную матрицу  A только в операциях умножения матрицы на вектор, а матрицу  B - в операциях умножения матрицы на вектор и решения системы с матрицей  B, поэтому всю необходимую информацию о матрицах  A и  B удобно задавать в виде трех соответствующих подпрограмм, в которых пользователь может максимально учесть специфику матриц. Особенно удобен метод Ланцоша для разреженных или имеющих регулярную структуру матриц.

В подпрограмме AGS0R гарантированной точностью по невязке для вычисленных собственных пар является величина  k √ε  || B^- 1 A ||, где  k ~ 1,  ε - относительная машинная точность, но часто возможно получение лучшей точности. Под невязкой для собственной пары (θ, y) здесь понимается величина  || B^- 1 Ay - θy ||_B / || y ||_B.

Выход из подпрограммы происходит, когда все требуемые собственные векторы будут вычислены с заданной точностью по невязке или если продолжение процесса (2) окажется невозможным или бесполезным.

Б.Парлетт. Симметричная проблема собственных значений. Численные методы. М.: Мир, 1983.

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

procedure AGS0R(N :Integer; KL :Integer; KR :Integer;
                FA :Proc_F5A; FB :Proc_F5A; SB :Proc_F5A; EPS :Real;
                var M :Integer; var VJ :Array of Real; K :Integer;
                L :Integer; var R :Array of Real;
                var R1 :Array of Real; var R2 :Array of Real;
                var EV :Array of Real; var ER :Array of Real;
                var V :Array of Real; var Y :Array of Real;
                var IP :Integer);

   Параметры 

   Версии 

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

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

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

Unit TAGS0R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FFAAGS0R_p, FFBAGS0R_p, FSBAGS0R_p, AGS0R_p;

function TAGS0R: String;

implementation

function TAGS0R: String;
var
M,IP,_i,J,I :Integer;
EV :Array [0..2] of Real;
ER :Array [0..2] of Real;
R :АRray [0..81] of Real;
R1 :Array [0..2] of Real;
R2 :Array [0..2] of Real;
V :Array [0..8] of Real;
Y :Array [0..14] of Real;
const
VJ :Array [0..2] of Real = ( 1.0,1.0,1.0 );
label
_55;
begin
for _i:=0 to High(R) do  //обнуление рабочего массива
 R[_i] := 0.0;
Result := '';  { результат функции }
M := 3;
IP := 0;
AGS0R(3,3,0,FFAAGS0R,FFBAGS0R,FSBAGS0R,1.E-3,M,VJ,5,82,R,R1,R2,EV,ER,V,Y,IP);
ResuLT := RЕSULT + FОRMАT('%S',['  IР=']);
ResuLT := RЕSult + Format('%5d ',[IР]);
ResULT := ResULT + FОrmat('%s',['   M=']);
Result := RЕSULT + Format('%5D ',[М]) + #$0D#$0A;
if ( IP < 0 )
 then goto _55;
Result := Result + Format('%s',['  EV=' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[EV[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  Y =' + #$0D#$0A]);
for I:=1 to 3 do
 begin
  for J:=1 to 3 do
   begin
    Result := Result + Format('  %20.16f ',
 [Y[(I-1)+(J-1)*3]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['  ER=' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 2 do
 begin
  Result := Result + Format('%20.16f ',[ER[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
_55:
UtRes('TAGS0R',Result);  { вывод результатов в файл TAGS0R.res }
exit;
end;

end.

Unit ffaags0r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

procedure ffaags0r(var X :Array of Real; var Z :Array of Real; N :Integer);

implementation

procedure ffaags0r(var X :Array of Real; var Z :Array of Real; N :Integer);
var
I :Integer;
const
RAB :Array [0..2] of Real = ( 1.0,9.0,14.0 );
label
_10;
begin
for I:=1 to N do
 begin
  Z[I-1] := RAB[I-1]*X[I-1]-Z[I-1];
_10:
 end;
end;

end.

Unit ffbags0r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

procedure ffbags0r(var X :Array of Real; var Z :Array of Real; N :Integer);

implementation

procedure ffbags0r(var X :Array of Real; var Z :Array of Real; N :Integer);
var
I :Integer;
const
RAB :Array [0..2] of Real = ( 1.0,1.0,2.0 );
label
_10;
begin
for I:=1 to N do
 begin
  Z[I-1] := RAB[I-1]*X[I-1];
_10:
 end;
end;

end.

Unit fsbags0r_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc;

procedure fsbags0r(var X :Array of Real; var Z :Array of Real;
                N :Integer);

implementation

procedure fsbags0r(var X :Array of Real; var Z :Array of Real;
                N :Integer);
var
I :Integer;
const
RAB :Array [0..2] of Real = ( 1.e0,1.e0,5.0e-1 );
label
_10;
begin
for I:=1 to N do
 begin
  Z[I-1] := X[I-1]*RAB[I-1];
_10:
 end;
end;

end.

Результаты:

       EV = (1., 7., 9.) ,
       ER = (1.E - 15, 4.E - 15, 9.E - 16) ,

               |  - 1.                5.E - 17                    - 9.E - 18 |
       Y =  |    2.E - 16      2.E - 15                       1.           |
               |    2.E - 16      0.70710678118655  - 3.E - 16 |

       IP = 0
       M = 3