function [f,ind,p] = frob(a,tol)

%       function [f,ind,p] = frob(a,tol)
%       Uses dig.m to find the FROB-enius Normal Form 'f' of A
%
%       ind =  equivalence classes (strongly connected components)
%              sorted out by accessibility  (left ---> right)
%              and separated by zeros
%        p  = permutation matrix such that f=p*a*p'

    [m,n]=size(a); if m~=n 'Not a Square Matrix',break, end
    if (nargin==1) tol=n*eps*norm(a,'inf'); end

[c,v]=dig(a,tol);
if v(1,1) == 67 v, break, end

% Use c to put the classes in 'decreasing' order
d=v(:,1);   % d is a vector of class representatives
num=max(size(d));
    i=1;j=2;
    while(i<num)
        while(j<=num)
             % if class j has access to class i, put j first (swap)
             if ( c(d(j),d(i)) ~= 0 )
                 note=v(i,:); v(i,:)=v(j,:); v(j,:)=note;
                 % update and start again
                 j=i+1; d=v(:,1);
             else
                 j=j+1;
             end
        end
        i=i+1; j=i+1;
    end

% Create a sequence of indices with classes separated by a zero
l=1;
for i=1:num
        flag=0; % I haven't encountered a zero yet
        for j=1:n
              if v(i,j)~=0
                  ind(l)=v(i,j); l=l+1;
              else
                  if  flag == 0
                      ind(l)=0; l=l+1;
                      flag=1; % zero taken care of
                  end
              end
        end
end
ind=ind(1:max(size(ind))-1); % get rid of annoying zero

%       Finally we will create the permutation matrix P, such that
%                          PAP' = F
p=zeros(size(a)); i=1;
for j=1:max(size(ind))
      if ind(j) ~= 0
          p(i,ind(j))=1; i=i+1;   % Well,ind has n nonzero entries !
      end
end


% Here is the Frobenius Normal form:
f=p*a*p';