function L=ic3_3d(A,n);
% Makes an incomplete Cholesky factorization of A,
% A is the 3D discrete Laplacian matrix on a cube,
% three extra subdiagonal filled-in in L
nnn=n*n*n;;nn=n*n;nnnmn=nnn-n+1;nnnmnn=nnn-nn+1;
bb=zeros(nnn,1);cc=zeros(nnn,1);aa=zeros(nnn,1);dd=zeros(nnn,1);
ee=zeros(nnn,1);ff=zeros(nnn,1);gg=zeros(nnn,1);
for i=1:nnn,
 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
 if i>nn as=as-dd(i-nn)*dd(i-nn);end
 if i>n-1 as =as-ee(i-n+1)*ee(i-n+1);end
 if i>nn-n as=as-ff(i-nn+n)*ff(i-nn+n);end
 if i>nn-1 as=as-gg(i-nn+1)*gg(i-nn+1);end
 aa(i)=sqrt(as);
 if i<nnn, bs=A(i,i+1);
   if i>n-1 bs=bs-ee(i-n+1)*cc(i-n+1);end
   if i>nn-1 bs=bs-gg(i-nn+1)*dd(i-nn+1);end
 bb(i)=bs/aa(i); end
 if i<nnn-n+2 es=A(i,i+n-1);
   if i>1 es=es-bb(i-1)*cc(i-1);end
   if i>nn-n es=es-ff(i-nn+n)*gg(i-nn+n);end
 ee(i)=es/aa(i);end
 if i<nnnmn, cs=A(i,i+n); 
   if i>nn-n cs=cs-dd(i-nn+n)*ff(i-nn+n);end
 cc(i)=cs/aa(i); end
 if i<nnn-nn+n+1 fs=A(i,i+nn-n);
   if i>n-1 fs=fs-ee(i-n+1)*gg(i-n+1);end
   if i>n fs=fs-dd(i-n)*cc(i-n);end
 ff(i)=fs/aa(i);end
 if i<nnnmnn+1 gs=A(i,i+nn-1);
   if i>1 gs=gs-bb(i-1)*dd(i-1);end
 gg(i)=gs/aa(i);end
 if i<nnnmnn, dd(i)=A(i,i+nn)/aa(i);end
end
L=spdiags([dd gg ff cc ee bb aa], [-nn,-nn+1,-nn+n,-n,-n+1,-1,0],nnn,nnn);