%  [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.
%
%
 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,
   spphase1;
   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.9;
 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;
   spphase2;
   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.