if ~exist('OCTAVE_VERSION')
	error('Function definitions in scripts available only in Octave. Edit the source code, extract functions to separate files and rerun');
end

clear;
disp('Problem 2 from midterm 2014');

function [M,Msum] = radnmarkov(N,density)
M = sprand(N,N,density)+speye(N,N);
M = M*spdiags(1./(sum(M)'),0,N,N);
Msum = norm(sum(M)-1,Inf)
end

function [v,l] = eigmarkov(M, Q)
% solves QMQ^-1 v = l*v for lambda - dominating eigenvalue

% it is equivalent to MQ^-1v = l*Q^-1 v or My = l * y with v = Q*y
% if, of course, Q is nonsingular; 
[V, L] = eigs(M, 6, 'lm');
[L, i] = sort(abs(diag(L)));
v = V(:,i(end)); l = L(end);
if norm(M*v-l*v,Inf) > sqrt(eps)
	warning('Inaccurate eigenpair');
end
if abs(l - L(end-1)) < sqrt(eps)*abs(l)
	warning('Possible multiple eigenvalues');
end
v = Q*v;
end


tic;
for N = [987,9876]
 d = 5/N;
 M = radnmarkov(N,d); Q = sprand(M);
 % note that usually Q will be singular so our reformulation makes the problem well posed even for singular Q
 [v,l] = eigmarkov(M, Q);
 % check if Q is singular... (if someone insists, because this takes both time and memory!)
 [~,R] = lu(Q);
 R = abs(full(diag(R))); % extract diagonal only
 if min(R)/max(R) < sqrt(eps)
 	warning('Q may be singular');
 end
 % note that computing the norm of an eigenvector is a NONSENSE!
 norm(v,Inf)
end
toc
