function powersR(rho,tau);
%
% This function illustrates the behaviour of norm(R^k)
% for k = 1:300.  
%
% Inputs:  
%         rho - A scalar specifying the intial spectral radius.
%               It is assumed that 0 < rho < 1.
%
%         tau - A non-negative number used to control the degree of
%               non-normality.  
%
% An upper triangular matrix  R = rho*D + tau*U  is generated
% and then norm(R^k,1)   is recorded for each power k and the
% sequence is plotted.  The diagonal matrix D is generated 
% to have spectral radius 1, and U is generated to be 
% strictly upper triangular.   So the spectral radius of R for
% the three plots will be  rho, rho/10, rho/100.
%
%
% Recommended Values  rho = .1,  tau = .001, .1, 1, 10
%
% 
%----------------------------------------------------------
%
%  D.C. Sorensen 
%  10  Oct 03

   nplots = 3;
   n = 100; m = 300;

   D = diag(randn(n,1));
   [Q,R] = qr(randn(n));

   for j = 1:n, R(j,j) = 0; end
   U = R;
%
%  U is strict upper triangular.
%

   dmax = max(abs(diag(D)));
   D = D/dmax;
%
%  D is diagonal with spec. radius 1 
%


   t = zeros(m,nplots);
   for j = 1:nplots,
       R =  rho*D + tau*U;
%
%      R is upper triangular with spec. radius  rho
%
       Rk = eye(n);
       for k = 1:m,
           Rk = R*Rk;
           t(k,j) = norm(Rk,1);
       end
       SpecRad =  rho
       rho = rho/10;
   end
%
%  Plot the three convergence histories 
%  to compare
%
   figure(1), clf
   semilogy([1:m],t(:,1), 'r-', 'linewidth',2)
   ylabel('power norm, ||R^k||', 'fontsize', 18)
   xlabel('iteration k', 'fontsize', 18)
   hold
   disp('Strike a Key')
   pause
   semilogy([1:m],t(:,2), 'k-', 'linewidth',2)
   ylabel('power norm, ||R^k||', 'fontsize', 18)
   xlabel('iteration k', 'fontsize', 18)
   disp('Strike a Key')
   pause
   semilogy([1:m],t(:,3), 'b-', 'linewidth',2)
   xlabel('iteration k', 'fontsize', 18)
   ylabel('power norm, ||R^k||', 'fontsize', 18)
   legend('\rho_1 = rho ','\rho_2 = rho/10','\rho_3 = rho/100')
   hold