function pcg_demo(m,problem,RHSfun,precond,tol,itmax,ic_opts)
%
% Demonstrates the performance of PCG on the problem 
%
%             -(Laplacian) u = f  or (Laplacian)^2 u = 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, PCG is applied to the problem 
% with an mxm mesh. In the second case, PCG 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. 
%
% The inputs are 
% Required:
%   m         = integer scalar or integer vector determining the 
%               number(s) of interior mesh points per side. 
%   problem   = string specifying name of the problem; acceptable values 
%               are 'Poisson' and 'biharmonic'. 
%   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', 'Jacobi', 'block Jacobi', 'SSOR', 
%		'incomplete Cholesky', and -- for the Poisson 
%		problem ONLY -- 'ADI'. 
%
% Optional:
%   tol       = stopping tolerance on the residual norm. 
%   itmax  = maximum allowable number of PCG iterations. 
%   ic_opts = options structure for MATLAB's incomplete Cholesky routine
%                   ichol called in demo_pcg; type "doc ichol" for information on 
%                  possible fields. If precond is 'incomplete Cholesky' and ichol 
%                  is not provided, then the preconditioner is classical incomplete 
%                  Cholesky with no fill-in. If precond is not 'incomplete Cholesky', 
%                  then ic_opts is ignored even if it is provided. 
%
% 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 PCG 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, PCG is allowed to iterate
% for no more than  itmax iterations; the default is itmax = 3000. 
% 
% M-files required: 
%   biharmonic.m          (sets up the discretized biharmonic operator)
%   laplacian.m           (sets up the discretized Laplacian)
%   loadrhs.m		  (loads the right-hand side f)
%   demo_pcg.m		  (performs PCG and other ancillary tasks)
%   seesol.m		  (optional: plots the solution surface, requires uncommenting lines below) 
% 

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

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

if ~(strcmp(precond,'none')|strcmp(precond,'Jacobi')|...
     strcmp(precond,'block Jacobi')|strcmp(precond,'SSOR')|...
     strcmp(precond,'incomplete Cholesky')|strcmp(precond,'ADI')) 
  error('pcg_demo: Preconditioner name is not an acceptable string; see comments.');
end
if ~(strcmp(problem,'Poisson')|strcmp(problem,'biharmonic'))
  error('pcg_demo: Problem name is not an acceptable string; see comments.');
end
if strcmp(problem,'biharmonic')&strcmp(precond,'ADI') 
  error ('pcg_demo: Can-t use ADI preconditioning on the biharmonic problem.'); 
end

if nargin < 7, 
  ic_opts = [ ]; 
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 

fig = gcf;
set(fig,'Units','normal','Position',[.4,.55,.5,.45]); clf;
figure(fig); hold on 

for j = length(m):-1:1 
    x = zeros(m(j)^2,1);
    if strcmp(problem,'Poisson') 
      [A,H_GL,V_GL] = laplacian(m(j));
      A = -A; H_GL = -H_GL; V_GL = -V_GL; 
    end
    if strcmp(problem,'biharmonic')
      A = biharmonic_op(m(j));
    end
    b = loadrhs(RHSfun,m(j)); 
    tic; 
    [x,itno,lrrnrm] = demo_pcg(x, A, b, precond, m(j), tol, itmax,ic_opts); 
    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(['PCG, preconditioner = ',precond,', ',problem, ' problem']); 
xlabel('iteration number'); 
ylabel('log_{10} relative residual norm'); 
% Uncomment the following to view the solution: 
%  figsol = figure(fig.Number+1);
%  set(figsol,'Units','normal','Position',[.4,.04,.5,.45]); clf;
%  figsol;seesol(y,m(length(m)));set(figsol,'Units','pixels');

set(fig,'Units','pixels');
hold off
