矩阵LU及PA=LU分解
clear all;
%************************************************************************************
%A=LU分解,非部分主元法
A=input('请输入一个矩阵\n');
U=double(A);
[m,n]=size(U);
L=eye(m);
EPSINON=10^(-15);
%************************************************************************************
%判断是不是方阵，不是方阵则不能进行A=L*U或P*A=L*U分解
if m~=n
fprintf('A不是方阵，不能进行LU分解;\n');
else fprintf('A是方阵;\n');
%判断A的行列式是否为0，为0则不能分解;不为0则进行下一步判断，如果各阶顺序主
%子式均为0，则可直接进行A=L*U分解，否则进行P*A=L*U分解
if det(A)==0
fprintf('A的行列式为零，不能进行分解.\n');
%***********************************************************************************
%判断能不能直接进行A=L*U分解，能分解则用户输入1可得到A=L*U,输入2可得P*A=L*U.
else total=0;
fprintf('A的行列式不为零，可以进行分解;\n');% total用来计数A的顺序主子式中不为0的阶数，如果为n，则能直接进行A=L*U分解。
for j=1:m
if det(A(1:j,1:j))~=0;
total=total+1;
end
end
if total==n
fprintf('A的各阶顺序主子式均不为0,能直接进行A=L*U分解；\n直接进行A=L*U分解请输入1，进行P*A=L*U分解请输入2:\n');
x=input('x=');
else x=2;
fprintf('A存在为0的顺序主子式，不能直接进行A=L*U分解；以下直接进行P*A=L*U分解：\n\n');
end
%**************************************************************************
switch(x)
%选择输入1进行A=L*U分解，或输入2进行P*A=L*U分解；case 1,表示进行A=L*U分解
case 1
for i=1:m %i表示行
L(i+1:m,i)=U(i+1:m,i)/U(i,i);
%U(i+1:m,i:n)=U(i+1:m,i:n)-U(i+1:m,i)*U(i,i:n)/U(i,i);
U(i+1:m,i:n)=U(i+1:m,i:n)-L(i+1:m,i)*U(i,i:n);
end
U(abs(U)
if sum(sum(L==L_lu))==m^2&&sum(sum(U==U_lu))==m^2
fprintf('与Matlab自带公式[L,U,P]=lu(A)所得结果一致。\n');
else fprintf('与[L,U]=lu(A)所得结果不同；因为本A=L*U程序为非部分主元法，而[L_lu,U_lu]=lu(A)是采用\n部分主元法得到；即在计算过程中[L_lu,U_lu]=lu(A)公式中的U_lu也进行了行变换\n');
%用difference_value来表示P*A与L*U的完全相等的程度
difference_value=A-L*U;
end
L,U
%case 2表示P*A=L*U分解
case 2
P=eye(m);
L1=zeros(m);
U=double(A);
for i=1:m
row_number1=find(abs(U(i:m,i))==max(max(abs(U(i:m,i)))));
%i每增加一次找出i到m行元素绝对值中的最大值的位置
row_number2=min(row_number1)+i-1; %如果两行首位元素的值一样，取row_number1中位置较小的元素的位置
%此处加（i-1）因为abs((A(i:m,i))取的是i到m行绝对值最大元素的位置，该位置剔除前i行
temp1=U(row_number2,:); %该列中绝对值最大元素所在的行与该列第一行（i行）互换
U(row_number2,:)=U(i,:);
U(i,:)=temp1;
temp2=P(row_number2,:); %单位矩阵做同样的行变换
P(row_number2,:)=P(i,:);
P(i,:)=temp2;
temp3=L1(row_number2,:); %构成L的下三角矩阵非对角线元素做同样的行变换
L1(row_number2,:)=L1(i,:);
L1(i,:)=temp3;
L=eye(m)+L1;
L1(i+1:m,i)=U(i+1:m,i)/U(i,i);
% U(i+1:m,i:n)=U(i+1:m,i:n)-U(i+1:m,i)*U(i,i:n)/U(i,i);
U(i+1:m,i:n)=U(i+1:m,i:n)-L1(i+1:m,i)*U(i,i:n);
%求出L下三角矩阵中非零元素的值，同时将i+1行第一个元素化为0
end
%将所得P、L、U与Matlab自身所带的公式[L_lu,U_lu,P_lu]=lu(A)进行比较
[L_lu,U_lu,P_lu]=lu(A);
if sum(sum(P-P_lu<=EPSINON))==m^2&&sum(sum(L-L_lu<=EPSINON))==m^2&&sum(sum(U-U_lu<=EPSINON))==m^2
fprintf('与Matlab自带公式[L,U,P]=lu(A)所得结果一致；difference_value表示本程序与Matlab公式\n之间的差值。当差值小于10^(-15)时，矩阵中的元素赋值为0.\n');
else fprintf('该程序错误。\n');
%用difference_value来表示P*A与L*U的完全相等的程度
end
difference_value=P*A-L*U;
U(abs(U)
end
end
end