function [L,U,y,x]=c2017210364(A,b)
% 判断是否为方阵
[m,n]=size(A);
if m~=n
disp('矩阵必须是方阵');
return;
end
% 判断各阶顺序主子式是否为0，其中的任何一个顺序主子式为0，则不能进行Doolittle分解,f=1
n=length(A);
f=0;
A1=A;
for i=1:n
if(det(A1(1:i,1:i))==0)
f=1;
end
end
if(f==1)
disp('输入的方阵不能进行Doolittle分解，请重新输入一个方阵')
return;
end
% % 第一种对矩阵进行L U 分解
% U=zeros(m);
% L=zeros(m);
% for j=1:m
% L(j,j)=1;
% end
% for j=1:m
% U(1,j)=A(1,j);
% end
% for i=2:m
% for j=1:m
% for k=1:i-1
% s1=0;
% if k==1
% s1=0;
% else
% for p=1:k-1
% s1=s1+L(i,p)*U(p,k);
% end
% end
% L(i,k)=(A(i,k)-s1)/U(k,k);
% end
% for k=i:m
% s2=0;
% for p=1:i-1
% s2=s2+L(i,p)*U(p,k);
% end
% U(i,k)=A(i,k)-s2;
% end
% end
% end
% disp('下三角矩阵')
% L
% disp('上三角矩阵')
% U
%对LU分解并通过LU求出解
n=length(b);%方程个数n
x=zeros(n,1);%未知向量
A(2:n,1)=A(2:n,1)./A(1,1);
for i=2:n-1
A(i,i)=A(i,i)-sum(A(i,1:i-1)'.*A(1:i-1,i));
for j=i+1:n
A(i,j)=A(i,j)-sum(A(i,1:i-1)'.*A(1:i-1,j));
A(j,i)=(A(j,i)-sum(A(j,1:i-1)'.*A(1:i-1,i)))/A(i,i);
end
end
A(n,n)=A(n,n)-sum(A(n,1:n-1)'.*A(1:n-1,n));
% A;
U=A;
L=A;
for i=1:n
L(i,i)=1;
end
for i=1:n-1
for j=i+1:n
L(i,j)=0;
end
end
for i=2:n
for j=1:i-1
U(i,j)=0;
end
end
%用LU分解求解线性方程组
y=zeros(n,1);
y(1)=b(1);
for i=2:n
y(i)=b(i)-L(i,1:i-1)*y(1:i-1);
end
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);
end
