%g_j.m ___________ Gaussian-Jordan Elimination (see also gj_all.m)
%@ K Chen (Dec 2002)
%__________________________________________________________________________ 
help g_j
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);  %%  Mfile=g_j.m          %% Matrix of all multipliers or m=spone(n);
   for k=1:n                           %% Step k
           eval(['E' num2str(k) '=eye(n);'])
           eval(['L' num2str(k) '=eye(n);'])
           eval(['U' num2str(k) '=eye(n);'])
       for i=[1:k-1 k+1:n]              %% Row i 
         m(i,k) = A(i,k)/A(k,k);        %% Multiplier between rows i and k
           eval(['E' num2str(k) '(i,k)=-m(i,k);'])
          if i<k
           eval(['U' num2str(k) '(i,k)=-m(i,k);'])
          else
           eval(['L' num2str(k) '(i,k)=-m(i,k);'])
          end
           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
         m(k,k) = 1/A(k,k);              %% unit diagonals
           for j=k:n                    %% Column j
             A(k,j) = m(k,k)*A(k,j);  %% Update row k (scaling)
           end                      %% j
         b(k) = m(k,k)*b(k);  %% Update row k 
    eval(['E' num2str(k) '(k,k)=m(k,k);'])
    eval(['U' num2str(k) '(k,k)=m(k,k);'])
 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; clear A0
who A x L* U* E* D M
fprintf('\nThe GJ-factorisation of A is:  \n   ')
fprintf('\nwhere one may check this overall decomposition:\n\t')

for k=n:-1:1
 fprintf('E%d*',k)
end
 fprintf('A = I\n',n)
 
disp('Other variables:  A L U D M '),    clear i j k  A0

