```function [flag,x,y,h]=verquadprog(A,b,c,D)
%
%    This is an INTLAB file. It requires to have INTLAB installed under
%    MATLAB to function properly.
%
%    For real data A, b, c, D, with D positive definite (A, D full or sparse),
%    either computes a verified optimal solution x, verified Lagrangian
%    multipliers y and a verified optimal value h of the quadratic
%    programming problem
%        min 0.5*x'*D*x+c'*x   subject to   A*x>=b, x>=0,
%    or verifies (in)feasibility, or yields no verified result. The
%    respective outcome is always described verbally in the variable flag.
%
%    Possible outputs:
%
%    flag='verified optimum   ' : x is verified to enclose an optimal solution,
%                                 y is verified to enclose a vector of Lagrangian multipliers,
%                                 h is verified to enclose the optimal value,
%    flag='verified feasible  ' : x is verified to enclose a feasible solution
%                                   (optimality nor unboundedness could be verified),
%                                 y, h are NaN's,
%    flag='verified infeasible' : y is a Farkas vector verified to satisfy y>=0, A'*y>=0, b'*y<0
%                                   (whose existence proves infeasibility),
%                                 x, h are NaN's,
%    flag='verified not PD    ' : D verified not to be positive definite,
%                                   x is a vector verified to satisfy x'*D*x<=0,
%                                   y, h are NaN's,
%    flag='sizes do not match ' : x, y, h are NaN's
%                                   (sizes of A, b, c, D do not match, or D is nonsymmetric),
%    flag='no verified result ' : x, y, h are NaN's
%                                   (no verified result could be found).
%
%    Complexity: The algorithm solves one quadratic programming problem and
%    at most three linear programming problems (independently of the size
%    of the original problem), and employs two verification procedures
%    which both run approximately in O(m^3) time, where m=size(A,1).
%

%
%    This work was supported by the Czech Republic National Research
%    Program "Information Society", project 1ET400300415.
%
%    WARRANTY
%
%    Because the program is licensed free of charge, there is
%    no warranty for the program, to the extent permitted by applicable
%    law. Except when otherwise stated in writing the copyright holder
%    and/or other parties provide the program "as is" without warranty
%    of any kind, either expressed or implied, including, but not
%    limited to, the implied warranties of merchantability and fitness
%    for a particular purpose. The entire risk as to the quality and
%    performance of the program is with you. Should the program prove
%    defective, you assume the cost of all necessary servicing, repair
%    or correction.
%
gr=getround;
setround(0);
warning('off');
% checking data
b=b(:); c=c(:); [m,n]=size(A);
flag='no verified result ';
x=repmat(infsup(NaN,NaN),n,1); y=repmat(infsup(NaN,NaN),m,1); h=infsup(NaN,NaN);
if ~(length(b)==m&&length(c)==n&&isequal(size(D),[n,n])&&isequal(D,D'))
flag='sizes do not match ';
setround(gr); return
end
if isintval(A)||isintval(b)||isintval(c)||isintval(D)           % data not thin
setround(gr); return
end
% verifying positive definiteness
[pd,xpd]=verpd(D);
if pd==0                                                        % D verified not positive definite
flag='verified not PD    '; x=xpd;
setround(gr); return
end
if pd==-1||pd==-2
setround(gr); return
end
Ao=[A -eye(m,m)];
% verifying infeasibility
yi=verinfeas(Ao,b);
if ~isnan(yi.inf(1))                                            % verified Farkas vector found
y=yi; flag='verified infeasible';
setround(gr); return
end
% verifying feasibility
xf=veropt(Ao,b,ones(n+m,1));
flag='no verified result ';
setround(gr); return
end
xf=xf(1:n);
% computing unverified optimum
options=optimset(options,'Display','off');
% solves min 0.5*x'*D*x + c'*x  subject to  A*x>=b, x>=0
if exitflag<=0
x=xf; flag='verified feasible  ';                           % previous feasible solution outputed
setround(gr); return
end
p=lambda.ineqlin;
% verifying optimality
M=[D -A'; A zeros(m,m)];                                        % data of the respective LCP
q=[c' -b']';
z=[x' p']';                                                     % y, z computed in this way satisfy (in infinite precision)
y=M*z+q;                                                        % y=M*z+q, y>=0, z>=0, y'*z=0
x=y-z;
I=eye(m+n,m+n);
I=infsup(I,I);                                                  % x satisfies (in infinite precision)
M=infsup(M,M);                                                  % (I+M)*x+(I-M)*abs(x)=2*q
A=I+M; B=I-M; b=infsup(2,2)*infsup(q,q);
xx=absvalverifn(A,B,b,x);
if isnan(xx.inf(1))
x=xf; flag='verified feasible  ';                           % previous feasible solution outputed
setround(gr); return                                        % no verified solution of LCP found
end
% extracting x^-=max(-x,0)                                      % xx is verified solution
for i=1:m+n
if 0<xx.inf(i)
xx(i)=infsup(0,0);
end
if xx.sup(i)<0
xx(i)=-xx(i);
end
if (xx.inf(i)<=0)&&(xx.sup(i)>=0)
xx(i)=infsup(0,-xx.inf(i));
end
end
% verified quantities
x=xx(1:n); p=xx(n+1:n+m); y=p;                                  % verified optimal solution and Lagrangian multipliers
D=infsup(D,D); c=infsup(c,c);
h=0.5*(x'*D*x)+c'*x;                                            % verified optimal value
flag='verified optimum   ';
setround(gr);
%
% Subfunctions
%
function [pd,x]=verpd(A)
%    VERPD    Verified positive definiteness of a real matrix.
%
%    For a symmetric real matrix A,
%    [pd,x]=verpd(A)
%    verifies positive definiteness or not-positive-definiteness of A,
%    or yields no verified result:
%
%    pd= 1 :  A verified positive definite,
%    pd= 0 :  A verified not positive definite, and
%             x is verified to satisfy x'*A*x<=0,
%    pd=-1 :  no verified result,
%    pd=-2 :  A not real or not symmetric.
%
%    If pd~=0, then x is a vector of NaN's.
%

%
gr=getround;
setround(0);
% checking data
n=size(A,1);
pd=-1; x=repmat(NaN,n,1);
if isintval(A)||~isequal(A,A')
pd=-2;
setround(gr); return
end
% verifying positive definiteness
if isspd(A)==1                                                  % routine isspd by S. M. Rump
pd=1;                                                       % verified positive definite
setround(gr); return
end
% verifying not-positive-definiteness
% first, checking the diagonal and the 2x2 condition
A=infsup(A,A);
for i=1:n
if A.sup(i,i)<=0                                            % diagonal entry nonpositive
pd=0;                                                   % verified not positive definite
x=zeros(n,1); x(i)=1;
setround(gr); return
end
for j=i+1:n
a=A(i,i)+A(j,j)-2*A(i,j);
if a.sup<=0                                             % 2x2 condition not satisfied
pd=0;                                               % verified not positive definite
x=zeros(n,1); x(i)=1; x(j)=-1;
setround(gr); return
end
end
end
% second, finding a random x satisfying x'*A*x<=0
for i=1:max(n^2,1e03)
x=2*rand(n,1)-1;
while isequal(x,zeros(n,1))
x=2*rand(n,1)-1;
end
a=x'*A*x;
if a.sup<=0                                                 % x'*A*x verified nonpositive
pd=0;                                                   % verified not positive definite
setround(gr); return
end
end
x=repmat(NaN,n,1);                                              % no verified result
setround(gr);
%
function [x,B,N]=veropt(A,b,c)                                  % unchanged against verlinprog
% B is the "basis index set" of an optimal solution of the LP problem
% min c'*x  subject to  A*x=b, x>=0,
% x is a verified basic feasible solution with this basis
% N is the set of nonbasic indices
[m,n]=size(A);
x=repmat(infsup(NaN,NaN),n,1);
B=repmat(NaN,m,1); N=repmat(NaN,n,1);
options=optimset('linprog');
options=optimset(options,'Display','off');                     % suppresses displaying linprog's error messages
[xopt,fval,exitflag]=linprog(c,-eye(n),zeros(n,1),A,b,[],[],[],options);
if exitflag<=0
return
end
[xx,J]=sort(xopt); B=J(n-m+1:n); N=J(1:n-m);                    % B is set of "basic" indices, N of "nonbasic" ones
% xB=verifylss(A(:,B),b);                                       % old version using full A; replaced by next 5 lines
AB=A(:,B);
if issparse(AB)
AB=full(AB);                                                % only the square submatrix taken full
end
xB=verifylss(AB,b);
return
end
x=infsup(zeros(n,1),zeros(n,1)); x(B)=xB;                       % verified "optimal" solution found
%
function y=verinfeas(A,b)                                       % unchanged against verlinprog
% y is a Farkas vector (i.e., satisfying A'*y>=0, b'*y<0)
% its existence implies infeasibility of A*x=b
[m,n]=size(A);
y=repmat(infsup(NaN,NaN),m,1);
ep=max(1e-08,max([m n 100])*max([norm(A,inf) norm(b,inf)])*eps);
Afv=[A' -A'    -eye(n) zeros(n,1);                              % Afv is (n+1)x(2*m+n+1)
b' -b' zeros(1,n)         1];
bfv=[zeros(n,1)' -1]';                                          % bfv is (n+1)x1
bfv=bfv+ep*[ones(1,n) -1]';                                     % perturbation to compensate roundoff errors (so that A'*y>=0)
yf=veropt(Afv,bfv,ones(2*m+n+1,1));                             % system: A'*y>=0, b'*y<=-1, y written as y=y1-y2
if ~isnan(yf.inf(1))
yf=mid(yf);
y1=yf(1:m);
y2=yf(m+1:2*m);
yf=y1-y2;                                                   % would-be Farkas vector
A=infsup(A,A); b=infsup(b,b); yf=infsup(yf,yf);             % (i.e., should satisfy A'*y>=0, b'*y<0)
alpha=A'*yf; beta=b'*yf;
if (all(alpha.inf>=0))&&(beta.sup<0)                        % infeasibility verified
y=yf;                                                   % Farkas vector outputed
return
else
return
end
else
return
end
%
function xx=absvalverifn(A,B,b,x)                     % absvalverifn of 2007-03-26
%    ABSVALVERIFN     Verification of a solution of the equation A*x+B*abs(x)=b.
%
%    For an approximate real solution x of the equation
%        A*x+B*abs(x)=b,     (1)
%    xx=absvalverifn(A,B,b,x)
%    either produces a tight interval vector xx verified to contain a solution
%    of (1), or fails (yields an interval vector xx of NaN's).

%
n=size(A,1);
xx=repmat(infsup(NaN,NaN),n,1);
if ~isintval(A), A=infsup(A,A); end                   % converting data from real to intval
if ~isintval(B), B=infsup(B,B); end
if ~isintval(b), b=infsup(b,b); end
if  isintval(x), return, end
if issparse(A)                                        % A, B made full because of verifylss
A=full(A);
end
if issparse(B)
B=full(B);
end
I=eye(n,n);
z=ones(n,1);
for j=1:n
if x(j)<0, z(j)=-1; end                           % sign vector of x
end
infinf=infsup(repmat(-Inf,n,1),repmat(Inf,n,1));
x1=verifylss(A+B*diag(z),b);                          % enclosure via verifylss
if ~(~isnan(x1.inf(1))&&all(z.*inf(x1)>=0)&&all(z.*sup(x1)>=0))
x1=infinf;                                        % failure to produce verified output
end
M=inv(I-abs(inv(A)*B));
if (~isnan(M.inf(1,1)))&&(all(all(M.inf>=0)))         % M verified nonnegative