% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   GAUSSELIM solves a system of n linear
%       equations in n unknowns
%
%       Math 225B, Winter 2011, HW1, #4
%       Prof. Q. Nie
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Solves ay = b
%           where
%               a is an nxn matrix
%               b is an nx1 column vector
%       Returns the solutions in
%       	y, an nx1 column vector
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    function y = GaussElim(a, b)
%
%   Check dimensions for consistency
%
    dm = size(a); dn = size(b);
    if (dm(1) ~= dm(2)) || (dn(1) ~= dm(1)) || dn(2) ~= 1
        error('Inconsistent dimensions');
    end
%
%   Set dimensions and elimination matrix
%
    n = dm(1);              % Number of equations and unknowns
    x = [a b];              % Concatenate b to a, x is n x (n+1)
    w = 1:n;                % Indices for column pivoting
    clear dm dn;            % Clear temporary variables
%
%   Loop for iterations with (n-k)xn submatrix
%
    for k = 1:n-1
        [c dm] = max(abs(x(k:n,k:n)));   % Find max. entry of submatrix
        [d dn] = max(abs(c));
        if d < eps
            error('Matrix is singular within machine accuracy');
        end
%
%       Perform full pivoting
%
        k1 = dm(dn) + k - 1;        % Get indices of max. entry
        k2 = dn + k - 1;
        z = x(k,:);                 % Swap rows
        x(k,:) = x(k1,:);
        x(k1,:) = z;
        u = x(:,k);                 % Swap columns
        x(:,k) = x(:,k2);
        x(:,k2) = u;
        t = w(k);                   % Swap column indices
        w(k) = w(k2);
        w(k2) = t;
%
%       Perform Gauss elimination
%
        t = x(k,k);                 % Get pivot element
        x(k,:) = x(k,:) / t;        % Divide row by pivot element
        for m = k+1:n
            t = x(m,k);             % Get first element this row
            z = x(m,:);             % Get row
            x(m,:) = z - t * x(k,:);% Subtract row to zero m,k element
        end
    end
    t = x(n,n);                     % Normalize final row
    x(n,:) = x(n,:) / t; 
%
%   Solve equations
%
    u = zeros(n,1);                 % Initialize solution vector
    for k = n:-1:1
        v = u;                      % Move backward through equations
        u(k) = x(k,n+1) - x(k,1:n) * v;
    end
%
%   Reverse pivoting of variables
%
    y = zeros(1,n);                 % Initialize output vector
    for k = 1:n
        y(w(k)) = u(k);             % Get kth variable solution
    end
%
%   Clear temporary variables
%
    clear t c d u w x z;
%
    end
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
