%run8.m ====> Approximate Inverse Preconditioners' GMRES test
%_______________________________________________
%HERE spai2  => SPAI preocnditioner
%     waz_t  => Factorized SPAI
%     def_no => Deflation
%e.g.  for k=0:3, figure;eval(['run8(' num2str(k) ')']); end

function run8(iprob)
if ~exist('matrix1.mat'); error('Get matrix1.mat FIRST'); end
if ~exist('fwt.m'); error('Get fwt.m FIRST'); end
if ~exist('gmres.m'); error('Get gmres.m FIRST'); end
   help run8
warning off MATLAB:flops:UnavailableFunction
    
if nargin<1
iprob=input('ch8 : which test matrix (0-3) ? ')
end
 eval(['load matrix' num2str(iprob)])
   Matrix, n=size(A,1); nnzA=nnz(A); 

rand('state',4210)
rhs = rand(n,1); if iprob==3;rhs = rhs+i*rand(n,1); end

if iprob==2
   rhs=diag(diag(A))\rhs; A=diag(diag(A))\A; % scaling
end

disp('_______________________________No Prec: Ax=b___________')
[x0,flag,relres,iter,res0]=gmres(A,rhs,50,1.0e-6,4);
     l0=length(res0);

 lev=3; dwt=4;
 mu=min(max(floor(2*nnzA/n),0), (dwt-1)*2^(lev-1));
 mu=max(mu,5),
 if iprob==3, mu=11; end
 A1=fwts(A,dwt,1,1, lev); b1=fwt(rhs,dwt,1,0, lev);
disp('_______________________________DWT 4____________________')
 r=n:-1:1; 
 B=banda(A1(r,r),mu,mu); [L U]=lu(B);
[xA,flag,relres,iter,resA]=gmres(A1(r,r),b1,50,1.0e-6,4,L,U);
      lA=length(resA);
disp('_______________________________DWTPer 4_________________')
 Ap = perm2(A1,1,1,lev); bp = perm2(b1,1,0,lev);
 r=n:-1:1; Ap=Ap(r,r);
 B=banda(Ap,mu,mu,mu,mu); [L U]=lu(B);
[xp,flag,relres,iter,resp]=gmres(Ap,bp,50,1.0e-6,4,L,U);
      lp=length(resp);

disp('_______________________________DWTPerMod 4_________________')
 Am = perm0(A1,1,1,lev,0); bm = perm0(b1,1,0,lev,0);
 r=n:-1:1; Am=Am(r,r);
 B=banda(Am,mu,mu,mu,mu); [L U]=lu(B);
[xm,flag,relres,iter,resm]=gmres(Am,bm,50,1.0e-6,4,L,U);
      lm=length(resm);

disp('_______________________________DWTPer1 4_________________')
 Am = perm0(A1,1,1,lev,1); bm = perm0(b1,1,0,lev,1);
 r=n:-1:1; Am=Am(r,r);
 B=banda(Am,mu,mu,mu,mu); [L U]=lu(B);
[x1,flag,relres,iter,res1]=gmres(Am,bm,50,1.0e-6,4,L,U);
      l1=length(res1);

hold off
  res0=res0/res0(1); res1=res1/res1(1); resA=resA/resA(1); 
  resp=resp/resp(1); resm=resm/resm(1); 
  
semilogy(l0,res0(end),'r*',1:l0,res0,'r-','markersize',12), hold on 
semilogy(lA,resA(end),'gv',1:lA,resA,'g-','markersize',12), hold on 
semilogy(lp,resp(end),'bd',1:lp,resp,'b-','markersize',12), hold on 
semilogy(lm,resm(end),'ms',1:lm,resm,'m-','markersize',12), hold on
semilogy(l1,res1(end),'c^',1:l1,res1,'c-','markersize',12)
 legend( '1','No precond', ...
         '2','DWT', ...
         '3','DWTPer', ...
         '4','DWTPerMod', ...
         '5','DWTPer1', 0)
 hold off, axis tight
 xlabel(['Test matrix' num2str(iprob) ': Cgce steps (mu=' num2str(mu)...
           ')']); ylabel('Residuals')
 s=text(60,0.1,Matrix);set(s,'Fontsize',12,'Color','b')
