%The Scaled Partial Pivoting Gaussian part
% (only update of 'A', i.e. the LU factorization)

function [A,l]  = Gauss(A)

% Given as input: 
%   a square coefficient matrix 'A'
% we will have as output:
%   the 'A=LU' factorization and the pivoting index vector 'l'

[m,n] = size(A);
if abs(m-n)>0
   error('The Matrix A is not square');
end
fprintf('This is the partial scaled Gaussian elimination part: \n\n')
% start with index vector 'l' being l=[1,2,...,n]
% and generate a scale vector 's'
for i=1:n
    l(i) = i;
    smax =0;
    for j=1:n
        smax = max(smax,abs(A(i,j)));
    end
    s(i)=smax;
end    
fprintf('The (global) scale vector is %1.f \n')
s
fprintf('The initial vector ell of indexes is %1.f \n')
l
%pause
% This is the scaled partial pivoting forward Gaussian elimination
for k=1:n-1
    rmax=0;
    for i=k:n
        r=abs(A(l(i),k)/s(l(i)));
        if (r>rmax)
            rmax = r;
            j=i;
        end
    end
    alpha = l(j);
    beta = l(k);
    l(j) = beta;
    l(k) = alpha;
%    l(j) = l(k);
fprintf('\n The iteration(column no.) in partial scaled Gauss elimination is %1.f\n',k)
fprintf('the vector ell of indexes in the partial scaled Gaussian  elimination is %1.f \n ')
l
    for i=k+1:n
        xmult=A(l(i),k)/A(l(k),k);
        xmultipliers(l(i),k)=xmult;
        %fprintf('the multipliers for column %1.f are %1.2f\n',[k xmult])
        A(l(i),k) = xmult;
        for j = k+1:n
            A(l(i),j) = A(l(i),j) - xmult*A(l(k),j);
        end
    end 
    fprintf('the multipliers for column %1.f \n',k)
    xmultipliers
% pause   
end 
