function specrad= sorspecrad(omega,D,L,U)
%sorspecrad Spectral radius of the SOR matrix.
%
%   specrad = sorspecrad(omega,D,L,U) computes the SOR matrix
%   inv(D+omega*L)*((1-omega)*D-omega*U), where L = diag(A,-1),
%   U = diag(A,1), and D = diag(diag(A)) and returns its
%   spectral radius.
   
BSOR = (D+omega*L)\((1-omega)*D-omega*U);
specrad = max(abs(eig(BSOR)));
