function sub=greedy(ien,nsub)
% sub=greedy(ien,nsub)
% decompose elements into substructures by a simple
% greedy algorithm
% try to have the same number of elements per subsctructure
% input
%   ien(a,e)    global node number for node a on element e
%   nsub        number of substructure for every element
% output
%   sub(e)      numer of substructure element e is in

% get element-node adjacency
[mnodes,nel]=size(ien); % max nodes per element, number of elements
nnod=max(max(ien));     % number of nodes
nod_el_adj=spalloc(nnod,nel,nnz(ien));
for e=1:nel 
    nodes=ien(:,e);             % nodes on this element
    nodes=nodes(find(nodes));   % get rid of zeros if any
    nod_el_adj(:,e)=sparse(nodes,1,1,nnod,1); % create column
end
% el_el_adj(i,j) is the number of nodes elements have in common
el_el_adj=nod_el_adj'*nod_el_adj;

sub=zeros(1,nel); % initialize what substructure element goes to
isub=0; % no substructures created
taken=0; % no elements assigned
for iel=1:nel
    if sub(iel)==0   % if element not assigned to substructure yet
        isub=isub+1; % start a new substructure
        cand=iel;    % start list of candidates to be added
    end
    while taken < nel*(isub/nsub) % keep adding until too big
        if length(cand)>0, % have adjacent candidate?
            e=cand(1);  % take first one (should be smarter)
            % IN MATLAB, COMPLEXITY QUADRATIC IN SUBDOMAIN SIZE
            cand(1)=[]; % delete from candidate list
            sub(e)=isub; % add to substructure
            taken=taken+1;
            neig=find(el_el_adj(:,e)); % all neighbors
            neig=neig(find(sub(neig)==0)); % pick only unassigned
            cand=union(cand,neig);
        else    % run out of candidates
            % does not happen often otherwise IN MATLAB COMPLEXITY QUADRATIC
            unassigned=find(sub==0);
            cand=unassigned(1);
        end
    end
end