function P2()
% Checks out NinePoint
clc
n1 = 150;
n2 = 100;
a = -2; b = -4 ; c = -4*(a+b);
f = randn(n1*n2,1);
Tol = .00000001;
MaxIt = 300;

for DropTol = [.001  .01  .1  1]
   fprintf('\n********************************************\nDropTol = %6.3f\n\n',DropTol)
   fprintf('CVL: \n')
   tic
   uCVL = NinePointCVL(n1,n2,a,b,c,f,Tol,MaxIt,DropTol);
   tCVL = toc;
   fprintf('tCVL = %5.3f\n\n',tCVL)
   fprintf('You: \n')
   tic
   uYou = NinePoint(n1,n2,a,b,c,f,Tol,MaxIt,DropTol);
   tYou = toc;
   fprintf('tYou = %5.3f\n\n',tYou)
   fprintf('norm(uCVL-uYou)/norm(uCVL) = %10.3e\n',norm(uCVL-uYou)/norm(uCVL))
end
   

function u = NinePointCVL(n1,n2,a,b,c,f,Tol,MaxIt,DropTol)
% f is a n1*n2-by-1 column vector
% c > 0, 4(a+b) + c = 0.
% Solves the system Au = f where A = A(n1,n2,a,b,c)
% Uses pcg with tolerance Tol and iteration maximum maxIt.
% Uses the incomplete cholesky factor produced by cholinc with drop
% tolerance DropTol

%  Get A in sparse format and compute its incomplete upper Cholesky factor R
% so A approx R'R
A = MakeSparseA(n1,n2,a,b,c);
R = cholinc(A,DropTol);
% Call preconditioned cg..
u = pcg(A,f,Tol,MaxIt,R',R);

function A = MakeSparseA(n1,n2,a,b,c)
% Set up A(n1,n2,a,b,c) in sparse format
In1 = speye(n1,n1);
In2 = speye(n2,n2);
En1 = spdiags(ones(n1,2),[-1,1],n1,n1);
En2 = spdiags(ones(n2,2),[-1,1],n2,n2);
A =  kron(In1,c*In2+b*En2) + kron(En1,b*In2+a*En2);