%WAZ_all.m - Bi-Conjugation Algorithm by Fox/Benzi/Tuma (SISC 98) work
%            W'AZ=I   (inverse LU decomposition)
%---------- Here we set up A as dense matrix for illustration only
% USAGE:  Input A on line before typing:   waz_all
%
% @ K Chen (Nov 2002)
%__________________________________________________________________________ 
help WAZ_ALL
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)

%========================================================================

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

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

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

   pq=w(1:i,i)'*A(1:i,i);  % Orig L Fox method works
   if i<n
   eval(['w' num2str(i) '=eye(n);'])
   end
   eval(['z' num2str(i) '=eye(n);'])

 for j=i+1:n
   qj=w(1:j,j)'*A(1:j,i);  % Orig L Fox method works
   pj=z(1:j,j)'*A(i,1:j)';   p_ji=pj/pq; q_ji=qj/pq;
   eval(['w' num2str(i) '(j,i)=-q_ji;'])
   eval(['z' num2str(i) '(i,j)=-p_ji;'])

   for k=1:j-1
     w(k,j)=w(k,j) - q_ji*w(k,i);
     z(k,j)=z(k,j) - p_ji*z(k,i);
   end
 end

  for k=1:i
    z(k,i)=z(k,i)/pq;
  end

   eval(['z' num2str(i) '(i,i)=1/pq;'])

end %========================Global end

%|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_ Step 3
disp('------------ Exact factorisation of A done by WAZ_all.m ---')
disp('Check: w''*A*z-eye(n), inv(A)-z*w'', also [L D U]=lu(A) etc')
disp('--_--------- System Ax=b becomes w''*A*z*y = w''*b --------')
fprintf('\nThe bi-factorisation of A  is:    w''*A*z = eye(n)\n   ')
fprintf('\nwhere one may check this overall decomposition:\n\t')

for k=n-1:-1:1
 fprintf('w%d*',k)
end
 fprintf('A')
for k=1:n
 fprintf('*z%d',k)
end
 fprintf(' = eye(%d)\n',n)
 
disp(' '),    clear i j k  q_ji p_ji pq
