unit LU_Utils;
interface
uses Classes;
type
 TDblArr = array of Double;
 TIntArr = array of Integer;
 TDblMatr = array of array of Double;
procedure CopyMatr(Src, Dest: TDblMatr);
//исходная матрица разрушается, при необходимости нужно сделать копию
procedure LUDecomp(Matr: TDblMatr; Indx: TIntArr; var DSign: Integer);
procedure LUBackSubst(Matr: TDblMatr; Indx: TIntArr; Reslt: TDblArr);
//решение системы линейных уравнений. XY на входе - правые части,
//на выходе - вектор решений
//LU-разложение матрицы может использоваться многократно с
//разными правыми частями (при модификации кода)
procedure SLAUSolver(Matr: TDblMatr; XY: TDblArr);
//исходная матрица заменяется обратной
procedure MatrInverse(Matr: TDblMatr);
//исходная матрица сохраняется
function Determinant(Matr: TDblMatr): Double;
implementation
uses SysUtils;
procedure CopyMatr(Src, Dest: TDblMatr);
var
 n, i: Integer;
begin
 n := Length(Src);
 for i := 0 to n - 1 do
  Move(Src[i, 0], Dest[i, 0], n * SizeOf(Double));
end;
procedure LUDecomp(Matr: TDblMatr; Indx: TIntArr; var DSign: Integer);
const
 Tiny = 1.0E-40;
var
 n, k, j, imax, i: Integer;
 sum, dum, big: Double;
 VV: TDblArr;
 DumPtr: Pointer;
begin
 n := Length(Matr);
 SetLength(VV, n);
 DSign := 1;
 for i := 0 to n - 1 do begin
  big := 0;
  for j := 0 to n - 1 do
  if (abs(Matr[i, j]) > big) then
  big := Abs(Matr[i, j]);
  if (big = 0.0) then
  raise Exception.Create('Stopped. Matrix is singular!!!');
  VV[i] := 1 / big;
 end;
 for j := 0 to n - 1 do begin
  for i := 0 to j - 1 do begin
  sum := Matr[i, j];
  for k := 0 to i - 1 do
  sum := sum - Matr[i, k] * Matr[k, j];
  Matr[i, j] := sum;
  end;
  big := 0;
  for i := j to n - 1 do begin
  sum := Matr[i, j];
  for k := 0 to j - 1 do
  sum := sum - Matr[i, k] * Matr[k, j];
  Matr[i, j] := sum;
  dum := VV[i] * abs(sum);
  if (dum > big) then begin
  big := dum;
  imax := i;
  end;
  end;
  if (j <> imax) then begin
  DumPtr := Pointer(Matr[imax]);
  Pointer(Matr[imax]) := Pointer(Matr[j]);
  Pointer(Matr[j]) := DumPtr;
  DSign := -DSign;
  VV[imax] := VV[j]
  end;
  Indx[j] := imax;
  if (Matr[j, j] = 0) then
  Matr[j, j] := Tiny;
  if (j <> n - 1) then begin
  dum := 1 / Matr[j, j];
  for i := j + 1 to n - 1 do
  Matr[i, j] := Matr[i, j] * dum;
  end;
 end;
end;
procedure LUBackSubst(Matr: TDblMatr; Indx: TIntArr; Reslt: TDblArr);
var
 n, j, ip, ii, i: Integer;
 sum: Double;
begin
 n := Length(Matr);
 ii := -1;
 for i := 0 to n - 1 do begin
  ip := Indx[i];
  sum := Reslt[ip];
  Reslt[ip] := Reslt[i];
  if (ii <> -1) then
  for j := ii to i - 1 do
  sum := sum - Matr[i, j] * Reslt[j]
  else
  if (sum <> 0.0) then
  ii := i;
  Reslt[i] := sum;
 end;
 for i := n - 1 downto 0 do begin
  sum := Reslt[i];
  if (i < n - 1) then
  for j := i + 1 to n - 1 do
  sum := sum - Matr[i, j] * Reslt[j];
  Reslt[i] := sum / Matr[i, i];
 end
end;
procedure SLAUSolver(Matr: TDblMatr; XY: TDblArr);
var
 Indx: TIntArr;
 DSign: Integer;
begin
 SetLength(Indx, Length(XY));
 LUDecomp(Matr, Indx, DSign);
 LUBackSubst(Matr, Indx, XY);
end;
procedure MatrInverse(Matr: TDblMatr);
var
 n, i, j: Integer;
 Indx: TIntArr;
 Col: TDblArr;
 DSign: Integer;
 Temp: TDblMatr;
begin
 n := Length(Matr);
 SetLength(Indx, n);
 SetLength(Col, n);
 SetLength(Temp, n, n);
 LUDecomp(Matr, Indx, DSign);
 for j := 0 to n - 1 do begin
  for i := 0 to n - 1 do
  Col[i] := 0;
  Col[j] := 1;
  LUBackSubst(Matr, Indx, Col);
  for i := 0 to n - 1 do
  Temp[i, j] := Col[i];
 end;
 CopyMatr(Temp, Matr);
 Temp := nil;
end;
function Determinant(Matr: TDblMatr): Double;
var
 n, j: Integer;
 Indx: TIntArr;
 DSign: Integer;
 Temp: TDblMatr;
begin
 n := Length(Matr);
 SetLength(Indx, n);
 SetLength(Temp, n, n);
 CopyMatr(Matr, Temp);
 LUDecomp(Temp, Indx, DSign);
 Result := DSign;
 for j := 0 to n - 1 do
  Result := Result * Temp[j, j];
end;
end.

unit Unit1;
interface
uses
 Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
 Dialogs, LU_Utils, StdCtrls;
type
 TForm1 = class(TForm)
  Memo1: TMemo;
  Button1: TButton;
  Button2: TButton;
  Button3: TButton;
  procedure Button1Click(Sender: TObject);
  procedure Button2Click(Sender: TObject);
  procedure Button3Click(Sender: TObject);
 private
  { Private declarations }
 public
  { Public declarations }
 end;
var
 Form1: TForm1;
implementation
{$R *.dfm}
function ShowArr(Arr: TDblArr): string;
var
 i: Integer;
begin
 Result := '';
 for i := 0 to Length(Arr) - 1 do
  Result := Result + FormatFloat('0.0000', Arr[i]) + ' ';
end;
procedure ShowMatr(Matr: TDblMatr; Strings: TStrings);
var
 n, i, j: Integer;
 s: string;
begin
 n := Length(Matr);
 for i := 0 to n - 1 do begin
  s := '';
  for j := 0 to n - 1 do
  s := s + FormatFloat('0.0000', Matr[i, j]) + ' ';
  Strings.Add(s);
 end;
end;
procedure TForm1.Button1Click(Sender: TObject);
var
 M:TDblMatr;
 B:TDblArr;
begin
SetLength(M,3,3);
SetLength(B,3);
 M[0, 0] := 0.5;
 M[0, 1] := 2;
 M[0, 2] := -1;
 M[1, 0] := -15;
 M[1, 1] := -0.1;
 M[1, 2] := 3.2;
 M[2, 0] := 0.7;
 M[2, 1] := 4.1;
 M[2, 2] := -0.2;
 B[0] := 11.4;
 B[1] := -196.02;
 B[2] := 10.22;
 Memo1.Clear;
 ShowMatr(M, Memo1.Lines);
 Memo1.Lines.Add('');
 SLAUSolver(M, B);
 ShowMatr(M, Memo1.Lines);
 Memo1.Lines.Add('');
 Memo1.Lines.Add(ShowArr(B));
end;
procedure TForm1.Button2Click(Sender: TObject);
var
 M:TDblMatr;
begin
 SetLength(M,3,3);
 M[0, 0] := 1;
 M[0, 1] := 1;
 M[0, 2] := 1;
 M[1, 0] := 1;
 M[1, 1] := 2;
 M[1, 2] := 1;
 M[2, 0] := 1;
 M[2, 1] := 1;
 M[2, 2] := 3;
 Memo1.Clear;
 ShowMatr(M, Memo1.Lines);
 Memo1.Lines.Add('');
 MatrInverse(M);
 ShowMatr(M, Memo1.Lines);
end;
procedure TForm1.Button3Click(Sender: TObject);
var
 M:TDblMatr;
begin
 SetLength(M,3,3);
 M[0, 0] := 1;
 M[0, 1] := 0;
 M[0, 2] := 1;
 M[1, 0] := 0;
 M[1, 1] := 2;
 M[1, 2] := 0;
 M[2, 0] := 1;
 M[2, 1] := 0;
 M[2, 2] := 3;
 Memo1.Clear;
 ShowMatr(M, Memo1.Lines);
 Memo1.Lines.Add('');
 Memo1.Lines.Add(FormatFloat('0.0000', Determinant(M)));
end;
end.