function [val,primals,duals] = iqp(H,c,A,b,verbose)

% [val,primals,duals] = iqp(H,c,A,b,verbose)
%
% A simple interior point algorithm for solving the quadratic program
%   maximize - x'Hx/2 + c'x st Ax+b >= 0
% when H is positive semidefinite.  Use H=[] for a linear program.
% Use verbose=[] to turn off messages.
% Detection of infeasibility/unboundedness is not very good: we just
% return NaN if it appears that the optimal solution to either the
% primal or the dual is very large.
% General robustness could also be improved.

% Copyright 1997 Geoff Gordon

% a few parameters
epsilon = 1e-8;    % accuracy requirement for solution
toobig = 1e+8;     % give up if solution gets too large
backoff = .99995;  % don't go quite all the way to boundary
maxiter = 75;      % stop after this many iterations
gamma = .5;        % how fast to try to decrease barrier
decay = .5;        %  -but effective gamma can be 1-(1-gamma)/decay
mingamma = 1e-8;   %  -but bounded below by mingamma
canreorder = 1;    % set to zero to prevent reordering

[m,n] = size(A);

% default args
if (nargin < 5) verbose = 1; end
if (H == []) H = sparse(n,n); end

% symmetric part of H is all that matters
H = (H + H') / 2;

% reorder for sparsity
if (canreorder & issparse(A) & issparse(H))
  order = symmmd(H + A' * A);
  H = H(order,order);
  A = A(:,order);
  c = c(order);
else
  order = [];
end

% simple stupid initialization
primals = zeros(n,1);
duals = ones(m,1);
slacks = ones(m,1);
complementarity = (slacks' * duals) / m;

% do until complementarity is small
iter = 0; pstep = 1; dstep = 1; step = 0; barrier = 100;
while ((complementarity > epsilon) & (iter < maxiter))

  % target barrier is a fraction of complementarity determined
  % by how recently we've run against the edge of feasibility
  % (exact code here is arbitrary, anything reasonable should work)
  step = decay * min(pstep, dstep) * (step + 1/decay);
  effgamma = max(mingamma, (1 - step * (1 - gamma)));
  barrier = min(barrier, complementarity * effgamma);

  % compute rhs for Newton equations
  invslacks = 1 ./ slacks;
  d = spdiags(duals .* invslacks, 0, m, m);
  rhs = A' * (barrier .* invslacks + duals - d * b) + c;

  % solve Newton equations -- matlab is smart enough to use
  % Cholesky decomposition here
  newprimals = (H + (A' * d * A)) \ rhs;

  % from newprimals, compute newslacks
  newslacks = A * newprimals + b;

  % from newprimals and newslacks, compute newduals
  newduals = barrier .* invslacks - d * (newslacks - slacks);

  % compute primal and dual stepsizes: longest possible w/o hitting zero
  pstep = 1 ./ max(max(slacks - newslacks, slacks) ./ slacks);
  dstep = 1 ./ max(max(duals - newduals, duals) ./ duals);
  if (min(newslacks) <= 0) pstep = backoff * pstep; end
  if (min(newduals) <= 0) dstep = backoff * dstep; end

  % do the update
  primals = primals + pstep * (newprimals - primals);
  slacks = slacks + pstep * (newslacks - slacks);
  duals = duals + dstep * (newduals - duals);
  complementarity = (slacks' * duals) / m;
  iter = iter + 1;

  % give up if numbers get too big
  if max(max(slacks),max(duals)) > toobig
    val = NaN;
    if (order ~= []) primals(order) = primals; end
    if (verbose ~= [])
      disp('Warning: program appears to be ill-posed');
    end
    return;
  end

  % report progress
  if (verbose ~= [])
    fprintf('iter %d, pstep %g, dstep %g, target %g, actual %g\n', ...
      iter, pstep, dstep, barrier, complementarity);
  end

end

val = primals' * c - primals' * H * primals / 2;
if (order ~= []) primals(order) = primals; end