function [x,itno,lrrnrm] = demo_pcg(x, A, b, precond, blk_size, tol, itmax, ic_opts)
%
% Implements the conjugate gradient method with selected optional
% preconditioners to produce x such that ||b-Ax|| <= atol, where 
% atol = tol*(initial residual norm). This routine is especially set up 
% for use with pcg_demo.m but can also be used independently to a 
% limited extent.
%
% The inputs are 
% Required:
%   x         = initial approximate solution 
%   A         = coefficient matrix
%   b         = right-hand side
% Optional:
%   precond   = string specifying name of selected preconditioner; 
%               acceptable values are 'none', 'Jacobi', 'block Jacobi', 
%               'SSOR', and 'incomplete Cholesky'. If 'block Jacobi' is 
%               specified, then the next argument m = block size must also 
%               be given. For some problems, e.g., Poisson problems, 
%	  	 'ADI' is also acceptable; in this case, the user must 
%	 	 define "horizontal" and "vertical" operators H_GL 
%		 and V_GL and declare them to be global variables. 
%   blk_size  = block size, necessary (and referenced) if and only if 
%		precond = 'block Jacobi'; if another preconditioner is 
%		used and values of tol and itmax are specified, this 
%		can be []. 
%   tol       = stopping tolerance; terminate when the residual norm 
%                 is <= tol*(initial residual norm). 
%   itmax  = maximum allowable number of iterations routine. 
%   ic_opts = options structure for MATLAB's incomplete Cholesky ichol. 
%                  If precond is 'incomplete Cholesky' and ichol is empty, 
%                  then the preconditioner is classical incomplete Cholesky 
%                  with no fill-in. If precond is not 'incomplete Cholesky', 
%                  then ic_opts is ignored. 
%
% The returns are 
%   x       = final iterate
%   itno    = vector of iteration numbers
%   lrrnrm  = vector of log relative residual norms
%

global H_GL V_GL; % (for use in ADI preconditioning)
%
% Check inputs and initialize.
%
n = length(x);
if (n ~= length(b) | [n,n] ~= size(A))
  error('demo_pcg: A, x, and b must have consistent dimensions.')
end
if nargin < 3, 
  error('demo_pcg: The number of arguments must be => 3; see comments.'); 
elseif nargin > 8, 
  error('demo_pcg: The number of arguments must be <= 8; see comments.'); 
end 
%
% Check preconditioner inputs and initialize. 
%
if nargin < 4, 
  precond = 'none'; 
elseif isempty(precond), 
  precond = 'none'; 
end 
if ~(strcmp(precond,'none')|strcmp(precond,'Jacobi')|...
     strcmp(precond,'block Jacobi')|strcmp(precond,'SSOR')|...
     strcmp(precond,'incomplete Cholesky')|strcmp(precond,'ADI')) 
   error('demo_pcg: Preconditioner name is not an acceptable string. See comments.');
else
  if strcmp(precond,'Jacobi')
    M = diag(diag(A)); 
  end
  if strcmp(precond,'block Jacobi') 
    if nargin < 5 
      error('demo_pcg: Block size must be specified with block Jacobi.');
    elseif isempty(blk_size) 
      error('demo_pcg: Block size must be specified with block Jacobi.');
    end
    m = blk_size; 
    k = fix(n/m); 
    if k*m ~= n
      error('demo_pcg: The order of A must be divisible by blk_size.');
    end
    if issparse(A)
      M = sparse(n,n); 
    else
      M = zeros(n,n);
    end
    for j = 1:k 
      M((j-1)*m+1:j*m,(j-1)*m+1:j*m) = A((j-1)*m+1:j*m,(j-1)*m+1:j*m); 
    end
    R = chol(M); 
  end
  if strcmp(precond,'SSOR') 
    % Sets up SSOR preconditioner with omega = 1.
    L = tril(A);
    D = diag(diag(A)); 
  end
  if strcmp(precond,'incomplete Cholesky')
%fprintf('\n m = %d. Begin incomplete Cholesky  ', sqrt(n)); tic 
if isempty(ic_opts), 
  L = ichol(A);
else
  L = ichol(A,ic_opts);
end
%fprintf('... end incomplete Cholesky. Time = %f\n', toc); 
  end
  if strcmp(precond,'ADI')
% NOTE: The following are pinhead choices of rho_adi. 
     rho_adi = .05; 
%     rho_adi = .02; 
%     rho_adi = 1/sqrt(size(H_GL,1));
     rhoI = rho_adi*speye(size(H_GL)); 
     MH = H_GL + rhoI;
     MV = V_GL + rhoI;
  end
end

if nargin < 6 
  tol = 1.e-6;
elseif isempty(tol), 
  tol = 1.e-6;
end
if nargin < 7, 
  itmax = max(200,sqrt(n)); 
elseif isempty(itmax), 
  itmax = max(200,sqrt(n)); 
end

%
% Initialize for the iterations. 
%
r = b - A*x;
if strcmp(precond,'none') 
  w = r;
else
  if strcmp(precond,'Jacobi')
    w = M\r;
  end
  if strcmp(precond,'block Jacobi')
    w = R\(R'\r);
  end
  if strcmp(precond,'SSOR')
    w = L'\(D*(L\r));  
  end
  if strcmp(precond,'incomplete Cholesky')
    w = L'\(L\r); 
  end
  if strcmp(precond,'ADI')
    w = rho_adi*(MV\(MH\r)); 
  end
end
rhosq = r'*w;
rho = sqrt(rhosq);
z = zeros(size(x));
p = zeros(size(x));
beta = 0;
itno = [0];
rsnrm = norm(r);
rsnrm0 = rsnrm;
atol = tol*rsnrm;
lrrnrm = [0];
j = 0;
%
% Iterate.
%
while rsnrm > atol & j < itmax 
   j = j+1;
   p = w + beta*p;
   Ap = A*p;
   alpha = rhosq/(p'*Ap); 
   z = z + alpha*p;
   r = r - alpha*Ap;
   if strcmp(precond,'none') 
     w = r;
   else
     if strcmp(precond,'Jacobi')
       w = M\r;
     end
     if strcmp(precond,'block Jacobi')
       w = R\(R'\r);
     end
     if strcmp(precond,'SSOR')
       w = L'\(D*(L\r));
     end
     if strcmp(precond,'incomplete Cholesky')
       w = L'\(L\r); 
     end
     if strcmp(precond,'ADI')
       w = rho_adi*(MV\(MH\r)); 
     end
   end
   rtw = r'*w; 
   beta = rtw/rhosq;
   rhosq = rtw; 
   itno = [itno,j];
   rsnrm = norm(r);
   lrrnrm = [lrrnrm,log10(rsnrm/rsnrm0)];
end
%
% Update x.
%
x = x + z;
