Решение системы уравнений методом Гаусса.

Любят у нас преподаватели давать первокурсникам эту задачку: решение системы линейных алгебраических уравнений Ax=b методом Гаусса. И метод Гаусса можно освоить и полезное дело сделать. Если будет у вас такая программа, любую задачу методом Гаусса сможете решить легко. Для тех кто не в курсе - программа на Паскале.

program Slau;
uses crt;
const size=30; {максимально допустимая размерность}
type matrix=array [1..size,1..size+1] 
 of real;
type vector=array [1..size] of real;

function GetNumber (s:string;
 a,b:real):real;
{Ввод числа из интервала a,b. 
 Если a=b, то число любое}
var n:real;
begin
 repeat
 write (s);
 {$I-}readln (n);{$I+}
 if (IoResult<>0) then 
 writeln ('Введено не число!')
 else if (ab)) then
 writeln ('Число не в интервале от ',
 a,' до ',b)
 else break;
 until false;
 GetNumber:=n;
end;

procedure GetMatrix (n,m:integer; 
 var a:matrix); {ввод матрицы}
var i,j:integer; si,sj: string [3];
begin
 for i:=1 to n do begin
 str (i,si);
 for j:=1 to m do begin
 str (j,sj);
 a[i,j]:=GetNumber ('a['+ si+ ','+ sj+
 ']=', 0,0);
 end;
 end;
end;

procedure GetVector (n:integer; 
 var a:vector); {ввод вектора}
var i:integer; si:string [3];
begin
 for i:=1 to n do begin
 str (i,si);
 a[i]:=GetNumber ('b['+si+']=',0,0);
 end;
end;

procedure PutVector (n:integer; 
 var a:vector); {вывод вектора}
var i:integer;
begin
 writeln;
 for i:=1 to n do writeln (a[i]:10:3);
end;

procedure MV_Mult (n,m:integer;
 var a:matrix;var x,b:vector);
{умножение матрицы на вектор}
var i,j:integer;
begin
 for i:=1 to n do begin
 b[i]:=0;
 for j:=1 to m do b[i]:=b[i]+a[i,j]*x[j];
 end;
end;

function Gauss (n:integer; var a:matrix;
 var x:vector):boolean;
{метод Гаусса решения СЛАУ}
{a - расширенная матрица системы}
const eps=1e-6; {точность расчетов}
var i,j,k:integer;
 r,s:real;
begin
 for k:=1 to n do begin {перестановка 
 для диагонального преобладания}
 s:=a[k,k];
 j:=k;
 for i:=k+1 to n do begin
 r:=a[i,k];
 if abs(r)>abs(s) then begin
 s:=r;
 j:=i;
 end;
 end;
 if abs(s)k then
 for i:=k to n+1 do begin
 r:=a[k,i];
 a[k,i]:=a[j,i];
 a[j,i]:=r;
 end; {прямой ход метода}
 for j:=k+1 to n+1 do a[k,j]:=a[k,j]/s;
 for i:=k+1 to n do begin
 r:=a[i,k];
 for j:=k+1 to n+1 do 
 a[i,j]:=a[i,j]-a[k,j]*r;
 end;
 end;
 if abs(s)>eps then begin {обратный ход}
 for i:=n downto 1 do begin
 s:=a[i,n+1];
 for j:=i+1 to n do s:=s-a[i,j]*x[j];
 x[i]:=s;
 end;
 Gauss:=true;
 end
 else Gauss:=false;
end;

var a,a1:matrix;
 x,b,b1:vector;
 n,i,j:integer;

begin
 n:=trunc(GetNumber 
 ('Введите размерность матрицы: ',2,size));
 GetMatrix (n,n,a);
 writeln ('Ввод правой части:');
 GetVector (n,b);
 for i:=1 to n do begin 
 {делаем расширенную матрицу}
 for j:=1 to n do a1[i,j]:=a[i,j];
 a1[i,n+1]:=b[i];
 end;
 if Gauss (n,a1,x)=true then begin
 write ('Решение:');
 PutVector (n,x);
 write ('Проверка:');
 MV_Mult (n,n,a,x,b1);
 PutVector (n,b1);
 end
 else write ('Решения нет');
 reset (input); readln;
end.