  function [L,U] = JKI1(A) 
%--------------------
% function [L,U] = JKI1(A) 
%--------------------
 n = size(A,1) ;
 L = zeros(n,n);   %% 
 U = zeros(n,n);   %% 
%%->L and U can occupy a single sparse matrix called LU
%% instead start as LU = sparse(n,n) - then follow notes
%%-------------------- main loop - proceeds by columns
 for j=1:n 
   w = A(:,j); 
%%NOTE: w is sparse and the next loop implicitly
%% solves a lower triangular solve with w as RHS (notes)
%% needs to be processed in Top. order. 
   for k=1:j-1 
     w = w - w(k)*L(:,k) ;
   end
   L(j+1:n,j) =   w(j+1:n) / w(j); 
   U(1:j,j)   =   w(1:j); 
 end
 %% instead of this:
 L = eye(n,n)+L;
 %% extract L and U from matrix LU -
 
