function [ L,U ] = LUdecomp(A)
%LUdecomp: decompose square matrix A as A=LU
% where L is lower triag and U is upper triag
[m, n]=size(A);
if m ~= n, error('Input must be a square matrix.'), end
L=zeros(n); U=zeros(n);
% remember A^(0) is A
AofK = A;
for k = 1:n
% at this point AofK is A^(k-1)
for j = k:n
U(k,j) = AofK(k,j);
end
% check that we don't divide by zero(!)
if U(k,k) == 0
error('** A^(k-1)_{k,k}==0 in LU decomp')
end
for i = k:n
L(i,k) = AofK(i,k)/U(k,k);
end
% now modify AofK so that we can use it in the
% next iteration
for i = k:n
for j = k:n
AofK(i,j) = AofK(i,j) - L(i,k)*U(k,j);
end
end
end
end