БЧА НИВЦ МГУ. ASS1R_P. Решение систем с матрицами специального вида (невырожденными)

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

Тексты тестовых примеров
tass1r_p.zip , tass1e_p.zip

Подпрограмма: ASS1R (модуль ASS1R_p)

Назначение

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

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

ASS1R вычисляет решение системы линейных алгебраических уравнений Ax = b (где A - квадратная матрица общего вида, в том числе разреженная, порядка N) с заданной точностью EPS методом сопряженных градиентов, а также сумму квадратов компонент вектора невязки Ax - b .

Н.С.Бахвалов. Численные методы. Изд - во "Наука", 1973.

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

procedure ASS1R(var B :Array of Real; N :Integer; SUBA :Proc_FA;
                SUBAT :Proc_FA; var X :Array of Real; var RSQ :Real;
                ITMAX :Integer; EPS :Real; var G :Array of Real;
                var H :Array of Real; var XI :Array of Real;
                var XJ :Array of Real; var IERR :Integer);

   Параметры

B -	вещественный вектор длины N, содержащий компоненты вектора правой части системы;
N -	порядок матрицы системы (тип: целый);
SUBA -	подпрограмма вычисления вектора Ax, первый оператор которой имеет вид: procedure SUBA (var X :Array of Real; var AX :Array of Real); где X - вещественный вектор длины N, содержащий текущий вектор решения X; AX - вещественный вектор длины N, содержащий результирующий вектор Ax;
SUBAT -	подпрограмма вычисления вектора A^Tx, первый оператор которой имеет вид: procedure SUBAT (var X :Array of Real; var ATX :Array of Real); где X - вещественный вектор длины N, содержащий текущий вектор решения X; ATX - вещественный вектор длины N, содержащий результирующий вектор A^Tx;
X -	вещественный вектор длины N, содержащий на входе в подпрограмму начальное приближение к решению системы, а на выходе - вычисленное решение системы;
RSQ -	вещественная переменная, содержащая вычисленную сумму квадратов компонент вектора невязки Ax - b;
ITMAX -	заданное максимальное допустимое количество итераций метода сопряженных градиентов (тип: целый);
EPS -	заданная точность, с которой необходимо вычислить решение системы (тип: вещественный);
G, H, - XI, XJ	вещественные векторы длины N, используемые в подпрограмме в качестве рабочих;
IERR -	целая переменная, служащая для сообщения об ошибках, обнаруженных в ходе работы подпрограммы; при этом

IERR=65 -	когда матрица системы вырождена;
IERR=66 -	когда заданная точность не может быть достигнута за максимальное количество итераций.

   Версии

ASS1E -

решение систем линейных алгебраических уравнений с разреженными матрицами методом сопряженных градиентов в режиме расширенной (Extended) точности. При этом параметры B, X, RSQ, EPS, G, H, XI, XJ должны иметь тип Extended.

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

UTAS10 -	подпрограмма выдачи диагностических сообщений при работе подпрограммы ASS1R.
UTAS11 -	подпрограмма выдачи диагностических сообщений при работе подпрограммы ASS1D.

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

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

Если значение RSQ велико, то это означает, что матрица системы близка к вырожденной и полученный вектор x представляет собой наилучшее приближение к решению в смысле наименьших квадратов.

В подпрограммах ASS1R и ASS1E имеется глобальная запись (структура данных) с именем _ASS1RR и элементом ITER. По окончании работы этих подпрограмм значение переменной ITER полагается равной количеству выполненных итераций метода сопряженных градиентов.

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

Unit TASS1R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, F1ASS1R_p, F2ASS1R_p, ASS1R_p;

function TASS1R: String;

implementation

function TASS1R: String;
var
N,I,J,ITMAX,_i,IERR :Integer;
R,EPS,RSQ :Real;
B :Array [0..9] of Real;
X :Array [0..9] of Real;
G :Array [0..9] of Real;
H :Array [0..9] of Real;
XI :Array [0..9] of Real;
XJ :Array [0..9] of Real;
label
_1,_2,_3;
begin
Result := '';  { результат функции }
N := 10;
for I:=1 to N do
 begin
  X[I-1] := 0.0;
_1:
  B[I-1] := 1.0;
 end;
R := 0.1;
for I:=1 to N do
 begin
  for J:=1 to N do
   begin
_2:
    _ASR.elm1[(I-1),(J-1)] := 0.0;
   end;
 end;
for I:=1 to N do
 begin
  _ASR.elm1[(I-1),(I-1)] := R;
_3:
  R := R+0.1;
 end;
EPS := 1.E-5;
IТМАХ := 10*N;
ASS1R(B,N,F1ASS1R,F2ASS1R,X,RSQ,ITMAX,EPS,
     G,H,XI,XJ,IERR);
Result := Result + #$0D#$0A;
for _i:=0 to 9 do
 begin
  Result := Result + Format('%20.16f ',[X[_i]]);
  if ( ((_i+1) mod 4)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%20.16f ',[RSQ]) + #$0D#$0A;
Result := Result + Format('%5d ',[_ASS1RR.elm1]) + #$0D#$0A;
UtRes('TASS1R',Result);  { вывод результатов в файл TASS1R.res }
exit;
end;

end.

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

procedure f1ass1r(var X :Array of Real; var V :Array of Real);

implementation

procedure f1ass1r(var X :Array of Real; var V :Array of Real);
var
N,I,J :Integer;
label
_1;
begin
N := 10;
for I:=1 to N do
 begin
  V[I-1] := 0.0;
  for J:=1 to N do
   begin
_1:
    V[I-1] := V[I-1]+_ASR.elm1[(I-1),(J-1)]*X[J-1];
   end;
 end;
end;

end.

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

procedure f2ass1r(var X :Array of Real; var V :Array of Real);

implementation

procedure f2ass1r(var X :Array of Real; var V :Array of Real);
var
N,I,J :Integer;
label
_1;
begin
N := 10;
for I:=1 to N do
 begin
  V[I-1] := 0.0;
  for J:=1 to N do
   begin
_1:
    V[I-1] := V[I-1]+_ASR.elm1[(J-1),(I-1)]*X[J-1];
   end;
 end;
end;

end.

Результаты: 

      X  =  ( 10 , 5 , 10/3 , 2.5 , 2 , 5/3 , 10/7 , 1.25 , 10/9 , 1 )

      RSQ  = 0.19785E - 8 
      ITER = 11