% This program constructs a stiffness matrix for a five-point-stencil
% on a simple rectangular lattice of size
%    nside1 x nside2
% 
% The figure illustrates the ordering of points in the grid
% for the case nside1=5 and nside2=8.
%
%        5-----10----15----20----25----30----35----40
%        |     |     |     |     |     |     |     |
%   n    4-----9-----14----19----24----29----34----39
%   s    |     |     |     |     |     |     |     |
%   i    3-----8-----13----18----23----28----33----38
%   d    |     |     |     |     |     |     |     |
%   e    2-----7-----12----17----22----27----32----37
%   1    |     |     |     |     |     |     |     |
%        1-----6-----11----16----21----26----31----36
%
%                        nside2
%
% To demonstrate how things work, the program solves a Dirichlet and a
% Neumann problem on the grid.
%
% Note that the function that constructs the stiffness matrix has several
% different calling options:
%
%   flag_matrix='orig'   The usual 5-point discrete Laplacian 4/-1/-1/-1/-1
%
%   flag_matrix='rcon'   A grid Laplacian where the conductivities of all
%                        bars are drawn at random.
%
%   flag_matrix='rcut'   A grid Laplacian for a grid where some bars have
%                        been cut at random.
%
%   flag_matrix='peri'   Conductivities vary periodically.
%
%   flag_matrix='crac'   Two long cracks in the grid.

function main_compressible

%LOCAL_orig
%LOCAL_rcon
LOCAL_helm

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_orig

rng('default')
rng(0)

nside       = 200;
nside1      = 2*nside+1;
nside2      = nside;
h           = 1/(nside-1);
N           = nside1*nside2;
flag_matrix = 'orig';

fprintf(1,'=== Domain size: %d x %d = %d\n',nside1,nside2,N)

%%% Construct the geometry.
[XX1,XX2] = meshgrid(linspace(-1-h,1+h,nside1),linspace(0,1,nside2));
xx        = [reshape(XX1,1,N);reshape(XX2,1,N)];

%%% Construct the stiffness matrix (this is a sparse matrix).
flag_anchor = 0;
A           = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor);

%%% Find the interior
I1 = find(xx(1,:) < -0.5*h);
I2 = find(0.5*h < xx(1,:));
I3 = find( (-0.5*h < xx(1,:)) & (xx(1,:) < 0.5*h) );
%plot(xx(1,I1),xx(2,I1),'b.',...
%     xx(1,I2),xx(2,I2),'g.',...
%     xx(1,I3),xx(2,I3),'r.')
%axis equal

tic
S3 = full(A(I3,I3) - ...
          A(I3,I1)*(A(I1,I1)\A(I1,I3)) - ...
          A(I3,I2)*(A(I2,I2)\A(I2,I3)));
fprintf(1,'Time for computing Schur complement = %f (sec)\n',toc)
ss = svd(S3);

nbox_max = 25;
ntot     = nside;
acc      = 1e-5;
NODES    = LOCAL_get_tree_and_lists(linspace(0,1,nside),nbox_max,0,1);
textshift = 3;

figure(1)
hold off
plot(ntot*[0,1,1,0,0],ntot*[0,0,1,1,0],'k')
hold on

for pbox = 1:size(NODES,2);
  if (NODES{5,pbox} == 2)
    ibox = NODES{4,pbox}(1);
    jbox = NODES{4,pbox}(2);
    indi = NODES{6,ibox} - 1 + (1:NODES{7,ibox});
    i1   = NODES{6,ibox};
    i2   = NODES{6,ibox} - 1 + NODES{7,ibox};
    indj = NODES{6,jbox} - 1 + (1:NODES{7,jbox});
    j1   = NODES{6,jbox};
    j2   = NODES{6,jbox} - 1 + NODES{7,jbox};
    B    = S3(indi,indj);
    k    = sum(svd(B) > (acc*ss(1)));
    plot([i1-1,i2,i2,i1-1,i1-1],...
         [j1-1,j1-1,j2,j2,j1-1],'k',...
         [j1-1,j1-1,j2,j2,j1-1],...
         [i1-1,i2,i2,i1-1,i1-1],'k')
    text(0.5*(i1+i2)-0.5-textshift,0.5*(j1+j2)-0.5,sprintf('%d',k),'Color','b')
    text(0.5*(j1+j2)-0.5-textshift,0.5*(i1+i2)-0.5,sprintf('%d',k),'Color','b')
  end
end
title(sprintf('Ranks structure of S_{3}. acc=%8.2e  nside=%d   type=%s',acc,nside,flag_matrix))
axis ij
axis equal
axis off
hold off
print('-depsc',sprintf('fig_S3_ranks_n%d_acc%d_%s.eps',ntot,round(-log10(acc)),flag_matrix))
export_fig(sprintf('fig_S3_ranks_n%d_acc%d_%s.pdf',ntot,round(-log10(acc)),flag_matrix),'-pdf','-transparent')
keyboard

figure(1)
acc = 1e-10;
hold off
plot(ntot*[0,1,1,0,0],ntot*[0,0,1,1,0],'k')
hold on

for pbox = 1:size(NODES,2);
  if (NODES{5,pbox} == 2)
    ibox = NODES{4,pbox}(1);
    jbox = NODES{4,pbox}(2);
    indi = NODES{6,ibox} - 1 + (1:NODES{7,ibox});
    i1   = NODES{6,ibox};
    i2   = NODES{6,ibox} - 1 + NODES{7,ibox};
    indj = NODES{6,jbox} - 1 + (1:NODES{7,jbox});
    j1   = NODES{6,jbox};
    j2   = NODES{6,jbox} - 1 + NODES{7,jbox};
    B    = S3(indi,indj);
    k    = sum(svd(B) > (acc*ss(1)));
    plot([i1-1,i2,i2,i1-1,i1-1],...
         [j1-1,j1-1,j2,j2,j1-1],'k',...
         [j1-1,j1-1,j2,j2,j1-1],...
         [i1-1,i2,i2,i1-1,i1-1],'k')
    text(0.5*(i1+i2)-0.5-textshift,0.5*(j1+j2)-0.5,sprintf('%d',k),'Color','b')
    text(0.5*(j1+j2)-0.5-textshift,0.5*(i1+i2)-0.5,sprintf('%d',k),'Color','b')
  end
end
title(sprintf('Ranks structure of S_{3}. acc=%8.2e  nside=%d   type=%s',acc,nside,flag_matrix))
axis ij
axis equal
axis off
hold off
print('-depsc',sprintf('fig_S3_ranks_n%d_acc%d_%s.eps',ntot,round(-log10(acc)),flag_matrix))
export_fig(sprintf('fig_S3_ranks_n%d_acc%d_%s.pdf',ntot,round(-log10(acc)),flag_matrix),'-pdf','-transparent')

keyboard
        
figure(1)
hold off
nhalf = round(nside/2);
semilogy((1/ss(1))*svd(S3(1:nhalf,(nhalf+1):end)),'x-')
title('Singular values of the (2,3) block of S_3')
ylabel('ss(j)/||S_3||')
xlabel('j')
print('-depsc',sprintf('fig_S3_svd_n%d_%s.eps',ntot,flag_matrix))
export_fig(sprintf('fig_S3_svd_n%d_%s.pdf',ntot,flag_matrix),'-pdf','-transparent')

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_rcon

rng('default')
rng(0)

nside       = 200;
nside1      = 2*nside+1;
nside2      = nside;
h           = 1/(nside-1);
N           = nside1*nside2;
flag_matrix = 'rcom';

fprintf(1,'=== Domain size: %d x %d = %d\n',nside1,nside2,N)

%%% Construct the geometry.
[XX1,XX2] = meshgrid(linspace(-1-h,1+h,nside1),linspace(0,1,nside2));
xx        = [reshape(XX1,1,N);reshape(XX2,1,N)];

%%% Construct the stiffness matrix (this is a sparse matrix).
flag_anchor = 0;
A           = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor);

%%% Find the interior
I1 = find(xx(1,:) < -0.5*h);
I2 = find(0.5*h < xx(1,:));
I3 = find( (-0.5*h < xx(1,:)) & (xx(1,:) < 0.5*h) );
%plot(xx(1,I1),xx(2,I1),'b.',...
%     xx(1,I2),xx(2,I2),'g.',...
%     xx(1,I3),xx(2,I3),'r.')
%axis equal

tic
S3 = full(A(I3,I3) - ...
          A(I3,I1)*(A(I1,I1)\A(I1,I3)) - ...
          A(I3,I2)*(A(I2,I2)\A(I2,I3)));
fprintf(1,'Time for computing Schur complement = %f (sec)\n',toc)
ss = svd(S3);

nbox_max = 25;
ntot     = nside;
acc      = 1e-10;
NODES    = LOCAL_get_tree_and_lists(linspace(0,1,nside),nbox_max,0,1);
textshift = 3;

figure(1)
hold off
plot(ntot*[0,1,1,0,0],ntot*[0,0,1,1,0],'k')
hold on

for pbox = 1:size(NODES,2);
  if (NODES{5,pbox} == 2)
    ibox = NODES{4,pbox}(1);
    jbox = NODES{4,pbox}(2);
    indi = NODES{6,ibox} - 1 + (1:NODES{7,ibox});
    i1   = NODES{6,ibox};
    i2   = NODES{6,ibox} - 1 + NODES{7,ibox};
    indj = NODES{6,jbox} - 1 + (1:NODES{7,jbox});
    j1   = NODES{6,jbox};
    j2   = NODES{6,jbox} - 1 + NODES{7,jbox};
    B    = S3(indi,indj);
    k    = sum(svd(B) > (acc*ss(1)));
    plot([i1-1,i2,i2,i1-1,i1-1],...
         [j1-1,j1-1,j2,j2,j1-1],'k',...
         [j1-1,j1-1,j2,j2,j1-1],...
         [i1-1,i2,i2,i1-1,i1-1],'k')
    text(0.5*(i1+i2)-0.5-textshift,0.5*(j1+j2)-0.5,sprintf('%d',k),'Color','b')
    text(0.5*(j1+j2)-0.5-textshift,0.5*(i1+i2)-0.5,sprintf('%d',k),'Color','b')
  end
end
title(sprintf('Ranks structure of S_{3}. acc=%8.2e  nside=%d   type=%s',acc,nside,flag_matrix))
axis ij
axis equal
axis off
hold off
print('-depsc',sprintf('fig_S3_ranks_n%d_acc%d_%s.eps',ntot,round(-log10(acc)),flag_matrix))
export_fig(sprintf('fig_S3_ranks_n%d_acc%d_%s.pdf',ntot,round(-log10(acc)),flag_matrix),'-pdf','-transparent')

figure(2)
hold off
nhalf = round(nside/2);
semilogy((1/ss(1))*svd(S3(1:nhalf,(nhalf+1):end)),'x-')
title('Singular values of the (2,3) block of S_3')
ylabel('ss(j)/||S_3||')
xlabel('j')
print('-depsc',sprintf('fig_S3_svd_n%d_%s.eps',ntot,flag_matrix))
export_fig(sprintf('fig_S3_svd_n%d_%s.pdf',ntot,flag_matrix),'-pdf','-transparent')
keyboard

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_helm

rng('default')
rng(0)

nside       = 40;
nside1      = 2*nside+1;
nside2      = nside;
kh          = 10;
lambda      = 2*pi/kh;
h           = 1/(nside-1);
N           = nside1*nside2;
flag_matrix = 'orig';

fprintf(1,'=== Domain size: %d x %d = %d\n',nside1,nside2,N)
fprintf(1,'    HELMHOLTZ problem. Divider = %0.2f lambda (%0.1f ppw)\n',...
        1/lambda,lambda*nside)

%%% Construct the geometry.
[XX1,XX2] = meshgrid(linspace(-1-h,1+h,nside1),linspace(0,1,nside2));
xx        = [reshape(XX1,1,N);reshape(XX2,1,N)];

%%% Construct the stiffness matrix (this is a sparse matrix).
flag_anchor = 0;
A           = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor);
keyboard

%%% Find the interior
I1 = find(xx(1,:) < -0.5*h);
I2 = find(0.5*h < xx(1,:));
I3 = find( (-0.5*h < xx(1,:)) & (xx(1,:) < 0.5*h) );
%plot(xx(1,I1),xx(2,I1),'b.',...
%     xx(1,I2),xx(2,I2),'g.',...
%     xx(1,I3),xx(2,I3),'r.')
%axis equal

tic
S3 = full(A(I3,I3) - ...
          A(I3,I1)*(A(I1,I1)\A(I1,I3)) - ...
          A(I3,I2)*(A(I2,I2)\A(I2,I3)));
fprintf(1,'Time for computing Schur complement = %f (sec)\n',toc)
ss = svd(S3);

nbox_max = 25;
ntot     = nside;
acc      = 1e-5;
NODES    = LOCAL_get_tree_and_lists(linspace(0,1,nside),nbox_max,0,1);
textshift = 3;

figure(1)
hold off
plot(ntot*[0,1,1,0,0],ntot*[0,0,1,1,0],'k')
hold on

for pbox = 1:size(NODES,2);
  if (NODES{5,pbox} == 2)
    ibox = NODES{4,pbox}(1);
    jbox = NODES{4,pbox}(2);
    indi = NODES{6,ibox} - 1 + (1:NODES{7,ibox});
    i1   = NODES{6,ibox};
    i2   = NODES{6,ibox} - 1 + NODES{7,ibox};
    indj = NODES{6,jbox} - 1 + (1:NODES{7,jbox});
    j1   = NODES{6,jbox};
    j2   = NODES{6,jbox} - 1 + NODES{7,jbox};
    B    = S3(indi,indj);
    k    = sum(svd(B) > (acc*ss(1)));
    plot([i1-1,i2,i2,i1-1,i1-1],...
         [j1-1,j1-1,j2,j2,j1-1],'k',...
         [j1-1,j1-1,j2,j2,j1-1],...
         [i1-1,i2,i2,i1-1,i1-1],'k')
    text(0.5*(i1+i2)-0.5-textshift,0.5*(j1+j2)-0.5,sprintf('%d',k),'Color','b')
    text(0.5*(j1+j2)-0.5-textshift,0.5*(i1+i2)-0.5,sprintf('%d',k),'Color','b')
  end
end
title(sprintf('Ranks structure of S_{3}. acc=%8.2e  nside=%d   type=%s',acc,nside,flag_matrix))
axis ij
axis equal
axis off
hold off
print('-depsc',sprintf('fig_S3_ranks_n%d_acc%d_%s.eps',ntot,round(-log10(acc)),flag_matrix))
export_fig(sprintf('fig_S3_ranks_n%d_acc%d_%s.pdf',ntot,round(-log10(acc)),flag_matrix),'-pdf','-transparent')
keyboard

figure(1)
acc = 1e-10;
hold off
plot(ntot*[0,1,1,0,0],ntot*[0,0,1,1,0],'k')
hold on

for pbox = 1:size(NODES,2);
  if (NODES{5,pbox} == 2)
    ibox = NODES{4,pbox}(1);
    jbox = NODES{4,pbox}(2);
    indi = NODES{6,ibox} - 1 + (1:NODES{7,ibox});
    i1   = NODES{6,ibox};
    i2   = NODES{6,ibox} - 1 + NODES{7,ibox};
    indj = NODES{6,jbox} - 1 + (1:NODES{7,jbox});
    j1   = NODES{6,jbox};
    j2   = NODES{6,jbox} - 1 + NODES{7,jbox};
    B    = S3(indi,indj);
    k    = sum(svd(B) > (acc*ss(1)));
    plot([i1-1,i2,i2,i1-1,i1-1],...
         [j1-1,j1-1,j2,j2,j1-1],'k',...
         [j1-1,j1-1,j2,j2,j1-1],...
         [i1-1,i2,i2,i1-1,i1-1],'k')
    text(0.5*(i1+i2)-0.5-textshift,0.5*(j1+j2)-0.5,sprintf('%d',k),'Color','b')
    text(0.5*(j1+j2)-0.5-textshift,0.5*(i1+i2)-0.5,sprintf('%d',k),'Color','b')
  end
end
title(sprintf('Ranks structure of S_{3}. acc=%8.2e  nside=%d   type=%s',acc,nside,flag_matrix))
axis ij
axis equal
axis off
hold off
print('-depsc',sprintf('fig_S3_ranks_n%d_acc%d_%s.eps',ntot,round(-log10(acc)),flag_matrix))
export_fig(sprintf('fig_S3_ranks_n%d_acc%d_%s.pdf',ntot,round(-log10(acc)),flag_matrix),'-pdf','-transparent')

keyboard
        
figure(1)
hold off
nhalf = round(nside/2);
semilogy((1/ss(1))*svd(S3(1:nhalf,(nhalf+1):end)),'x-')
title('Singular values of the (2,3) block of S_3')
ylabel('ss(j)/||S_3||')
xlabel('j')
print('-depsc',sprintf('fig_S3_svd_n%d_%s.eps',ntot,flag_matrix))
export_fig(sprintf('fig_S3_svd_n%d_%s.pdf',ntot,flag_matrix),'-pdf','-transparent')

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This is the key function in this program. 
% For documention on the call see the top of the file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function A = LOCAL_get_A(nside1,nside2,flag_matrix,flag_anchor)

if (nargin < 4)
  flag_anchor = 0;
end

if (nargin == 1)
  nside2 = nside1;
end

BARS  = LOCAL_get_links(nside1,nside2);
nbars = size(BARS,2);

if strcmp(flag_matrix,'orig')
    
  BARS(3,:) = ones(1,nbars);
  
elseif strcmp(flag_matrix,'rcon')
    
  BARS(3,:) = 1 + rand(1,nbars);
  
elseif strcmp(flag_matrix,'rcom')
    
  BARS(3,:) = 10.^(-10*rand(1,nbars));
  
elseif strcmp(flag_matrix,'rcut')
    
  p = 0.04;
  BARS(3,:) = (rand(1,nbars) > p);
  
elseif strcmp(flag_matrix,'peri')
  
  [XX1,XX2] = meshgrid(0:(nside1-1),0:(nside2-1));
  xx        = [reshape(XX1,1,nside1*nside2);reshape(XX2,1,nside1*nside2)];
  h         = 1/(max([nside1,nside2])-1);
  Ncell     = round(min(nside1,nside2)/25);
  XX1sc     = h*xx(1,BARS(1,:));
  XX2sc     = h*xx(2,BARS(1,:));
  YY1sc     = h*xx(1,BARS(2,:));
  YY2sc     = h*xx(2,BARS(2,:));
  ZZ1sc     = 0.5*(XX1sc+YY1sc);
  ZZ2sc     = 0.5*(XX2sc+YY2sc);
  BARS(3,:) = 1 - 0.9*((cos(Ncell*pi*ZZ1sc).*cos(Ncell*pi*ZZ2sc)).^2);
    
elseif strcmp(flag_matrix,'crac')

  % We create a list of bars to be removed:
  n1a = round(0.4*nside1);
  n1b = round(0.6*nside1);
  n2a = round(0.4*nside2);
  n2b = round(0.7*nside2);
  gap = 4;
  
  ii1 = n1a*nside2 + (1:n2a);
  jj1 = ii1 + nside2;
  ii2 = n1a*nside2 + n2a + 1 + nside2*(1:(n1b-n1a));
  jj2 = ii2 - 1;
  ii3 = n1b*nside2 + 1 + (n2a:n2b);
  jj3 = ii3 + nside2;
  ii4 = (n1b-gap)*nside2 + (nside2:(-1):(n2a+gap+2));
  jj4 = ii4-nside2;
  ii5 = (n1a-1)*nside2 + (n2a+gap+1) + nside2*(1:(n1b-n1a-gap));
  jj5 = ii5+1;
  
  ii = [ii1,ii2,ii3,ii4,ii5];
  jj = [jj1,jj2,jj3,jj4,jj5];
  
  ni   = length(ii);
  BARS(3,:) = ones(1,nbars);
  BARS = [BARS,[ii;jj;-ones(1,ni)]];
      
end

if (flag_anchor == 0)
    
  A = sparse([BARS(1,:), BARS(1,:), BARS(2,:), BARS(2,:)],...
             [BARS(1,:), BARS(2,:), BARS(1,:), BARS(2,:)],...
             [BARS(3,:),-BARS(3,:),-BARS(3,:), BARS(3,:)]);
         
else
    
  ii_rim = [1:nside1,...
            (nside1+1):nside1:(nside1*(nside2-2)+1),...
            (2*nside1):nside1:(nside1*(nside2-1)),...
             nside1*(nside2-1) + (1:nside1)];
     
  A = sparse([BARS(1,:), BARS(1,:), BARS(2,:), BARS(2,:),ii_rim],...
             [BARS(1,:), BARS(2,:), BARS(1,:), BARS(2,:),ii_rim],...
             [BARS(3,:),-BARS(3,:),-BARS(3,:), BARS(3,:),ones(1,2*(nside1+nside2-2))]);
    
end

if strcmp(flag_matrix,'crac')
  [ii,jj,aa] = find(A);
  ind = find(abs(aa) < 1e-12);
  aa(ind) = zeros(size(ind));
  A = sparse(ii,jj,aa);
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This supplementary functions returns a list of all links in a
% rectangular grid of size n1 x n2.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function BARS = LOCAL_get_links(n1,n2)

if (nargin == 1)
  n2 = n1;
end

BARS  = zeros(3,2*n1*n2 - n1 - n2);

ndone = 0;
for i1 = 1:(n1-1)
  for i2 = 1:(n2-1)
    j_c = n2*(i1-1) + i2;
    j_n = n2*(i1-1) + i2 + 1;
    j_e = n2*i1 + i2;
    BARS(1:2,ndone + (1:2)) = [j_c, j_c; j_n, j_e];
    ndone = ndone + 2;
  end
  j_c = n2*i1;
  j_e = n2*(i1+1);
  BARS(1:2,ndone + 1) = [j_c; j_e];
  ndone = ndone + 1;
end
for i2 = 1:(n2-1)
  j_c = n2*(n1-1) + i2;
  j_n = n2*(n1-1) + i2 + 1;
  BARS(1:2,ndone + 1) = [j_c; j_n];
  ndone = ndone + 1;
end

BARS(3,:) = ones(1,ndone);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function draws the grid corresponding to the sparse matrix A.
% Node "i" has coordinates xx(:,i) in the plane.
% Two nodes i and j are connected iff A(i,j) is non-zero.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function LOCAL_draw_lattice(xx,A)

[ii,jj,aa] = find(A);

hold off
plot([xx(1,ii);xx(1,jj)],[xx(2,ii);xx(2,jj)],'k')
hold on
plot(xx(1,:),xx(2,:),'b.')
hold off

bord = 2;
axis([min(xx(1,:))-bord,...
      max(xx(1,:))+bord,...
      min(xx(2,:))-bord,...
      max(xx(2,:))+bord])
axis equal

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function constructs the complement of the index vector 
% ind1 in the set [1,2,3,...,ntot].
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ind2 = LOCAL_index_complement(ind1,ntot)

ind       = 1:ntot;
ind(ind1) = 2*ntot*ones(1,length(ind1));
ind       = sort(ind);
ind2      = ind(1:(ntot-length(ind1)));

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function NODES = LOCAL_get_tree_and_lists(xx,nbox_max,x_left,x_right)

ntot = length(xx);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STEP 1: Set up a tree structure.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

TREE = cell(7,ntot);

% Construct the top node.
TREE{1,1} = [x_left,x_right]; % Set the geometry.
TREE{2,1} = 0;                % Set the level.
TREE{3,1} = NaN;              % Set the father.
TREE{4,1} = [];               % Set the list of children.
TREE{5,1} = 0;                % Set the number of children.
TREE{6,1} = 1;                % Set the index of the first node.
TREE{7,1} = ntot;             % Set the index of the last node.  
ibox_last = 0;
ncreated  = 1;
% Create smaller nodes via hierarchical subdivision.
% We sweep one level at a time.
while (ncreated > 0)
  ibox_first = ibox_last + 1;
  ibox_last  = ibox_last + ncreated;
  ncreated   = 0;
  for ibox = ibox_first:ibox_last
    nbox = TREE{7,ibox};
    if (nbox > nbox_max)
      xleft              = TREE{1,ibox}(1);
      xright             = TREE{1,ibox}(2);
      xmid               = 0.5*(xleft + xright);
      ind                = TREE{6,ibox} - 1 + (1:TREE{7,ibox});
      nhalf              = sum(xx(ind) <= xmid);
      if (nhalf > 0) % Create the left child.
        ison         = ibox_last + ncreated + 1;
        TREE{4,ibox} = [TREE{4,ibox},ison]; % Update the child list for ibox.
        TREE{5,ibox} = TREE{5,ibox} + 1;    % Update the child list for ibox.
        TREE{1,ison} = [xleft,xmid];
        TREE{2,ison} = TREE{2,ibox} + 1;
        TREE{3,ison} = ibox;
        TREE{4,ison} = [];
        TREE{5,ison} = 0;
        TREE{6,ison} = TREE{6,ibox};
        TREE{7,ison} = nhalf;
        ncreated     = ncreated + 1;
      end
      if (nhalf < TREE{7,ibox})
        ison         = ibox_last + ncreated + 1;
        TREE{4,ibox} = [TREE{4,ibox},ison]; % Update the child list for ibox.
        TREE{5,ibox} = TREE{5,ibox} + 1;    % Update the child list for ibox.
        TREE{1,ison} = [xmid,xright];
        TREE{2,ison} = TREE{2,ibox} + 1;
        TREE{3,ison} = ibox;
        TREE{4,ison} = [];
        TREE{5,ison} = 0;
        TREE{6,ison} = TREE{6,ibox} + nhalf;
        TREE{7,ison} = TREE{7,ibox} - nhalf;
        ncreated     = ncreated + 1;
      end
    end
  end
end

% This whole section is devoted to the simple task of
% decreasing the size of TREE.
% The operation TREE = TREE{:,1:nboxes} apparently
% does not work for "cell" structures.
nboxes = ibox_last;
NODES = cell(50,nboxes);
for ibox = 1:nboxes
  for j = 1:size(TREE,1)
    NODES{j,ibox} = TREE{j,ibox};
  end
end
clear TREE

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 2: Create the interaction lists.
%         These require the construction of two auxiliary lists:
%           leaf_list = a single vector hold all leaves, sorted.
%           nei_lists{ibox} = all neighbors of ibox
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Construct "leaf_list".
leaf_list  = zeros(1,nboxes);
leftpoints = zeros(1,nboxes);
ileaf      = 0;
for ibox = 1:nboxes;
  if (NODES{5,ibox}==0)
    ileaf             = ileaf+1;
    leaf_list(ileaf)  = ibox;
    leftpoints(ileaf) = NODES{1,ibox}(1);
  end
end
nleaves     = ileaf;
[trash,ind] = sort(leftpoints(1:nleaves));
leaf_list   = leaf_list(ind);

% Construct "nei_lists".
nei_lists = cell(1,nboxes);
nei_lists{1} = []; % The root has no neighbors;
for ibox = 2:nboxes
  nei_lists{ibox} = [];
  ifath = NODES{3,ibox};
  for kbox = sort([ifath,nei_lists{ifath}]) % Loop over the father and his neighbors.
    for jbox = NODES{4,kbox} % jbox is a POTENTIAL neighbor of ibox
      if LOCAL_isneighbor(NODES,ibox,jbox)
        nei_lists{ibox} = [nei_lists{ibox},jbox];
      end
    end
  end
end

% Create lists 1 and 5
for j = 1:(nleaves-1)
  ibox_left  = leaf_list(j  );
  ibox_right = leaf_list(j+1);
  if (NODES{1,ibox_left}(2) == NODES{1,ibox_right}(1))
    % ibox_left and ibox_right are indeed neighbors.
    if (NODES{2,ibox_left} > NODES{2,ibox_right})
      % ibox_right is on a COARSER level.
      NODES{18,ibox_left } = [NODES{18,ibox_left },ibox_right];
      NODES{10,ibox_right} = [NODES{10,ibox_right},ibox_left ];
    elseif (NODES{2,ibox_left} < NODES{2,ibox_right})
      % ibox_right is on a FINER level.
      NODES{10,ibox_left } = [NODES{10,ibox_left },ibox_right];
      NODES{18,ibox_right} = [NODES{18,ibox_right},ibox_left ];
    else
      % ibox_right is on the SAME level.
      NODES{10,ibox_left } = [NODES{10,ibox_left },ibox_right];
      NODES{10,ibox_right} = [NODES{10,ibox_right},ibox_left ];
    end
  end
end

% Create list 2
for ibox = 2:nboxes
  NODES{12,ibox} = [];
  ifath = NODES{3,ibox};
  for ifathnei = nei_lists{ifath}
    for jbox = NODES{4,ifathnei}
      if LOCAL_iswellsep(NODES,ibox,jbox)
        NODES{12,ibox} = [NODES{12,ibox},jbox];
      end
    end
  end
end

% Create lists 3 and 4
% Specifically, for every leaf "ibox", find list3(ibox).
% Creation of list4 is piggybacked as it is the dual of list3.
for ibox = 1:nboxes
  NODES{14,ibox} = [];
  NODES{16,ibox} = [];
end
for ibox = leaf_list
  boxes_to_check = zeros(1,50);
  ntocheck       = 0;
  for jbox = nei_lists{ibox}
    boxes_to_check(ntocheck + (1:NODES{5,jbox})) = NODES{4,jbox};
    ntocheck = ntocheck + NODES{5,jbox};
  end
  while (ntocheck > 0)
    jbox = boxes_to_check(1);
    if LOCAL_iswellsep(NODES,jbox,ibox)
      NODES{14,ibox} = [NODES{14,ibox},jbox];
      NODES{16,jbox} = [NODES{16,jbox},ibox];
    else
      boxes_to_check(ntocheck + (1:NODES{5,jbox})) = NODES{4,jbox};
      ntocheck = ntocheck + NODES{5,jbox};
    end
    boxes_to_check(1) = boxes_to_check(ntocheck);
    ntocheck = ntocheck-1;
  end
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function output = LOCAL_isneighbor(NODES,ibox,jbox)

if (NODES{2,ibox} == NODES{2,jbox}) % ibox and jbox are on the same level
  if (NODES{1,ibox}(2) == NODES{1,jbox}(1)) 
    output = 1;
  elseif (NODES{1,ibox}(1) == NODES{1,jbox}(2)) 
    output = 1;
  else
    output = 0;
  end
else % ibox and jbox are not on the same level
  output = 0;
end
    
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function checks if jbox is well-separated from ibox

function output = LOCAL_iswellsep(NODES,ibox,jbox)

xleft  = NODES{1,ibox}(1);
xright = NODES{1,ibox}(2);
width  = xright - xleft;
xmid   = 0.5*(xleft+xright);
dist   = min(abs(NODES{1,jbox} - xmid));
if (dist >= (1.499999999999*width))
  output = 1;
else
  output = 0;
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

