% Script File: ShowGe1
% Illustrates reduction to reduced row echelon form for
% arbitrary system (and arbitrarily many RHS without pivoting.
clear
clc
format short
disp('Steps through the reduction of an augmented matrix')
disp('to upper echelon form')
input('Strike Any Key to Continue.');
clc
% Cases for problem 1.2.5; follow the model to enter other problems
% Uncoment the relevant lines by removing the '%', 
% or create a new augmented matrix A = [A | b]
% (5a)
% A = [1 -2;2 -1]; b = [3 ; 9];
% (5b)
% A = [2 -3;-4 6]; b=[5;8];
% (5c)
% A = [1 1;2 3;3 -2]; b = [0;0;0];
% (5d)
% A = [3 2 -1;1 -2 2;11 2 3]; b = [4;1;14];
% (5e)
% A = [2 3 1;1 1 1;3 4 2]; b = [1;3;4];
% (5f) 
  A = [1 -1 2;2 3 -1;7 3 4]; b = [4;1;7];
% (5g) 
% A = [1 1 1 1;2 3 -1 -1;3 2 1 1;3 6 -1 -1]; b = [0;2;5;4];
% (5h)
% A = [1 -2;2 1;-5 8]; b = [3;1;4];
% (5i) 
% A = [-1 2 -1;-2 2 1;3 2 2;-3 8 5]; b = [2;4;5;17];
% (5j) 
% A = [1 2 -3 1;-1 -1 4 -1;-2 -4 7 -1]; b = [1;6;1];
% (5k) 
% A = [1 3 1 1;2 -2 1 2;1 -5 0 1]; b = [3;8;5];
% (5l)
% A = [1 -3 1;2 1 -1;1 4 -2;5 -8 2]; b = [1;2;1;5];
[m,n] = size(A);
[m1,r] = size(b);
n1 = n+r;
if m1 ~= m 
    disp('error; row dimensions of A, b must agree')
else
    A = [A b]
    v = zeros(m,1);
    for k=1:m-1
        if A(k,k) ~= 0
            v(k+1:m) = A(k+1:m,k)/A(k,k);
            for i=k+1:m
                clc
                Before = A     
disp(sprintf('Multiply row %1.0f by %7.4f / %7.4f ',k,A(i,k),A(k,k)))
disp(sprintf('and subtract from row %1.0f:',i))
                A(i,k:n1) = A(i,k:n1) - v(i)*A(k,k:n1);
                After = A
                input('Strike Any Key to Continue.');
            end
            disp(' ')
        else
            kk = 1;
            while A(k+kk,k) == 0 & k+kk < m
                kk = kk+1;
            end
            if kk <=m & A(kk,k) ~= 0
disp('zero pivot encountered but rest of column contains nonzero element(s)')
disp('simple method fails')
            end
        end
    end
    clc
disp('Forward elimination is complete: upper triangular form')
    U = triu(A)
    pause
    flag1 = 0; flag = 0;
% disp('Back substitution:')
    for i = 1:m
        if flag1 == 0
            j = i;
            while A(i,j) == 0 & j < n1
                j = j+1;
            end
            if j > n 
                if j<=n1 & A(i,j) ~= 0
disp('System is inconsistent: zero row but nonzero RHS')
                    flag = 1;
                else
disp(sprintf('rows from %2i-th on are zero; system is consistent',i))
                   ifinal = i-1; 
                   flag1 = 1;
                end
            else
                p(i) = j; %keep track of pivot location for each row
                for k = j+1:n1
                    A(i,k) = A(i,k)/A(i,j);
                end
                A(i,j) = 1;
                ifinal = i;
            end
        end
    end

    if flag == 0
        input('Strike Any Key to Continue.');
disp('Scale each row by its pivot (first nonzero element):')
        Scaled = A
        ifinal
disp('pivots (and basic variables) found on columns:')
        pivots = p
disp('remaining columns of original matrix corresp. to free variables')
        input('Strike Any Key to Continue.');
        for k=ifinal:-1:2
            v(1:k-1) = A(1:k-1,p(k));
            j = p(k);
            for i = 1:k-1
                clc
disp('Carry out backward elimination:')
                Before = A 
disp(sprintf('Multiply row %1.0f by %7.4f / %7.4f ',k,A(i,k),A(k,k)))
disp(sprintf('and subtract from row %1.0f:',i))
                A(i,j:n1) = A(i,j:n1) - v(i)*A(k,j:n1);
                After = A
                input('Strike Any Key to Continue.')
            end
disp('Final Reduced row echelon form; ');
disp('variables corresponding to nonpivotal columns ("free variables")')
disp('can be set to arbitrary values, and variables corresponding to')
disp('pivotal columns (basic variables) can be solved for') 
        end
    end
end
