Текст подпрограммы и версий
aeh3r_p.zip , aeh3e_p.zip 		Тексты тестовых примеров
taeh3r_p.zip , taeh3e_p.zip

Подпрограмма:  AEH3R (модуль AEH3R_p)

Назначение

Вычисление собственных значений, принадлежащих заданному интервалу, их номеpов и соответствующих собственных векторов вещественной симметрической матрицы.

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

Исходная задача Ax = λx для вещественной симметрической матрицы  A решается путем приведения ее с помощью ортогональных преобразований  P к симметрической трехдиагональной матрице A1 = PTAP. Вычисление собственных значений, принадлежащих заданному интервалу, и их номеpов задачи А1y = λy осуществляется с помощью метода бисекций; при этом собственные значения обеих задач совпадают. Для вычисления собственных вектоpов используется метод обратных итераций. Собственные векторы обеих задач связаны соотношением  x = Py.

Уилкинсон, Райнш. "Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра", M., "Машиностроение", 1976.

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

procedure AEH3R(N :Integer; var A :Array of Real;
                var V :Array of Real; var EV :Array of Real;
                var RLB :Real; var RUB :Real; MM :Integer;
                var M :Integer; var IRAB1 :Array of Integer;
                var RAB1 :Array of Real; var IERR :Integer);


   Параметры
N - порядок исходной матрицы (тип: целый);
A - вещественный двумерный массив размера N на N, содержащий исходную матрицу; в подпрограмме используется только полный нижний треугольник матрицы  A;
V - вещественный двумерный массив размера N на MM, содержащий в первых  M столбцах вычисленные ортонормированные собственные векторы, соответствующие  M собственным значениям;
EV - вещественный одномерный массив длины MM, содержащий вычисленные в возрастающем порядке собственные значения, принадлежащие заданному интервалу;
      RLB -
      RUB   		заданные нижняя и верхняя границы интервала собственных значений (тип: вещественный); если RLB > RUB, то собственные значения не вычисляются;
MM - оценка свеpху числа собственных значений матрицы  A, принадлежащих заданному интервалу, MM ≤ N (тип: целый); если вычисленное число  M собственных значений, принадлежащих заданному интервалу, больше, чем MM, то собственные значения и соответствующие собственные векторы не вычисляются;
M - целая переменная, в которой запоминается вычисленное число собственных значений, принадлежащих заданному интервалу;
IRAB - целый вектоp длины MM, содержащий индексы расположенных в возрастающем порядке  M собственных значений;
RAB - вещественный вектоp длины 8 на N, используемый как рабочий;
IERR - целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом:
IERR=3*N+1 - если значение MM меньше числа вычисленных собственных значений  M на интервале; при этом собственные значения и соответствующие собственные векторы не вычисляются;
IERR= - k - если для вычисления собственного вектоpа с индексом  k потребовалось более 5 итераций; при этом компоненты этого вектоpа полагаются равными нулю. Если таких собственных вектоpов несколько, то значение IERR полагается равным индексу последнего из них. 
   Версии
AEH3E - вычисление собственных значений, принадлежащих заданному интервалу, их номеpов и соответствующих собственных вектоpов симметрической матрицы, заданной с расширенной (Extended) точностью. Массивы A, V, EV, RAB и переменные RLB, RUB имеют тип Extended. 
   Вызываемые подпрограммы
UTAE10 - подпрограмма выдачи диагностических сообщений при работе подпрограмм AEH3R и AEH3E. 
   Замечания по использованию
  Подпрограмма AEH3R сохраняет полный верхний треугольник исходной матрицы, остальные элементы используются для запоминания ортогональных преобразований.

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

Unit TAEH3R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, AEH3R_p;

function TAEH3R: String;

implementation

function TAEH3R: String;
var
J,I,N,MM,_i,M,IERR :Integer;
RLB,RUB :Real;
V :Array [0..11] of Real;
EV :Array [0..2] of Real;
IRAB1 :Array [0..2] of Integer;
RAB1 :Array [0..31] of Real;
const
A :Array [0..15] of Real = ( 1.0,0.42,0.54,0.66,0.42,1.0,0.32,0.44,0.54,0.32,
1.0,0.22,0.66,0.44,0.22,1.0 );
begin
Result := '';  { результат функции }
Result := Result + #$0D#$0A;
Result := Result + Format('%s',
 [' ВЫЧИСЛЕНИЕ НЕКОТОРЫХ СОБСТВЕННЫХ ЗНАЧЕНИЙ И COOTBETCTBYЮЩИХ COБCTBEHHЫX' + #$0D#$0A +
 ' ВЕКТОРОВ ВЕЩЕСТВЕННОЙ СИММЕТРИЧЕСКОЙ MATPИЦЫ' + #$0D#$0A +
 ' A' + #$0D#$0A]);
for I:=1 to 4 do
 begin
  for J:=1 to 4 do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*4]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
N := 4;
ММ := 3;
RLB := 0.0;
RUB := 1.0;
AEH3R(N,A,V,EV,RLB,RUB,MM,M,IRAB1,RAB1,IERR);
Result := Result + Format('%s',[' PEЗYЛЬTAT']) + #$0D#$0A; 
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' A' + #$0D#$0A]);
for I:=1 to 4 do
 begin
  for J:=1 to 4 do
   begin
    Result := Result + Format('%20.16f ',[A[(I-1)+(J-1)*4]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
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',[' V' + #$0D#$0A]);
for I:=1 to 4 do
 begin
  for J:=1 to 3 do
   begin
    Result := Result + Format('%20.16f ',[V[(I-1)+(J-1)*4]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',[' IERR' + #$0D#$0A]);
Result := Result + Format('%8d ',[IERR]) + #$0D#$0A;
UtRes('TAEH3R',Result);  { вывод результатов в файл TAEH3R.res }
exit;
end;

end.

Результаты:

Собственные значения в интервале (0., 1.):

       EV  =  (0.2422607,  0.6382838,  0.7967067) , 

Собственные векторы, соответствующие вычисленным
в интервале (0., 1.) собственным значениям:

       V1  =  (0.7188460,  0.0956990,  - 0.3874355,  - 0.5692064) , 
       V2  =  (0.3804499,  - 0.8502755,  - 0.0358896,  0.3619412) , 
       V3  =  (0.0503284,  - 0.2372265,  0.8128462,  - 0.5295958) , 

       IERR  =  0