clear all; 

G = [0 0 0 0 1/2 0;
    1/4 0 0 0 0 0;
    0 1/2 0 0 0 0;
    1/4 1/2 0 0 1/2 0;
    1/4 0 0 1 0 1; 
    1/4 0 0 0 0 0];

% G = [0 0 0 0 1/2 0;
%     1/3  0 0 0 0 0; 
%     0 1/2 0 0 0 0;
%     1/3 1/2  0 0 1/2 0;
%     1/3 0 1 1 0 1;
%     0 0 0 0 0 0; 
%      ]

[EV_G, ev_G] = eig(G); 

eigenvalues_G = diag(ev_G);

v0_G = [1/6 1/6 1/6 1/6 1/6 1/6]'; % initial guess for EV of G
tol=0.00001; max_iter = 10000; EV_dom = [] ; 
%%%%%%%%%%%% Power method to G %%%%%%%%%%%
for i=[1:max_iter]
    v_now = G*v0_G; % the Power method iteration v_k = G*v_{k-1}
    v_next = G*v_now; % the Power method iteration v_{k+1} = G*v_k
    if abs(v_next - v_now) < tol % when successive iterates are close
        EV_dom = v_next; 
        break;
    end
    v0_G = v_now; 
    ev_dom = dot(G*v_now, v_now)/dot(v_now,v_now); % Rayleigh quotient 
    if i==max_iter-1
        fprintf('Power method failed to converge!! OOPS!!\n')
    end
end
if isempty( null(G - eye(6)) ) ~= 1 % when: isempty( null(G - eye(6)) ) == 1, 
    % we have only the "O-vector" that satisfies the eq. Gv = v, i.e. no EV
    % of G corresponding to ev = 1
    ev_dom
    EV_dom % corresponds to soln (R1 R2 R3 R4 R5 R6)
    [temp page_rank] = sort(EV_dom,'descend');
    for j=1:6
        fprintf('Page %d \n',page_rank(j));
    end
end