function testSOR(matrix)
% ------------------------------------------------------
%              test function for SOR
% input: matrix = 1 or 2 for two test matrices
% ------------------------------------------------------

if ~nargin; matrix = 1; end;

fprintf('\n  ====== CAAM 454/554 SOR Testing ======\n')
format1 = 'Omega = %5.2f  Iter = %5i  rerr = %8.4e\n'; 
format2a = ': Elapsed time = %g\n\n'; 

N = input(' N = ');
Range = input(' omega range = ');
Range = unique(max(0.0001,min(1.999,Range)));
do_plot = numel(Range) > 2;

switch matrix
    case 1; A = gallery('poisson',N);
    case 2; A = gallery('wathen',N,N);
    otherwise; error('matrix must be 1 or 2');
end
n = size(A,1);
if matrix == 1, xs = min(1,numgrid('H',N)); xs = xs(:); end
if matrix == 2, xs = ones(n,1); end
b = A*xs;
tol = 1e-6;

Str = [' Poisson '; ' Wathen  '];
fprintf(['\n\t---' Str(matrix,:) 'Matrix, n = %i ---\n'],n)

if exist('zSOR','file')
    fprintf('\nRun zSOR:\n')
    solver = 'zSOR';
    format2 = [solver format2a];
    if do_plot; figure(1); end
    run_solver;
end

if exist('mySOR','file')
    fprintf('\nRun mySOR:\n')  
    solver = 'mySOR';
    format2 = [solver format2a]; 
    if do_plot; figure(2); end
    run_solver;
end


%%%%%%%%%%%%%%%%%%%
function run_solver

lR = length(Range);
IT = zeros(lR,1);
k = 1; RE = IT; t0 = tic;
for omega = Range
   [x, iter] = feval(solver,A,b,omega,5*n,tol/lR);
   rerr = norm(x-xs)/norm(xs);
   IT(k) = iter; RE(k) = rerr;
   fprintf(format1,omega,iter,rerr);
   k = k + 1;
end
fprintf(format2,toc(t0));

% ------ plotting ------
if do_plot
    plot(Range, IT, Range, IT, '*');
    title([solver ': Omega vs. Iterarion'])
    ylabel('Iteration');  xlabel('Omega')
    [mint, at] = min(IT); bestO = Range(at);
    str = ['\omega = ' num2str(bestO)];
    str = [str ' requires ' num2str(mint) ' iterations'];
    a = axis; t = 0.9;
    xp = t*a(1) + (1-t)*a(2);
    yp = (1-t)*a(3) + t*a(4);
    text(xp,yp,str,'fontsize',14);
    set(gca,'fontsize',14)
end

if matrix == 1 && not(do_plot) && N < 100
    figure(3)
    surfc(reshape(x,[N N]))
    title('Poisson Equaton Solution in 2D')
    set(gca,'fontsize',16)
end

end %run_solver

end %main
