```function [x,all,E]=verabsvaleqnall(A,B,b,k,t) % uses performwithy, xysolution, absvalverifn, sgn, timeprint
%    VERABSVALEQNALL     Verified all solutions of the equation A*x+B*abs(x)=b.
%
%    This is an INTLAB file. It requires to have INTLAB installed under
%    MATLAB to function properly.
%
%    For real square matrices A, B of the same size and a matching real vector b,
%        [X,all,E]=verabsvaleqnall(A,B,b)
%    produces a tight interval matrix X verified to contain as columns all
%    or some of the verified solutions of the equation
%        A*x+B*abs(x)=b   (1)
%    (explained below).
%
%    Possible outputs:
%
%    all== 1:   it is verified that X contains as columns all solutions
%               of (1), all of them verified; thus if X is empty, then
%               it is verified that (1) has no solution,
%
%    all==-1:   columns of X are verified solutions of (1), but it is not
%               verified that X contains all the solutions.
%
%    The command
%        [X,all,E]=verabsvaleqnall(A,B,b,k)
%    (i.e., with additional integer input argument k) enforces the computation
%    to be stopped after a prescribed number of solutions k has been found;
%    this can essentially reduce the computation time. In this case it is
%    always all==-1. In particular, use
%        [X,all,E]=verabsvaleqnall(A,B,b,1)
%    if you are interested in one solution only.
%
%    The command
%        [X,all,E]=verabsvaleqnall(A,B,b,k,1)
%    (i.e., with additional input argument "1") performs the same task, but
%    it also produces screen output of the form
%        Expected remaining time: ... sec.
%        Number of solutions found: ...
%    This feature, however, slows down the actual computation. Use k=Inf if
%    you do not wish to bound the number of solutions to be found.
%
%    The structured array E explains details of the output. It has three
%    fields: E.error, E.where, E.value.
%
%    EXAMPLE 1 (matrices 7-by-7: 10 solutions).
%    >> tic, n=7; rand('state',671); A=2*rand(n,n)-1, B=2*rand(n,n)-1, b=2*rand(n,1)-1, [x,all]=verabsvaleqnall(A,B,b), toc
%    A =
%       -0.1479   -0.5985   -0.2265   -0.2292   -0.2426   -0.4978    0.4772
%        0.3503    0.7914   -0.8554    0.2560   -0.4149   -0.3221   -0.5674
%       -0.8144    0.8176   -0.9111   -0.9181    0.1953   -0.9376    0.0201
%        0.1143   -0.8706   -0.1203    0.5198   -0.6242   -0.7633   -0.1536
%        0.7850   -0.7964    0.6195   -0.5218    0.9041    0.7736    0.9708
%       -0.4198   -0.5983    0.9180   -0.5057   -0.6677    0.1967    0.0734
%       -0.1962    0.6255   -0.3860    0.1035    0.4396   -0.7893   -0.9860
%    B =
%       -0.8464   -0.5703   -0.9208   -0.0867    0.2831    0.9318    0.8203
%       -0.7984    0.3861   -0.1074   -0.1288    0.8478    0.8475    0.8466
%        0.3445    0.4156    0.7606   -0.4585    0.9195    0.0428    0.0485
%       -0.1394    0.8962   -0.2990   -0.2622   -0.6214   -0.5709   -0.1978
%        0.8221    0.1798   -0.2713    0.9308   -0.9663    0.9149   -0.0731
%        0.8508   -0.2720   -0.7906   -0.8783    0.5006   -0.9402    0.6437
%        0.7253    0.0865    0.5792   -0.1374   -0.0348    0.4932   -0.2036
%    b =
%       -0.6525
%        0.3719
%        0.6019
%       -0.3199
%        0.2327
%       -0.3168
%        0.5135
%    intval x =
%      Columns 1 through 5
%    [    0.1483,    0.1484] [   -0.6616,   -0.6615] [   -4.3204,   -4.3203] [   -1.8898,   -1.8897] [    0.2118,    0.2119]
%    [    0.4041,    0.4042] [    0.5318,    0.5319] [   -0.2405,   -0.2404] [    0.4361,    0.4362] [    0.3699,    0.3700]
%    [    0.7863,    0.7864] [    0.5816,    0.5817] [   -1.2115,   -1.2114] [   -0.3516,   -0.3515] [    0.1696,    0.1697]
%    [    0.1354,    0.1355] [    0.9473,    0.9474] [    6.0073,    6.0074] [    2.5083,    2.5084] [    0.0233,    0.0234]
%    [   -0.1988,   -0.1987] [   -0.3511,   -0.3510] [   -2.2161,   -2.2160] [   -0.7775,   -0.7774] [    0.2023,    0.2024]
%    [   -0.3220,   -0.3219] [   -0.2792,   -0.2791] [   -0.1361,   -0.1360] [   -0.0891,   -0.0890] [   -0.0745,   -0.0744]
%    [    0.2679,    0.2680] [    0.6275,    0.6276] [    2.8806,    2.8807] [    1.2929,    1.2930] [    0.0600,    0.0601]
%   Columns 6 through 10
%    [    0.2842,    0.2843] [   -1.9018,   -1.9017] [    0.2798,    0.2799] [    0.1584,    0.1585] [    0.1048,    0.1049]
%    [    0.2851,    0.2852] [   -0.3675,   -0.3674] [    0.2884,    0.2885] [    0.3711,    0.3712] [    0.3815,    0.3816]
%    [   -0.0842,   -0.0841] [   -0.7374,   -0.7373] [   -0.0792,   -0.0791] [    0.1642,    0.1643] [    0.2708,    0.2709]
%    [   -0.0107,   -0.0106] [    2.2569,    2.2570] [   -0.0202,   -0.0201] [   -0.1676,   -0.1675] [   -0.2814,   -0.2813]
%    [    0.2235,    0.2236] [   -1.0900,   -1.0899] [    0.2207,    0.2208] [    0.1049,    0.1050] [   -0.0258,   -0.0257]
%    [    0.0125,    0.0126] [    0.4787,    0.4788] [    0.0113,    0.0114] [   -0.0478,   -0.0477] [   -0.0704,   -0.0703]
%    [    0.0044,    0.0045] [    0.8552,    0.8553] [   -0.0032,   -0.0031] [   -0.0900,   -0.0899] [   -0.1584,   -0.1583]
%    all =
%         1
%    Elapsed time is 1.017689 seconds.
%    Comment: All 10 verified solutions have been found.
%
%    WARNING. The algorithm may generally take up to 2^n steps, where n=size(A,1);
%    it is therefore recommended to run the file for moderate values of n
%    (say, n<=17) only. The exception is the case of a unique solution
%    which, in contrast, can be usually found very quickly. Also, the
%    equation (1) may have up to 2^n solutions, as shown by the simple example
%        zeros(n,n)*x+eye(n,n)*abs(x)=ones(n,1),
%    whose solutions are all the plus/minus vectors in R^n, i.e., 2^n of
%    them.
%
%    EXAMPLE 2 (matrices 10-by-10: 1024 solutions).
%    >> tic, n=10; rand('state',1); A=0.1*(2*rand(n,n)-1), B=rand(n,n); B=inv(B), b=rand(n,1), [x,all]=verabsvaleqnall(A,B,b); sols=size(x,2), all, toc
%    A =
%      Columns 1 through 7
%        0.0906    0.0797    0.0537   -0.0211   -0.0119    0.0939    0.0799
%        0.0408   -0.0142   -0.0881    0.0383    0.0994   -0.0393    0.0878
%        0.0908   -0.0601    0.0254    0.0139    0.0125   -0.0252   -0.0452
%        0.0196   -0.0394   -0.0470   -0.0314   -0.0068   -0.0303   -0.0046
%        0.0681    0.0077   -0.0375    0.0190   -0.0162    0.0402   -0.0493
%       -0.0114    0.0820    0.0045   -0.0452   -0.0577    0.0229   -0.0276
%        0.0674    0.0051   -0.0183   -0.0904   -0.0528    0.0518    0.0054
%        0.0037   -0.0386    0.0786    0.0676   -0.0372    0.0598   -0.0238
%       -0.0956   -0.0931    0.0148   -0.0795    0.0740   -0.0569    0.0647
%       -0.0248    0.0431    0.0136    0.0457   -0.0496    0.0787   -0.0371
%      Columns 8 through 10
%        0.0326   -0.0732    0.0081
%        0.0295    0.0870    0.0821
%       -0.0059   -0.0965    0.0162
%        0.0436   -0.0878   -0.0590
%       -0.0962    0.0015    0.0768
%       -0.0043   -0.0246   -0.0429
%        0.0979   -0.0795   -0.0703
%        0.0955    0.0010    0.0952
%       -0.0004   -0.0540   -0.0464
%       -0.0918   -0.0612   -0.0713
%    B =
%      Columns 1 through 7
%       -0.2352   -0.5766    1.8918   -0.4172   -0.0450   -1.4132    0.1931
%       -0.2286    0.1210   -0.1436   -0.6798   -0.0777    0.1745   -0.3110
%       -0.7952   -0.0076    0.8798   -0.4007    0.3096   -0.1172    0.8636
%        0.3066   -0.1892   -0.0645    0.5488    0.9801   -0.5211   -0.3190
%        0.6403    1.4887    0.3371   -0.2049   -0.0406    0.0974    0.1317
%       -0.0278    0.8313   -0.5236    0.7833   -0.5226    0.1913   -0.0099
%        0.5736    0.2986    0.1326   -0.5930   -0.3257   -1.0060    0.6667
%       -0.0876   -1.5485   -1.2575    1.3963    0.2402    1.0701   -1.1621
%       -0.1125   -0.5930   -0.4584    0.3201   -0.3340    0.6841   -0.4610
%        0.3664   -0.4505   -1.2536    0.1229    0.2388    1.2614   -0.0462
%      Columns 8 through 10
%       -0.4376   -0.1736    0.9970
%        0.2743    0.3365    0.5909
%        0.1284   -0.8769    0.1664
%       -0.5032    0.7230   -0.4455
%       -0.6928   -0.7540   -0.5845
%       -0.2025   -0.0502    0.0523
%       -0.1276    0.3920    0.0121
%        1.5229    0.5034   -0.5145
%        0.4155    1.3115   -0.6138
%        0.3281   -0.4432    0.1292
%    b =
%        0.1951
%        0.9186
%        0.7427
%        0.0069
%        0.8383
%        0.3371
%        0.0645
%        0.7455
%        0.3865
%        0.2190
%    sols =
%            1024
%    all =
%         1
%    Elapsed time is 8.318330 seconds.
%    Comment: The equation has exactly 1024 (=2^(10)) solutions, all of them found verified.
%

%
%    Looks for verified solutions of A*x+B*abs(x)=b through full search
%    over all systems (A+B*diag(y))*x=b, diag(y)*x>=0, y a plus/minus-one
%    vector (ynset algorithm).
%
%    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.
%
%    History
%
%    2007-10 to 2007-12-10  original version (unpublished)
%    2008-11-02   changed to admit interval data (because of verlcp)
%    2008-11-04   output variable E added, case of unique solution handled
%                 via verabsvaleqn, input variable k added, help written
%    2008-11-05   special call of verabsvaleqn in case of k==1 added
%    2008-11-07   the same elaborated for the case of a call by verlcpall;
%                 print of the number of solutions currently found
%    2008-11-23   examples added; final version (html)
%
gr=getround;
setround(0);
% setting default outputs
x=[];   % x empty; later to be incremented
all=-1; % all==1: verified that all solutions found
E.error='verabsvaleqnall: none';
E.where='NaN';
E.value='NaN';
if (nargin<5)||((nargin==5)&&~isequal(t,1))
t=0; % no time display
end
if nargin<4
k=Inf; % default: all solutions
end
if nargin<3
E.error='verabsvaleqnall: few input data';
setround(gr); return
end
b=b(:);
[m,n]=size(A);
kk=k; % in order input k not to interfere with k in the ynset loop
if ~(isequal(m,n)&&isequal(size(B),[n,n])&&isequal(length(b),n))
E.error='verabsvaleqnall: sizes do not match';
setround(gr); return
end
if nargin>=4&&~((kk>0&&round(kk)==kk)||kk==Inf) % kk nonpositive or not integer, not infinity
E.error='verabsvaleqnall: wrong bound on the number of solutions';
setround(gr); return
end
I=eye(n,n);
% next lines changed against verabsvaleqnall of 2007-12-10 (to admit interval data A, B, b)
% A1=infsup(A,A); B1=infsup(B,B); % old version
if ~isintval(A)
A1=infsup(A,A);
else
A1=A;
end
if ~isintval(B)
B1=infsup(B,B);
else
B1=B;
end
N=inv(A1);
M=inv(I-abs(N*B1));
if isa(M,'double') % bridging the current gap in verifylss
M=infsup(M,M);
end
if ~isnan(M.inf(1,1))&&~any(any(M.inf<0))
reg=1; % verified regular (unique solution)
else
reg=-1;
end
tic, sh=1; told=-sh; % put separately because of performwithy under t~=1, and of premature stop
if isequal(t,1) % display of remaining time required
tic, sh=1;  % time shift (used in "timeprint")
told=-sh;
end
if reg==1||kk==1 % one solution to be outputted
if ~isintval(A)&&~isintval(B)&&~isintval(b) % all data real
xr=verabsvaleqn(A,B,b); % verabsvaleqn requires noninterval data
if ~isnan(xr.inf(1)) % verabsvaleqn successful
x=xr;
if reg==1
all=1; % solution unique
E.error='verabsvaleqnall: none; full search stopped prematurely, unique solution found';
else % kk==1
all=-1;
E.error='verabsvaleqnall: none; full search stopped prematurely, prescribed number of solutions reached';
end
if isequal(t,1)
clc
% disp(['Expected remaining time: ', num2str(0), ' sec.'])
disp(['Elapsed time: ', num2str(round(toc)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
end
setround(gr); return % and exit verabsvaleqnall
end
else % not all data real (case of verlcpall)
% making real input A2, B2, b2
if ~isintval(A), A2=A; else A2=mid(A); end
if ~isintval(B), B2=B; else B2=mid(B); end
if ~isintval(b), b2=b; else b2=mid(b); end
xr0=verabsvaleqn(A2,B2,b2); % verabsvaleqn requires noninterval data
if ~isnan(xr0.inf(1)) % solution with noninterval data found
xr=absvalverifn(A,B,b,xr0,M,N); % verification of the solution xr0 based on noninterval data
if ~isnan(xr.inf(1)) % verification successful
x=xr;
if reg==1
all=1; % solution unique
E.error='verabsvaleqnall: none; full search stopped prematurely, unique solution found';
else % kk==1
all=-1;
E.error='verabsvaleqnall: none; full search stopped prematurely, prescribed number of solutions reached';
end
if isequal(t,1)
clc
% disp(['Expected remaining time: ', num2str(0), ' sec.'])
disp(['Elapsed time: ', num2str(round(toc)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
end
setround(gr); return % and exit verabsvaleqn
end
end
end
end
iter=0; % iteration counter
allfound=1; % allfound will remain 1 if it is verified that all solutions have been found
% ynset algorithm starts
% initial y for ynset
z=zeros(1,n);
x0=verifylss(A,b);
if ~isnan(x0.inf(1))
x0=mid(x0); % made real to get sign
y= sgn(x0); % initial y taken as the sign vector of the solution of Ax=b (supported by the result in Fiedler et al.)
else % default setting if verifylss fails
y=ones(1,n);
end
% perform the task with initial y
[x,iter,all,allfound,told,exit,Eperformwithy]=performwithy(A,B,b,y,M,N,x,reg,t,all,allfound,n,iter,sh,told);
if exit==1
all=1; % to be sure
E.error='verabsvaleqnall: none; full search stopped prematurely, unique solution found';
clc
disp(['Elapsed time: ', num2str(round(toc)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
setround(gr); return
end
if size(x,2)==kk % prescribed bound reached
all=-1;
E.error='verabsvaleqnall: none; full search stopped prematurely, prescribed number of solutions reached';
clc
disp(['Elapsed time: ', num2str(round(toc)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
setround(gr); return
end
if ~isequal(Eperformwithy.where,'NaN')
E=Eperformwithy;
E.error=['verabsvaleqnall:' E.error];
end
% ynset algorithm
while any(z~=ones(1,n))
% determining new y
k=find(z==0,1);
z(1:(k-1))=zeros(1,k-1);
z(k)=1;
y(k)=-y(k);
% perform the task with y
[x,iter,all,allfound,told,exit,Eperformwithy]=performwithy(A,B,b,y,M,N,x,reg,t,all,allfound,n,iter,sh,told);
if exit==1
all=1; % to be sure
E.error='verabsvaleqnall: none; full search stopped prematurely, unique solution found';
clc
disp(['Elapsed time: ', num2str(round(toc)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
setround(gr); return
end
if size(x,2)==kk % prescribed bound reached
all=-1;
E.error='verabsvaleqnall: none; full search stopped prematurely, prescribed number of solutions reached';
clc
disp(['Elapsed time: ', num2str(round(toc)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
setround(gr); return
end
if ~isequal(Eperformwithy.where,'NaN')
E=Eperformwithy;
E.error=['verabsvaleqnall:' E.error];
end
end % of ynset
if isequal(allfound,1)
all=1; % verified that all solutions found
end
if isequal(t,1) % remaining time==0
clc
% disp(['Expected remaining time: ', num2str(0), ' sec.'])
disp(['Elapsed time: ', num2str(round(toc)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
end
setround(gr);
% full search performed
% end of verabsvaleqnall
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Subfunctions: performwithy, xysolution, absvalverifn, sgn, timeprint
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function [x,iter,all,allfound,told,exit,E]=performwithy(A,B,b,y,M,N,x,reg,t,all,allfound,n,iter,sh,told) % uses xysolution
% performs the main task of the ynset algorithm (called twice)
% N=inv(A);
% M=inv(I-abs(N*B));
% x is the matrix whose columns are solutions already found
% only exit and E are new variables
E.error='performwithy: none';
E.where='NaN';
E.value='NaN';
iter=iter+1;
[xy,eexist,Exysolution]=xysolution(A,B,b,y,M,N);
if ~isnan(xy.inf(1))
if isempty(x)
x=xy;
if reg==1 % spectral condition satisfied, unique solution found
all=1;
if isequal(t,1)
% clc
% disp(['Expected remaining time: ', num2str(0), ' sec.'])
disp(['Elapsed time: ', num2str(round(toc)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
end
exit=1; return % and exit verabsvaleqnall
end
else
if ~ismember(xy',x','rows') % add only if not there yet
x=[x xy];
end
end
else
E=Exysolution;
E.error=['performwithy:' E.error];
end
if ~(isequal(eexist,1)||isequal(eexist,0)) % verified (non)existence of xy-solution not established
allfound=allfound-1;
end
if isequal(t,1) % "Expected remaining time" printed only here
told=timeprint(2^n-iter,iter,sh,told,x);
end
exit=-1;
% end of performwithy
%
function [xy,exist,E]=xysolution(A,B,b,y,M,N) % uses absvalverifn
% solution of (A+B*diag(y))*x=b, diag(y)*x>=0 for particular y
% M, N passed to absvalverifn
% N=inv(A);
% M=inv(I-abs(N*B));
% xy interval vector (verified or of NaN's)
% exist=1 for verified solution, =0 for verified nonexistence, =-1 otherwise
n=length(b);
xy=repmat(infsup(NaN,NaN),n,1); exist=-1;
E.error='xysolution: none';
E.where='NaN';
E.value='NaN';
xcurr=verifylss(A+B*diag(y),b); % current solution, interval quantity
if isnan(xcurr.inf(1))&&~isintval(A)&&~isintval(B)&&~isintval(b) % verifylss fails
xcurr=pinv(A+B*diag(y))*b; % replaced by lsq solution; real quantity
xcurr=absvalverifn(A,B,b,xcurr,M,N); % verification; interval quantity
end
if ~isnan(xcurr.inf(1))
yx=diag(y)*xcurr;
if all(yx.inf>=0)
xy=xcurr; exist=1; % verified solution
return
end
if all(yx.sup>=0) % possibility of solution
xx=absvalverifn(A,B,b,mid(xcurr),M,N); % verification
if ~isnan(xx.inf(1))
xy=xx; exist=1; % verified solution
return
end
end
if any(yx.sup<0) % solution verified not to exist
exist=0; % verified no solution for y
return
end
end
if isnan(xcurr.inf(1))
E.error='xysolution: verifylss failed';
E.where=y;
E.value=xcurr;
else % ~isnan(xcurr.inf(1))
E.error='xysolution: verification failed';
E.where=y;
E.value=xcurr;
end
% end of xysolution
%
function xx=absvalverifn(A,B,b,x,M,N) % uses only verifylss
%    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).
%    M, N passed in order not to be computed anew.
%    N=inv(A);
%    M=inv(I-abs(N*B));

%
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
z=sgn(x);
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
if (~isnan(M.inf(1,1)))&&(all(all(M.inf>=0)))         % M verified nonnegative
else
x2=infinf;                                        % failure to produce verified output
end
x3=intersect(x1,x2);                                  % intersection of the two enclosures
if x3~=infinf                                         % at least one enclosure computed
xx=x3;                                            % verified enclosure of the solution of Ax+B|x|=b
end
% end of absvalverifn
%
function z=sgn(x)
% signum of x for real or intval x
% x column or row, z column
n=length(x);
z=zeros(n,1);
if ~isintval(x) % real vector
for i=1:n
if x(i)>=0
z(i)=1;
else
z(i)=-1;
end
end
else % intval vector
for i=1:n
if x.inf(i)>0
z(i)=1;
end
if x.sup(i)<0
z(i)=-1;
end
% otherwise z(i)=0
end
end
% end of sgn
%
function told=timeprint(z,d,sh,told,x)
% prints remaining time
if toc-told>=sh
clc
disp(['Expected remaining time: ', num2str(round((toc/d)*z)), ' sec.'])
disp(['Number of solutions found: ', num2str(size(x,2))])
told=toc;
end
% end of timeprint
%
```