function [xout, lambdaout, howout, success] = pdipm6_qp(H,c,A,b,VLB,VUB,x0,N,verbose,opt)
% PDIPM6_QP  Quadratic program solver.
%   Wrapper function for a primal-dual interior point QP solver
%   using the calling syntax of qp() from version 1 of the
%   Optimization Toolbox, with the exception that verbose == 0
%   means no output. 
%
%   From the Optimization Toolbox v.1 docs ...
%     X=QP(H,f,A,b) solves the quadratic programming problem:
%  
%              min 0.5*x'Hx + f'x   subject to:  Ax <= b 
%               x    
%  
%         X=QP(H,f,A,b,VLB,VUB) defines a set of lower and upper
%         bounds on the design variables, X, so that the solution  
%         is always in the range VLB <= X <= VUB.
%  
%         X=QP(H,f,A,b,VLB,VUB,X0) sets the initial starting point to X0.
%  
%         X=QP(H,f,A,b,VLB,VUB,X0,N) indicates that the first N constraints 
%         defined by A and b are equality constraints.
%  
%         X=QP(H,f,A,b,VLB,VUB,X0,N,DISPLAY) controls the level of warning
%         messages displayed.  Warning messages can be turned off with
%         DISPLAY = -1.
%  
%         [x,LAMBDA]=QP(H,f,A,b) returns the set of Lagrangian multipliers,
%         LAMBDA, at the solution.
%  
%         [X,LAMBDA,HOW] = QP(H,f,A,b) also returns a string HOW that 
%         indicates error conditions at the final iteration.
%  
%         QP produces warning messages when the solution is either unbounded
%         or infeasible. 

%   MATPOWER
%   $Id: pdipm6_qp.m,v 1.6 2009/11/04 15:58:43 ray Exp $
%   by Ray Zimmerman, PSERC Cornell
%   Copyright (c) 2009 by Power System Engineering Research Center (PSERC)
%   See http://www.pserc.cornell.edu/matpower/ for more info.

%% default args
nx = length(c);
if nargin < 9
    verbose = 0;
    if nargin < 8
        N = 0;
        if nargin < 7
            x0 = [];
            if nargin < 6
                VUB = [];
                if nargin < 5
                    VLB = [];
                end
            end
        end
    end
end
if isempty(H)
    H = sparse(nx, nx);
end
if isempty(VUB)
    VUB = Inf * ones(nx, 1);
end
if isempty(VLB)
    VLB = -Inf * ones(nx, 1);
end
if isempty(x0)
    x0 = zeros(nx, 1);
    k = find(VUB < 1e10 & VLB > -1e10);
    x0(k) = (VUB(k) + VLB(k)) / 2;
    k = find(VUB < 1e10 & VLB <= -1e10);
    x0(k) = VUB(k) - 1;
    k = find(VUB >= 1e10 & VLB > -1e10);
    x0(k) = VLB(k) + 1;
end

%% set options and function args
opt.verbose = verbose;
if isfield(opt, 'cost_mult')
    cost_mult = opt.cost_mult;
else
    cost_mult = 1;
end
ipm_f = @qp_f;
ipm_gh = @qp_gh;
ipm_hess = @qp_hessian;

%% set up lower bound
l = -Inf * ones(size(b));
l(1:N) = b(1:N);

%% run it
[xout, f, info, Output, Lambda] = ...
  pdipm6(ipm_f, ipm_gh, ipm_hess, x0, VLB, VUB, A, l, b, opt, H, c, opt.cost_mult);

success = (info > 0);
if success
    howout = 'success';
else
    howout = 'failure';
end
lambdaout = [   -Lambda.mu_l + Lambda.mu_u;
                Lambda.lower;
                Lambda.upper    ];


%%-----  objective function  -----
function [f, df, d2f] = qp_f(x, H, c, varargin)
f = 0.5 * x' * H * x + c' * x;
if nargout > 1
    df = H * x + c;
    if nargin > 2
        d2f = H;
    end
end

%%-----  constraints  -----
function [g, h, dg, dh] = qp_gh(x, varargin)
g = []; h = []; dg = []; dh = [];

%%-----  hessian  -----
function Lxx = qp_hessian(x, lambda, H, c, cost_mult)
Lxx = H * cost_mult;
