function sub=greedy(ien,nsub)
% sub=greedy(ien,nsub)
% decompose elements into substructures
% 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

% 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);
isub=0;
new_sub=true;
% start
sub=zeros(1,nel); % what substructure element goes to
isub=1; % start substructure 1
taken=0; % no elements assigned yet
for iel=1:nel
    if sub(iel)==0   % element not assigned to substructure yet
        isub=isub+1; % start a new substructure
        ssub=[];     % as empty
        cand=iel;    % start list of candidates to be added
    end
    while taken <= nel*(isub/nsub) % keep adding
        if length(cand)>0, % have adjacent candidate?
            e=cand(1);  % take first; should be smarter
            % IN MATLAB, COMPLEXITY QUADRATIC IN SUBDOMAIN SIZE
            cand(1)=[]; % delete from candidate list
            sub(e)=isub; % add to substructure
            neig=find(el_el_adj(:,e)); % all neighbors
            neig=neig(find(isub(neig)==0)); % pick only unassigned
            cand=union(cand,neig)
        if 
        if  % not too big?
            sub(iel)=isub; % add element to substructure
            siz=siz+1;
        end
    if new_sub
        % start new substructure
        % find first available element
        while sub(e)==0
            if imask==nel,
                
        
        % add all neighbors
    
    

end