function conjgrad(L,dim)
% L is the number of discretization points in each direction (1,2,3)
% dim is the dimension (1,2,3)
I = speye(L,L);
E = sparse(2:L,1:L-1,1,L,L);
T = 2*I-E-E';
if (dim==1)
 b=ones([L,1]);
 A=T;
 cpu0 = cputime;
 pcg(A,b,1.e-10,5000);
 cpu1 = cputime;
elseif (dim==2)
 A = kron(T,I)+kron(I,T);
 b=ones([L*L,1]);
 cpu0 = cputime;
 pcg(A,b,1.e-6,5000);
 cpu1 = cputime;
elseif (dim==3)
 A = kron(T,kron(I,I))+kron(I,kron(I,T))+kron(I,kron(T,I));
 b=ones([L*L*L,1]);
 cpu0 = cputime;
 pcg(A,b,1.e-6,5000);
 cpu1 = cputime;
end
fprintf(" cpu time pcg %e  \n",cpu1-cpu0);
spy(A);
R=chol(A);
cpu2=cputime;
fprintf(" cpu time choleski %e \n",cpu2-cpu1);
whos
end