function gmres_demo(m,coeffs,RHSfun,precond,tol,itmax)
%
% Demonstrates the performance of GMRES with optional left
% preconditioning on the problem 
%
%             -(Laplacian) u + cu + du_x + eu_y = f  
%
% with Dirichlet boundary conditions on the unit square in R^2 
% discretized on an evenly spaced mxm mesh of interior points in
% the unit square [0,1]x[0,1]. 
%
% m can be either an integer or a vector of integers with components 
% increasing in size. In the first case, GMRES is applied to the problem 
% with an mxm mesh. In the second case, GMRES is applied to
% problems with  m(j)xm(j) meshes,  j = 1, ..., length(m). 
%
% NOTE: If m is an integer vector, make sure the components are 
% increasing in size. 
%
% Currently, the preconditioning options are 'ILU' (incomplete LU
% using MATLAB's ilu command with the default settings) and 'RAS'
% (restricted additive Schwarz). Restricted additive Schwarz is set
% up to use a 4x4 array of 16 subdomains with 3 grid lines of
% overlap.   
%
% The inputs are 
% Required:
%   m         = integer scalar or integer vector determining the 
%               number(s) of interior mesh points per side. 
%   coeffs    = row vector of coefficients [c d e].
%   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). 
%   precond   = string specifying the preconditioner; acceptable 
%               values are 'none', 'ILU' and 'RAS'. 
% Optional:
%   tol       = stopping tolerance on the residual norm. 
%   itmax  = maximum allowable number of GMRES 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 GMRES run is x = 0. The iterations terminate 
% when the residual norm is <= tol*(initial residual norm); the
% default is tol = 1.e-6. For safety, GMRES is allowed to iterate
% for no more than itmax iterations; the default is itmax = 3000. 
% 
% M-files required: 
%   condif.m      (sets up the discretized convection-diffusion operator) 
%   gmres514.m   (performs GMRES and other ancillary tasks)
%   loadrhs.m      (loads the right-hand side)
%   ras.m      (performs restricted additive Schwarz)
% 

global H_GL V_GL; % (for use in ADI preconditioning)

if nargin < 4, 
  error('gmres_demo: The number of arguments must be => 4; see comments.'); 
elseif nargin > 6, 
  error('gmres_demo: The number of arguments must be <= 6; see comments.');
end

if ~(strcmp(precond,'none')|strcmp(precond,'ILU')|strcmp(precond,'RAS'))
  error('gmres_demo: Preconditioner name is not an acceptable string; see comments.');
end

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

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

% Uncomment the following to view the solution: 
%set(figure(2),'Units','normal','Position',[.4,.04,.5,.45]); clf;
%set(figure(1),'Units','normal','Position',[.4,.55,.5,.45]); clf;
fig = gcf; clf; hold on  

% Set the restart value for gmres514. 
restart = 50;

fprintf(['\n GMRES(',num2str(restart),'), preconditioner = ',precond,...
       ',\n Convection-Diffusion Problem: c = ',num2str(coeffs(1)),...
       ', d = ', num2str(coeffs(2)),', e = ', num2str(coeffs(3)),'\n']); 

for j = length(m):-1:1 
  x = zeros(m(j)^2,1);
  A =  condif(m(j),coeffs); 
  b = loadrhs(RHSfun,m(j)); 
  tic; 
    if strcmp(precond,'none')
      M = [];
    end
    if strcmp(precond,'ILU')
      [L,U] = ilu(A);
      M = @(v)U\(L\v);
    end
    if strcmp(precond,'RAS')
      nsd = 4; msd = nsd\m(j); lap = 3;
      M = @(v)reshape(ras(zeros(m(j),m(j)),coeffs,nsd,msd,lap,reshape(v,m(j),m(j))),m(j)^2,1); 
    end
    [x,itno,lrrnrm] = gmres514(A,b,x,restart,M,tol,itmax);
    time = toc;
    plot(itno,lrrnrm);  
    text(itno(size(itno,2)),lrrnrm(size(lrrnrm,2))+.1,['m=' num2str(m(j))]);
%    text(itno(size(itno,2)),lrrnrm(size(lrrnrm,2))+.4,['m=' num2str(m(j))]);
%    text(itno(size(itno,2)),lrrnrm(size(lrrnrm,2)),[num2str(time) 'secs']);
    drawnow 
    fprintf('\n For m = %d:  %d iterations; time = %g seconds', m(j), length(itno), time);
    if j == length(m)
       y = x; 
    end
end 
fprintf('\n\n');
title(['GMRES(',num2str(restart),'), preconditioner = ',precond,...
       ', Convection-Diffusion: c = ',num2str(coeffs(1)),...
       ', d = ', num2str(coeffs(2)),', e = ', num2str(coeffs(3))]); 

xlabel('iteration number'); 
ylabel('log_{10} relative residual norm'); 
% Uncomment the following to view the solution: 
%figure(2);seesol(y,m(length(m)));set(figure(2),'Units','pixels');
%set(figure(1),'Units','pixels');
drawnow
hold off
