%% Anything on a row following a % is ignored by Matlab, 
%% so that's how one can add `comments' on the program
%%
%% Program GElim.m
%%
%% a program to do GAUSSIAN ELIMINATION of a MATRIX A, 
%% as shown and explained in Math 320 lecture on Fri Sept 27, 2002
%%
%% Fabian Waleffe, 9/27/2002
%%
%%

N=5;            % This is the number of equations, you can change it
A=randn(N,N+1); % This defines a random matrix called A
                % with N rows and N+1 columns
                % This AUGMENTED MATRIX defines a linear system of 
                % N equations in N unknowns
                % see Edwards and Penney Chap 3.2.

%% we deduced and experimented with the Gaussian elimination algorithm in B107.
%% Here's what we ended with

for j=1:N, for i=j+1:N, L=A(i,j)/A(j,j); A(i,:)=A(i,:)-L*A(j,:),pause,end,end

%% That one-liner is nice for cut-and-paste and recalls 
%% but the algorithm comes out better written as follows:

for j=1:N,
   for i=j+1:N,
       L=A(i,j)/A(j,j);
       A(i,:)=A(i,:)-L*A(j,:);
   end
end

%%------------------------------------------------------------------
%% If you play with this enough, you'll see that the numbers
%% can get pretty big and this algorithm can break down. In fact
%% it basically always breaks down if N is big.
%% The ultimate breakdown is when any A(j,j) becomes =0, obviously.
%% To help with this undesirable thing we can add some ROW SWAPPING
%% (also called PIVOTING):
%% Before we compute L, we first reorder all the equations
%% in such a way that abs(A(j,j)) is the biggest among the A(j:N,j)
%% i.e. biggest among all A(k,j) elements with k=j, j+1, ..., N-1, N
%%

%% Here's the beast with pivoting:

N=5;            
A=randn(N,N+1);  % Let's redefine a matrix to play with

for j=1:N,
   [am,ind]=max(abs(A(j:N,j))); jm=j-1+ind;   %% jm is the index of the max
   A([j jm],:)=A([jm j],:);                   %% swap Row j and Row jm
   for i=j+1:N,
       L=A(i,j)/A(j,j);                       %% Multiplier
       A(i,:)=A(i,:)-L*A(j,:), pause,         %% Row i becomes Row i - L*Row j
   end,
end

%% TRY IT!!
