function D = howardn(P,R)



% Algoritmus Howardova iteračního postupu pro libovolny pocet stavu a

% strategii

% Vstupni parametry:

%   P...matice prechodu vsech strategii naskladane do jedne matice

%   R...matice vynosu vsech strategii naskladane do jedne matice

%   Mame-li jednotlive matice P1,P2,...,Pk a R1,R2,...,Rk, doporucuje se

%   uvadet funkci ve tvaru: howardn([P1,P2,...,Pk],[R1,R2,...,Rk])

% Vystup:

%   D...matice obsahujici optimalni strategii pro jednotlive obdobi,

%   jednotliva obdobi jsou oznacena cislem radku



[a,b]=size(P);



% vypocet matice Q, ktera obsahuje v radcich vektory q pro jednotlive

% strategie

for k=1:(b/a)

    Q(k,:)=diag(P(1:a,(1+a*(k-1)):(a*k))*R(1:a,(1+a*(k-1)):(a*k))');

end;



% matice D, obsahujici vektory strategii pro jednotlive obdobi

[x,d]=max(Q);

D(1,:)=d;



% nastaveni pocatecnich podminek pro cyklus

konec=false;

j=1;



while konec==false



% indexujeme strategii z minuleho obdobi    

    d=D(j,:);

       

% sestrojime rovnici podle vektoru strategii    

    for k=1:a

        PP(k,:)=P(k,(a*d(k)-a+1):(a*d(k)));

        QQ(k)=Q(d(k),k);

    end;

    

% vyresime rovnici a ziskame vektor v, pro pripad pozdejsiho vypisu jej

% take ulozime do matice V, kterou vsak nebudeme ve vystupu zobrazovat

    A=eye(a)-PP;

    v=inv([ones(a,1),A(:,1:(a-1))])*QQ';

    v(a+1)=0;

    V(j,:)=v;

    

% nahrazeni prvni slozky ve vektoru v za 1, zbavujeme se tak g, ktere vsak

% v dalsim vypoctu nepotrebujeme

    v(1)=1;



% pocitame optimalni strategii pro soucasne obdobi

    for k=1:(b/a)

        M(k,:)= [Q(k,:)',P(1:a,(1+a*(k-1)):(a*k))]*v;

    end;

       

    [x,m]=max(M);

        

% a vysledek ulozime do nasi matice D    

    D(j+1,:)=m;

    

% kontrola, zdali je splnena podminka pro ukonceni algoritmu

    if sum(all(D(j,:)==D(j+1,:)))==1

        konec=true;

    end;



% navysujeme krok v obdobi

    j=j+1;    

    

end;    



end