program ludec;     { LU Decomposition in Place }
{                    Jeff Crawford     TRW     3-21-86

References:   "Computer Methods for Circuit Analysis and Design", Vlach
              Singhal;  Van Nostrand Reinhold Company, 1983.

              "Applied Numerical Analysis", Curtis F. Gerald, Patrick
              O. Wheatley, Addison-Wesley Publishing Co., 1984.

This program is unique compared to other LU Decompositions seen in the
references above.  Due to the nature of the algorithm, it is possible to
do the reduction "in-place" rather than identifying two additional matrices
L & U, thus using more memory.

This program will readily allow computations on complex systems by
structuring the input matrix as shown below:

          ! a11 a12  -b11 -b12 !  ! c11 c12 !  =  ! x11 x12 !
          ! a21 a22  -b21 -b22 !  ! c21 c22 !  =  ! x21 x22 !
          ! b11 b12   a11  a12 !  ! d11 d12 !  =  ! y11 y12 !
          ! b21 b22   a21  a22 !  ! d21 d22 !  =  ! y21 y22 !

where in reality the equations are of the form

! a11 +jb11  a12 + jb12 ! ! c11 + jd11  c12 + jd12 ! = ! x11 +jy11  x12 +jy12 !
! a21 +jb21  a22 + jb22 ! ! c21 + jd21  c22 + jd22 !   ! x21 +jy21  x22 +jy22 !
                                                                              }
type     mat1 = array[1..20,1..20] of real;

var      A         : mat1;
         i,j,m,k   : integer;
         n         : integer;
         Sum       : real;
         B,X       : array[1..10] of real;

begin

clrscr;
textcolor(10);
gotoxy(10,10);
write('Enter Order of System  ');
read(n);
for i:= 1 to n do
    begin
    for j:= 1 to n do
         begin
         gotoxy(8,10);
         clreol;
         write('A[',i:2,' ',j:2,' ] = ');
         read(A[i,j]);
         end;
    gotoxy(8,10);
    clreol;
    write('         B[',i,'] = ');
    read(B[i]);
    end;
clrscr;
for j:= 2 to n do  A[1,j]:= A[1,j]/A[1,1];
for i:= 2 to n do                      { I is Main Outer Loop Counter }
    begin
    for k:= i to n do                  { K is Which Row of L Working on }
         begin
         Sum:= 0.0;
         for j:= 1 to i-1 do           { J is # of Terms to Sum }
              begin
              Sum:= Sum + A[k,j]*A[j,i];
              end;
         A[k,i]:= A[k,i] - Sum;
    end;

    for m:= i+1 to n do                { M = Which Column Row Element is in }
         begin
         Sum:= 0.0;
         for j:= 1 to i-1 do  Sum:= Sum + A[i,j]*A[j,m];
         A[i,m]:= (A[i,m] - Sum) / A[i,i];
         end;
    end;

B[1]:= B[1]/A[1,1];                    { Forward Substitution }
for i:= 2 to n do
    begin
    Sum:= 0.0;
    for k:= 1 to i-1 do      Sum:= Sum + A[i,k]*B[k];
    B[i]:= ( B[i] - Sum ) / A[i,i];
    end;
X[n]:= B[n];
j:= n-1;
while j >= 0 do
    begin
    sum:= 0.0;
    for k:= j+1 to n do      Sum:= Sum + A[j,k]*X[k];
    X[j]:= B[j] - Sum;
    j:= j - 1;
    end;
clrscr;
gotoxy(1,4);
writeln('                           Solution Vector   ');
writeln('                       LU Decomposition Method ');
writeln;
for i:= 1 to n do writeln('                         X[ ',i,' ] = ',X[i]:10:6);
end.