function rho = powerIter(n,OP,tol)
% Compute approximate spectral radius of linear operator OP.
%
% This version can handle the case when the dominant eigenvalue is complex.
% It is susceptible to slow convergence when there are several eigenvalues
% with modulus close to RHO. Accordingly, this should only be used to
% estimate the dynamic range.

%  Copyright 2023 The MathWorks, Inc.
RNG = RandStream('mt19937ar','Seed',0);
v0 = RNG.randn(n,1);
v0 = v0/norm(v0);
w0 = OP(v0);
rho = norm(w0);
if rho==0
   return
elseif ~isfinite(rho)
   rho = Inf;  return
end

% Power iterations
% Maintain OP(V) = W where V has normalized columns.
v = [v0 w0/rho];
w = [w0 zeros(n,1)];
j = [2 1];
MAXITER = 25;
log_rho = zeros(1,MAXITER);  % stores log(rho)
for ct=1:MAXITER
   [q,r] = qr(v,0);
   % Apply OP
   u = v(:,j(1));
   OPu = OP(u);
   g = norm(OPu);
   if g==0
      rho = 0;  return
   end
   w(:,j(1)) = OPu;
   v(:,j(2)) = OPu/g;
   % Update spectral radius estimate (Rayleigh quotient)
   if abs(r(2,2))>0.1
      % This computes max(eig(Q'*OP(Q))) using OP(V) = W.
      % Needed for convergence when dominant eigenvalue is complex.
      %eig((q'*w)/r)
      log_rho(ct) = log(max(abs(eig((q'*w)/r))));
   else
      log_rho(ct) = log(abs(u'*OPu));
   end
   % Check for convergence (max variation < TOL over last 3 iterations)
   if ct>5 && max(abs(log_rho(ct)-log_rho(ct-2:ct-1)))<tol
      rho = exp(log_rho(ct));   return
   end
   % Circulate buffer
   j = j([2 1]);
end

% LOG averaging tends to smoothen oscillations when failed to converge
NAVG = 6;
log_avg = sum(log_rho(ct-NAVG+1:ct))/NAVG;
rho = exp(log_avg);
%fprintf('POWER IT = %d\n',ct)
