function [silist,si_size, s]=blockspectral(msanum, racenum, ndist)

global avgdeg avgdegsr

warning off MATLAB:eigs:SigmaChangedToLA

if racenum==5
    racenum=4;
    race=1;
else
    race=racenum;
    racenum=1;
end

[aij1,aij2,aij3,aij4,values,sizeA,blockids]=neighbors(msanum,ndist);

%loop through all races, or just one if one race is specified
s_index=1;        %index for individual si's

for j=1:racenum
    
    %selects appropriate indices based on race
    if race==1
        ai=aij1(1,:);
        aj=aij1(2,:);
    elseif race==2
        if isempty(aij2)
            ai=[];
            aj=[];
        else
            ai=aij2(1,:);
            aj=aij2(2,:);
        end
    elseif race==3
        if isempty(aij3)
            ai=[];
            aj=[];
        else
            ai=aij3(1,:);
            aj=aij3(2,:);
        end
    elseif race==4
        if isempty(aij4)
            ai=[];
            aj=[];
        else
            ai=aij4(1,:);
            aj=aij4(2,:);
        end
    end

    
    
    if isempty(ai)==0       %requirement for neighbors in order to attempt running this
        values_curr=values(race,1:size(ai,2));      %values=neigh in "neighbors.m"
    
        for ii=1:sizeA(race)
            %finds connected components.  "done_rows" are rows whose columns
            %have already been added to the to-do list stored in "matrix_rows"
            % matrix_rows == "in" list
            % done_rows == "done" list
            done_rows=ai(1);
            matrix_rows=aj(find(ai==ai(1)));

            for i=1:size(ai,2)
                matrix_rows=setdiff(matrix_rows,done_rows);
                if isempty(matrix_rows)
                    break
                end
                done_rows=union(done_rows,matrix_rows(1));
                next_ai=aj(find(ai==matrix_rows(1)));
                matrix_rows=union(matrix_rows,next_ai);
            end
            
            %gets block ids
            sub(ii).id=sort(blockids(race,done_rows));
            
            %finds indices of ai/aj that are contained within the connected set
            %"done_rows"
            index=find(ai==done_rows(1));
            for jj=2:length(done_rows)
                index=union(index,find(ai==done_rows(jj)));
            end
        
            %condenses, removing zero rows and columns, but preserving index
            %pairs
            aitemp=ai(index)';
            ajtemp=aj(index)';
            valtemp=values_curr(index)';
            aiaj=cat(2,aitemp,ajtemp,valtemp);      %this looks like a stata dataset (same format)
            for ij=1:2
                aiaj=sortrows(aiaj,ij);
                iaiaj=aiaj(1,ij);
                ivalue=1;
                for kk=1:size(aiaj,1)           
                    if aiaj(kk,ij)==iaiaj           % counting unique rows and overwriting old numbers
                        aiaj(kk,ij)=ivalue;
                    else
                        iaiaj=aiaj(kk,ij);
                        ivalue=ivalue+1;
                        aiaj(kk,ij)=ivalue;
                    end
                end
            end
        
            %saves indices and values as entry ii in structure array, then
            %deletes those values from ai,aj, values_curr before finding next
            %submatrix
            sub(ii).ai=aiaj(:,1);
            sub(ii).aj=aiaj(:,2);
            sub(ii).val=aiaj(:,3);
        
            ai(index)=[];
            aj(index)=[];
            values_curr(index)=[];
        
            if isempty(ai)
                break
            end
            clear aitemp ajtemp valtemp;
        end
    
        %calculated largest eigenvalue for each submatrix, stored as "largest",
        %and also stores dimension of each submatrix
        
        for k=1:size(sub,2)
            dim(k)=max(sub(k).ai);
            A=sparse(sub(k).ai,sub(k).aj,sub(k).val,dim(k),dim(k));
            [U,L]=eigs(A,1,'LR');
            largest(k)=abs(L(1,1));
            u = abs(U);
            suma = ones(1, size(A,1)) * u;
            u=u ./ suma;
            
            clear A U L;
            
            %individual si's stored in s
            s(s_index:(s_index+dim(k)-1),1) = k;                 %which connected component
            s(s_index:(s_index+dim(k)-1),2) = (sub(k).id)';      %block id
            s(s_index:(s_index+dim(k)-1),3) = u * largest(k);   %individual si's
            s(s_index:(s_index+dim(k)-1),4) = largest(k);                    %connected component si
            
            si_size=size(s,1);
            
            s_index=s_index+dim(k);

        end

    else        %if there are no neighbors, si=0
        largest=0;
        dim=0;
        sub=[];
    end
        
    singletons=sizeA(race)-sum(dim);
    %calculates spectral index and saves for output along with race & # of
    %connected components
    silist(j,1)=race;
    silist(j,2)=sum(largest.*dim)/sizeA(race);
    silist(j,3)=size(sub,2);
    silist(j,4)=singletons;
    silist(j,5)=median(dim);
    silist(j,6)=max(dim);
    silist(j,7)=min(dim);
    silist(j,8)=sizeA(race);
    silist(j,9)=sizeA(race)/sum(sizeA);         %race percentage
    silist(j,10)=avgdeg(race);
    silist(j,11)=avgdegsr(race);
    
    clear largest dim sub;
    
    race=race+1;
end