function r=Rho(omega,mu);
% RHO computes spectral radius for SOR
%   r=Rho(omega,mu) computes the spectral radius of the iteration
%   matrix for SOR, given the spectral radius mu of the Jacobi
%   iteration matrix and the relaxation parameter omega

if omega<2/(1+sqrt(1-mu^2))
  r=((omega*mu)^2-2*(omega-1)+sqrt((omega*mu)^2*...
    ((omega*mu)^2-4*(omega-1))))/2;
else
  r=abs(omega-1);
end;  