% Triangular Decomposition:  A=LU
%
% version 1.1  November 21, 2011
%
% The input matrix must be square.  It may be unsymmetric.
%
% Matlab already has an lu decomposition function: [L,U] = lu(A)
%
% This script demonstrates the underlying steps for educational purposes.
%
disp(' Use Pivot?')
ip = input(' 1=yes 2=no ');
%
%
clear AA;
clear A;
clear G;
clear L;
clear U;
% 
disp(' ');
disp(' Select file input method ');
disp('   1=file preloaded into Matlab ');
disp('   2=Excel file ');
file_choice = input('');
%
if(file_choice==1)
        AA = input(' Enter the matrix name:  ');
end
if(file_choice==2)
        [filename, pathname] = uigetfile('*.*');
        xfile = fullfile(pathname, filename);
%        
        AA = xlsread(xfile);
%         
end
%
A=AA;
nn=size(A);
n=nn(1,1);
%
L=zeros(n,n);
U=zeros(n,n);
%
for i=1:n
   L(i,i)=1.; 
end
%
i=1;
if(ip==1)
    [L,A]=pivot(L,A,i,n);
end
%
U(1,1)=A(1,1);
for i=2:n
    L(i,1)=A(i,1)/U(1,1);
    U(1,i)=A(1,i);
end
%
for i=2:n
%
    if(ip==1)
        [L,A]=pivot(L,A,i,n);         
    end
%
    sum=0.;
    for k=1:i-1
        sum=sum+U(k,i)*L(i,k);
    end
    U(i,i)=A(i,i)-sum;
%
    for j=i+1:n
%        
        sum=0.;
        for k=1:j-1
            sum=sum+U(k,i)*L(j,k);
        end
        L(j,i)=(A(j,i)-sum)/U(i,i);
%
        sum=0.;
        for k=1:j-1
            sum=sum+U(k,j)*L(i,k);
        end
        U(i,j)=(A(i,j)-sum);
    end
end
A
L
U
L*U