function [optx,zmin] = quad_prog(G,d,A,b,toler,maxiter)
% usage: [xopt,zmin] = quad_prog(G,d,A,b,toler,maxiter)
% description: This code solves the quadratic program
%                minimize z = x'*G*x + d'*x
%                 sbj to A*x >= b, x >= 0.
% It assumes that A is mxn and that G is
% nxn positive definite, d is nx1 and b is mx1.  Arguments
% tol (termination criterion tolerance) and maxit (maximum
% number of iterations) are optional. It uses a
% gradient projection method whose code is taken from
% the program gradproj.m due to C. T. Kelley.

persistent Hp dp

if (nargin >1 && nargin < 4)
  error('quad_prog: arguments G,d,A,b are required')
elseif (nargin == 1) % this is a subfunction call
  optx = 0.5*G'*Hp*G + G'*dp;
  zmin = Hp*G + dp;
else % this is the main call 
% set stopping parms
if nargin < 5
  tol=1e-6;
else
  tol = toler;
end
if nargin < 6
  maxit=10;
else
  maxit = maxiter;
end
% set some constants and minimal bulletproofing
[m,n] = size(A);
if (size(b) ~= [m,1])
  error('quad_solv: b must be mx1, where A is mxn');
end
if (size(G) ~= [n,n])
  error('quad_solv: G must be nxn, where A is mxn');
end
if (size(d) ~= [n,1])
  error('quad_solv: d must be nx1, where A is mxn');
end
% test for consistency of equations
[optx,zmin,is_bounded,sol] = lp(ones(n+m,1),[A,-eye(m)],b);
if (~sol)
  error('quad_prog: system Ax >= b, x >= 0 has no solution');
end
% symmetrize G and account for 1/2 factor
Hp = G+G';
invG = inv(Hp);
A = [A; eye(n)]; % adjust to account for x >= 0
b = [b; zeros(n,1)];
Hp = A*invG*A';
dp = -(b + A*invG*d);
% set up for Kelley iteration
x0 = zeros(m+n,1); % initial guess
f = 'quad_prog';
low = zeros(m+n,1);
up = Inf*zeros(m+n,1);
% The gradient projection method commences here
xc=x0; ndim=length(up); kku=zeros(ndim,1); kkl=zeros(ndim,1);
for i=1:ndim
        kku(i)=up(i); kkl(i)=low(i);
  if kkl(i) > kku(i)
        error(' lower bound exceeds upper bound')
        end
end
%
% put initial iterate in feasible set
%
if norm(xc - max(kkl,min(kku,xc))) > 0
      disp(' initial iterate not feasibile ');
      xc=max(kkl,min(kku,xc));
end
alp=1.d-4;
if nargin < 6
maxit=1000; 
end
if nargin < 5
tol=1.d-6;
end
itc=1; 
[fc,gc]=feval(f,xc);
numf=1; numg=1; numh=0;
ithist=zeros(maxit,5);
xt=max(kkl,min(kku,xc-gc));
pgc=xc - max(kkl,min(kku,xc-gc));
ia=0; for i=1:ndim; if(xc(i)==kku(i) | xc(i)==kkl(i)) ia=ia+1; end; end;
ithist(1,5)=ia/ndim;
ithist(1,1)=norm(pgc); ithist(1,2) = fc; ithist(1,4)=itc-1; ithist(1,3)=0; 
while(norm(pgc) > tol & itc <= maxit)
        lambda=1;
        xt=max(kkl,min(kku,xc-lambda*gc)); ft=feval(f,xt);
        numf=numf+1;
  iarm=0; itc=itc+1;
        pl=xc - xt; fgoal=fc-(pl'*pl)*(alp/lambda);
%
%       simple line search
%
        q0=fc; qp0=-gc'*gc; qc=ft;
  while(ft > fgoal)
                lambda=lambda*.1;
    iarm=iarm+1;
    xt=max(kkl,min(kku,xc-lambda*gc));
                pl=xc-xt;
    ft=feval(f,xt); numf = numf+1;
    if(iarm > 10) 
    disp(' Armijo error in gradient projection')
                histout=ithist(1:itc,:); costdata=[numf, numg, numh];
    return; end
                fgoal=fc-(pl'*pl)*(alp/lambda);
  end
  xc=xt; [fc,gc]=feval(f,xc); numf=numf+1; numg=numg+1;
        pgc=xc-max(kkl,min(kku,xc-gc)); 
  ithist(itc,1)=norm(pgc); ithist(itc,2) = fc; 
  ithist(itc,4)=itc-1; ithist(itc,3)=iarm;
ia=0; for i=1:ndim; if(xc(i)==kku(i) | xc(i)==kkl(i)) ia=ia+1; end; end;
ithist(itc,5)=ia/ndim;
end
x=xc; 
histout=ithist(1:itc,:); costdata=[numf, numg, numh];
% projected gradient ends here
% now make corrections for fact that x is really lambda
optx = invG*(A'*x-d);
zmin = optx'*G*optx + d'*optx;
end