Подпрограмма:  AES1R (модуль AES1R_p)

Назначение

Вычисление всех различных собственных значений вещественной симметричной разреженной или имеющей регулярную структуру матрицы методом Ланцоша.

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

Метод Ланцоша для заданной симметричной матрицы А порядка N и начального вектора v ∈ Е^N ,  v ≠ 0 ,  строит последовательность векторов Ланцоша  v₁, v₂, ..., v_j и последовательность трехдиагональных матриц  Т₁, Т₂, ..., Т_j по следующим формулам процесса (1):

      v1 = v / (vTv)1/2  ,     β1 = 0  , 

      ui = Avi - βivi-1  ,     αi = uiTvi  , 

  T1 =  α1 ,                            i = 1 ,

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

   wi+1 = ui - αivi ,    βi+1 = (wTi+1wi+1)1/2 ,    vi+1 = wi+1 / βi+1  ,
    i = 1, ..., j 

Подпрограмма AES1R использует следующее свойство метода Ланцоша: в спектре трехдиагональной матрицы Т_j , полученной при выполнении достаточно большого числа j (j ~ 2N,...,6N) шагов алгоритма (1), появляются приближения ко всем различным собственным значениям исходной матрицы А.

В подпрограмме AES1R на шагах с номерами

       p = mi = M1 + (i - 1)*MH ,   i = 1, 2, ...,   (М1, МН - заданные числа) 

процесс (1) прерывается и вычисляются собственные значения матриц Т_p. Для разделения собственных значений матриц Т_p на "истинные" и "лишние" применяется критерий, использующий собственные значения матрицы Т_p^{^} , получаемой из Т_p вычеркиванием первой строки и первого столбца. Для истинных собственных значений матриц Т_p вычисляются оценки точности.

Пусть К_j - число истинных собственных значений матрицы Т_j , К'_j - число истинных значений Т_j , являющихся приближениями к собственным значениям матрицы А с точностью max { ЕРS, 10*j*ε }*|| А ||₂, где ЕРS - заданная пользователем величина, ε - относительная машинная точность, || А ||₂ - спектральная норма матрицы А.

Выход из подпрограммы AES1R происходит

1) при выполнении одного из соотношений

        K'p = Kp = Kp1 ,   K'p = N ;   ( где p1 = mi-1 ) 

2) по заданному ограничению на число шагов процесса (1);

3) если очередное вычисленное значение β_j+1 очень мало;

4) если при работе подпрограммы AEE2R, используемой для вычисления собственных значений трехдиагональных матриц, произошла фатальная ошибка.

Метод Ланцоша использует исходную матрицу только в операциях умножения матрицы на вектор, поэтому информацию о матрице А удобно задавать в виде подпрограммы умножения матрицы на вектор, в которой пользователь может максимально учесть специфику матрицы.

Сullum J., Donath W.Е. Lanczos and thе computation in specified_p intervals of thе spectrum of large, sparse real symmetric_p matrices. Sparse Мatrix Рroc., Кnoxville, Тen., 1978. Рhiladelphia, Pа, 1979, 220-255.

Х.Д. Икрамов. Разреженные матрицы. Математический анализ. Т. 20. (Итоги науки и техн. ВИНИТИ АН СССР). М., 1982.

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

procedure AES1R(N :Integer; var M0 :Integer; M1 :Integer;
                MH :Integer; M2 :Integer; var EPS0 :Real; FF :Proc_F4A;
                var VJ :Array of Real; var W :Array of Real;
                var B :Array of Real; var A :Array of Real;
                var NAC :Integer; var TETA :Array of Real;
                var TETA1 :Array of Real; var IR :Array of Integer;
                var IP :Integer);

   Параметры 

   Версии :   нет 

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

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

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

Unit TAES1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, FAES1R_p, AES1R_p;

function TAES1R: String;

implementation

function TAES1R: String;
var
M0,I,K,IP :Integer;
EPS0 :Real;
W :Array [0..5] of Real;
B :Array [0..56] of Real;
ТЕТА :Array [0..27] of Real;
RR :Array [0..27] of Real;
A :Array [0..27] of Real;
IR :Array [0..27] of Integer;
const
VJ :Array [0..5] of Real = ( 1.0,1.0,1.0,1.0,1.0,1.0 );
begin
Result := '';  { результат функции }
EPS0 := 1.E-6;
M0 := 0;
AES1R(6,M0,4,4,28,EPS0,FAES1R,VJ,W,B,A,K,TETA,RR,IR,IP);
Result := Result + Format('%s',[' K,IP,M0 =']);
Result := Result + Format(' %4d  %4d %4d ',[K,IP,M0]) + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' TETA' + #$0D#$0A]);
for I:=1 to K do
 begin
  Result := Result + Format('%20.16f ',[TETA[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' RR' + #$0D#$0A]);
for I:=1 to K do
 begin
  Result := Result + Format('%20.16f ',[RR[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' IR' + #$0D#$0A]);
for I:=1 to K do
 begin
  Result := Result + Format('%5d ',[IR[I-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TAES1R',Result);  { вывод результатов в файл TAES1R.res }
exit;
end;

end.

Результаты:

      K  =  6,  IP  =  0
      EV  =   (1., 2., 3., 4., 5., 6.)
      ER  =   (1., 1., 1., 1., 1., 1.)
      IR  =   (-12, -12, -12, -12, -12, -12)