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