function L=ic0_2d(A,n);
% Makes an incomplete Cholesky factorization of A,
% A is the 2D discrete Laplacian matrix on a square,
% no extra subdiagonal filled-in in L
nm1=n-1;nn=n*n;nnmn=nn-n+1;
bb=zeros(nn,1);cc=zeros(nn,1);aa=zeros(nn,1);
for i=1:nn,
 as=A(i,i);
 if i>1, as=as-bb(i-1)*bb(i-1); end
 if i>n, as=as-cc(i-n)*cc(i-n); end
 aa(i)=sqrt(as);
 if i<nn, bb(i)=A(i,i+1)/aa(i); end
 if i<nnmn, cc(i)=A(i,i+n)/aa(i); end
end
L=spdiags([cc bb aa], [-n,-1,0],nn, nn);