function x=prdirect(alpha,Pcc,v)
if abs(alpha-1)<eps
    [V,D] = eigs(Pcc',3,'LR');
    if abs(D(2,2)-D(1,1)) < 1000*eps
        error('limit vector not computed');
    end
    % TODO handle multiple components in 
    % a symmetric walk
    x=V(:,1);
    x = x./sum(x);
    return
end

n = size(Pcc,1);

if v==1,
    v = ones(n,1)./n;
else
    v = v./sum(v);
end
    
A = speye(n) - alpha*Pcc';
b = (1-alpha)*v;
x = A\b;
