function resid=s8a3(n,pcflag)

% generate A (ignore factor of -1/h^2, we only want to test cgpc...)
M  = 2*n-1; N  = n-1;	% interior grid is of size M x N
NB = N;     NA = M*N;	% sizes of B and A
h  = 1/(2*n);           % meshwidth
d  = 4*ones(NA,1);                      % diagonal
du = -ones(NA-1,1); du(NB:NB:NA-n+1)=0; % upper and lower diagonal
do = -ones((M-1)*NB,1);         % off blocks
A  = sparse(diag(d)+diag(du,1)+diag(du,-1)+diag(do,n-1)+diag(do,-(n-1)));
b  = ones(NA,1);		% rhs

% solve with preconditioned cg
[u,resid] = cgpc(A,-b);

% plot the convergence

plot(0:length(resid)-1,log10(resid),'*');
xlabel('iteration k'); ylabel('norm of residual');

% plot of voltage u
figure;
a = 0.5; b = 1;x=0:h:a; y=0:h:b;
[X,Y] = meshgrid(0:h:a,0:h:b);
Z = reshape(u,N,M)';
Z = [zeros(1,N+2); zeros(M,1) Z zeros(M,1);zeros(1,N+2)];
surf(X,Y,Z); xlabel('x'); ylabel('y');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [x,Rsave]=cgpc(A,b)

% conjugate gradients with incomplete 
% Cholesky preconditioning to solve Ax + b = 0
% in:   A       the (n x n)-Matrix
%       b       the negative rhs
% out:  x       solution
%       Rsave   norm(residuals,2) 

tol = 1e-25; tol2=tol*tol;      % stop if norm(residual)<tol
Rsave = [];                     % save norm of resdiduals 

% init
C=cholinc(A,'0');
n=length(b); x = zeros(n,1); r = b; p = x; Rp = r'*r; Rsave = Rp;
z=x;
%  loop
for k = 0:10*n
  zold = z; 
  z     = C\(C'\r);
  if k==0 sigma = 0; else sigma = z'*r/(zold'*rold); end
  p     = sigma*p - z;
  u     = A*p;
  rho   = z'*r/(p'*u);
  x     = x + rho*p;
  rold  = r;          r     = r + rho*u;
  Rp    = r'*r;
  Rsave = [Rsave sqrt(Rp)]
  if (Rp < tol2) break; end
end
