%g_h.m ___________ Gauss-Huard Elimination 
%@ K Chen (Dec 2002)
%__________________________________________________________________________ 
help g_h
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_h.m         %% Matrix of all multipliers or m=spone(n);
   for k=1:n                           %% Step k
           eval(['H' num2str(k) '=eye(n);'])
       for i=1:k-1                      %% Row i 
         m(i,k) = A(i,k)/A(k,k);        %% Multiplier between rows i and k
           eval(['H' num2str(k) '(i,k)=-m(i,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
         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(['H' num2str(k) '(k,k)=m(k,k);'])
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k+1
         i=k+1;                      %% Row i 
      if i<=n
         for k1=1:k    % k Rows to elimin k+1
           m(i,k1) = A(i,k1);        %% Multiplier between rows i and k
           eval(['H' num2str(k) '(i,k1)=-m(i,k1);'])
           b(i) = b(i)  - m(i,k1)*b(k1);       %% Update b
           for j=i:n                          %% Column j
             A(i,j) = A(i,j) - m(i,k1)*A(k1,j);  %% Update A
           end                               %% j
         end                               %% Next row k1
           for j=1:i-1                  %% Column j
             A(i,j) = 0;
           end                        %% j
 eval(['H' num2str(k) '(i,k)=H' num2str(k) '(i,1:k)*H'... % Special position
                                   num2str(k) '(1:k,k);'])  
      end     % if
 end                              %% k    ==>   Check solution x=b
         A=A0; clear A0
who A x L* U* E* D M
fprintf('\nOne may check this overall ')
fprintf(' GH-factorisation of A :  \n   ')

for k=n:-1:1
 fprintf('H%d*',k)
end
 fprintf('A = I\n',n)
 
 clear i j k  A0
