%WAZ_T.m ---- Bi-Conjugation Algorithm by Fox/Benzi/Tuma (SISC 98) work
%waz_t.m ---- WAZT(p,tau) method for Z*W' \approx A^{-1}
%             W'AZ=I  or inv(A)=ZW'  (Incomplete inverse LU decomposition)
%            
%Usage:  [W Z]=waz_t(A, p,tau);                        % assume A=sparse(A)
%
%               where   p = nnz for each column
%                     tau = threshold             { just like ILUT(p,tau) }
%
%e.g.    n = 16, A=sprand(n,n,.6)+10*speye(n);
%   [W Z]=waz_t(A,n/4,0.05); erz=norm(Z*W'-inv(A),1); nzz=nnz(W'+Z);
%   disp('____________________Comparison with the ILU____________________')
%   [L U]=ilu_0(A);          er0=norm(inv(U)*inv(L)-inv(A),1); nz0=nnz(L+U);
%   [L U]=ilu_t(A,n/4,0.05); ert=norm(inv(U)*inv(L)-inv(A),1); nzt=nnz(L+U);
%                            ert_LU=norm(L*U-A,1);
%   fprintf('\nW''AZ:  |Z*W''-A^{-1}|_1   \t= %f with NNZ = %d\n',erz,nzz)
%   fprintf(' ILU:  |U^{-1}L{-1}-A^{-1}|\t=%f + ILU0 er= %f\n',ert,er0)
%   fprintf(' ILU:  |L*U - A|_1         \t=%f with NNZ = %d\n\n',ert_LU,nzt)
%  
%                %See also waz.m : diag(W)=I like L and Z like U  (in LDU=A)

function [w,z]=waz_t(A, pp, TOL)
if nargin<3, help waz_t; return; end
n=size(A,1);

%========================================================================
%|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_ Step 1   SISC(1998,v.19,p.9)
   w = eye(n); z=w;  p=[];q=[];

%|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_ Step 2   Algorithm BC

for i=1:n %==================Global start

   pq=w(1:i,i)'*A(1:i,i);  % Orig L Fox method works
   v_i = TOL*norm(A(i,:),inf);

 for j=i+1:n
   q(j)=w(1:j,j)'*A(1:j,i)/pq;  % Benzi/Tuma (sisc98 + JCP02)
   p(j)=z(1:j,j)'*A(i,1:j)'; 
 end
   for k=1:i, z(k,i)=z(k,i)/pq; end  % Unit diagonals in Z

   for j=i+1:n
    for k=1:i
     wk_tmp=w(k,j) - q(j)*w(k,i);
       if abs(wk_tmp)<=v_i, wk_tmp=0; end
     w(k,j)=wk_tmp; 
     wk_tmp=z(k,j) - p(j)*z(k,i);
       if abs(wk_tmp)<=v_i, wk_tmp=0; end
      z(k,j)=wk_tmp; 
    end
   end %j
%___________________________________________________________ Enforce pp
    J = find( abs(z(1:i-1,i))>eps  );  %--------- Extract NZ's Z
    if length(J)>pp-1,
     [d IU]=sort(-abs(z(J,i))); J=J(IU(1:pp-1)); 
     Filter = sparse(i-1,1); Filter(J)=z(J,i); z(1:i-1,i)=Filter;
    end
    J = find( abs(w(1:i-1,i))>eps  );  %--------- Extract NZ's W
    if length(J)>pp-1
    [d IL]=sort(-abs(w(J,i))); J=J(IL(1:pp-1)); 
    Filter = sparse(i-1,1); Filter(J)=w(J,i); w(1:i-1,i)=Filter;
    end
if rem(i-1,round(n/5))==0
 fprintf('\t  %3d%%  Done  <%s.m>\n',round(i/n*100),mfilename)
end
end %i ========================Global end
