Jak widać z powyższego — i tak właśnie bardzo często zdarza się w najrozmaitszych zadaniach obliczeniowych matematyki stosowanej — nasz probleme skaluje się, to znaczy można go postawić zarówno dla niewielkiej liczby niewiadomych (na przykład dla n=2 i K=M=4), jak też dla bardzo dużej liczby parametrów (wystarczy wziąć dostatecznie duże n: już dla n=64 liczba samych stanów źródeł, 264, przekroczy możliwości adresowania współczesnych procesorów!).

Możemy więc na próbę, dla testów i wstępnych obserwacji spróbować wykorzystać środowisko Octave dla małych n,K,M. Utworzone przez nas narzędzia dla zadania małej skali mogą przydać się np, do pre- lub postprocessingu danych lub wyników zadania w dużej skali.

Naszą macierz W będziemy tworzyć wierszami, według schematu:

for r=0:M-1
	for q = 0:K-1
		for s = 0:2^n-1
			% wyznacz numer stanu i odpowiadający trójce (s,q,r)
			% wyznacz częstości w_ij przejścia od stanu i=(s,q,r) do innego stanu j
		end
	end
end

Alternatywnie, moglibyśmy użyć jednej pętli po wszystkich numerach stanów:

for i=1:N
	% wyznacz trójkę (s,q,r) odpowiadającą stanowi i
	% wyznacz częstości przejścia od stanu (s,q,r) do innego stanu j
end

Będzie więc nam potrzebna funkcja, która danemu stanowi, reprezentowanemu przez trójkę s,q,r, przyporządkowuje jego numer kolejny, na przykład

function i = triple2state(x,z,p)
% nadaje numer "i" stanowi reprezentowanemu trojka [x(1),x(2),x(3)]
% zakladamy, ze x(1) = 0..2^n-1, x(2) = 0..K-1, x(3) = 0..M-1,
% 	gdzie z = [2^n K M]
% p jest sposobem numerowania, np. aby dostac sposob "sqr", bierzemy p = [1 2 3],
% 	dla sposobu "qrs" kladziemy p = [2 3 1], itp.

% schemat "p(1)-p(2)-p(3)"
 
i = 1 + x(p(1)) + z(p(1)) * ( x(p(2)) + z(p(2)) * x(p(3)) ); 	
 	% stany powinnismy numerowac od jedynki, bo numer stanu
 	% bedzie indeksem w macierzy!

end

Najwyższy czas, by sięgnąć do literatury. W monografii [28], dotyczącej numerycznego wyznaczania różnych cech — także stanów stacjonarnych — łańcuchów Markowa jest omówiona metoda, w której poprzez usunięcie jednego wiersza i jednej kolumny z macierzy W można zadanie sprowadzić do zwykłego zadania rozwiązywania układu równań liniowych z macierzą nieosobliwą. Zapewne dla średnich rozmiarów macierzy W to byłaby całkiem dobra alternatywa — wszak W jest macierzą rozrzedzoną i taką pozostanie po usunięciu wiersza i kolumny — ale musielibyśmy nieco obawiać się postępującego wypełnienia bloków czynników rozkładu LU macierzy.

Inny sposób podejścia do zadania wyznaczenia stanu stacjonarnego łańcucha Markowa (por. [27]) to wykonanie tzw. odwrotnej metody potęgowej na macierzy (raz tylko sfaktoryzowanej!) WT-ϵ⁢I. Jest to metoda iteracyjna rozwiązywania zagadnienia własnego, której zasadniczym składnikiem jest wielokrotne rozwiązywanie układu równań postaci WT-ϵ⁢I. Aby więc metoda była efektywna, musimy umieć wyznaczyć rozkład macierzy WT-ϵ⁢I na prostsze składniki, najlepiej — wykorzystując fakt, że macierz WT jest mocno rozrzedzona. Więcej na temat metody potęgowej i odwrotnej metody potęgowej można przeczytać w wykładzie z Matematyki Obliczeniowej~II. Jednak w naszym przypadku domyślnie stosowana w Octave metoda rozkładu macierzy rozrzedzonych zabiera zbyt wiele czasu i jest mało efektywna, gdy N jest duże.

Na szczęście, jest jeszcze jedna droga, wykorzystująca dodatkowe właściwości naszego zadania. (Skąd my to znamy?…) Skoro mamy do czynienia z łańcuchem Markowa, to możemy zapisać zadanie w równoważnej postaci z macierzą stochastyczną³⁰Czyli taką, której elementy są nieujemne i w każdej kolumnie sumują się do jedności. P,

gdzie P=δ⁢WT+I, przy czym δ=0.99/maxi⁡wi⁢i (dowód równoważności i stochastyczności P zostawiamy Czytelnikowi jako niezobowiązujące ćwiczonko, zob. [27]). Zatem poszukiwany wektor π jest wektorem własnym macierzy (kolumnami) stochastycznej P, odpowiadającym wartości własnej równej 1. Jest to znacząca informacja, gdyż to pociąga za sobą, że:

• istnieje wektor własny odpowiadający wartości własnej 1 (rzeczywiście, lewostronny wektor własny dla 1 to wektor 1,…,1.

• wartość własna 1 jest dominującą wartością własną (rzeczywiście, λ⁢P≤P1=1) i nie ma innych wartości własnych o module równym 1 (wynika to z twierdzenia Gerszgorina).

Ponadto nasza macierz P jest nieredukowalna (w języku łańcuchów Markowa, łańcuch jest nieprzywiedlny: każdy stan jest osiągalny z każdego innego, w skończonej liczbie kroków), więc na mocy twierdzenia Perrona–Frobeniusa [3, 28] wartość własna 1 macierzy P jest pojedyncza, a odpowiada jej wektor własny π o dodatnich współrzędnych.

Te wszystkie fakty oznaczają, że — zamiast kosztownego SVD — możemy spróbować użyć na przykład metody potęgowej wyznaczania dominującej wartości własnej macierzy P!

Należy pamiętać, że metoda potęgowa, oprócz ewidentnej zalety (jest metodą iteracyjną, więc można obliczenia przerwać w bardziej dogodnym momencie, a iteracja opiera się na bardzo tanim mnożeniu, πk+1=P⁢πk, przez macierz rozrzedzoną P; ma też minimalne wymagania pamięciowe, co istotne przy dużych n), ma też wady: jeśli druga co do wielkości wartość własna macierzy, λ2, jest bliska dominującej — λ2≈1 — zbieżność metody może być bardzo wolna.

W Octave możemy wyznaczyć dominującą wartość własną na dwa sposoby: zaimplementować wprost metodę potęgową, albo wykorzystać funkcję eigs, wykorzystującą bardziej zaawansowane techniki wyznaczania ekstremalnych wartości własnych. Poniżej wybieramy właśnie ten drugi wariant:

opts.maxit = MAXIT;
opts.tol = TOL;
tic; [V, D, info] = eigs(spW, 6, 'lm', opts); toc
D = diag(D); [m i] = max(abs(D));
if (abs(m-1) > 1e1*eps)
	warning(['Dominant eigenvalue = ', num2str(m, '%e'), ' not equal to 1. Difference: ',  num2str(abs(m-1))]);
end
Pi = V(:,i);
Pi = abs(Pi); % ustalamy znak Pi
Pi = Pi/sum(Pi);

Pozostaje jeszcze przerobić macierz W na macierz P; w naszym przypadku zastąpimy funkcję W = sptelecomtrans(..) funkcją [P, W] = sptelecom(..), wyznaczającą przede wszystkim macierz P, a macierz W — tylko opcjonalnie.

function [spP, spW] = sptelecom(n,K,M,K1,rate)
% create (sparse) transition probability matrix spP and (if requested)
% transition rate matrix spW
%
% queue length 0..K-1
% source states 0..2^n-1
% processing phases 0..M-1
% queue threshold K1

if nargin < 5
	rate = [1 1 1 1 1];
end

a = 	rate(1); 	% off -> on
b = 	rate(2);	% on -> off
alpha =	rate(3); 	% src -> queue
beta = 	rate(4);	% priority src -> queue
mu = 	rate(5);	% advance processing stage

N = 2^n;
z = [N, K, M];
ord = [1 2 3];

% fprintf(stderr,'Problem data:\nn=%d N=%d K=%d M=%d K1=%d\na=%g b=%g alpha=%g beta=%g mu=%g\n', n, N, K, M, K1, a, b, alpha, beta, mu);

dim = triple2state(z-1,z,ord); % total dimension of W, dim = N*M*K
nz = 0; % count nonzero elements of W

% disp(['Constructing sparse matrix of size ', num2str(dim), ' and at most ', num2str(dim*(n+3)),' zeros']);

IW = nan(dim*(n+3),1); % oszacowanie liczby elementow niezerowych
JW = nan(size(IW));
VW = zeros(size(IW));
RW = zeros(dim,1);

for s = 0:N-1 % stan zrodla
	for q = 0:K-1 % stan kolejki
		for r = 0:M-1 % stan procesora
			S = triple2state([s,q,r],z,ord); % numer biezacego stanu

			%%%%%%%% zmiana stanu z powodu zmiany stanu zrodla

			for i = 1:n
				bit = bitget(s,i);
				sn = bitset(s,i,~bit); % zmienia i-ty bit na przeciwny: to nowy numer stanu zrodla
				SN = triple2state([sn,q,r],z,ord); % numer nowego stanu zrodla
				
				nz = nz+1;
				IW(nz) = S; JW(nz) = SN;
				
				if bit == 0 
					% przejscie od wylaczonego do wlaczonego
					VW(nz) = a;
					RW(S) = RW(S) + a;
				else
					% przejscie od wlaczonego do wylaczonego
					VW(nz) = b;
					RW(S) = RW(S) + b;
				end
			end
			
			%%%%%%%% zmiana stanu z powodu zmiany stanu kolejki

			if q < (K-1) % tylko wtedy mozna wydluzyc kolejke, w przeciwnym razie wszystkie pakiety sa odrzucane
				SN = triple2state([s,q+1,r],z,ord); % nowy stan: przyjeto pakiet do kolejki
				bs = bitsum(s,n); % ile jest wlaczonych sposrod n mozliwych
				b0 = bitget(s,1); % czy zerowy jest wlaczony
				% fprintf(stderr,'s = %d bit = %d sum = %d\n',s,b0, bs); 
				if bs > 0 % cokolwiek nadaje
					
					nz = nz+1;
					IW(nz) = S; JW(nz) = SN;
					VW(nz) = 0;
					
					if b0 > 0 % specjalne traktowanie pierwszego zrodla
						VW(nz) = beta;
						RW(S) = RW(S) + beta;
					end
					if q < (K1-1) % pakiety z innych sa odrzucane jesli kolejka ma dlugosc K1 lub wieksza
						VW(nz) = VW(nz) + (bs-b0)*alpha;
						RW(S) = RW(S) + (bs-b0)*alpha;
					end
				end
			end


			%%%%%%%% zmiana stanu z powodu zmiany stanu procesora
			
			switch r
			case (M-1)
				SN = triple2state([s,q,0],z,ord);
				nz = nz+1;
				IW(nz) = S; JW(nz) = SN;
				VW(nz) = mu;
				RW(S) = RW(S) + mu;
			case 0
				if (q > 0)
					SN = triple2state([s,q-1,1],z,ord);
					nz = nz+1;
					IW(nz) = S; JW(nz) = SN;
					VW(nz) = mu;
					RW(S) = RW(S) + mu;
				end
			otherwise
				SN = triple2state([s,q,r+1],z,ord);
				nz = nz+1;
				IW(nz) = S; JW(nz) = SN;
				VW(nz) = mu;
				RW(S) = RW(S) + mu;
			end	
		end % for...
	end % for...
end % for...

% suma wyrazow w wierszu: RW(i) = pstwo pozostania w stanie i
% W + diag(RW) to macierz czestosci
% konstrukcja macierzy prawdopodobienstw

deltaW = 0.99/max(abs(RW));
VW = VW*deltaW;
IW(nz+[1:dim]) = [1:dim];
JW(nz+[1:dim]) = [1:dim];
VW(nz+[1:dim]) = 1-RW*deltaW;
nz = nz + dim;

disp([num2str(nz), ' nonzeros in the probability matrix P']);

file=fopen('mout.dat','w');
fprintf(file, '%d %d %g\n', [JW(1:nz),IW(1:nz),VW(1:nz)]');
fclose(file);

spP = sparse(JW(1:nz),IW(1:nz),VW(1:nz),dim,dim);

% sprawdzmy, czy teraz W jest kolumnami stochastyczna
if (norm(spsum(spP,1)-1,Inf) > 10*eps) || any(spP(:) < 0) || any(spP(:) > 1)
	warning('P may be not a Markov transition probability matrix');
end

if nargout > 1
	VW = VW/deltaW; % niezbyt eleganckie
	VW(nz-dim+1:nz) = -RW;
	spW = sparse(IW(1:nz),JW(1:nz),VW(1:nz),dim,dim);
end

end % function

function bs = bitsum(s,n)
% wyznacza liczbe zapalonych bitow liczby s
bs = sum(bitget(s,1:n));
end % function

function info = isprobability(p)
% sprawdza, czy 0 <= p <= 1
info = all(p >= 0) & all(p <= 1);
if (~info) && (nargout < 1)
	warning([inputname(1), ' = [', num2str(p,' %g'), '] is not a probability']);
end
end % function

Parametry zadania zapiszemy w prostym skrypcie Octave, o nazwie telecomparam.m:

a 	= 3.33;
b 	= 3.52;
alpha 	= 1;
beta 	= 2;
mu = 34e-1; % uwaga

n = 13;
K = 8;
M = 5;
K1 = 7;
MAXIT = 100;
TOL = 1e-6;

Aby wygenerować macierz P, wystarczy więc wydać polecenia:

telecomparam;
P = sptelecom(n,K,M,K1,[a,b,alpha,beta,mu]);

W zasadzie wystarczy wprost przetłumaczyć tekst procedury z Octave na język C, pamiętając wszakże o następujących pułapkach, które czyhają na nas po drodze:

• tablice w C indeksujemy od zera, dlatego

  • licznik nz będziemy inicjalizować na -1,

  • funkcja triple2state będzie zwracać s+N*(q+K*r), a nie jak w Octave, 1+s+N*(q+K*r),

• bity w C numerujemy od zera, dlatego pętla po bitach źródeł powinna być dla i=0,…,n-1: for (i = 0; i < n; i++), a bitem źródła uprzywilejowanego źródła jest bit zerowy b0 = bitget(s,0); funkcja bitget(s,i) powinna akceptować wartości i=0,…,n-1;

• w kodzie zakładamy, że bitget() zwraca zero lub 1, dlatego nie możemy użyć prostego s & (1<<i), które może zwracać wartości większe od 1;

Ponadto, nie będziemy już generować (niepotrzebnej w tej chwili) macierzy częstości W, a tylko samą macierz P.

Warto zwrócić uwagę na fragment pozwalający prosto odczytać wszystkie parametry zadania zapisane w pliku Octave telecomparam.m.

Głównym celem programu sptelecom.c jest wyznaczenie macierzy rozrzedzonej P w formacie współrzędnych i zapisanie jej do pliku (aby można było następnie wczytać ją do Octave i potraktować eigs).

Tu kryje się jeszcze jedna pułapka. Zapisując macierz w formacie tekstowym, ryzykujemy utratę dokładności reprezentacji wartości elementów macierzy P. Dlatego, lepiej będzie zapisywać dane w formacie binarnym — i tak faktycznie zrobimy: patrz kod źródłowy programu w języku C.

/*
	gcc -o sptelecom sptelecom.c -lm
	time ./sptelecom <  telecomparam.m
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define NPARAM 11 /* liczba parametrow do wczytania z pliku */

/* zmienne, do ktorych wczytamy parametry z pliku */
int n, N, K, M, K1;
double a, b, alpha, beta, mu;
int MAXIT; double TOL;

int triple2state(int s, int q, int r)
{
	return(s+N*(q+K*r)); /* zero-based */
}

int bitget(int s, int i)
{
	return( (s & (1<<i)) > 0 ); /* zwraca zero lub 1 */
}

int bitflip(int s, int i)
{
	return( s ^ (1<<i) );
}

int bitsum(int s)
{
int bs = 0, i;

if (s <= 1)
	return(s);
for (i=0; i < n; i++)
{
	bs += bitget(s,i);
}
return(bs);
}

void multW(int *IW, int *JW, double *VW, double *X, double *Y, int dim, int nz)
{
/*
 Mnozy macierz W (w formacie wspolrzednych IW,JW,VW) przez wektor X i wynik zapisuje do Y
*/
int i;	
	for(i = 0; i < dim; i++)
		Y[i] = 0.0;
	for(i = 0; i < nz; i++)
		Y[IW[i]] += VW[i]*X[JW[i]];
}

int main(int argc, char **argv)
{
int s,q,r, S, SN, dim, nz, nnz, sn, bit, i;
int *IW, *JW;
double *VW, *RW;
double *X, *Y, nrm;
FILE *fp;

/* prymitywny lecz skuteczny sposob wczytania parametrow */
/*	
	#include "telecomparam.m"
*/

{
/* standardowe wczytanie parametrow */
char name[32];
double val;
double params[NPARAM];
char names[NPARAM][15] = {"n", "K", "M", "K1", 
	"a", "b", "alpha", "beta", "mu", "MAXIT", "TOL" };
int i,k;

	for (i=0; i < NPARAM; i++) 
		params[i] = -1;

	for (i=0; i < NPARAM; i++) 
	{
		scanf("%s = %lf;%*[^\n]\n", name, &val);
		fprintf(stderr, "\tReading: %g for %s\n", val, name);
		for (k = 0; k < NPARAM; k++)
			if (strcmp(name,names[k]) == 0)
			{
				params[k] = val;
				break;
				
			}
	}
	
	k = 1;
	for (i=0; i < NPARAM; i++)
	{ 
		if (params[i] < 0)
		{
			k = -1;
			fprintf(stderr, "Parameter '%s' not set\n", names[i]);
		}
	}
	if (k == -1)
	{
		fprintf(stderr, "Exiting!\n");
		exit(-1);
	}

	n 	= params[0];
	K 	= params[1];
	M 	= params[2];
	K1 	= params[3];
	a 	= params[4];
	b 	= params[5];
	alpha 	= params[6];
	beta 	= params[7];
	mu 	= params[8];
	MAXIT 	= params[9];
	TOL 	= params[10];
}

if (argc > 1)
{
	MAXIT = atoi(argv[1]);
	fprintf(stderr, "\tChanging MAXIT to %d\n", MAXIT);
}

N = 1<<n;
fprintf(stderr,"Problem data:\nn=%d N=%d K=%d M=%d K1=%d\na=%g b=%g alpha=%g beta=%g mu=%g\n", n, N, K, M, K1, a, b, alpha, beta, mu);

dim = triple2state(N-1,K-1,M-1)+1;
nz = -1; /* indeks ostatnio dodanego niezerowego elementu macierzy W */
nnz = dim*(n+3); /* nieco przeszacowana maksymalna liczba niezerowych elementow w macierzy W */
fprintf(stderr,"Constructing sparse matrix of size %d and at most %d nonzeros.\n", dim, nnz);

IW = (int *)calloc(nnz, sizeof(int)); 
JW = (int *)calloc(nnz, sizeof(int)); 
VW = (double *)calloc(nnz, sizeof(double)); 
RW = (double *)calloc(dim, sizeof(double)); /* tablice wyzerowane! */

for (s = 0; s < N; s++) /* stan zrodla */
	for (q = 0; q < K; q++) /* stan kolejki */
		for (r = 0; r < M; r++) /* stan procesora */
		{
		S = triple2state(s,q,r); /* numer biezacego stanu */
		
		/*
		 * zmiana stanu z powodu zmiany stanu zrodla
		 */
		 
			for (i = 0; i < n; i++) 
			{
				bit = bitget(s,i);
				sn = bitflip(s,i); /* zmienia i-ty bit na przeciwny: to nowy numer stanu zrodla */
				SN = triple2state(sn,q,r); /* numer nowego stanu zrodla */
				
				nz++;
				IW[nz] = S; JW[nz] = SN;
				
				if (bit == 0)
				{ 
					/* przejscie od wylaczonego do wlaczonego */
					VW[nz] = a;
					RW[S] += a;
				}
				else
				{
					/* przejscie od wlaczonego do wylaczonego */
					VW[nz] = b;
					RW[S] += b;
				}
			}
			
		/*
		 *  zmiana stanu z powodu zmiany stanu kolejki
		 */
		 
			if (q < (K-1)) /* tylko wtedy mozna wydluzyc kolejke, w przeciwnym razie wszystkie pakiety sa odrzucane */
			{
				int bs, b0;
				
				SN = triple2state(s,q+1,r); /* nowy stan: przyjeto pakiet do kolejki */
				bs = bitsum(s); /* ile jest wlaczonych */
				b0 = bitget(s,0); /* czy zerowy jest wlaczony */
				/* fprintf(stderr,"s = %d bit = %d sum = %d\n",s,b0, bs); */
				if (bs > 0) /* cokolwiek nadaje */
				{	
					nz++;
					IW[nz] = S; JW[nz] = SN;
					VW[nz] = 0;
					
					if (b0 > 0) /* specjalne traktowanie pierwszego zrodla */
					{
						VW[nz] = beta;
						RW[S] += beta;
					}
					if (q < (K1-1)) /* pakiety z innych sa odrzucane jesli kolejka ma dlugosc K1 lub wieksza */
					{
						VW[nz] += (bs-b0)*alpha;
						RW[S] += (bs-b0)*alpha;
					}
				}
			}
		 
		/*
		 * zmiana stanu z powodu zmiany stanu procesora
		 */
			if (r == (M-1))
			{
				SN = triple2state(s,q,0);
				nz++;
				IW[nz] = S; JW[nz] = SN;
				VW[nz] = mu;
				RW[S] += mu;
			}
			else if (r == 0)
			{
				if (q > 0)
				{
					SN = triple2state(s,q-1,1);
					nz++;
					IW[nz] = S; JW[nz] = SN;
					VW[nz] = mu;
					RW[S] += mu;
				}
			}
			else 
			{
				SN = triple2state(s,q,r+1);
				nz++;
				IW[nz] = S; JW[nz] = SN;
				VW[nz] = mu;
				RW[S] += mu;
			}	
		 
		}
nz++;

/* 
 suma wyrazow w wierszu: RW(i) = pstwo pozostania w stanie i
 W + diag(RW) to macierz czestosci
 konstrukcja macierzy prawdopodobienstw
*/

{
double deltaW = 0.0;

// max of RW
for (i = 0; i < dim; i++)
	deltaW = (RW[i] > deltaW ? RW[i] : deltaW);
	  
deltaW = 0.99/deltaW;
for (i = 0; i < nz; i++)
	VW[i] *= deltaW;
for (i = 0; i < dim; i++) 
{
	IW[nz+i] = JW[nz+i] = i;
	VW[nz+i] = 1.0-RW[i]*deltaW;
}
nz += dim;
}

free(RW);
fprintf(stderr, "%d nonzeros in the probability matrix\n", nz);

fp = fopen("cout.dat","wb"); /* zapis macierzy W w formacie binarnym, to wazne! */
fwrite(&nz, sizeof(int), 1, fp);
fwrite(JW, sizeof(int), nz, fp);
fwrite(IW, sizeof(int), nz, fp);
fwrite(VW, sizeof(double), nz, fp);
fclose(fp);

/*
 metoda potegowa, tylko dla porownania z eigs z Octave
*/

X = (double *)calloc(dim, sizeof(double));
Y = (double *)calloc(dim, sizeof(double));

for (i = 0; i < dim; i++)
	X[i] = 1.0/dim;

for (s = 0; s < MAXIT; s++)
{
	multW(JW,IW,VW,X,Y,dim,nz);
	multW(JW,IW,VW,Y,X,dim,nz);
	nrm = 0.0;
	for (i = 0; i < dim; i++)
	{
		double diff = fabs(Y[i]-X[i]);
		nrm = (diff > nrm) ? diff : nrm;
	}
	if (nrm < TOL)
		break;
}

fprintf(stderr,"Total %d iterations\n", s);

fp = fopen("coutx.dat","w"); /* wynik: wektor stanu stacjonarnego */ 
for(i = 0; i < dim; i++)
	fprintf(fp, "%g\n", X[i]);
fclose(fp);

return(0);
}

Mając macierz P w formacie współrzędnych, w postaci binarnej, musimy najpierw zmusić Octave do wczytania jej właśnie w tej postaci:

function [spW, dim] = sptelecomload(filename)
	fp = fopen(filename,'rb');
	nz = fread(fp, 1, 'int');
	IJ = fread(fp, [nz,2], 'int');
	IJ = IJ+1; % indeksy nie od zera, ale od 1
	dim = max(IJ(:,1));
	V = fread(fp, nz, 'double');
	fclose(fp);

	spW = sparse(IJ(:,1), IJ(:,2), V, dim, dim);
end

Teraz możemy już rozwiązać zadanie:

1. Wpisać poprawne wartości parametrów do pliku telecomparam.m

2. Wygenerować plik z macierzą P:

   make sptelecom
   ./sptelecom < telecomparam.m
   

3. Wyznaczyć wektor π:

   telecomparam;
   system ('make cout.dat');
   
   [P, dim] = sptelecomload('cout.dat');
   
   opts.maxit = MAXIT;
   opts.tol = TOL;
   tic; [V, D, info] = eigs(P, 2, 'lm', opts); toc
   D = diag(D); [m i] = max(abs(D));
   if (abs(m-1) > 1e1*eps)
   	warning(['Dominant eigenvalue = ', num2str(m, '%e'), ' not equal to 1. Difference: ',  num2str(abs(m-1))]);
   end
   Pi = V(:,i);
   Pi = abs(Pi);
   Pi = Pi/sum(Pi); 
   residP = norm(Pi-P*Pi,Inf)
   if (abs(residP) > 1e1*eps)
   	warning(['Residual ||(P-I)*Pi|| = ', num2str(residP)]);
   end
   
   semilogy([0:K-1],sum(sum(reshape(Pi,[2^n,K,M]),3),1),'o-');
   % semilogy([0:2^n-1],sum(sum(reshape(Pi,[2^n,K,M]),3),2),'o-');
   
   file=fopen('moutx.dat','w');
   fprintf(file, '%23.18e\n', Pi');
   fclose(file);
   
   

Ponieważ Octave potrafi uruchamiać z poziomu sesji inne programy systemowe, wykorzystaliśmy to do automatycznego wygenerowania nowych danych, jak tylko zmieniły się parametry zapisane w pliku, lub zmianie uległ sam kod źródłowy.