% ilu_t.m ----- ILUT(p,tau) method for L*U \approx A
% Usage:  [L U]=ilu_t(A, p,tau);
%                assume A=sparse(A)
% e.g. rand('seed',0); A=rand(8); A=A.*(A>.4) + diag(rand(8,1))  % sparse A
%      [L U]=ilu_0(A);          er1=norm(L*U-A,1), nz1=nnz(L+U)
%      [L U]=ilu_t(A,4,0.05);   er2=norm(L*U-A,1), nz2=nnz(L+U)
function [L,U]=ILU_t_method(A, p,tau)
    if nargin<3, help ilu_t; return;end; 
n = size(A,1); 
    TOL=1.0E-12;   v=[]; B=sparse(n,n);
for k=1:n,  v(k)=tau*norm(A(k,:),inf);  end

for i=1:n
    w=A(i,:);
  for k=1:(i-1)
  if abs(w(k))>TOL %-------------------- pattern nz(A)
    w(k) = w(k) / B(k,k);
    if abs(w(k))<v(i), w(k)=0; end     % L threshold

    if abs(w(k))>=v(i) % Step (7)
     J=(k+1):n;
     tmp=w(J) - w(k)*B(k,J); if abs(tmp)<=v(i),tmp=0; end
     w(J) = tmp;
    end %if

  end %if TOL
  end %k 
    J = find( abs(w)>eps  );  %--------- Extract NZ's
   JL = intersect(1:i-1, J);  JU = intersect((i+1):n, J);
    if length(JL)>p, [d IL]=sort(-abs(w(JL))); JL=JL(IL(1:p)); end
  if length(JU)>p-1, [d IU]=sort(-abs(w(JU))); JU=JU(IU(1:p-1)); end

    J = [JL i JU]; % indexes to be updated by new row (diag included)
    Filter = sparse(1,n); Filter(J)=w(J);
    B(i,:) = Filter ;     % row i updated (L/U)

end %i 
  L = full( tril(B,-1) ) + speye(n);
  U = full( triu(B) );
