function [G, vres] = lr(A, tol, maxiter)
% Logarithmic reduction
m = size(A,1);
% Inizializzazione
step = 0;
H=inv(eye(m)-A(:,:,2));
V1= H*A(:,:,1);
V3= H*A(:,:,3);
W = V1*V3+V3*V1;
W = eye(m) - W;
G = V1;
U = eye (m);
err = 1-min(sum(G,2));
% Processo iterativo
% while abs(err) > tol 
while (err>tol & step<maxiter)
    step = step + 1;
    U = U*V3;
    H = inv(W);
    V1 = H * V1 * V1;
    V3 = H * V3 * V3;
    W = V1*V3 + V3*V1;
    W = eye(m) - W;
    G = G + U* V1;
    err = 1- min(sum(G,2))
    res=norm( (A(:,:,3)*G+A(:,:,2))*G+A(:,:,1)-G, inf );
    if step== 1
        vres = res;
    else
        vres=[vres,res];
    end
    step, res
end