```function [flag,x,y,h]=verlinprog(A,b,c)
%    VERLINPROG       Verified linear programming.
%
%    This is an INTLAB file. It requires to have INTLAB installed under
%    MATLAB to function properly.
%
%    For a real matrix A (full or sparse) and matching real vectors b, c,
%    [flag,x,y,h]=verlinprog(A,b,c)
%    either computes verified optimal solution x, verified dual optimal solution y
%    and verified optimal value h of the linear programming problem
%        min c'*x   subject to   A*x=b, x>=0,
%    or verifies (in)feasibility, or verifies unboundedness, 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 a primal optimal solution,
%                                 y is verified to enclose a dual optimal solution,
%                                 h is verified to enclose the optimal value,
%    flag='verified unbounded ' : x is verified to enclose a primal feasible solution xo, and
%                                 y is verified to enclose a vector yo such that the objective
%                                   tends to -Inf along the feasible half-line {xo+t*yo | t>=0},
%                                 h is NaN,
%    flag='verified feasible  ' : x is verified to enclose a primal feasible solution
%                                   (optimality nor unboundedness could be verified),
%                                 y, h are NaN's,
%    flag='verified infeasible' : y is verified to enclose a Farkas vector yo satisfying A'*yo>=0, b'*yo<0
%                                   (whose existence proves primal infeasibility),
%                                 x, h are NaN's,
%    flag='no verified result ' : x, y, h are NaN's
%                                   (no verified result could be found).
%    flag='sizes do not match ' : x, y, h are NaN's
%                                   (sizes of A, b, c do not match).
%
%    Complexity: The algorithm solves at most four linear programming
%    problems (independently of the size of the original problem) and uses
%    a verification procedure which runs 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);
b=b(:); c=c(:); [m,n]=size(A); p=length(b); q=length(c);
flag='no verified result ';
x=repmat(infsup(NaN,NaN),n,1); y=repmat(infsup(NaN,NaN),m,1); h=infsup(NaN,NaN);
if ~(m==p&&n==q)||(m>n),
flag='sizes do not match ';
setround(gr); return
end
if isintval(A)||isintval(b)||isintval(c),                       % error('data not real')
setround(gr); return
end
if issparse(b)
b=full(b);
end
if issparse(c)
c=full(c);
end
% verifying infeasibility
yi=verinfeas(A,b);
if ~isnan(yi.inf(1))                                            % verified Farkas vector found
y=yi; flag='verified infeasible';
setround(gr); return
end
% verifying feasibility
xf=veropt(A,b,ones(n,1));
flag='no verified result ';
setround(gr); return
end
% verifying unboundedness
yu=verunbound(A,c);
if ~isnan(yu.inf(1))                                            % verified descent direction found
x=xf; y=yu; flag='verified unbounded ';
setround(gr); return
end
% verifying optimality
[xo,B,N]=veropt(A,b,c);
if isnan(xo.inf(1))                                             % verified feasible primal solution with basis B not found
x=xf; flag='verified feasible  ';                           % previous feasible solution outputed
setround(gr); return
end
% yB=verifylss(A(:,B)',c(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
yB=verifylss(AB',c(B));
x=xo; flag='verified feasible  ';                           % candidate for optimum outputed as feasible solution
setround(gr); return
end
c=infsup(c,c); A=infsup(A,A);
crit=c'-yB'*A;                                                  % criterial row (dual feasibility)
crit=crit(N);                                                   % nonbasic part of it
x=xo; flag='verified feasible  ';                           % candidate for optimum outputed as feasible solution
setround(gr); return
end
% verified quantities                                           % verified primal and dual feasible solutions found
x=xo;                                                           % x is a verified primal optimal solution
y=yB;                                                           % y is a verified dual optimal solution
h1=c'*x; h2=b'*y;
h=intersect(h1,h2);                                             % h is a verified optimal value (duality theorem)
if isnan(h.inf)
h=h1;
end
flag='verified optimum   ';
setround(gr);
%
% Subfunctions
%
function [x,B,N]=veropt(A,b,c)
% 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 (because of verifylss)
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)
% y verified to enclose a Farkas vector yo (i.e., satisfying A'*yo>=0, b'*yo<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 y=verunbound(A,c)
% y verified to enclose a vector yo satisfying A*yo=0, yo>=0, c'*yo<=-1
% under feasibility its existence implies unboundedness
[m,n]=size(A);
y=repmat(infsup(NaN,NaN),n,1);
Aunb=[A zeros(m,1);                                             % Aunb is (m+1)x(n+1)
c'        1];
bunb=[zeros(1,m) -1]';                                          % bunb is (m+1)x1
yunb=veropt(Aunb,bunb,ones(n+1,1));                             % yunb is (n+1)x1
if ~isnan(yunb.inf(1))
y=yunb(1:n);                                                % y satisfies A*y=0, y>=0, c'*y=-1
return
end
```