%Gauss-Jordan Elimination
%Description: This procedure is used to calculate the matrix inverse of a
%given square matrix using Gauss-Jordan Elimination algorithm.
%This algorithm first augments the square matrix with an identity matrix of
%the same dimension, then reduces the square matrix to reduced row
%echelon form and obtain the inverse of the given matrix by undoing the
%augmented matrix.
%Input: a non-singular square matrix A
%Output: matrix inverse of A and also the debug information if
%necessary.
%Note: This procedure is numerically unstable. Besides, zero pivot
%element of the matrix could cause this procedure to terminate without
%returning the correct solution.
%This program is written by Xiaoke Yang @ School of Automation Science and
%Electrical Engineering, Beihang University.You can copy, modify or 
%redistribute it freely.Welcome reports of bugs, suggestions, etc. 
%yxkmlstar@gmail.com
%Last modified by Xiaoke Yang, Oct.10,2007

% A=[2 -1 0;-1 2 -1;0 -1 2]
A=rand(10)*10;
N=length(A(:,1));
b=eye(N);
Au=[A b];                   %augmented matrix
%forward elimination to obtain the row echelon form
for i=1:N-1
    if(abs(A(i,i))<1e-15)
        disp(sprintf('A(%d,%d)=0, program terminated!\n',i,i))
        return;
    end
    m=-1/Au(i,i);
    for j=i+1:N             %j indicating the row
        for k=i+1:2*N         %k indicating the column
            Au(j,k)=Au(j,k)+Au(j,i)*m*Au(i,k);
        end
        Au(j,i)=0;
    end
    disp(sprintf('after %d elimination, \n[A b]=',i))
    disp(Au)
end
%backward elimination to obtain the reduced row echelon form
for i=N:-1:1                    %i indicating the column of the original matrix
    for j=2*N:-1:i              %j indicating the column of the augmented matrix
        Au(i,j)=Au(i,j)/Au(i,i);
        for k=i-1:-1:1          %k indicating the row
            Au(k,j)=Au(k,j)-Au(i,j)*Au(k,i);
        end
    end
    disp(sprintf('after %d backward elimination, \n[A b]=',i))
    disp(Au)
end
invA=Au(:,N+1:2*N);
disp(sprintf('The inverse matrix is \n invA='));
disp(invA);