function testnormdecay(varargin)
% Constructs sequence of matrices with maximal spectral radius below a given threshold.
%
% Options:
%   'RNG',val                       random number generator seed
%   'dim',val                       dimension
%   'delta',val                     threshold of spectral radius. Default=1
%
% Output:
%   A nice picture
%
% Written by: tommsch, 2018

RNG=parsem('RNG',varargin,0); if(RNG); rng(RNG); end;
dim=parsem('dim',varargin,3);
delta=parsem('delta',varargin,1);

A=randn(dim); A=2*A/rho(A);
B=randn(dim); B=0.5*B/rho(B);

clf; hold on;
MAXNORM_NORMALIZED=zeros(1,0);
MAXNORM=zeros(1,0);
RHO=zeros(1,0);

clf; axis equal;
P=[''];
X=eye(dim);
REFERENCE=1;

for i=1:20000
    if(rho(A*X)^(1/(length(P)+1))>delta);
        if(rho(B*X)^(1/(length(P)+1))>delta);
            fprintf('ERROR\n'); break;
        else
            X=B*X; P=['B' P];
        end
    else
        X=A*X; P=['A' P];
    end

    RHO(i)=rho(X)^(1/length(P));
    if(RHO(i)>delta);
        fprintf('ERROR');
    end;
	REFERENCE(i)=REFERENCE(end)*delta;
    MAXNORM_NORMALIZED(i)=norm(X)^(1/length(P));
    MAX(i)=norm(X);
end
semilogy(MAXNORM_NORMALIZED,'or'); hold on;
semilogy(MAX,'xk'); hold on;
semilogy(REFERENCE,'*g'); hold on;
%semilogy(RHO,'bo'); hold on;
vprintf('Product: %v\n', P);
%vprintf('max(RHO) = %i\n', max(RHO));
vprintf('max(normalized norm) = %i\n', max(MAXNORM_NORMALIZED ,[],2));

vprintf('red: normalized norm of matrix-product.\n','cpr',[.8 0 0]);
vprintf('black: norm of matrix-product.\n');
vprintf('green: geometric decay with speed delta.\n','cpr',[0 .8 0]);


%vprintf('blue: normalized spectral radius of matrix-product.\n','cpr',[0 0 .8]);
%tjsr({A,B},'verbose',1);

end


function dummy; end %#ok<DEFNU> %Generates an error, if the 'end' of a function is missing. 
