annotate toolboxes/FullBNT-1.0.7/graph/scc.m @ 0:e9a9cd732c1e tip

first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [c,v] = scc(a,tol)
wolffd@0 2
wolffd@0 3 % Finds the strongly connected sets of vertices
wolffd@0 4 % in the DI-rected G-raph of A
wolffd@0 5 % c = 0-1 matrix displaying accessibility
wolffd@0 6 % v = displays the equivalent classes
wolffd@0 7 %
wolffd@0 8 % v(i,j) is the j'th member of the i'th equiv class (0 padded)
wolffd@0 9 %
wolffd@0 10 % http://www.math.wsu.edu/math/faculty/tsat/matlab.html
wolffd@0 11
wolffd@0 12 [m,n] = size(a);
wolffd@0 13 if m~=n 'Not a Square Matrix', return, end
wolffd@0 14 b=abs(a); o=ones(size(a)); x=zeros(1,n);
wolffd@0 15 msg='The Matrix is Irreducible !'; v='Connected Directed Graph !';
wolffd@0 16 if (nargin==1) tol=n*eps*norm(a,'inf'); end
wolffd@0 17
wolffd@0 18 % Create a companion matrix
wolffd@0 19 c = b>tol*o;
wolffd@0 20 if (c==o)
wolffd@0 21 % msg, return
wolffd@0 22 v = 1:length(a);
wolffd@0 23 return
wolffd@0 24 end
wolffd@0 25
wolffd@0 26
wolffd@0 27 % Compute accessibility in at most n-step paths
wolffd@0 28 for k=1:n
wolffd@0 29 for j=1:n
wolffd@0 30 for i=1:n
wolffd@0 31 % If index i accesses j, where can you go ?
wolffd@0 32 if c(i,j) > 0 c(i,:) = c(i,:)+c(j,:); end
wolffd@0 33 end
wolffd@0 34 end
wolffd@0 35 end
wolffd@0 36 % Create a 0-1 matrix with the above information
wolffd@0 37 c>zeros(size(a)); c=ans; if (c==o) msg, return, end
wolffd@0 38
wolffd@0 39 % Identify equivalence classes
wolffd@0 40 d=c.*c'+eye(size(a)); d>zeros(size(a)); d=ans;
wolffd@0 41 v=zeros(size(a));
wolffd@0 42 for i=1:n find(d(i,:)); ans(n)=0; v(i,:)=ans; end
wolffd@0 43
wolffd@0 44 % Eliminate displaying of identical rows
wolffd@0 45 i=1;
wolffd@0 46 while(i<n)
wolffd@0 47 for k=i+1:n
wolffd@0 48 if v(k,1) == v(i,1)
wolffd@0 49 v(k,:)=x;
wolffd@0 50 end
wolffd@0 51 end
wolffd@0 52 i=i+1;
wolffd@0 53 end
wolffd@0 54 j=1;
wolffd@0 55 for i=1:n
wolffd@0 56 if v(i,1)>0
wolffd@0 57 h(j,:)=v(i,:);
wolffd@0 58 j=j+1;
wolffd@0 59 end
wolffd@0 60 end
wolffd@0 61 v=h;
wolffd@0 62
wolffd@0 63
wolffd@0 64