function thetaVals = MakeSingular(A)
% A is an n-by-n nonsingular stochastic matrix.
% thetaVals is a column n-vector with the property that
% if theta = thetaVals(k), then A(theta,k) is singular.

[L,U,p] = lu(A,'vector');
[n,n] = size(A);
for k=1:n
% Determine thetaVals(k)...
% Define the kth column of A(theta,k) by

% w = [mu*A(1:k-1,k); theta*A(k,k); mu*A(k+1:n,k)]

% where mu satisfies

% mu*sum(A(1:k-1,k)) + theta*A(k,k) + mu*sum(A(k+1:n,k)) = 1
% i.e.,
% mu = (1 - theta*A(k,k))/(1 - A(k,k))

% Note that A(theta,k) = A + u*ek' where u = w - A(:,k) and ek = I(:,k)
%
% A(theta,k) will be singular if

% -1 = ek'*inv(A)u = ek'*inv(A)*(w - A(:,k)) = ek'inv(A)*w -1

% Thus, we must have

% 0 = ek'*inv(A)*w
% = ek'*inv(A)*( mu*A(:,k) + (theta*A(k,k) - mu*A(k,k))*ek )
% = mu + (theta - mu)A(k,k)*v(k)

% where the vector v solves Av = ek. Solving this for theta we get

% theta = (1-A(k,k)*v(k))/(A(k,k)*(1 - v(k)))

% So all we need is the diagonal of inv(A), i.e., the kth component
% of v = U\(L\(P*ek)). Can cut some corners both in the L-solve and the
% U-solve.
% Notice that P*ek = I(:,q(k)) where
q(p) = 1:n;
% Thus, the rhs in Ly = P*ek is a vector with zeros in components 1:q(k)-1
y = [zeros(q(k)-1,1);L(q(k):n,q(k):n)\[1;zeros(n-q(k),1)]];
% Now we want just the k-th component of the solution to Uv = y.
v(k:n) = U(k:n,k:n)\y(k:n);
theta = (1-A(k,k)*v(k))/(A(k,k)*(1 -v(k)));
thetaVals(k) = theta;
% Thus, v(k) costs 2(n-k)^2 flops instead of 2n^2 flops
end

% So 2/3n^3 for LU plus 2/3n^3 for all the v(k) = 4/3 n^3 total
% If you just do v = U\(L\(P*ek)) its 8/3 n^3 total.

% -1 if you ignore the special features of the triangular solves
% -2 if you just did B =inv(A) then v(k) = B(k,k) in the loop