%IPM1MEX  Interior point method for aX>0 and no equality constraints 
%
% This is the MATLAB equivalent implementation of libmeq/ipm1.c.
% A more detailed help is available in IPM1MEX.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.
function [x,z,stat] = ipm1mexm(uH,c,a,x,z,tol,nit)
 stat = 0;
 n  = size(c,1);
 H = zeros(n);
 ku = triu(true(n,n));
 H(ku) = uH;
 H = H + triu(H,1)';
 for kit=1:nit
  stat = ~any(x==0);
  if ~stat, return;end
  d = z./x;
  U = H + diag(d);
  [U,p] = chol(U);
  if p>0, stat = false; return, end
  dx = H*x + c;
  dx = U\(U'\dx);
  dz = z - d.*dx;
  alpha = 1;
  k = (a*dx > 0); alpha = min([alpha ; x(k)./dx(k)]);
  k = (a*dz > 0); alpha = min([alpha ; z(k)./dz(k)]);
  tau = sum((x - alpha*dx).*(z - alpha*dz)); % Since r0 and r1 can be both very small,
  tau = (tau/sum(x.*z))^2*(tau/n);           % this is more precise and less prone to NaNs
  e = (dx.*dz-tau)./x;
  dz = U\(U'\e);
  dx = dx + dz;
  dz = e + z - d.*dx;
  alpha = 1;
  k = (a*dx > 0); alpha = min([alpha ; x(k)./dx(k)]);
  k = (a*dz > 0); alpha = min([alpha ; z(k)./dz(k)]);
  x = x - 0.99*alpha*dx;
  z = z - 0.99*alpha*dz;
  stat = max(abs(dx)) < tol;
  if stat, break, end
end
