% illustrate nested dissection on n x n mesh
%  input m where n=2^m-1,
n = 2^m-1; N = n*n;
T = 2*eye(n) + diag(-ones(n-1,1),1) + diag(-ones(n-1,1),-1);
A = kron(T,eye(n)) + kron(eye(n),T);
LA = (chol(A))'; 
f=1;
figure(f), clf, spy(LA,'r'), hold on, spy(tril(A,0),'k'), spy(triu(A,1),'b')
title(['Natural order: N=',int2str(N),', A black and blue, L black and red, nnz(L)=',int2str(nnz(LA))])
pause
% Do nested dissection once
p1 =[];
for i=1:n, p1=[p1,[(i-1)*n+1:(i-1)*n+(n-1)/2]]; end
for i=1:n, p1=[p1,[(i-1)*n+(n+3)/2:i*n]]; end
for i=1:n, p1=[p1,(i-1)*n+(n+1)/2]; end
% norm(sort(p1)-(1:N))
Ap1=A(p1,p1);LAp1=(chol(Ap1))';
f=f+1;
figure(f), clf, spy(LAp1,'r'), hold on, spy(tril(Ap1,0),'k'), spy(triu(Ap1,1),'b')
title(['Nested Dissection 1x: N=',int2str(N),', A black and blue, L black and red, nnz(L)=',int2str(nnz(LAp1))])
pause
% Do nested dissection twice
nm1o2 = (n-1)/2;
p12 = [1:nm1o2^2,nm1o2^2+nm1o2+1:nm1o2*n,nm1o2^2+1:nm1o2^2+nm1o2];
p12 = [p12,p12+nm1o2*n,(N-n+1:N)];
% norm(sort(p12)-(1:N))
p2 = p1(p12);
% norm(sort(p2)-(1:N))
Ap2=A(p2,p2);LAp2=(chol(Ap2))';
f=f+1;
figure(f), clf, spy(LAp2,'r'), hold on, spy(tril(Ap2,0),'k'), spy(triu(Ap2,1),'b')
title(['Nested Dissection 2x: N=',int2str(N),', A black and blue, L black and red, nnz(L)=',int2str(nnz(LAp2))])
pause
% Do nested dissection 3x
p23 = [];
nn=(n-1)/2;
for i=1:nn, p23=[p23,[(i-1)*nn+1:(i-1)*nn+(nn-1)/2]]; end
for i=1:nn, p23=[p23,[(i-1)*nn+(nn+3)/2:i*nn]]; end
for i=1:nn, p23=[p23,(i-1)*nn+(nn+1)/2]; end
p3=p2;
p3(1:nn^2)=p3(p23);
p3(1+nn^2:nn^2+nn^2)=p3(p23+nn^2);
p3(1+nn*n:nn*n+nn^2)=p3(p23+nn*n);
p3(1+nn*n+nn^2:nn*n+nn^2+nn^2)=p3(p23+nn*n+nn^2);
% norm(sort(p3)-(1:N))
Ap3=A(p3,p3);LAp3=(chol(Ap3))';
f=f+1;
figure(f), clf, spy(LAp3,'r'), hold on, spy(tril(Ap3,0),'k'), spy(triu(Ap3,1),'b')
title(['Nested Dissection 3x: N=',int2str(N),', A black and blue, L black and red, nnz(L)=',int2str(nnz(LAp3))])
pause
% complete nested dissection
ndperm = nd(1:N);
And = A(ndperm,ndperm);LAnd=(chol(And))';
f=f+1;
figure(f), clf, spy(LAnd,'r'), hold on, spy(tril(And,0),'k'), spy(triu(And,1),'b')
title(['Full Nested Dissection: N=',int2str(N),', A black and blue, L black and red, nnz(L)=',int2str(nnz(LAnd))])

