unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, unit1,Forms, Dialogs;

procedure LUDCMP(var A:matrx2;N:integer;var INDX:array of real;var D:real);
procedure LUBKSB(var A:matrx2; N:integer; var INDX,B:array of real);
procedure MPROVE(A,ALUD:matrx2; N:integer;
            INDX:array of real; B:array of real; var X:array of real);

implementation
procedure LUDCMP(var A:matrx2;N:integer;var INDX:array of real;var D:real);
const
    NMAX = 100;    TINY = 1E-20;
var
    VV:array[0..100] of real;
    I,J,K,IMAX:Integer;
    AAMAX,SUM,DUM:Real;
begin
    D:=1;
    For I:=1 To N do
    begin
        AAMAX:=0;
        For J:=1 To N do
            If (Abs(A[I, J]) > AAMAX) Then AAMAX:=Abs(A[I, J]);
        If AAMAX = 0 Then
        begin
            ShowMessage('Singular matrix');
            exit;
        end;
        VV[I]:= 1 / AAMAX;
    end;
    For J:=1 To N do
    begin
        If J > 1 Then
        begin
            For I:=1 To J - 1 do
            begin
                Sum:=A[I, J];
                If I > 1 Then
                begin
                    For K:=1 To I - 1 do
                        Sum:=Sum - A[I, K] * A[K, J];
                    A[I, J]:=Sum;
                end;
            end;
        end;
        AAMAX:=0;
        For I:=J To N do
        begin
            Sum:=A[I, J];
            If J > 1 Then
            begin
                For K:=1 To J - 1 do
                    Sum:=Sum - A[I, K] * A[K, J];
                A[I, J]:=Sum;
            end;
            DUM:=VV[I] * Abs(Sum);
            If DUM >= AAMAX Then
            begin
                IMAX:=I;
                AAMAX:=DUM;
            end;
        end;
        If J <> IMAX Then
        begin
            For K:=1 To N do
            begin
                DUM:=A[IMAX, K];
                A[IMAX, K]:=A[J, K];
                A[J, K]:=DUM;
            end;
            D:=-D;
            VV[IMAX]:=VV[J];
        end;
        INDX[J]:=IMAX;
        If J <> N Then
        begin
            If A[J, J] = 0 Then A[J, J]:=TINY;
            DUM:=1 / A[J, J];
            For I:=J + 1 To N do
                A[I, J]:=A[I, J] * DUM;
        end;
    end;
    If A[N, N]=0 Then A[N, N]:=TINY;
end;

procedure LUBKSB(var A:matrx2; N:integer; var INDX,B:array of real);
var
    I,LL,J,II:integer;
    SUM:real;
BEGIN
    II:=0;
    For I:=1 To N DO
    begin
        LL:=Trunc(INDX[I]);
        Sum:=B[LL];
        B[LL]:=B[I];
        If II <> 0 Then
        begin
            For J:=II To I - 1 do
                Sum:=Sum - A[I, J] * B[J];
        end
        Else If Sum <> 0 Then
            II:=I;
        B[I]:=Sum;
    end;
    For I:=N DownTo 1 do
    begin
        Sum:=B[I];
        If I < N Then
        begin
            For J:=I + 1 To N do
                Sum:=Sum - A[I, J] * B[J];
        end;
        B[I]:=Sum / A[I, I];
    end;
End;

procedure MPROVE(A,ALUD:matrx2; N:integer;
            INDX:array of real; B:array of real;var X:array of real);
var
    R:array [0..100] of real;
    I,J:integer;
    SDP:real;
begin
    For I:=1 To N do
    begin
        SDP:=-B[I];
        For J:=1 To N do
            SDP:=SDP + A[I, J] * X[J];
        R[I]:=SDP;
    end;
    LUBKSB(ALUD, N, INDX, R);
    For I:=1 To N do
        X[I]:=X[I] - R[I];
end;
end.