unit InvertMatrix;

interface InvertMatrixGauss;

implementation

(****************************************************************
Обращение матрицы алгоритмом Гаусса.

Принимаемые параметры:
    N :        - размер матрицы
    A :        - матрица [1..N, 1..N]
    Epsilon    - Погрешность   расчётов.   Если   число по модулю
                 меньше Epsilon,  то оно считается нолем.

Результат:
    Значение функции  - True при успехе.
    A                 - обращенная матрица при успехе
****************************************************************)
function InvertMatrixGauss(
             const N       : Integer;
             var   A       : array of array of real;
             const Epsilon : Real):Boolean;
var
  b : array of Real;
  c : array of Real;
  d : Real;
  i : Integer;
  j : Integer;
  k : Integer;
  p : Integer;
  w : Real;
  y : Real;
  z : array of Integer;
  
begin
  d:=1;
  Result := True;

  SetBounds( b, [1,n] );
  SetBounds( c, [1,n] );
  SetBounds( z, [1,n] );

  for i:=1 to n do
  begin
    z[i]:=i;
  end;

  i:=1;
  repeat
    k:=i;
    y:=a[i,i];
    if i+1<=n then
    begin
      j:=i+1;
      repeat
        w:=a[i,j];
        if AbsReal(w)>AbsReal(y)
        then
        begin
          k:=j;
          y:=w;
        end;
        j:=j+1
      until not(j<=n);
    end;
    d:=d*y;
    if AbsReal(y)<=Epsilon then
    begin
      Result := False;
      i:=n+1;
    end
    else
    begin
      y:=1/y;
      for j:=1 to n do
      begin
        c[j]:=a[j,k];
        a[j,k]:=a[j,i];
        a[j,i]:=-c[j]*y;
        b[j]:=a[i,j]*y;
        a[i,j]:=b[j];
      end;
      j:=z[i];
      z[i]:=z[k];
      z[k]:=j;
      a[i,i]:=y;
      for k:=1 to n do
      begin
        if k<>i
        then
        begin
          for j:=1 to n do
          begin
            if j<>i then
            begin
              a[k,j]:=a[k,j]-b[j]*c[k];
            end;
          end;
        end;
      end;
    end;
    i:=i+1
  until not(i<=n);


  if Result then
  begin
    for i:=1 to n do
    begin
      k:=z[i];
      while k<>i do
      begin
        for j:=1 to n do
        begin
          w:=a[i,j];
          a[i,j]:=a[k,j];
          a[k,j]:=w;
        end;
        p:=z[i];
        z[i]:=z[k];
        z[k]:=p;
        d:=-d;
        k:=z[i];
      end;
    end;
  end;
end;

end.