function [xopt,lambda,how] = qpdantz(H,f,A,b,xmin,maxiter)
%QPDANTZ command will be removed in a future release.  If you have
%   Optimization Toolbox, replace "qpdantz" with "quadprog" command. 
%
% QPDANTZ solves a quadratic program using an active set method.
%
%   [x,lambda,status]=QPDANTZ(H,f,A,b,xmin,maxiter) solves the quadratic
%   programming problem:
%
%     min 0.5*x'Hx + f'x   subject to:  Ax <= b, x >= xmin
%      x
%
%   where:
%        x = the optimal solution vector.
%   lambda = the dual variables (Lagrange multipliers) at the solution.
%   status = 'feasible' if x satisfies the constraints
%          = 'infeasible' if no feasible x could be found
%          = 'unreliable' if maximum number of iterations is exceeded.

%   Author: A. Bemporad
%   Copyright 1986-2011 The MathWorks, Inc.

narginchk(4,6);

if nargin<5 || isempty(xmin)
    xmin=-1e3*ones(size(f(:)));
    ctrlMsgUtils.warning('MPC:computation:AssumeXMin');    
end

if nargin<6,
    maxiter=200;
end
if maxiter<1,
    ctrlMsgUtils.error('MPC:computation:InvalidMAXITER','qpdantz');
end

mnu=length(f);
nc=length(b);

if ~isequal(size(H),[mnu mnu])
    ctrlMsgUtils.error('MPC:computation:InvalidHessianOrLinCost','qpdantz');            
end

if size(A,2)~=mnu,
    ctrlMsgUtils.error('MPC:computation:InvalidConsMatrixOrLinCost','qpdantz');                
end

if size(A,1)~=nc,
    ctrlMsgUtils.error('MPC:computation:InvalidConsMatrixOrRHS','qpdantz');                
end

% H must be symmetric. Otherwise set H=(H+H')/2
if norm(H-H') > eps
    ctrlMsgUtils.warning('MPC:computation:HessianNonsymmetric');    
    H = (H+H')*0.5;
end

% Force double precision
xmin=double(xmin);
H=double(H);
f=double(f);
A=double(A);
b=double(b);
maxiter=double(maxiter);

% This is a constant term that adds to the initial basis in each QP.
sw=ctrlMsgUtils.SuspendWarnings;
a=-H*xmin(:);    
% Inverted hessian
H=H\eye(mnu);
delete(sw);
rhsc=b(:)-A*xmin(:);
rhsa=a-f(:);
TAB=[-H H*A';A*H -A*H*A'];
basisi=[H*rhsa;
    rhsc-A*H*rhsa];
ibi=-(1:mnu+nc)';
ili=-ibi;
[basis,ib,il,iter] = qpdantzg(TAB,basisi,ibi,ili,maxiter);

if iter>maxiter,
    how='unreliable';
elseif iter <=0
    how='infeasible';    
else
    how='feasible';
end

xopt=zeros(mnu,1);

for j=1:mnu
    if il(j) <= 0
        xopt(j)=xmin(j);
    else
        xopt(j)=basis(il(j))+xmin(j);
    end
end

% Add additonal check to avoild numerical issues in QPDANTZG
if any((A*xopt-b)>sqrt(eps)*abs(b))
    how='infeasible';    
end

if nargout>=2,
    lambda=zeros(nc,1);
    for j=1:nc,
        if ib(mnu+j) <= 0
            lambda(j)=0;
        else
            lambda(j)=basis(ib(mnu+j));
        end
    end
end

ctrlMsgUtils.warning('MPC:utility:ObsoleteQPdantz');    