%=======================================================================
% ch08_sor.m
%
% Solve A*x = b by successive over-relaxation (SOR).
%
% Origin:
%   MATLAB translation of ch08_sor.f90 by Alexander Godunov, based on
%   GS_sor.f90. The sequential SOR update and completed-sweep residual
%   test are retained.
%
% Example:
%   A = [8 2 4; 2 6 1; 1 1 8], b = [2; 6; 4].
%   Exact solution: x = [-0.2; 1; 0.4].
%
% Change omega and the initial guess to study convergence. SOR does
% not converge for every matrix or every choice of omega.
%=======================================================================

clear; clc;

A = [8.0, 2.0, 4.0;
     2.0, 6.0, 1.0;
     1.0, 1.0, 8.0];
b = [2.0; 6.0; 4.0];
x = zeros(3,1);      % Initial guess
omega = 1.2;        % Relaxation parameter
eps = 1.0e-8;       % Absolute residual tolerance
max_iter = 100;
n = size(A,1);

fprintf('Successive over-relaxation (SOR)\n');
fprintf('\nMatrix A and right-hand side b:\n');
for i = 1:n
    fprintf('%12.5f%12.5f%12.5f   %12.5f\n', A(i,:), b(i));
end
fprintf('\nInitial guess x:');
fprintf('%12.6f', x);
fprintf('\nRelaxation parameter omega = %.4f\n', omega);
fprintf('Residual tolerance         = %.4e\n', eps);

[x, iterations, residual, info] = sor(A, b, x, omega, eps, max_iter);
if info == 1
    error('Invalid input: check dimensions, omega, tolerance, and diagonal.');
end

fprintf('\nIterations = %d\nSolution x:\n', iterations);
for i = 1:n
    fprintf('x(%d) = %12.6f\n', i, x(i));
end
fprintf('Maximum residual |A*x-b| = %.4e\n', residual);
if info == 2
    error('Iteration limit reached before the tolerance was met.');
end


function [x, iterations, residual, info] = sor(A, b, x, omega, eps, max_iter)
%=======================================================================
% Solve A*x = b by sequential successive over-relaxation sweeps.
%
% Input:
%   A        - square coefficient matrix, left unchanged
%   b        - right-hand side, left unchanged
%   x        - initial guess, overwritten by the final estimate
%   omega    - relaxation parameter, 0 < omega < 2
%   eps      - requested maximum absolute residual (positive)
%   max_iter - maximum number of complete sweeps (positive)
%
% Output:
%   iterations - number of complete sweeps
%   residual   - max_i |sum_j A(i,j)*x(j) - b(i)|
%   info       - 0 if residual <= eps; 1 for invalid input;
%                2 if the iteration limit was reached
%
% During a sweep, x(j) for j < i has its new value, while x(j) for
% j > i still has its value from the preceding sweep. At omega = 1,
% the update reduces to Gauss-Seidel. The residual is checked only
% after a complete sweep, using the same criterion as the Fortran.
%=======================================================================

n = size(A,1);
iterations = 0;
residual = realmax;
info = 1;
if n < 1 || size(A,2) ~= n || numel(b) ~= n || numel(x) ~= n || ...
        max_iter < 1 || eps <= 0 || omega <= 0 || omega >= 2
    return
end
b = b(:);
x = x(:);

% A zero diagonal makes the update undefined. Strict diagonal
% dominance is checked as a warning, not imposed as a requirement.
for i = 1:n
    if A(i,i) == 0.0
        return
    end
end
for i = 1:n
    off_diagonal = 0.0;
    for j = 1:n
        if j ~= i
            off_diagonal = off_diagonal + abs(A(i,j));
        end
    end
    if abs(A(i,i)) <= off_diagonal
        fprintf('Warning: A is not strictly diagonally dominant.\n');
        break
    end
end

% The supplied guess may already solve the system.
residual = maximum_residual(A, b, x);
if residual <= eps
    info = 0;
    return
end

% Update x(i) immediately so later rows use the new values.
for k = 1:max_iter
    for i = 1:n
        delta = b(i);
        for j = 1:n
            delta = delta - A(i,j)*x(j);
        end
        x(i) = x(i) + omega*delta/A(i,i);
    end
    iterations = k;
    residual = maximum_residual(A, b, x);
    if residual <= eps
        info = 0;
        return
    end
end
info = 2;
end


function residual = maximum_residual(A, b, x)
%-----------------------------------------------------------------------
% Evaluate the original system at the completed iterate.
%-----------------------------------------------------------------------
n = size(A,1);
residual = 0.0;
for i = 1:n
    row_sum = -b(i);
    for j = 1:n
        row_sum = row_sum + A(i,j)*x(j);
    end
    residual = max(residual, abs(row_sum));
end
end