% cf_split.m for C/F splitting for an AMG
%
% Usage:  [C F] = cf_split(A) % where A is an input matrix
%                             %  C + F = 1:n
% e.g. (1) --------------------------------------------------------------
% n=13; G=numgrid('S',n);A=delsq(G); figure
%    lap_lab; n2=size(A,1); xy=[x(:) y(:)]; [C F]=cf_split(A); figure 
%    gplot(A,xy,'go-'), hold on, gplot(A(C,C),xy(C,:),'rs'), hold off
% figure, plot([0 1],[0 1],'o')
% for k=1:n2, text(x(k),y(k),num2str(k));hold on, end, hold off
%%     (2) --------------------------------------------------------------
% n=16; G=numgrid('A',n);A=delsq(G); figure
%    lap_lab; n2=size(A,1); xy=[x(:) y(:)]; [C F]=cf_split(A); figure 
%    gplot(A,xy,'go-'), hold on, gplot(A(C,C),xy(C,:),'rs')
%    for k=1:length(C), text(xy(C(k),1)+1/n^2,xy(C(k),2),num2str(k)); end

function [C,F] = cf_split(A)
		% Assume A is sparse so storage is less a problem 
if nargin<1, help cf_split, return, end
n=length(A); C=[]; F=[]; U=1:n; lambda=[]; tau=0.1;

%%%%%%%%%%%%%%%%%%%%%%% INITIAL (1-2) %%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:n  % Node list -- Initial Measure
   Count = find_S_T(A, i, tau); lambda(i) = Count;
end

%%%%%%%%%%%%%%%%%%%%%%% MAIN STEPS 3-7 %%%%%%%%%%%%%%%%%%%%%%%%%%%
  fprintf('For n=%d, ''%s.m'' has done:',n,mfilename)
  iter=1; Fadd=[];
while length(U)>0
  [val k]=max(lambda(U));
  k=U(k);
  C = [C k];  % Found one C-variable(4)
  U = setdiff(U, k);               % Remove k on U list  (4) Done

  [Count Index]=find_S_T(A,k,tau); 
    Index = setdiff(Index,F); 
    F=[F Index];                   % Add to F-variables' list
    U=setdiff(U, Index);           % Remove from U-list  (5) Done

    [I J K]=find(A(Index,:)); J = unique(J); % Work out N_F \ F
    Fadd=[Fadd(:); J(:)];   J=Fadd; 
    Index = setdiff(J,F);         % Work out N_F \ F
    Index = setdiff(Index,C);     % Work out N_F \ C (boundary of F)

    %...................Work on i \n S_j^T..........................
     for j=1:length(Index)
          p=Index(j);
          [Count Set_j]=find_S_T(A,p,tau);
        for L=1:length(Set_j)
          [Count Set_k]=find_S_T(A,Set_j(L),tau);
          Set_U = intersect(Set_k, U);
          Set_F = intersect(Set_k, F);
          lambda(Set_j(L))=length(Set_U) + 2*length(Set_F);
        end
     end %j------------------------------------ (6) done
    %...................Work on i \n S_j^T..........................
    if rem(iter,10)==0,
      k=n-length(U); fprintf(' %d%%',round(k/n*100));
    end
      iter=iter+1;
end %%while length(U)>0
    fprintf('\n Obtained coarse nodes = %d and fine nodes = %d\n',...
    length(C), length(F))
%===========================================END=====================

function [Counti,index_list] = find_S_T(A, k, tau)
      Ni = A(:,k); index_list=[]; Counti = 0;
      Col = find(Ni<0); % Get a row index list / Candidates for S_i^T
     for k=1:length(Col)
      Row = find( A(Col(k),:)<0 ); % Get a col index list 
      Row = setdiff(Row, Col(k)) ; % Take away diag position
      Val = abs(A(Col(k),Row)); 
      if length(Val)>0, 
       Val = Val > tau*mean(Val); Val = find(Val); % n-coupled guys if any
      end
      if length(Val)>0
        Counti = Counti + 1; 
        index_list = [ index_list Col(k)]; % Accumulate members in S^T
      end
     end %k
        index_list = setdiff( index_list, k);
