%t_band.m ====> Banded Preconditioned GMRES test
%_______________________________________________
%HERE  banda => banded arrow extraction
%      bandb => banded block extraction
%      bandg => banded matrix with a gap
%      ilu_t => threshold ILU

if ~exist('matrix1.mat'); error('Get matrix1.mat FIRST'); end
if ~exist('banda.m'); error('Get banda.m FIRST'); end
if ~exist('gmres.m'); error('Get gmres.m FIRST'); end
   help t_band
    
iprob=input('ch4 : which test matrix (0-3) ? ')
 eval(['load matrix' num2str(iprob)])
     Matrix

rand('state',4210)
rhs = rand(size(A,1),1); if iprob==3;rhs = rhs+i*rand(size(A,1),1); end
disp('-------------------------------Original Ax=b-----------')
[x0,flag,relres,iter,res0]= gmres(A,rhs,50,1.0e-6,4);
flag,l0=length(res0);

disp('-------------------size 1------Banded  PAx=b-----------')
 B=banda(A, 1, 1); [L U]=lu(B); % Won't work due to infinite
[xa,flag,relres,iter,resa]=gmres(A,rhs,50,1.0e-6,4,L,U);
 flag,la=length(resa);

disp('-------------------size far 32-Banded  PAx=b again-----')
 B=bandb(A, -32, 1, 1); [L U]=lu(B); % Should work for infinite
[xb,flag,relres,iter,resb]=gmres(A,rhs,50,1.0e-6,4,L,U);
 flag,lb=length(resb);

disp('-------------------far permuted banda  PAx=b again-----')
 [B p p2]=bandb(A, -32, 1, 1); A=A(p2,p2);
 B=banda(A, 1, 1); [L U]=lu(B); % Should work for infinite
[x2,flag,relres,iter,res2]=gmres(A,rhs,50,1.0e-6,4,L,U);
 flag,l2=length(res2);

disp('-------------------------------ilut - ILU threshold----')
tau=0.1; [L U]=ilu_t(A, 31,tau);
[xt,flag,relres,iter,rest]=gmres(A,rhs,50,1.0e-6,4,L,U);
 flag,lt=length(rest);

hold off
  res0=res0/res0(1); resa=resa/resa(1);
  resb=resb/resb(1); res2=res2/res2(1);
  rest=rest/rest(1); 
semilogy(l0,res0(end),'r*',1:l0,res0,'r-','markersize',12), hold on 
semilogy(la,resa(end),'bd',1:la,resa,'b-','markersize',12), hold on 
semilogy(lb,resb(end),'ms',1:lb,resb,'m-','markersize',12), hold on 
semilogy(l2,res2(end),'go',1:l2,res2,'g-','markersize',12), hold on
semilogy(lt,rest(end),'cv',1:lt,rest,'c-','markersize',12), 
 legend( '1','No precond', ...
         '2','band(1,1)', ...
         '3','Off band(1,1)', ...
         '4','Off band(1,1) permuted', ...
         '5','ilu\_t', 0)
 hold off, axis tight
 xlabel(['Test matrix' num2str(iprob) ': Cgce steps']); ylabel('Residuals')
 s=text(100,0.01,Matrix);set(s,'Fontsize',12,'Color','b')
