function [G,vres]=cr(A, tol, maxiter)
% Computes the minimal solvent G of A_3 X^2 +  A_2 X + A_1= X
% by means of cyclic reduction
% eps: stop condition
m=size(A,1);
A1=-A(:,:,3); B1=eye(m)-A(:,:,2); C1=-A(:,:,1); HB=B1;
iter=0;
vres=[];
while (norm(A1,1)>tol & iter<maxiter)
    V=inv(B1);
    U=V*A1;
    W=V*C1;
    V=A1*W;
    B1=B1-C1*U-V;
    HB=HB-V;
    A1=-A1*U;
    C1=-C1*W;
    iter=iter+1
    norma_A1=norm(A1, 1)
    G = inv(HB)*A(:,:,1);
    res = norm( (A(:,:,3)*G+A(:,:,2))*G+A(:,:,1)-G, inf)
      vres=[vres,res];
end
% G=inv(HB)*A(:,:,1);
% res=norm( (A(:,:,3)*G+A(:,:,2))*G+A(:,:,1)-G, inf );
