%==========================================================%
%Full Pivoting Gauss Elimination Method====================%
%Author:Yao Fei  ==========================================%
%Rev:2.0   Date:2006.11.19=================================%
%==========================================================%

% x是方程的解
% amax中存放每一步的全主元
% P*A*Q=LU

function [amax1,P,Q,L,U,x]=geap(A,n,b)

B=A;
P=1:n;                              %行矢量，用于记录行交换的情况
Q=1:n;                              %列矢量，用于记录列交换的情况
U=zeros(n);
L=zeros(n);

%进行LU分解
amax1(1)=max(max(abs(B))); % 求出矩阵A中绝对值最大元素,作为第一个主元素
for k=1:n-1
    [s1,s2]=find(abs(B)==max(max(abs(B)))); %s1,s2分别为矩阵中绝对值最大元素的行和列下标。
                                            %因对B所有元素查找，故为全主元
                                            
    B([1,s1],:)=B([s1 1],:);       %进行行交换，将该最大值元素所在的行交换到第一行
    m=k+s1-1;   
    P([k m])=P([m k]);              %通过将行矢量P的第k个元素与第m个元素交换位置，记录行交换的情况
    U([m k],:)=U([k m],:);          %相应地进行行交换，调整U矩阵
    L([m k],:)=L([k m],:);          %相应地进行行交换，调整L矩阵
   
    B(:,[1,s2])=B(:,[s2 1]);       %进行列交换，将该最大值元素所在的列交换到第一列
    m=k+s2-1;                       %至此当前主元已经成为第一个元素
    Q([k m])=Q([m k]);              %通过将列矢量Q的第k个元素与第m个元素交换位置，记录列交换的情况
    U(:,[m k])=U(:,[k m]);          %相应地进行列交换，调整U矩阵
    L(:,[m k])=L(:,[k m]);          %相应地进行列交换，调整L矩阵
   
    piv=B(2:n-k+1,1)/B(1,1);        %piv是包含n-k个元素的归一化主元列矢量
    u=B(1,2:n-k+1);                 %u是n-k个元素的行矢量，等于B的第一行，从第二个元素开始，用于消元
    U(k,k:n)=B(1,:);                %用B的第一行替换U的第k行
    L(k+1:n,k)=piv;                 %用归一化主元列矢量替换L的第k列
    B(1,:)=[];                     
    B(:,1)=[];                      %得到原来B的第一个元素（b11）的余子式，现在B被降了一阶，为n-k阶
    B=B-piv*u;                      %消元过程
   
    amax1(k+1)=max(max(abs(B)));     %确定下一步的全主元
end
U(n,n)=B;                           %现在把最后一个元素补上
L=L+eye(n);                         %至此已经完成A的LU分解！！

%前向消去
y = zeros(n,1);
for k = 1:n
   j = 1:k-1;                       %当k为1时，j为空矩阵（1 x 0阶）
   y(k) = b(P(k)) - L(k,j)*y(j);
end

%后向替换，得到解x(次序未调整)
xtemp = zeros(n,1);
x= zeros(n,1);
for k = n:-1:1
   j = k+1:n;
   xtemp(k) = (y(k) - U(k,j)*xtemp(j))/U(k,k);
end

%根据列矢量Q，调整解的次序，得到最终的解x
for  k=1:n
    x(Q(k))=xtemp(k);
end