% g_e.m ___________ Gaussian Elimination (see also ge_all.m)
% @ K Chen (Nov 2002)
%__________________________________________________________________________ 
help g_e
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), A0=A;
if ~exist('b'), disp('New RHS b set for you'); x=(1:n)'; b=A*x; clear x;end

 m=eye(n);  %%  See mfile=g_e.m  %% Matrix of all multipliers or m=spone(n);
   for k=1:n                           %% Step k
       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
           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
   x=U\b;
A=A0; clear A0
who A x L U D M
