% t_dete.m ===> Banded Arrow via detect.m GMRES test
%___________________________________________________
%HERE  banda => banded arrow extraction
%      bandb => banded block extraction
%      bandg => banded matrix with a gap
%      ilu_t => threshold ILU

if ~exist('matrix3.mat'); error('Get matrix3.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_dete

iprob=input('ch4 : which test matrix (0-3) ? ')
 eval(['load matrix' num2str(iprob)])
   Matrix, n=size(A,1); nnzA=nnz(A); nnzAp=floor(nnzA/n);

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

  mu=min(floor( (2*nnzA-5*n)/2),5);
  if mu<=0, mu=1; end
disp(['---------Detection----Banded arrow A(r,c)x=b  col=' num2str(mu)])
[r,c,no]=detect(A,ratio,0);
if no==0, [r,c,no]=detect(A,ratio/2,0); end
if no==0, [r,c,no]=detect(A,ratio/4,0); end
B=banda(A(r,c), 2,2, mu,mu); [L U]=lu(B);
[x1,flag,relres,iter,res1]=gmres(A(r,c),rhs,50,1.0e-6,4,L,U);
flag,l1=length(res1); figure; hold off

  mu=min(floor( (2*nnzA-3*n)/2),5);
  if mu<=0, mu=1; end
disp(['---------Detection----Banded arrow A(r,c)x=b  col=' num2str(mu)])
[r,c,no]=detect(A,ratio,0);
if no==0, [r,c,no]=detect(A,ratio/2,0); end
if no==0, [r,c,no]=detect(A,ratio/4,0); end
B=banda(A(r,c), 1,1, 20,20); [L U]=lu(B);
[x2,flag,relres,iter,res2]=gmres(A(r,c),rhs,50,1.0e-6,4,L,U);
flag,l2=length(res2);

disp('-------------ilu0---+detection----------------------------')
[r,c,no]=detect(A,ratio,0);
if no==0, [r,c,no]=detect(A,ratio/2,0); end
if no==0, [r,c,no]=detect(A,ratio/4,0); end
[L U]=ilu_0(A(r,c));
[xt,flag,relres,iter,rest]=gmres(A(r,c),rhs,50,1.0e-6,4,L,U);
flag,lt=length(rest);

disp(['-----------ilut-------No Detection used nz per c=' num2str(nnzAp)...
      '-------'])
[L U]=ilu_t(A,2*nnzA,0.1);
[xto,flag,relres,iter,resto]=gmres(A,rhs,50,1.0e-6,4,L,U);
flag,lto=length(resto);

rest=rest/rest(1); res1=res1/res1(1); res2=res2/res2(1); 
 resto=resto/resto(1);
semilogy(lt,rest(end),'rv',1:lt,rest,'r-','markersize',12), hold on
semilogy(lto,resto(end),'mv',1:lto,resto,'m-','markersize',12), hold on
semilogy(l1,res1(end),'go',1:l1,res1,'g-','markersize',12), hold on
semilogy(l2,res2(end),'bd',1:l2,res2,'b-','markersize',12), 
 legend('1','ilu0',...
        '2 no detect','ilu\_t',...
        '3','fat arrow',...
        '4','slim arrow',0)
 hold off, axis tight
 xlabel(['Test matrix' num2str(iprob) ': Cgce steps']); ylabel('Residuals')
 s=text(100,0.05,Matrix);set(s,'Fontsize',12,'Color','b')
