function [X,all]=absvaleqnall(A,B,b) % ABSolute VALue EQuatioN, ALL solutions % VERSION FOR POSTING
%
%  A,B square nxn, b in R^n.
%  all= 1:  X contains all solutions of A*x+B*abs(x)=b as its columns,
%  all=-1:  each column of X is a solution of A*x+B*abs(x)=b, but it is
%           not guaranteed that all of them are included.
%
% EXAMPLE (n=7, 10 solutions).
%
% 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]=absvaleqnall(A,B,b)
% 
% 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
% 
% X =
% 
%     0.2842   -1.9018    0.1484   -0.6615   -4.3204   -1.8897    0.2118    0.1584    0.1048    0.2798
%     0.2852   -0.3674    0.4041    0.5318   -0.2405    0.4361    0.3700    0.3711    0.3815    0.2885
%    -0.0841   -0.7374    0.7863    0.5816   -1.2114   -0.3516    0.1697    0.1642    0.2708   -0.0792
%    -0.0106    2.2570    0.1354    0.9473    6.0074    2.5083    0.0233   -0.1676   -0.2813   -0.0201
%     0.2235   -1.0900   -0.1987   -0.3510   -2.2160   -0.7775    0.2024    0.1049   -0.0257    0.2208
%     0.0125    0.4788   -0.3220   -0.2792   -0.1360   -0.0891   -0.0745   -0.0478   -0.0703    0.0114
%     0.0045    0.8552    0.2679    0.6275    2.8807    1.2929    0.0600   -0.0899   -0.1583   -0.0032
% 
% all =
% 
%      1
%
%  MIT License for ABSVALEQNALL:
%
%  Copyright 2009-2017 Jiri Rohn (rohn@cs.cas.cz)
%
%  Permission is hereby granted, free of charge, to any person obtaining a copy
%  of this software and associated documentation files (the "Software"), to deal
%  in the Software without restriction, including without limitation the rights
%  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
%  copies of the Software, and to permit persons to whom the Software is
%  furnished to do so, subject to the following conditions:
% 
%  The above copyright notice and this permission notice shall be included in all
%  copies or substantial portions of the Software.
% 
%  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
%  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
%  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
%  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
%  SOFTWARE.
%
X=[]; all=1;
n=length(b);
ep=max(n,100)*max([norm(A,inf) norm(B,inf) norm(b,inf)])*eps; 
if ((rank(A)==n)&&(max(abs(eig(abs(inv(A))*abs(B))))<1))
    x=absvaleqn(A,B,b);
    X=x; return
end
y=zeros(1,n); 
z=ones(1,n); 
if rank(A+B*diag(z))==n
    x=(A+B*diag(z))\b;
    C=-inv(A+B*diag(z))*B;
    if diag(z)*x >= -ep*ones(n,1), X=[X x]; end 
else
    all=-1; return
end
while any(y~=ones(1,n))
    k=find(y==0,1);
    y(1:(k-1))=zeros(1,k-1);
    y(k)=1;
    z(k)=-z(k);
    if abs(1-2*z(k)*C(k,k))>ep
        alpha=(2*z(k))/(1-2*z(k)*C(k,k));
        x=x+alpha*x(k)*C(:,k);
        C=C+alpha*C(:,k)*C(k,:);
        if diag(z)*x >= -ep*ones(n,1)
            x=(A+B*diag(z))\b; 
            C=-inv(A+B*diag(z))*B;
            X=[X x]; 
        end
    else
        all=-1; return
    end
end
%
function [x,S,iter]=absvaleqn(A,B,b) % ABSolute VALue EQuatioN % VERSION FOR POSTING
%
% ~isempty(x): x solves A*x+B*abs(x)=b (A,B square), S is empty,
% ~isempty(S): S is a singular matrix satisfying abs(S-A)<=abs(B), x is empty.
%  iter: number of iterations (it may be zero).
%
%  MIT License for ABSVALEQN:
%
%  Copyright 2005-2016 Jiri Rohn (rohn@cs.cas.cz)
%
%  Permission is hereby granted, free of charge, to any person obtaining a copy
%  of this software and associated documentation files (the "Software"), to deal
%  in the Software without restriction, including without limitation the rights
%  to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
%  copies of the Software, and to permit persons to whom the Software is
%  furnished to do so, subject to the following conditions:
% 
%  The above copyright notice and this permission notice shall be included in all
%  copies or substantial portions of the Software.
% 
%  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
%  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
%  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
%  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
%  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
%  OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
%  SOFTWARE.
%
b=b(:); n=length(b); I=eye(n,n);
ep=n*(max([norm(A,inf) norm(B,inf) norm(b,inf)]))*eps;
x=[]; S=[]; iter=0;
if rank(A)<n, S=A; return, end
x=A\b; z=sgn(x);
if rank(A+B*diag(z))<n, S=A+B*diag(z); x=[]; return, end
x=(A+B*diag(z))\b;
C=-inv(A+B*diag(z))*B;
X=zeros(n,n);
r=zeros(1,n);
while any(z.*x<-ep)
    k=find(z.*x<-ep,1);
    iter=iter+1;
    if 1+2*z(k)*C(k,k)<=0
        S=A+B*(diag(z)+(1/C(k,k))*I(:,k)*I(k,:)); x=[]; return
    end
    if ((k<n)&&(all(r(k)>r(k+1:n))))||((k==n)&&(r(k)>0))
        x=x-X(:,k); z=sgn(x);
        ct=A*x;
        jm=abs(B)*abs(x);
        y=zeros(1,n);
        for i=1:n
            if jm(i)>ep, y(i)=ct(i)/jm(i); else y(i)=1; end
        end
        S=A-diag(y)*abs(B)*diag(z); x=[]; return
    end
    X(:,k)=x;
    r(k)=iter;
    z(k)=-z(k);
    alpha=2*z(k)/(1-2*z(k)*C(k,k));
    x=x+alpha*x(k)*C(:,k);
    C=C+alpha*C(:,k)*C(k,:);
end 
%
function z=sgn(x) % SiGN vector of x (column)
n=length(x); 
z=zeros(n,1);
for j=1:n
    if x(j)>=0, z(j)=1; else z(j)=-1; end 
end
%