function [xo,Ot,nS,lambda]=IPquad(H,c,x0,A,b,tol,mxit,ip)
%   Solves the linear programming problem
%      min S(x)=c'x + 1/2x'Hx  subject to: A x = b, x >= 0
%   using the Interior Point algorithm.
%
%   [xo,Ot,nS,lambda]=IPquad(H,c,x0,A,b,tol,mxit,ip)
%
%   H: hessian matrix of the objective function
%   c: coefficients of the objective function
%   x0: initial point (must be feasible)
%   A: constraint matrix
%   b: constraint vector
%   tol: tolerance (default = 1e-4)
%   mxit: maximum number of iterations (default = 100)
%   ip: (0): no plot (default), (>0) plot figure ip with pause, (<0) plot figure ip
%   xo: optimal point
%   Ot: optimal value of S
%   nS: number of objective function evaluations
%   lambda: Lagrange multipliers

%   Copyright (c) 2015 by LADES-PEQ-COPPE-UFRJ
%   $Revision: 1.0 $  $Date: 2015/08/27 09:45:15 $
%   Argimiro R. Secchi (arge@peq.coppe.ufrj.br)

 if nargin < 5,
   error('IPquad requires five input arguments');
 end
 if nargin < 6 | isempty(tol),
   tol=1e-4;
 end
 if nargin < 7 | isempty(mxit),
   mxit=100;
 end
 if nargin < 8 | isempty(mxit),
   ip=0;
 end

[m,n]=size(A);
erro=norm(c+H*x0);
I=eye(n);
um=ones(n,1);
%gama=0.9995;
gama=0.8;

nS=0;
Q=I;
x1=x0;
x0=um;
c0=c;
H0=H;

while erro > tol & nS < mxit,
 nS=nS+1;
 D=diag(x1);
 Q=Q*D;
 A=A*D;
 c=D*c;
 H=D*H*D;
 iH=inv(H);
 d=-c-H*x0;

 B=inv(A*iH*A')*A*iH;
 P=iH*(I-A'*B);
 dp=P*d;
 lambda=B*d;
 z=-x0./dp;
 j=find(z>0);
 alfa=gama*min(z(j));
 x1=x0+alfa*dp;
 xo=Q*x1
 x01=Q*x0;
 if ip,
   
   if n > 2,
   plot3([x01(1);xo(1)],[x01(2);xo(2)],[x01(3);xo(3)],'-ro'); hold on
   else
   plot([x01(1);xo(1)],[x01(2);xo(2)],'-ro'); hold on    
   end
   
   if ip > 0,
     disp('Pause: hit any key to continue...'); pause;
   end
 end
 erro=norm(dp);
end

Ot=c0'*xo+0.5*xo'*H0*xo;
