%==========================================================%
%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