%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   myLU.m
%
% If you put this file in the working directory for MatLab then you will be
% able to access the function myLU()
%    [L, U] = myLU( matrix )
%
% This function attempts to compute the LU factorization for a matrix
% WITHOUT doing any row swaps.  If this is not possible, it throws an
% error.  The input matrix does not need to be square or have pivots on the
% diagonal.  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [L,A] = myLU(A)
numrows = size(A,1);    % record number of rows of A
L = eye(numrows);       % eye( n ) makes an n by n matrix with 1 on diagonal
oldpivot = 0;           % this will be used to check if row swaps were needed

for ( row=1:numrows )   % this is a loop; first row = 1, then row = 2, etc
    pivot = find(A(row,:), 1);      % find( __ , 1) returns index of first nonzero element
    if (pivot < oldpivot)           % if new pivot is left of old pivot we have problem!
        error('Cannot be reduced without row swaps!')
    end
    
    for ( j=(row+1):numrows )       % now we will reduce all rows below the pivot row
        L(j,row) = A(j,pivot)/A(row,pivot);    % first compute the multiplier
        A(j,:) = A(j,:) - L(j,row)*A(row,:);   % then reduce the matrix
    end
    
    if (isempty(pivot))             % if the row was all zero
        oldpivot = size(A,2)+1;       % then the old pivot is off the end of matrix
    else
        oldpivot = pivot;             % otherwise the old pivot is this one
    end
end

end