function x = Gauss_demo(A,c)
%GAUSS_DEMO: Demonstrates the steps in the Gauss elimination method.
% Checks consistency and existence of a set of simultaneous linear 
% algebraic equations using rank(A) and rank(Augmented) and performs
% the Gauss elimination method. Row pivoting is carried out.
% The program writes out the step-by-step procedure.
%
%   GAUSS_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 JORDAN_DEMO, GAUSS, JORDAN, SEIDEL, JACOBI
%
%
%(c) A. Constantinides
%October 1, 2002

disp('**********************************************************************')
disp('*                                                                    *')
disp('*                THE GAUSS ELIMINATION METHOD                        *')
disp('*                                                                    *')
disp('*         FOR SIMULTANEOUS LINEAR ALGEBRAIC EQUATIONS                *')
disp('*                                                                    *')
disp('*                      (Gauss_demo.m)                                *')
disp('*                                                                    *')
disp('**********************************************************************')

c = (c(:).')';    % Make sure it is a column vector
tol=1e-6;
n = length(c);
[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_of_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 elimination procedure
disp(' ')
disp('START THE GAUSS ELIMINATION 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
   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
   lk = unit;
   for m = k + 1 : n
      fprintf('\nDIVIDE ROW %2i  BY  %5g ', k,  aug(k, k))
      fprintf('\nMULTIPLY ROW %2i  BY  %5g  AND SUBTRACT FROM ROW %2i \n',  k, aug(m, k), m)    
      lk(m , k) = - aug(m , k) / aug(k , k);
      lk*aug
      pause
   end
   aug = lk * aug;
 end 
 disp(' ')
 disp('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX')

end

x = zeros(n , 1);

% Back substitution
 
% Apply the back-substitution formulas
if rank_of_A == length(A) 
    X(n) = aug(n , n + 1) / aug(n , n);
    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(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
  U=[];
    for i=1:n
    U=[U aug(:,i)];
    end
disp('Diagonal matrix U and its determinant')
U
Determinant=det(U)
  disp(' ')
  disp('XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX')
end