wolffd@0: function [c,v] = scc(a,tol)
wolffd@0: 
wolffd@0: %       Finds the strongly connected sets of vertices
wolffd@0: %                in the DI-rected G-raph of A
wolffd@0: %          c = 0-1 matrix displaying accessibility
wolffd@0: %          v = displays the equivalent classes
wolffd@0: %
wolffd@0: % v(i,j) is the j'th member of the i'th equiv class (0 padded)
wolffd@0: %
wolffd@0: % http://www.math.wsu.edu/math/faculty/tsat/matlab.html
wolffd@0: 
wolffd@0: [m,n] = size(a);
wolffd@0: if m~=n 'Not a Square Matrix', return, end
wolffd@0: b=abs(a); o=ones(size(a)); x=zeros(1,n);
wolffd@0: msg='The Matrix is Irreducible !'; v='Connected Directed Graph !';
wolffd@0: if (nargin==1) tol=n*eps*norm(a,'inf'); end
wolffd@0: 
wolffd@0: % Create a companion matrix
wolffd@0: c = b>tol*o;
wolffd@0: if (c==o)
wolffd@0:   %  msg, return
wolffd@0:   v = 1:length(a);
wolffd@0:   return
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % Compute accessibility in at most n-step paths
wolffd@0: for k=1:n
wolffd@0:   for j=1:n
wolffd@0:     for i=1:n
wolffd@0:       % If index i accesses j, where can you go ?
wolffd@0:       if c(i,j) > 0  c(i,:) = c(i,:)+c(j,:); end
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: % Create a 0-1 matrix with the above information
wolffd@0: c>zeros(size(a)); c=ans; if (c==o) msg, return, end
wolffd@0: 
wolffd@0: % Identify equivalence classes
wolffd@0: d=c.*c'+eye(size(a)); d>zeros(size(a)); d=ans;
wolffd@0: v=zeros(size(a));
wolffd@0: for i=1:n find(d(i,:)); ans(n)=0; v(i,:)=ans; end
wolffd@0: 
wolffd@0: % Eliminate displaying of identical rows
wolffd@0: i=1;
wolffd@0: while(i<n)
wolffd@0:   for k=i+1:n
wolffd@0:     if v(k,1) == v(i,1)
wolffd@0:       v(k,:)=x;
wolffd@0:     end
wolffd@0:   end
wolffd@0:   i=i+1;
wolffd@0: end
wolffd@0: j=1;
wolffd@0: for i=1:n
wolffd@0:   if v(i,1)>0
wolffd@0:     h(j,:)=v(i,:);
wolffd@0:     j=j+1;
wolffd@0:   end
wolffd@0: end
wolffd@0: v=h;
wolffd@0: 
wolffd@0: 
wolffd@0: