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,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'); set(gca(),"auto_clear","off");
   else
   plot([x01(1);xo(1)],[x01(2);xo(2)],'-ro'); set(gca(),"auto_clear","off");    
   end
   
   if ip > 0,
     disp('Pause: hit any key to continue...'); halt;
   end
 end
 erro=norm(dp);
end

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