%run9.m ====> Recursive Wavelet Schur Preconditioners' GMRES test
%_______________________________________________
%HERE gmres_nr.m for Flexible GMRES / Wavelet ML Schur / Richardson (I)
%     richa_nr.m for Richardson / Wavelet Schur / RichardsonA       (II)
%
% Usage: run9(prob_no, levs, gamma, method, restart)
%   or : run9(prob_no, levs, gamma)
%   or : run9(prob_no, levs)
%   or : run9(prob_no)
%
% for  prob_no = Problem number (0-3) 
%          levs = number of wavelet levels (order is set inside) say 2
%           gamma = cycling pattern (ML) say 2 or 1
%           method = Method number 1 (I) or 2 (II)
%           restart = GMRES(restart)
%
%e.g  run9(1,5,2,2)
% or  for k=1:3, figure;eval(['run9(' num2str(k) ',4,2,1)']); end

function run9(iprob, lev, gamma, meth, m) % m for restart in gmres(m)
% Internal notation mu: cycling pattern
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
 global order levs mu prob
 clear R* M* A* B* C* b* z* x* i
    
if nargin<1, help run9
iprob=input('ch9 : which test matrix (0-3) ? ');
end
  eval(['load matrix' num2str(iprob)]), n=size(A,1); nnzA=nnz(A); 
if nargin<2, lev= floor(log(n)/log(2))-7; end
if nargin<3, gamma=1; end
if nargin<4, help run9
 meth=input('ch9 : which method (1-2) ? ');
end
if nargin<5, m=25; end
   Matrix
   levs=lev; mu=gamma; method=meth; steps=ceil(n/m); %Sent to Global/Use later
   prob=iprob; restart=m; order=4; % DWT order 

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

if ~exist('A'), disp('No A?'); return; end
%|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_|_| SOLVER_NOW |_| |_|
  n = length(A); 
  if iprob<=2
     b=diag(diag(A))\b; A=diag(diag(A))\A;
  end %% 2-stage work
  clear R* N* D*

if method==1  % Recommended CC(02) Algorithm I GMRES+Richardson(mu)
 fprintf('gmres_nr.m for prob = %d with N=%4d lev=%d  gamma=%d\n', ... \
            prob,n, levs,mu)
 [xo,res,git]=gmres_nr(A,x,b,restart, steps,1.0e-6);
 s=sprintf('gmres_nr.m for prob = %d with N=%4d lev=%d\n',prob,n, levs);

elseif method==2 % CC(02) Alg I (Richadson(mu) all)
steps,restart
    fprintf('richa_nr.m for prob = %d with N=%4d lev=%d  gamma=%d\n', ... \
            prob,n, levs,mu)
 [xo,res,git]=richa_nr(A,x,b,steps*restart,1.0e-6); 
 s=sprintf('richa_nr.m for prob = %d with N=%4d lev=%d\n',prob,n, levs);

else, disp('method out of range ?'); return
end
  semilogy(res,'-ro'); hold on;title(s)
    semilogy([1 max(length(res),40)], [0.1 1.0E-7],'g:'); hold off
 s = sprintf('%d Iterations of a %d-level method for N=%d (%s)', git,levs,n, Matrix);
 s = sprintf('Residuals for Prob %d',prob);  ylabel(s)
xlabel(['Test matrix' num2str(iprob) ': Cgce steps (mu=' num2str(mu) ', M'...
        num2str(method) ', levs=' num2str(levs) ')']); 
 s=text(20,0.5,Matrix);set(s,'Fontsize',12,'Color','b')
