% Show the red-black ordering sparsity profile of the matrix 
% corresponding to the Lapace operator on a square 
% discretized with the five points stencil and N odd!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (mod(N,2) ==0)
   disp('N should be odd')
   exit;
end
% Generate a grid with NxN internal nodes
N2=N*N; NNS2=N2/2+1;
G=numgrid('S',N+2);
if (N <= 9) 
  G 
end 
% plot the grid
spy(G)
title('Finite difference grid')
pause;
% Generate the matrix and plot it's sparsity profile
A=delsq(G); 
D=zeros(N2);
for i=1:N2
   for j=1:N2
     D(i,j)=A(i,j);
   end;
end;
if (N <= 5) 
D
end 
spy(D)
title('Sparsity of A with column numbering')
pause
for i=1:fix(N2/2)
   p(2*i-1)=i; p(2*i)=fix(NNS2)+i;
end;
p(N2)=fix(N2/2)+1;
GRB=zeros(N+2);
for i=2:N+1
   for j=2:N+1
     GRB(i,j)=p(G(i,j));
   end;
end;
if (N <= 7) 
  GRB 
end
pause
% Sparsity profile of the RB matrix
A=delsq(GRB);
DRB=zeros(N2);
for i=1:N2
   for j=1:N2
     DRB(i,j)=A(i,j);
   end;
end;
if (N <= 5) 
DRB
end 
spy(DRB) 
title('Sparsity of A with RB numbering')