function [F,P,E]=frobform(A)
% [F,P,E]=frobform(A) da la forma de Frobenius F de una matriz no negativa nxn
%           A irreducible no primitiva. Tambien devuelve la matriz de permutacion
%           P tal que P'AP=F. Asimismo se obtiene el cell array 1xd E de los 
%           SUBCONJUNTOS E{1},E{2},...,E{d} de {1,2,...,n} que la "cadena de
%           Markov", asociada a A, recorre periodicamente.
%
%
%    Ultima fecha de depuracion de errores: 13 de mayo de 2002
%    Autor: Juan-Miguel Gracia.

if irredble(A)==0, 
    disp('                    ============================================')    
    disp('                    ¡La matriz dada es reducible! Final de la ejecucion.')
    disp('                        Intente FROBFORM con una matriz irreducible.    ')
    disp('                    ============================================ '),
    break
end
n=length(A);

% Indice de imprimitividad d
%           primera forma de calculo.
p=poly(A);
q=invierte(p);
a=find(abs(q)>0.00001);
a=a-1;
a=diff(a);
d1=mcd(a);

%          segunda forma de calculo.
e=eig(A);
ea=abs(e);
rho=max(ea);
h=ea-rho;
h2=find(abs(h)<0.00001);
d2=length(h2);

if d1~=d2
    warning('¡Los valores de d son diferentes!')
else
    d=d2;
end



if d==1,
    F=A; P=eye(n); E={1:n};
    disp('                               ')
    disp('La matriz dada es primitiva; por lo que su forma')
    disp('de Frobenius coincide con ella misma'),
    disp('                               '),
    break
end


% Subconjuntos C{1},C{2},...,C{d}
S=sign(A);
b=1;
B=sign(S^(b*d)); 
v=1:n;
k=1; 

if B==eye(n) & n>1,
    C{1}=[1];
    for k=2:n
        C{k}=min(find(S(C{k-1}(1),:)==1));    
    end
    long=length(C);
    Gastados=[];
    for j=1:long
        Gastados=union(Gastados,C{j});      
    end
    Gastados=[Gastados,zeros(1,n-length(Gastados))];
    if all(Gastados==1:n), 
        % Permutacion P
        s=[];
        for j=1:length(C)
            s=[s,C{j}];
        end
        n=length(s);E=C;
        for i=1:n
            for j=1:n
                if i==s(j),
                    P(i,j)=1;
                else
                    P(i,j)=0;
                end
            end
        end
        % Forma de Frobenius F
        F=P.'*A*P; break
    end
end



%=========================================================
while v~=[] 
    if k==1     
        i=min(v);
        C{k}=find(B(i,:)==1);
        k=k+1; if k>=d+1, break, end
    else
        v=setdiff(v,C{k-1});
        if isempty(v), break, end
        i=min(v);
        C{k}=find(B(i,:)==1);
        k=k+1;
        if k>=d+1, 
            long=length(C);
            Gastados=[];
            for j=1:long
                Gastados=union(Gastados,C{j});      
            end
            Gastados=[Gastados,zeros(1,n-length(Gastados))];
            if all(Gastados==1:n),
                % reordenacion de los subconjuntos C{1}-->C{2}-->...-->C{d}-->C{1}
                E{1}=C{1};
                for u=1:d-1
                    f=min(find(S(E{u}(1),:)==1));
                    for l=1:d
                        if  ismember(f,C{l})==1, 
                            E{u+1}=C{l}; break,
                        end
                    end
                end
                % Permutacion P
                s=[];
                for j=1:length(E)
                    s=[s,E{j}];
                end
                n=length(s);
                for i=1:n
                    for j=1:n
                        if i==s(j),
                            P(i,j)=1;
                        else
                            P(i,j)=0;
                        end
                    end
                end
                
                % Forma de Frobenius F
                F=P.'*A*P;
            else
                b=b+1;
                B=sign(S^(b*d));    k=1; C=[];    v=1:n;
            end
        end
    end
end 
