%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%generate a random LP first   min c'x s.t.   Ax=b  x>=0
sigmaa=.1;    % centering parameter
toler=1e-6;   % stopping tolerance
nump=10;      % number of problems to solve
%%%%%%
for ii=1:nump,             %%%   solve nump problems
disp(' ')
m=45;               % constant number of constraints
n=ii*60;            % number of variables increasing with ii
%m=ii*10;n=ii*40;
disp(['starting new problem with m,n  ',num2str([m n])]);
A=randn(m,n) +[.1*eye(m)  zeros(m,n-m)];   % full row rank A
x=10*rand(n,1);
b=A*x;    % ensure primal feas.
z=10*rand(n,1);
c=A'*randn(m,1)+z;   % ensure dual feas.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%solve using p-d i-p method
x=ones(n,1);z=x;   y=randn(m,1); % starting points
muu=z'*x/n;         %initial log-barrier parameter estimate
rd= A'*y+z-c;       % dual feas. residual
rp=A*x-b;           % primal feas. residual
noiters=0;
fullstepp=0;        % indicator for taking a full Newton primal Step
fullstepd=0;        % indicator for taking a full Newton dual Step
while (muu +norm(rd)+norm(rp))/(norm(A)+norm(b)+norm(c)) > toler,
noiters=noiters+1;
D2=diag(x./z);      % D^2=Z^{-1}X
R=chol(A*D2*A');
rc=z.*x-muu*sigmaa*ones(n,1);   % compl. slack. residual with centering
%%% solve for search direction   [dx dy dz]
dy=(R')\(-rp-A*D2*rd+A*(rc./z));
dy=R\dy;
dz=-rd-A'*dy;
dx=-rc./z-D2*dz;
inds=find(dx<0);
xstep=.97*min(-x(inds)./dx(inds));
if fullstepp==0,
   if xstep>1,
      xstep=1;
      disp(['taking a full Newton primal step at iter  ',num2str(noiters)]);
      fullstepp=1;
   end
end
inds=find(dz<0);
zstep=.97*min(-z(inds)./dz(inds));
if fullstepd==0,
   if zstep>1,
      zstep=1;
      disp(['taking a full Newton dual step at iter  ',num2str(noiters)]);
      fullstepd=1;
   end
end
%%%%%UPDATES
x=x+xstep*dx;
z=z+zstep*dz;
y=y+zstep*dy;
muu=z'*x/n;
rd= A'*y+z-c;
rp=A*x-b;
rc=z.*x-muu*sigmaa*ones(n,1);
end
disp(['number of iterations is  ',num2str(noiters)]);
disp(['final duality gap +residuals is   ',num2str(...
               (muu +norm(rd)+norm(rp))/(norm(A)+norm(b)+norm(c)))]);
disp(['final residuals is   ',num2str(...
               (norm(rd)+norm(rp))/(norm(A)+norm(b)+norm(c)))]);
end
