% ge_all.m ___________ Gaussian Elimination with illustrations
% @ K Chen (Nov 2002)                         (see also g_e.m)
%__________________________________________________________________________ 
help ge_all
if ~exist('A'), disp('New matrix A set for you');
  n =  7; A=rand(n); %%A=3*rand(n).*(rand(n)>.6); % sparse
          A=A+9*eye(n); 
end,       n=size(A,1)
if ~exist('b'), disp('New RHS b set for you'); x=(1:n)'; b=A*x; clear x;end
A0=A;

 m=eye(n);  %%  See mfile=g_e.m  %% Matrix of all multipliers or m=spone(n);
   for k=1:n                           %% Step k
     eval(['L' num2str(k) '=eye(n);'])
       for i=k+1:n                       %% Row i 
         m(i,k) = A(i,k)/A(k,k);           %% Multiplier between rows i and k
           b(i) = b(i)  - m(i,k)*b(k);       %% Update b
           eval(['L' num2str(k) '(i,k)=m(i,k);'])
           for j=k:n                          %% Column j
             A(i,j) = A(i,j) - m(i,k)*A(k,j);  %% Update A
           end                               %% j
       end                                 %% i
 end                                     %% k    ==>   Check solution x=b
   L=m;  U=A;                          %% The A=LU factorisation completed
   L=m;  D=diag(diag(A));  M=D\A;    %% The A=LDM factorisation completed

   A = A0; %% Set to the original A

fprintf('\nThe GE-factorisation of A is:  A=LU or A=LDM \n   ')
fprintf('\nwhere one may check this overall decomposition:\n\t')

for k=n-1:-1:1
 fprintf('inv(L%d)*',k)
end
 fprintf('A = U\n',n)
fprintf('  and  this relationship\n\t')
   fprintf('L = ')
  for k=1:n-2
   fprintf('L%d*',k)
  end
   fprintf('L%d\n',n-1)
 
disp('Other variables:  A L U D M '),    clear i j k  A0

