profile on
clear all
%%% Generate random data for an LP in SEF
%%% with strictly feasible  starting points.
m=10;n=19;    % dimensions of problem
toler=1e-8;
disp(['dimension m,n of problem is: ',num2str([m n])])
disp(['tolerance for duality gap is: ',num2str(toler)])
A=[rand(m,n-m) eye(m)];b=10*(rand(m,1)+1);c=rand(n,1)+1;
x=ones(n,1);s=x;y=ones(m,1);
b=A*x;c=A'*y-s;   % ensure feasible start
iters=0;
while sum(x.*s) > toler,
   iters=iters+1;
   %% Find affine scaling search direction
   [dx,dy,ds]=intpteff(A,b,c,x,y,s);
   [dx2,dy2,ds2]=intpt(A,b,c,x,y,s);
   dir2=[dx2;dy2;ds2];
   dir=[dx;dy;ds];
   inds=find(dx<0);
   alphamin=min(-x(inds)./dx(inds));
   inds=find(ds<0);
   alphamin=min(alphamin,min(-s(inds)./ds(inds)));
   %%% Backtrack if necessary
   alpha=1;
   if alphamin < 1,
     while min( min(x+alpha*dx),min(s+alpha*ds) ) <= 0,
        alpha = .64*alphamin;
     end
   end
   normerror=norm(dir-dir2);
   disp(['iteration; gap, step , direction error :  ',...
        num2str([iters sum(x.*s) alpha normerror])]);
   %% update
   x=x+alpha*dx;s=s+alpha*ds; y=y+alpha*dy;
end
disp(' ')
disp(['number of iterations is: ',num2str(iters)]);
disp(' ')
disp(['optimal x s to duality gap: ',num2str(sum(x.*s))])
%[x s]
profile report
