clc,clear;
A=[10 -7 0 1;
    -3 2.099999 6 2;
    5 -1 5 -1;
    2 1 0 2];
b=[8 5.900001 5 1]';
vpa(A*[0 -1 1 1]'-b,100)
Realx=[0 -1 1 1]'
vpa(Realx,9)
[L,U,x,detA]=GaussElim(A,b)
vpa(x,9)
err1=vpa(max(abs(x'-Realx)),100)

%

[P,L,U,x,detA]=GaussElim2(A,b)
vpa(x,9)
err2=vpa(max(abs(x'-Realx)),100)
err2-err1

%

%乐于助人排除bug
clear,clc;
A=[10 -7 0 1;
    -3 2.099999 6 2;
    5 -1 5 -1;
    2 1 0 2];
b=[8 5.900001 5 1]';
A0=A; %用A0存储原始矩阵
n=size(A,1); %n为A的维度，因为A是方阵，所以指定了列数，这样n只会返回行数值
P=eye(n); %定义原始排列阵P
%下求P
if(abs(A(1,1))~=max(abs(A(:,1))))
    w=find(abs(A(:,1))==max(abs(A(:,1)))); %定位第一列中最大元素的所在行
    v=A(w,:);
    A(w,:)=A(1,:);
    A(1,:)=v; %将最大值所在的行与第一行交换
    v=P(w,:);
    P(w,:)=P(1,:);
    P(1,:)=v; %交换信息存储在P中
end
for i=2:1:n
    A(i,:)=A(i,:)-(A(i,1)/A(1,1)).*A(1,:); %高斯消元
end
for i=2:1:n-1
    if (abs(A(i,i))~=max(abs(A(i:n,i))))
        w=find(abs(A(i:n,i))==max(abs(A(i:n,i)))); %定位第i列中最大元素的所在行
        w=w+i-1; %将w转换为在A中的行数
        v=A(w,:);
        A(w,:)=A(i,:);
        A(i,:)=v; %将最大值所在的行与第i行交换
        v=P(w,:);
        P(w,:)=P(i,:);
        P(i,:)=v; %交换信息存储在P中
    end
    for k=i+1:1:n
        A(k,:)=A(k,:)-(A(k,i)/A(i,i)).*A(i,:); %高斯消元
    end
end


A=P*A0;b=P*b;




n=size(A,1); %n为A的维度，因为A是方阵，所以指定了列数，这样n只会返回行数值
L=zeros(n);
U=zeros(n);
det=A(1,1); %定义行列式计算初始值
for i=1:1:n
    U(1,i)=A(1,i); %赋予U第一行的值
end
for i=1:1:n
    L(i,1)=A(i,1)/U(1,1); %赋予L第一列的值
end
for i=2:1:n
    for j=i:1:n
        R=0;
        S=0;
        for k=1:1:i-1
            R=R+L(i,k)*U(k,j);
            S=S+L(j,k)*U(k,i); %计算求和数
        end
        U(i,j)=A(i,j)-R;
        L(j,i)=(A(j,i)-S)/U(i,i); %按行按列计算U与L
    end
    det=det*U(i,i); %每算一次行便记录对角元，最终返回矩阵的值
end
detA=det;
y(1)=b(1);
for i=2:1:n
    t=0;
    for j=1:1:i-1
        t=t+L(i,j)*y(j);
    end
    y(i)=b(i)-t;
end %求解Ly=b
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
    t=0;
    for j=i+1:1:n
        t=t+U(i,j)*x(j);
    end
    x(i)=(y(i)-t)/U(i,i);
end %求解Ux=y

%
%

function [L,U,x,detA]=GaussElim(A,b)
n=size(A,1);
L=zeros(n);
U=zeros(n);
det=A(1,1); 
for i=1:1:n
    U(1,i)=A(1,i); %赋予U第一行的值
end
for i=1:1:n
    L(i,1)=A(i,1)/U(1,1); %赋予L第一列的值
end
for i=2:1:n
    for j=i:1:n
        R=0;
        S=0;
        for k=1:1:i-1
            R=R+L(i,k)*U(k,j);
            S=S+L(j,k)*U(k,i); %计算求和数
        end
        U(i,j)=A(i,j)-R;
        L(j,i)=(A(j,i)-S)/U(i,i); %按行按列计算U与L
    end
    det=det*U(i,i); %每算一次行便记录对角元，最终返回矩阵的值
end
detA=det;
y(1)=b(1);
for i=2:1:n
    t=0;
    for j=1:1:i-1
        t=t+L(i,j)*y(j);
    end
    y(i)=b(i)-t;
end %求解Ly=b
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
    t=0;
    for j=i+1:1:n
        t=t+U(i,j)*x(j);
    end
    x(i)=(y(i)-t)/U(i,i);
end %求解Ux=y
end


function [P,L,U,x,detA]=GaussElim2(A,b)
A0=A; %用A0存储原始矩阵
n=size(A,1); %n为A的维度，因为A是方阵，所以指定了列数，这样n只会返回行数值
P=eye(n); %定义原始排列阵P
%下求P
if(abs(A(1,1))~=max(abs(A(:,1))))
    w=find(abs(A(:,1))==max(abs(A(:,1)))); %定位第一列中最大元素的所在行
    v=A(w,:);
    A(w,:)=A(1,:);
    A(1,:)=v; %将最大值所在的行与第一行交换
    v=P(w,:);
    P(w,:)=P(1,:);
    P(1,:)=v; %交换信息存储在P中
end
for i=2:1:n
    A(i,:)=A(i,:)-(A(i,1)/A(1,1)).*A(1,:); %高斯消元
end
for i=2:1:n-1
    if (abs(A(i,i))~=max(abs(A(i:n,i))))
        w=find(abs(A(i:n,i))==max(abs(A(i:n,i)))); %定位第i列中最大元素的所在行
        w=w+i-1; %将w转换为在A中的行数
        v=A(w,:);
        A(w,:)=A(i,:);
        A(i,:)=v; %将最大值所在的行与第i行交换
        v=P(w,:);
        P(w,:)=P(i,:);
        P(i,:)=v; %交换信息存储在P中
    end
    for k=i+1:1:n
        A(k,:)=A(k,:)-(A(k,i)/A(i,i)).*A(i,:); %高斯消元
    end
end
%求出交换阵P后交给高斯消元法得到L,U,x,detA的值
[L,U,x,detA]=GaussElim(P*A0,P*b);
end