function ks_method_demo(m,coeffs,method,precond,RHSfun,restart,tol,itmax)
%
% Demonstrates the performance of Krylov subspace methods on the 
% convection-diffusion problem 
%
%   -(Laplacian)u + cu + du_x + eu_y = f   in D
%                                  u = 0   on the boundary of D
%
% where D is the unit square in R^2 discretized on an evenly 
% spaced square mesh of interior points. 
%
% Specifically, GMRES, BCG, BiCGSTAB, CGS, or QMR is applied 
% to the problem for a given c,d,e,and f and square meshes with 
% size determined by m. 
%
% m can be either an integer or a vector of integers with components 
% increasing in size. In the first case, the method is applied to 
% the problem with an mXm mesh. In the second case, the method is 
% applied to problems with m(j)Xm(j) meshes, j = 1, ..., length(m). 
%
% NOTE: If m is an integer vector, the components must be increasing 
% in size. 
%
% NOTE: The demo uses the Matlab routines for the Krylov subspace 
% methods. 
%
% The inputs are 
% Required:
%   coeffs    = vector of coefficients [c d e].
%   m         = scalar integer or integer vector determining the 
%               number(s) of interior mesh points per side. 
%   method    = string specifying the method; acceptable values are 
%		'GMRES', 'BCG', 'BiCGSTAB', 'CGS', and 'QMR'.
%   precond   = string specifying the preconditioner; acceptable 
%               values are 'none', 'Jacobi', 'block Jacobi', 'ILU', 
%               and 'ADI'. 
%   RHSfun    = function handle (@functionname) for evaluating the 
%               right-hand side function f at given (x,y). The routine 
%               should accept scalar x and y as arguments and return f(x,y). 
% Optional:
%   restart   = number of GMRES steps before restarting. If GMRES is 
%		used and this is not specified, then the method does 
%		not restart, i.e., in effect uses restart = n; this 
%		may require undesirably large storage if m is large. 
%		If a method other than GMRES is used and values of tol 
%		and itmax are specified, then restart can be []. 
%   tol       = stopping tolerance on the residual norm; stop when 
%              residual norm <= tol*(*initial residual norm). 
%   itmax     = maximum allowable number of iterations. 
%
% The output is a graph of iteration numbers versus log residual norms 
% for each mesh specified by m. The program also prints out the run time 
% (in seconds) for each value of m. 
%
% The initial guess for each run is x = 0. The iterations terminate 
% when the residual norm is less than or equal to tol*norm(b), 
% where b is the vector of f-values; the default is tol = 1.e-4. 
% For safety, the methods are allowed to iterate for no more than 
% itmax iterations; the default is itmax = 3000. 
%
% M-files required: 
%   adi_pre.m 
%   condif.m
%   RHSfun.m (called by loadrhs.m) 
%   loadrhs.m 
%

% For ADI preconditioning: Declare global variables and 
% set adi_init_flag_GL to a nonzero value to initialize 
% the ADI precontioner. 
global H_GL V_GL adi_init_flag_GL; 
adi_init_flag_GL = 1; 

if nargin < 5, 
  error('ks_method_demo: The number of arguments must be => 5; see comments.'); 
elseif nargin > 8, 
  error('ks_method_demo: The number of arguments must be <= 8; see comments.');
end

if ~(norm(size(coeffs)-[1 3])==0|norm(size(coeffs)-[3 1])==0), 
  error('ks_method_demo: coeffs must be a vector of length 3; see comments.'); 
end

if ~(strcmp(method,'GMRES')|strcmp(method,'BCG')|...
     strcmp(method,'BiCGSTAB')|strcmp(method,'CGS')|strcmp(method,'QMR')), 
  error('ks_method_demo: Method name is not an acceptable string; see comments.');
end
if ~(strcmp(precond,'none')|strcmp(precond,'Jacobi')|...
     strcmp(precond,'block Jacobi')|...
     strcmp(precond,'ILU')|strcmp(precond,'ADI')) 
  error('ks_method_demo: Preconditioner name is not an acceptable string. See comments.');
end

if nargin < 7, 
  tol = 1.e-4;
elseif isempty(tol) 
  tol = 1.e-4;
end 

if nargin < 8, 
  itmax = 3000; 
elseif isempty(itmax) 
  itmax = 3000; 
end 

if strcmp(method,'GMRES'),
  if nargin < 6, 
    restart = m^2;
  elseif isempty(restart), 
    restart = m^2;
  end
  maxit = floor(restart\itmax); 
end

if strcmp(method,'GMRES'), 
  fprintf('\n\n GMRES(%d), preconditioner = %s',restart,precond);
else
  fprintf('\n\n %s, preconditioner = %s',method,precond);
end
fprintf('\n Convection-diffusion problem, c = %g, d = %g, e = %g', coeffs(1), coeffs(2), coeffs(3)); 

set(figure(1),'Units','normal','Position',[.4,.55,.5,.45]); clf;
%set(figure(2),'Units','normal','Position',[.4,.04,.5,.45]); clf;
figure(1); hold on  

for j = length(m):-1:1 
    x = zeros(m(j)^2,1);
    [A,H_GL,V_GL] = condif(m(j),coeffs);
    b = loadrhs(RHSfun,m(j));
    tic; 
    if strcmp(precond,'none'), 
      if strcmp(method,'GMRES'), 
        [x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit); 
      elseif strcmp(method,'BCG'), 
        [x,flag,relres,iter,resvec] = bicg(A,b,tol,itmax); 
      elseif strcmp(method,'BiCGSTAB'), 
        [x,flag,relres,iter,resvec] = bicgstab(A,b,tol,itmax); 
      elseif strcmp(method,'CGS'),
        [x,flag,relres,iter,resvec] = cgs(A,b,tol,itmax); 
      else
        [x,flag,relres,iter,resvec] = qmr(A,b,tol,itmax); 
      end
    elseif strcmp(precond,'Jacobi'), 
      M = diag(diag(A)); 
      if strcmp(method,'GMRES'), 
        [x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,M); 
      elseif strcmp(method,'BCG'), 
        [x,flag,relres,iter,resvec] = bicg(A,b,tol,itmax,M); 
      elseif strcmp(method,'BiCGSTAB'), 
        [x,flag,relres,iter,resvec] = bicgstab(A,b,tol,itmax,M); 
      elseif strcmp(method,'CGS'),
        [x,flag,relres,iter,resvec] = cgs(A,b,tol,itmax,M); 
      else
        [x,flag,relres,iter,resvec] = qmr(A,b,tol,itmax,M); 
      end
    elseif strcmp(precond,'block Jacobi')|strcmp(precond,'ILU'), 
      if strcmp(precond,'block Jacobi'), 
        n = m(j)^2; 
        k = fix(m(j)\n); 
        if k*m(j) ~= n
          error('ks_method_demo: The order of A must be divisible by the block size.');
        end
        if issparse(A)
          M = sparse(n,n); 
        else
          M = zeros(n,n);
        end
        for i = 1:k 
          M((i-1)*m(j)+1:i*m(j),(i-1)*m(j)+1:i*m(j)) = ... 
                            A((i-1)*m(j)+1:i*m(j),(i-1)*m(j)+1:i*m(j)); 
        end
        [M1,M2] = lu(M); 
      else 
        [M1,M2] = ilu(A);
      end
      if strcmp(method,'GMRES'), 
        [x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,M1,M2); 
      elseif strcmp(method,'BCG'), 
        [x,flag,relres,iter,resvec] = bicg(A,b,tol,itmax,M1,M2); 
      elseif strcmp(method,'BiCGSTAB'), 
        [x,flag,relres,iter,resvec] = bicgstab(A,b,tol,itmax,M1,M2); 
      elseif strcmp(method,'CGS'),
        [x,flag,relres,iter,resvec] = cgs(A,b,tol,itmax,M1,M2); 
      else
        [x,flag,relres,iter,resvec] = qmr(A,b,tol,itmax,M1,M2); 
      end
    elseif strcmp(precond,'ADI'), 
      if strcmp(method,'GMRES'), 
        [x,flag,relres,iter,resvec] = gmres(A,b,restart,tol,maxit,@adi_pre); 
      elseif strcmp(method,'BCG'), 
        [x,flag,relres,iter,resvec] = bicg(A,b,tol,itmax,@adi_pre); 
      elseif strcmp(method,'BiCGSTAB'), 
        [x,flag,relres,iter,resvec] = bicgstab(A,b,tol,itmax,@adi_pre); 
      elseif strcmp(method,'CGS'),
        [x,flag,relres,iter,resvec] = cgs(A,b,tol,itmax,@adi_pre); 
      else
        [x,flag,relres,iter,resvec] = qmr(A,b,tol,itmax,@adi_pre); 
      end
    end
    time = toc;
    itno = size(resvec,1)-1; 
    plot([0:1:itno]',log10(resvec(1,1)\resvec)); 
    %text(itno,log(resvec(size(resvec,1),1))+.01,['m=' num2str(m(j))]); 
    text(itno,log10(resvec(1,1)\resvec(size(resvec,1),1))+.1,... 
         ['m=' num2str(m(j))]); 
    drawnow; resvec = 0; 
    if strcmp(method,'GMRES'), 
      fprintf('\n For m = %d:  time = %g sec., its. = %d, relres = %g, flag = %d', ...
          m(j), time, restart*(iter(1)-1)+iter(2), relres, flag); 
    else
      fprintf('\n For m = %d:  time = %g sec., its. = %d, relres = %g, flag = %d', ...
          m(j), time, iter, relres, flag); 
    end
    if j == length(m)
       y = x; 
    end
end 

fprintf('\n\n');

if strcmp(method,'GMRES'), 
  title([method,'(',num2str(restart),')' ', preconditioner = ', precond,...
          ', convection-diffusion problem -{\Delta}u+ ', ...
         num2str(coeffs(1)),'u+',num2str(coeffs(2)),'u_x+',...
         num2str(coeffs(3)),'u_y = f.']); 
else 
  title([method, ', preconditioner = ', precond,...
          ', convection-diffusion problem -{\Delta}u+ ', ...
         num2str(coeffs(1)),'u+',num2str(coeffs(2)),'u_x+',...
         num2str(coeffs(3)),'u_y = f.']); 
end 

xlabel('iteration number'); 
ylabel('log_{10} relative residual norm'); 
set(figure(1),'Units','pixels');
hold off
