function [L,U] = lu_decompose(A)
%   lu decompose
%   L:下三角矩阵
%	U:上三角矩阵
%	A:输入矩阵
[m,n]=size(A);

L=eye(m);
L(:,1)=A(:,1)/A(1,1);%L第一列赋值

U=zeros(m,n);
U(1,:)=A(1,:);%U第一行赋值

for i=2:m
    for j=2:n
        if i<=j
            U(i,j)=A(i,j)-sum(L(i,1:i-1).*U(1:i-1,j)');%递推表达式(1-6)
        else
            if U(j,j)==0
                L(i,j)=0;
            else
                L(i,j)=(A(i,j)-sum(L(i,1:j-1).*U(1:j-1,j)'))/U(j,j);%递推表达式(1-7)
            end
        end
    end
end