function x_new=linprogip(n,eps)
% linear programming solved with interior point method
% min c^T x under the constraints A x = b and x>=0
% the KKT conditions are
%   c-A^T mu -lambda = 0
%   A x -b = 0
%   x_i lambda_i = 0 i=1,...,n
%   x>=0 and lambda>=0
% Newton method is used for the 3 first conditions
p=1;
A=ones(p,n);
b=sparse(p,1);
b(1)=1;
c=ones(n,1);
c(1)=-1;
x_old=max(eps,zeros(n,1));
lambda_old=max(eps,zeros(n,1));
mu_old=zeros(p,1);
for iter=1:10
  mat1 = horzcat(sparse(n,n),-speye(n,n),???);
  mat2 = horzcat(A,sparse(p,n),sparse(p,p));
  mat3 = horzcat(sparse(1:n,1:n,lambda_old,n,n),sparse(1:n,1:n,???,n,n),sparse(n,p));
  mat  = vertcat(mat1,mat2,mat3);
  rhs1 = c-A'*mu_old-???;
  rhs2 = A*x_old-b;
  rhs3 = sparse(n,1);
  for i=1:n
   rhs3(i)=x_old(i)*???;
  end
  rhs = vertcat(rhs1,rhs2,rhs3);
  sol=mat\rhs;
  x_new=max(eps,x_old-sol(1:n));
  lambda_new=max(eps,lambda_old-sol(n+1:???));
  mu_new=mu_old-sol(???:2*n+p);
  discrep=norm(x_new-x_old)/norm(x_new);
  printf ("iter: %d Discrepancy: %f \n",iter,discrep);
  %fprintf ("iter: %d Discrepancy: %f \n",iter,discrep); %matlab
  x_old=x_new;
  lambda_old=lambda_new;
  mu_old=mu_new;
  if (discrep<0.001)
   break
  end
end
end
