function x= Jordan_demo(A,c)
%JORDAN_DEMO: Demonstrates the steps in the Gauss-Jordan reduction method.
% Checks consistency and existence of a set of simultaneous linear 
% algebraic equations using rank(A) and rank(Augmented) and performs
% the Gauss-Jordan reduction method. No pivoting is carried out.
% The program writes out the step-by-step procedure.
%
%   JORDAN_DEMO(A,c) solves the set of linear algebraic equations Ax=c.
%
%   A = the matrix of coefficients
%   c = the vector of constants
%   x = vector of solution 
%
%   See also GAUSS_DEMO, GAUSS, JORDAN, SEIDEL, JACOBI
%
%(c) A. Constantinides
%October 1, 2002

clc
disp('**********************************************************************')
disp('*                                                                    *')
disp('*              THE GAUSS-JORDAN REDUCTION METHOD                     *')
disp('*                                                                    *')
disp('*       FOR SIMULTANEOUS LINEAR ALGEBRAIC EQUATIONS                  *')
disp('*                                                                    *')
disp('*                       (Jordan_demo.m)                              *')
disp('*                                                                    *')
disp('**********************************************************************')

c = (c(:).')';  % Make sure it is a column vector
tol=1e-6;
n= length(A);
[nr nc] = size(A);

% Check coefficient matrix and vector of constants
if nr ~= nc
   error('Coefficient matrix is not square.')
end
if nr ~= n
   error('Coefficient matrix and vector of constants do not have the same length.')
end

% Check consistency and existence
disp('Augmented matrix [A|c]:')
aug = [A c]						% Augmented matrix
det_of_A=det(A);
rank_of_A=rank(A);
rank_of_aug=rank(aug);
fprintf('\nDeterminant of matrix A = %g ', det_of_A)
fprintf('\nRank of matrix A = %2i ', rank_of_A)
fprintf('\nRank of augmented matrix A = %2i ', rank_of_aug)
if rank(A)==rank_of_aug
   fprintf('\n\nThe equations are consistent because rank(A) = rank(A augmented).')   
end
if rank_of_A ==length(A)
   fprintf('\nThere will be a unique solution because rank(A) = n.') 
   if c==0
   fprintf('\nBut the solution is trivial because the system is homogeneous.') 
   end
end
   
if rank_of_A < length(A)
    fprintf('\n\nMatrix A is singular of rank = %2i',rank_of_A)
    if rank_of_A == rank_of_aug
        fprintf('\nThere will be an infinite number of solutions because rank(A) < n.')
    end
    if rank_of_A < rank_of_aug
            fprintf('\nThere will be no solution because rank(A) < rank(A augmented).')
    end
end

unit = eye(n);          	    % Unit matrix
aug = [A c];					% Augmented matrix
disp(' ');disp(' ')
disp('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX')

% Beginning of the Gauss-Jordan procedure
disp(' ')
disp('START THE GAUSS-JORDAN REDUCTION METHOD:')
disp(' ')
disp('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX')

for k = 1 : n
   pivot = abs(aug(k , k));
   prow = k;
   pcol = k;
   
   % Locating the maximum pivot element
   for row = k : n
      if abs(aug(row , pcol)) > pivot
         pivot = abs(aug(row , pcol));
         prow = row;
      end
   end
   
   % Interchanging the rows
   disp(' ')
   pr = unit;
   tmp = pr(k , :);
   pr(k , : ) = pr(prow , : );
   pr(prow , : ) = tmp;
   aug = pr * aug;
   if prow ~= k
        disp('PERFORM PIVOTING OF ROWS:')
        fprintf('INTERCHANGE ROW %2i  WITH ROW %2i\n',k,prow)
        aug
        pause
        disp(' ')
   end
   
   % Reducing the elements below diagonal to zero in the column k
 if abs(aug(k,k)) > tol 
     fprintf('NORMALIZATION: DIVIDE ROW %2i  BY  %5g \n', k,  aug(k,k))
       for J = n+1:-1:k ; 
          aug(k, J) = aug(k, J)/aug(k, k);
       end
       aug
       pause
     for I = 1: n
       if I ~= k
         fprintf('\nREDUCTION: MULTIPLY ROW %2i  BY  %5g  AND SUBTRACT FROM ROW %2i ', k, aug(I, k), I)
         for J = n+1:-1:k ; 
           aug(I, J) = aug(I, J)  -  aug(I, k)* aug(k, J);
         end
         disp(' ')
         aug
         pause
       end
     end
 end
    disp(' ')
    disp('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX')
end

x = zeros(n , 1);
  
%end

% Apply the back-substitution formulas
if rank_of_A == length(A) 
    X(n) = aug(n, n + 1);
    NCOUNT = n - 1;
end 
if rank_of_A < length(A) & rank_of_A ==rank_of_aug
   NMR = n-rank_of_A;
   fprintf('\nTHE PROGRAM SETS %2i UNKNOWN(S) TO UNITY AND ', NMR)
   fprintf('\nREDUCES THE PROBLEM TO FINDING THE OTHER %2i UNKNOWN(S) ', rank_of_A)
   for JJ = 1:NMR
       X(n + 1 - JJ) = 1;
   end
   NCOUNT = rank_of_A;
end  
if rank_of_A < length(A) & rank_of_A < rank_of_aug
    fprintf('\n\nThere are no solutions to this problem\n\n')
else
  for I = NCOUNT:-1:1
      SUM = 0;
      for J = I + 1:n
      SUM = SUM + aug(I, J) * X(J);
      end 
  X(I) = (aug(I, n + 1) - SUM) / aug(I, I);
  end
  disp(' ');
  disp('The vector of results is:')
  x=X'
  % Verification
  disp(' ')
  disp('Displayed below is A*x which should correspond to the vector c:')
  Ax=A*x
  disp(' ')
  disp('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX')
end