m=100; n1=50;
A = [randn(m,n1), eye(m)];
b = rand(m,1);
c = [randn(n1,1); zeros(m,1)];
x0 = [zeros(n1,1);b];


options = optimset('LargeScale','off');
[m,n] = size(A);

disp('Running linprog')
[xlp,fval,exitflag,output,lamba] = ...
          linprog(c,[],[],A,b,zeros(n,1),inf*ones(n,1),x0,options);

disp(fval)
disp('Running lpfeasdir')
xfd = lpfeasdir(A,b,c,x0);

disp('Relative error in the solution vector computed by lpfeasdir')
norm(xlp-xfd)/norm(xlp)