%  [x,y,obhis]=solqp(Q,A,b,c,(toler,beta))
%
%  See the User Guide at the end of this file.
%
%
%  Start Phase 1: try to find an interior feasible point.
%
%
global x y obhis;
retval = 1;

 if exist('toler') ~= 1 
   toler=1.e-5; 
 end
 if exist('beta') ~= 1 
   beta=0.8;    
 end
 if exist('alpha') ~= 1 
   alpha=0.95;    
 end
 [m,n] = size(A);
 disp('Search for a feasible point:')
 a=b-A*ones(n,1);
 x=ones(n+1,1);
 z=0;
 ob=x(n+1); 
 obhis=[ob];
 gap = ob - z;
 while gap >= toler,
####### start   spphase1;
%   spphase1 
%
% solve the scaled least squares against two vectors 
%
% [i,j,sa]=find(A);
% AX = sparse(i,j,sa.*x(j),m,n);
% clear i j sa;
% AX'\[ones(n,1) x.*c1];
% 
 dx = ones(n,1)./x(1:n);
 DD = sparse(1:n,1:n,dx.*dx,n,n);
 [DD A';A sparse(m,m)]\[dx sparse(n,1); sparse(m,1) a];
%
 y1=ans(n+1:n+m,1);
 y2=ans(n+1:n+m,2);
 clear dx ans DD;
 w1=(1/ob - a'*y1)/(1/ob^2 - a'*y2);
 w2=1/(1/ob^2 - a'*y2);
 y1=y1-w1*y2;
 y2=-w2*y2;
%
 w1=b'*y1;
 w2=b'*y2;
 y1=y1/(1+w1);
 y2=y2-w2*y1;
 u=[x(1:n).*(-y2'*A)';x(n+1)*(1-y2'*a);w2/(1+w1)];
 v=[x(1:n).*(y1'*A)' ;x(n+1)*(y1'*a)  ; 1/(1+w1)];
%
%  update the dual and the objective lower bound
%
 if min(u-z*v)>=0, 
   y = y2+z*y1; 
   z=b'*y; 
 end;
 clear y1 y2 w1 w2;
%
%  find the descent direction
%  
 u=u-z*v-((ob-z)/(n+2))*ones(n+2,1);
 nora=max(u);
%
%  update the solution along the descent direction
%
 if nora==u(n+1),
   alpha=1.;
 end;
 v=ones(n+2,1)-(alpha/nora)*u;
 x=(x.*v(1:n+1))/v(n+2);
 clear u v
%
%return
%
%  This is the Phase 1 procedure called by SPSOLQP. 
#######################################################
   ob=x(n+1);
   obhis=[obhis ob];
   gap = ob - z,
   if z > 0,
     gap = -1;
     disp('The system has no feasible solution.'),
     return
   end;
 end;
 clear a
% 
% Start Phase 2
%
%
 disp('Search for an optimal solution:');
 alpha = 0.99;
 x=x(1:n);
 comp=rand(n,1);
 [speye(n) A';A sparse(m,m)]\[comp;sparse(m,1)];
 comp=ans(1:n);
 clear ans;
 nora=min(comp./x);
 if nora < 0,
   nora = -.01/nora;
 else
   nora = max(comp./x);
   if nora == 0,
     disp('The problem has a unique feasible point');
     return
   end;
   nora = .01/nora;
 end;
 x = x + nora*comp;
 obvalue=x'*(Q*x)/2+c'*x,
 obhis=[obvalue];
 lower =-inf;
 zhis=[lower];
 gap=1;
 lamda=max([1 abs(obvalue)/sqrt(sqrt(n))]);
 iter=0;
 while gap >= toler,
   iter=iter+1;
############################### start   spphase2;
%   spphase2 
%
 lamda=(1.-beta)*lamda;
% if gap <= 5*toler;
%   lamda = lamda/2;
% end;
 go=0;
 dx = ones(n,1)./x;
%
%  Repeatly solve an ellipsoid constrained QP problem by solving a linear
%  system equation until find a positive solution.
%
 while go <= 0,
   DD = sparse(1:n,1:n,(lamda*dx).*dx,n,n);
%
%   u=[Q+DD A';A sparse(m,m)]\[-(Q*x+c)+(lamda/n)*dx;sparse(m,1)];
%
   u=[Q+DD A';A sparse(m,m)]\[-(Q*x+c);sparse(m,1)];
   xx=x+u(1:n);
   go=min(xx);
   if go > 0,
     ob=xx'*Q*xx/2+c'*xx;
     go = min([go obvalue-ob+eps]);
   end;
   lamda=2*lamda;
   if lamda >= (1+abs(obvalue))/toler,
     disp('The problem seems unbounded.');
     return
   end; 
 end;
%
 y=-u(n+1:n+m);
 u=u(1:n);
 nora = min(u./x);
 if nora < 0,
   nora=-alpha/nora;
 elseif nora == 0,
   nora=alpha;
 else
   nora=inf;
 end
%
 w1 = u'*Q*u;
 w2 = -u'*(Q*x+c);
 if w1 > 0,
  nora=min([w2/w1,nora]);
 end;
 if nora == inf,
  ob = -inf;
 else
   x =x+nora*u;
   ob=x'*Q*x/2+c'*x;
 end;
 clear u dx xx DD w1 w2
%
%  This is the Phase 2 procedure called by SPSOLQP. 
####################################################
   if ob == -inf,
     gap = 0;
     disp('The problem is unbounded.');
     return
   else 
     obhis=[obhis ob]; 
     comp=Q*x+c-A'*y; 
     if min(comp)>=0, 
       zhis(iter+1)=ob-x'*comp;
       lower=zhis(iter+1),
       gap=(ob-lower)/(1+abs(ob));
       obvalue=ob,
     else
       zhis(iter+1)=zhis(iter);     
       lower=zhis(iter+1);
       gap=(obvalue-ob)/(1+abs(ob));
       obvalue=ob,       
     end;
   end;
 end;
 disp('A (local) optimal solution is found.');
 return
%
%  This program solves quadratic program in standard form:
%
%     minimize    0.5*(x'*Q*x)+c'*x
%     subject to  A*x=b, x>=0.
%
%  Input 
%      Q: Sparse symmetric objective matrix.
%      A: Sparse constraint left-hand matrix
%      b: constraint right-hand column vector
%      c: objective column vector
%      toler: relative stopping tolerance: the objective value close to 
%             the local optimal one in the range of tolerance. 
%             Default value: 1.e-5.
%      beta : step size: 0 < beta < 1. Default value: .95.
%
%  Output
%     x: (local) optimal solution
%     y: optimal dual solution (Lagrangien multiplier)
%     obhis : objective value history vs iterations
%
%  Subroutines called : spphase1 and spphase2
%
%  To run the program, just type: spsolqp
%
%    This program is the implementation of the interior ellipsoidal trust
%  region and barrier function algorithm with dual solution updating
%  technique in the standard QP form. Two phases are used: the first uses 
%  SPPHASE1 to find an interior feasible point and the second uses SPPHASE2
%  to find a local optimal solution.
%
%  Technical Reference
%  
%  Y. Ye, "An extension of Karmarkar's algorithm and the trust region method
%         for convex quadratic programming," in Progress in Mathematical
%         Programming (N. Megiddo ed.), Springer-Verlag, NY (1989) 49-63.
%
%  Y. Ye, "On affine-scaling algorithm for nonconvex quadratic programming,"
%         Math. Programming 56 (1992) 285-300.
%
%  Comment: Each iteration we solve a linear KKT system like
%
%  ( Q+mu X^{-2}   A^T )(dx) = c'
%  (     A          0  )(dy) = 0
%
%  where X = diag(x)  which is a positive diagonal matrix.