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
%