function gauss_elim_pivot        %                   F2012      J. M. Hill
% This program reads an nxn augmented coefficient matrix from a file or the
% keyboard, stores keyboard entry to a file, employs Gaussian Elimination 
% with partial pivoting to solve the set of equations.
clc;                             %clear command window
clear;                           %clear variables and functions from memory
answer=input('Type k for kbd input or f to read input from a file: ','s');
if(lower(answer) == 'k')
    [n,a,b] = get_input_kbd;     %get input from keyboard
    write_matrix_file(n,a,b);    %write input data to file
else
    [n,a,b] = read_matrix_file;  %read matrix from file 
end
disp(' ');                       %put a blank line in output
write_matrix_screen(n,a,b);      %write input data to screen
disp(' ');
[a,b] = gauss(n,a,b);       %produce an upper-triangular aug. coef. matrix      
write_matrix_screen(n,a,b);
disp(' ');
x = solve(n,a,b);           % solve upper triangular matrix
write_vector(n,x);          %write answers
%--------------------------------------------------------------------------     
function write_matrix_screen(n,a,b)                   
for i = 1:n
    for j = 1:n
        fprintf(1,'%3.2g ',a(i,j));  %write matrix to screen
    end
    fprintf(1,'%3.2g',b(i));         %write b(i) to screen
    fprintf(1,'\n');                 %new line on screen
end
%--------------------------------------------------------------------------     
function write_matrix_file(n,a,b)                       
fid = fopen('data.txt','w');         %open file for writing
for i = 1:n
    for j = 1:n
        fprintf(fid,'%3.2g ',a(i,j));%write matrix to file
    end
    fprintf(fid,'%3.2g',b(i));       %write b(i) to file
    fprintf(fid,'\n');               %new line in file
end
fclose(fid);
%--------------------------------------------------------------------------
function [n,a,b] = get_input_kbd                          
n =  input('input the number of equations: ');
for i = 1:n                          
    for j = 1:n
        fprintf('a(%d,%d)= ',i,j);  %prompt for each element of a on row i
        a(i,j) = input(' ');        %read each element
    end
    fprintf('b(%d)= ',i);           %prompt for b(i)
    b(i) = input(' ');              %read b
end
%--------------------------------------------------------------------------
function [n,a,b] = read_matrix_file                       
c = load('data.txt');   %load entire matrix
[n,m] = size(c);        %determine size, n rows, m columns
b = c(:,m);             %strip b vector out of c matrix  
for i = 1:n
    for j = 1:n
        a(i,j) = c(i,j);%strip a matrix out of c matrix
    end
end                     %alternatively: a = c(:,1:n) can replace the loop; 
%--------------------------------------------------------------------------
function write_vector(n,x)  % write vetcor to screen
for i = 1:n
    fprintf(1,'x(%d) = %g \n',i,x(i));
end
%--------------------------------------------------------------------------
function [a,b] = elim(n,a,b,i) % zero elements under current pivot element
for k = i+1:n
    q = -a(k,i)/a(i,i);
    a(k,i) = 0;
    b(k) = q*b(i)+b(k);
    for j = i+1:n
        a(k,j) = q*a(i,j) + a(k,j);
    end
end
%--------------------------------------------------------------------------
function x = solve(n,a,b)
x(n) = b(n)/a(n,n);
for k = 1:n-1
    q = 0;
    for j = 1:k
        q = q + a(n-k,n+1-j)*x(n+1-j);
        x(n-k) = (b(n-k) - q)/a(n-k,n-k);
    end
end
%--------------------------------------------------------------------------
function [a,b] = gauss(n,a,b)
for i = 1:n-1
    big = abs(a(i,i));                %
    ibig = i;                         %
    for k = i+1:n                     %
        aba = abs(a(k,i));            %
        if(aba > big)                 %% identify largest pivot element
            big = aba;                %%
            ibig = k;                 %
        end                           %
    end                               %
    if(abs(big) < 1.e-12)                %
        disp('matrix is singular');      %% test for sigular matrix
        return;                          %
    end                                  %
    if(ibig ~= i)                     %
        for k = 1:n                   %
            temp = a(i,k);            %
            a(i,k) = a(ibig,k);       %
            a(ibig,k) = temp;         %% swap rows to put largest coeff
        end                           %% on main diagonal
        temp = b(i);                  %
        b(i) = b(ibig);               %
        b(ibig) = temp;               %
    end                               %                      
    [a,b] = elim(n,a,b,i);        % zero elements under current pivot elem.
    write_matrix_screen(n,a,b);   %write intermediate steps to screen
    disp(' ');
end