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;