clf
pos=get(gcf,'pos');
pos(3)=pos(4)*2;
set(gcf,'pos',pos);

pausetime=0.01;

for method=1:5
  N=10;
  A=sparse(toeplitz([2,-1,zeros(1,N-2)]));
  I=speye(N,N);
  B=kron(A,I)+kron(I,A);
  [x,y]=ndgrid(1:N,1:N);
  B0=B;
  L0=chol(B0)';
  
  switch method
   case 1
    p=1:N^2;
    methodstr='Natural Ordering';
   case 2
    p=symrcm(B);
    methodstr='Reverse Cuthill-McKee';
   case 3
    p=realmmd(B);
    methodstr='Minimum Degree';
   case 4
    p=symamd(B);
    methodstr='SYMAMD (MATLAB Default)';
   case 5
    % Precomputed (for 10-by-10 only)
    p=[13 12 11 1 2 31 41 42 32 21 22 4 15 25 24 35 34 44 45 3 14 23 33 43 61 62 72 73 81 91 92 71 82 64 65 75 85 95 94 84 63 74 83 93 51 52 53 54 55 6 7 18 27 37 48 17 28 38 9 10 20 30 40 50 8 19 29 39 49 57 67 78 87 97 98 77 88 69 70 80 90 100 68 79 89 99 47 58 59 60 5 16 26 36 46 56 66 76 86 96];
    methodstr='Nested Dissection';
  end

  B=B(p,p);
  R=chol(B);
  xy=[x(p);y(p)]';
  Bplt=B;
  
  nflops=0;
  n=size(B,1);
  L=zeros(n);
  
  labelpars={'fontname','helvetica','fontsize',16,'horiz','center'};
  flashpars={'fontname','helvetica','fontsize',25,'horiz','center'};
  titlepars={'fontname','helvetica','fontsize',32,'horiz','center'};
  lpars={'marker','.','linestyle','none','markersize',64};
  tpars={'fontname','helvetica','fontsize',16,'horiz','center','col','g'};
  for i=1:n
    subplot(1,2,2)
    spy(Bplt,16), xlabel('');
    nzhandle=text(50,114.9,['nz = ',int2str(nnz(Bplt))],flashpars{:});

    subplot(1,2,1)
    cla
    gplot(B,xy);
    line(xy(i:N^2,1),xy(i:N^2,2),'color','k',lpars{:});
    for j=i:n
      degree=length(find(B(:,j)))-1;
      text(xy(j,1),xy(j,2),int2str(degree),tpars{:});
    end
    axis equal,axis off,axis([1,N,1,N])
    
    flophandle=text(5.5,-.25,['flops = ',int2str(nflops)],flashpars{:});
    titlehandle=text(11.5,11.2,methodstr,titlepars{:});
    
    pause(pausetime);
    if i==1
      pause
    end
    
    subplot(1,2,2)
    Btmp=sparse(N^2,N^2);
    Btmp(i,i:N^2)=Bplt(i,i:N^2);
    Btmp(i:N^2,i)=Bplt(i:N^2,i);
    hold on, spy(Btmp,'r',25), xlabel(''); hold off
    delete(nzhandle);
    if i==n
      nzhandle=text(50,114.9,['nz = ',int2str(nnz(Bplt))],flashpars{:}); 
    else
      nzhandle=text(50,113.25,['nz = ',int2str(nnz(Bplt))],labelpars{:});
    end
    
    subplot(1,2,1)
    line(xy(i,1),xy(i,2),'color','r',lpars{:});
    delete(flophandle);
    if i==n
      flophandle=text(5.5,-.25,['flops = ',int2str(nflops)],flashpars{:});
    else
      flophandle=text(5.5,-.1,['flops = ',int2str(nflops)],labelpars{:});
    end
    
    pause(pausetime);
    
    degree=length(find(B(:,i)))-1;
    nflops=nflops+2*degree^2;
    
    L(i,i)=sqrt(B(i,i));
    L(i+1:n,i)=B(i+1:n,i)/L(i,i);
    B(i+1:n,i+1:n)=B(i+1:n,i+1:n)-L(i+1:n,i)*L(i+1:n,i)';
    Bplt(i+1:n,i+1:n)=Bplt(i+1:n,i+1:n)-L(i+1:n,i)*L(i+1:n,i)';
    B(i, i:n)=0;
    B(i:n, i)=0;
    Bplt(i+1:n,i)=0;
  end
  pause
end