function woptimal = optomega(A)
%OPTOMEGA Estimate the optimal omega for use with sor by
%computing the iteration matrix.

D = diag(diag(A));
L = tril(A,-1);
U = triu(A,1);

%the optimum value of omega is in the interval 0 < omega < 2
omegavals = 0:.01:2.0;
specrad = zeros(length(omegavals),1);
for i=1:length(omegavals)
   %compute the spectral radius of BSOR(omega) and store it in specrad
   specrad(i) = sorspecrad(omegavals(i),D,L,U);
end

plot(omegavals, specrad);
title('Spectral radius as a function of omega');
xlabel('omega');
ylabel('Spectral radius');

%find the minimum of the spectral radius as a function of omega
[~, i] = min(specrad);
woptimal = omegavals(i);
