A=[1 2 3;1 3 5; 3 1 2];
[N,N]=size(A);
exc=zeros(N,2);
P=zeros(N,N);
A1=A;
B1=zeros(N,N);
INVB1=zeros(N,N);
B2=zeros(N,N);
INVB2=zeros(N,N);
I=eye(N);
T=zeros(N,N);           % 單位矩陣的分量
Ek=I;

if A(1,1)==0            % 交換使得第一行第一個元素不為0
    for  i=2:N
        if A(1,i)~=0
            r=i;
            break
        end
    end
    A1(:,1)=A(:,r);
    A1(:,r)=A(:,1);
end

for i=1:N
    P(:,i)=A1(:,i);
    T(:,i)=I(:,i);
end
B1=[P(:,1),T(:,2:N)];
t=A1(1,1);
INVB1=[P(:,1)/t,T(:,2:N)];

k=1;
for i=2:N-1
    Ps=INVB1*P(:,i);
    if Ps(k)==0
        break
    else
        Ps=-Ps/Ps(k);
        Ps(k)=-Ps(k);
        Ek(:,k)=Ps;
    end
    B2=[P(:,1:k+1),T(:,k+2:N)];
    INVB2=Ek*INVB1;