Текст подпрограммы и версий
asg8r_p.zip , asg8e_p.zip , asg8c_p.zip 		Тексты тестовых примеров
tasg8r_p.zip , tasg8e_p.zip , tasg8c_p.zip

Подпрограмма:  ASG8R (модуль ASG8R_p)

Назначение

Решение вещественной системы линейных алгебраических уравнений А*x = b или AT*x = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы A.

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

Для заданной вещественной квадратной матрицы А порядка N выполняется треугольная факторизация А = L*U, где L - нижняя треугольная матрица, а U - верхняя треугольная матрица, вычисляется величина, обратная числу обусловленности матрицы А : 

     RCOND = 1/(|| A ||1*|| A-1||1) ,
 гдE
     || A ||1 = maxj=1,...,N   ( | a1j | + | a2j | + ... + | aNj | ) .
 Затем решается система
       L*y = b
 (для системы AT*x = b решается  UT*y = b) ,
 и затем решение  x находится как решение системы 
       U*x = y
 (для AT*x = b решается LT*x= y).

Дж.Форсайт, М.Малькольм, К.Моулер. Машинные методы математических вычислений. М.: Мир, 1980.

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

procedure ASG8R(var A :Array of Real; M :Integer; N :Integer;
                var NLEAD :Array of Integer; var B :Array of Real;
                LTR :Integer; L :Integer; var RCOND :Real;
                var X :Array of Real; var IERR :Integer);


   Параметры
A - вещественный двумерный массив размера М на N, в котором задается исходная матрица; на выходе из подпрограммы на соответствующих местах находятся элементы матрицы U и поддиагональные элементы матриц исключения метода Гаусса  Li , i = 1, ..., N - 1;
M - первая размерность массива А в вызывающей программе (тип: целый);
N - порядок матрицы системы (тип: целый);
NLEAD - целый вектор длины N, содержащий на выходе информацию о выполненных в ходе факторизации перестановках (см. замечания по использованию);
B - вещественный вектор длины N, в котором задается правая часть системы;
LTR - задает вид решаемой системы (тип: целый):
LТR = 0 - если решается система А*x = b,
LТR ≠ 0 - если решается система АT*x = b;
L - признак решаемой системы (тип: целый);
L = 0 - если ситема с данной матрицей решается впервые,
L ≠ 0 - если система с заданной матрицей решается повторно;
RCOND - вещественная переменная, содержащая на выходе вычисленное значение величины, обратной числу обусловленности матрицы А (см. замечания по использованию);
X - вещественный вектор длины N, содержащий на выходе решение системы (см. замечания по использованию);
IERR - целая переменная, содержащая на выходе информацию о прохождении счета; при этом:
IЕRR=65 - если М ≤ 0 или N ≤ 0 ;
IЕRR=66 - если в процессе счета возникло переполнение (это говорит о том, что либо ||А||1, либо некоторые элементы матрицы U, либо некоторые компоненты решения системы превосходят по абсолютной величине максимальное представимое на данной машине число);
IЕRR=-К - если в результате факторизации в К - ой строке матрицы U диагональный элемент равен нулю (это свидетельствует о вырожденности матрицы А). Если таких строк у матрицы U несколько, то значение К полагается равным номеру последней из них;
IЕRR=67 - если система несовместна (см. замечания по использованию). 
   Версии
ASG8E - решение системы линейных алгебраических уравнений А*x = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы А для вещественных А и b, заданных с расширенной (Extended) точностью.
ASG8C - решение системы линейных алгебраических уравнений А*x = b методом Гаусса с выбором ведущего элемента по столбцу и оценка числа обусловленности матрицы А для комплексных А и b. 
   Вызываемые подпрограммы
AFG4R - подпрограмма треугольной факторизации и оценки числа обусловленности матрицы А.
UTAFSI - подпрограмма выдачи диагностических сообщений. 
   Замечания по использованию
  1. 

В подпрограмме АSG8E массивы А, В и Х и переменная RСОND имеют тип Extended, для треугольного разложения и оценки числа обусловленности матрицы А вызывается подпрограмма АFG4E.

  2. 

В подпрограмме АSG8С массивы А, В и Х имеют тип Complex, для треугольного разложения и оценки числа обусловленности матрицы А вызывается подпрограмма АFG4С.

  3. 

На выходе К - й элемент вектора NLЕАD равен номеру строки, переставленной на К - м шаге факторизации с К - ой строкой матрицы А. Так как факторизация Гаусса требует N - 1 шагов, то NLЕАD (N) = N.

  4. 

Если задано L ≠ 0, т.е. система решается повторно, то перед выполнением подпрограммы нужно запомнить, если требуется, вычисленную ранее оценку числа обусловленности, а в качестве матрицы А при обращении к подпрограмме нужно взять результат предыдущего обращения к АSG8R или к АFG4R, т.е. матрица А должна быть уже факторизована, т.к. при L ≠ 0 треугольная факторизация не выполняется, а переменной RСОND присваивается значение ноль.

  5.  Если вырабатывается значение переменной IЕRR, отличное от нуля, то выдается соответствующее диагностическое сообщение, и, если IЕRR ≥ 65, то происходит выход из подпрограммы. Если система совместна, но матрица А вырождена (IЕRR < 0), т.е. для некоторых номеров К  U (К, К) = 0, то Х (К) полагается равным 1.

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

Unit TASG8R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, ASG8R_p;

function TASG8R: String;

implementation

function TASG8R: String;
var
M,N,LTR,L,J,I,JI,JJ,II,IERR :Integer;
RCOND :Real;
NLEAD :Array [0..3] of Integer;
Х :Array [0..3] of Real;
const
A :Array [0..15] of Real = ( 7.9,8.5,4.3,3.2,5.6,-4.8,4.2,-1.4,5.7,0.8,-3.2,
-8.9,-7.2,3.5,9.3,3.3 );
B :Array [0..3] of Real = ( 7.0,7.0,7.0,7.0 );
begin
RЕSult := '';  { результат функции }
M := 4;
N := M;
LTR := 0;
L := 0;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   A=' + #$0D#$0A]);
for I:=1 to M do
 begin
  for J:=1 to N do
   begin
    Result := Result + Format('   %20.16f ',
 [A[(I-1)+(J-1)*4]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
ASG8R(A,M,N,NLEAD,B,LTR,L,RCOND,X,IERR);
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   NLEAD=' + #$0D#$0A]);
for JI:=1 to N do
 begin
  Result := Result + Format('   %3d ',[NLEAD[JI-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   A=' + #$0D#$0A]);
for JJ:=1 to N do
 begin
  for II:=1 to M do
   begin
    Result := Result + Format('   %20.16f ',
 [A[(II-1)+(JJ-1)*4]]) + #$0D#$0A;
   end;
 end;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   IERR=']);
Result := Result + Format('%3d ',[IERR]);
Result := Result + Format('%s',['     RCOND=']);
Result := Result + Format('%20.16f ',[RCOND]) + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   B=' + #$0D#$0A]);
for J:=1 to N do
 begin
  Result := Result + Format('   %20.16f ',[B[J-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
Result := Result + #$0D#$0A;
Result := Result + Format('%s',['   X=' + #$0D#$0A]);
for J:=1 to N do
 begin
  Result := Result + Format('   %20.16f ',[X[J-1]]) + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TASG8R',Result);  { вывод результатов в файл TASG8R.res }
exit;
end;

end.


Результат:

               |   8.5           -4.8            0.8               3.5        |
               | -0.92941    10.06118     4.95647     -10.45294 |
      A  =  | -0.50588    -0.65879    -9.40171       2.40526  |
               | -0.37647    -0.04046    -0.73072      12.65817 |

      NLEAD  =   (2, 2, 4, 4)

      RCOND  =  0.41764

               |  1.023054 |
               |  0.273777 |
      X  =  | -0.462957 |
               | -0.003275 |