function [x,ok,xx]=findvertex(A,b,x0,eps)
%FINDVERTEX Find starting vertex for the simplex algorithm.
%
%[x,ok,xx]=findvertex(A,b,x0,eps)
%A   - constraint coefficient matrix
%b   - constraint right-hand side
%x0  - starting point. Must be a feasible point.
%eps - numerical equality tolerance.
%x   - vertex, or last point visited if no vertex exist.
%ok  - nonzero if x is a vertex, otherwise zero.
%xx  - point sequence.

% v1.0  2000-09-17. Niclas Borlin, niclas@cs.umu.se.
% v1.1  2001-03-26. Added sorting of residuals.

x=x0;
xx=x;
ok=0;
n=length(x);

while (1)
	% Find active set.
	res=A*x-b;
	if (any(res<-eps))
		error('Supplied point x0 is not feasible');
	end
	active=abs(res)<eps;
	Aactive=A(active,:);
	if (rank(Aactive)==n)
		% We have found a corner. Return.
		ok=1;
		return;
	end
	
	% Find a linearly independent inactive constraint.
	% Try the smallest residuals first.
	inactive=find(~active);
	[dummy,i]=sort(res(inactive));
	% Sort row numbers by smallest residual.
	inactive=inactive(i);
	found=0;
	for j=inactive'
		newAactive=[Aactive;A(j,:)];
		if (rank(newAactive)>rank(Aactive))
			% Found a linearly independent constraint.
			found=1;
			break;
		end
	end
	
	if (~found)
		% No linearly independent inactive constraint. Return with error.
		ok=0;
		return;
	end
	
	% Calculate search direction towards the new constraint.
	p=pinv(newAactive)*[zeros(size(Aactive,1),1);-1];
	p(abs(p)<eps)=0;

	% Find constraints decreasing along p.
	Ap=A*p;
	
	D=find(Ap<-eps);
	
	% For each decreasing constraint, calculate possible step length.
	steps=+inf*ones(size(b));
	steps(D)=res(D)./(-Ap(D));
	sigma=min(steps);
	
	% Take step.
	x=x+sigma*p;
	xx=[xx,x];
end
