function Ainv= Jordan_inverse(A)
%JORDAN_INVERSE: Calculates the inverse of matrix A using the Gauss-Jordan 
% Reduction method and writes out the step-by-step procedure
%
%   JORDAN_INVERSE(A) Calculates the inverse of matrix A  
%
%   A = the matrix to be inverted
%   Ainv = the inverse of A
%
%
%(c) A. Constantinides
%October 1, 2002

clc
disp('**********************************************************************')
disp('*                                                                    *')
disp('*               THE GAUSS-JORDAN REDUCTION METHOD                    *')
disp('*                                                                    *')
disp('*                     FOR MATRIX INVERSION                           *')
disp('*                                                                    *')
disp('*                      (Jordan_inverse.m)                            *')
disp('*                                                                    *')
disp('**********************************************************************')

N= length(A);
Q=0;
disp('Augmented matrix [A|I]:')
Augm=[A,eye(N)]
fprintf('\nDeterminant of matrix A = %g ', det(A))
fprintf('\nRank of matrix A = %2i ', rank(A))
   
if rank(A) < length(A)
    fprintf('\n\nMatrix A is singular of rank = %2i',rank(A))
    fprintf('\nSingular matrices do not have inverse')
end
% Beginning of the Gauss-Jordan reduction procedure
for K = 1:N
    if Augm(K,K) ~=0
     disp(' '); disp(' ')
     fprintf('NORMALIZATION: DIVIDE ROW %2i  BY  %5g ', K,  Augm(K, K))
       for J = 2*N:-1:K ; 
          Augm(K, J) = Augm(K, J)/Augm(K, K);
       end
       for I = 1: N
       if I ~= K
         fprintf('\nREDUCTION: MULTIPLY ROW %2i  BY  %5g  AND SUBTRACT FROM ROW %2i ',  K, Augm(I, K), I)
       for J = 2*N:-1:K ; 
        Augm(I, J) = Augm(I, J)  -  Augm(I, K)* Augm(K, J);
       end
       disp(' ')
       Augm
       pause
       end
    end
  end
  
end


if rank(A) < length(A) 
    fprintf('\n\nThere is no inverse for this matrix\n\n')
else
  disp(' ');disp(' ')
  disp('The inverse of the matrix is:')
 
  Ainv=Augm(:,N+1:2*N)
end