function [x,iter,vres]=rs(P,eps,maxiter);
% Risolve (I-P')x=0 mediante iterazioni basate su regular splittings
% I-P'=M-N, che si ottengono (s)commentando parti di programma sotto.
% Input: P stocastica, eps per la condizione di arresto, maxiter=numero max
% di iterazione
% Output: approssimazione di x, iterazioni iter effettuate, vettore dei
% residui
n=size(P,1);
A=eye(n)-P';

%%%%%%%%%%%%%%%%%%%%%%%
% Metodo di Gauss-Seidel
% M=D-L, N=U
%%%%%%%%%%%%%%%%%%%%%%%
 M=tril(A);
 N=-triu(A,1);

%%%%%%%%%%%%%%%%%%%%%%%
% Metodo di Jacobi
% M=D, N=L+U
%%%%%%%%%%%%%%%%%%%%%%%
% M=diag(diag(A));
% N=-tril(A,-1)-triu(A,1);

%%%%%%%%%%%%%%%%%%%%%%%
% Metodo di Gauss-Seidel "rovesciato"
% M=D-U; N=L
%%%%%%%%%%%%%%%%%%%%%%%
% M=triu(A);
% N=-tril(A,-1);

% x=ones(n,1)/n;
x = rand(n,1); x=x/sum(x);
iter=0;
res=norm(P'*x-x,1);
vres=res;
while (res>eps & iter<maxiter)
iter=iter+1;
y=M\(N*x);
res=norm(y-x,1);
vres=[vres,res];
x=y/sum(y);
end
