%=======================================================================
% ch08_gaussian_elimination.m
%
% Solve A*x = b by Gaussian elimination with scaled row pivoting.
%
% Origin:
%   MATLAB translation of ch08_gaussian_elimination.f90 by Alexander
%   Godunov, based on Gauss_2.f90. The active-row scaling, pivot
%   selection, elimination, and back substitution are retained.
%
% Example:
%   A = [0 1 2; 2 1 4; 2 4 6], b = [4; 3; 7].
%   Exact solution: x = [-2.5; 0; 2]. A row interchange is necessary.
%
% The program prints the original system, the upper triangular system,
% the solution, and the residual checked with the original A and b.
% A small residual alone does not establish accuracy for an
% ill-conditioned system.
%=======================================================================

clear; clc;

A = [0.0, 1.0, 2.0;
     2.0, 1.0, 4.0;
     2.0, 4.0, 6.0];
b = [4.0; 3.0; 7.0];
n = size(A,1);

fprintf('Gaussian elimination with scaled row pivoting\n');
fprintf('\nOriginal matrix 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

% The routine changes its matrix and right-hand side. Keep A and b
% unchanged for the independent check of the solution.
[x, A_work, b_work, info] = gauss_scaled(A, b);
if info ~= 0
    if info < 0
        error('Invalid matrix or right-hand-side size.');
    end
    error('No unique solution: zero pivot at step %d.', info);
end

fprintf('\nMatrix A and right-hand side b after elimination:\n');
for i = 1:n
    fprintf('%12.5f%12.5f%12.5f   %12.5f\n', A_work(i,:), b_work(i));
end

fprintf('\nSolution x:\n');
for i = 1:n
    fprintf('x(%d) = %12.6f\n', i, x(i));
end

% Evaluate the original equations using explicit loops.
max_residual = 0.0;
for i = 1:n
    residual = -b(i);
    for j = 1:n
        residual = residual + A(i,j)*x(j);
    end
    max_residual = max(max_residual, abs(residual));
end
fprintf('\nMaximum residual |A*x-b| = %.4e\n', max_residual);


function [x, A, b, info] = gauss_scaled(A, b)
%=======================================================================
% Solve A*x = b by forward elimination and back substitution.
%
% Input:
%   A    - square coefficient matrix
%   b    - right-hand side
%
% Output:
%   x    - solution if info = 0
%   A    - upper triangular matrix after elimination
%   b    - transformed right-hand side
%   info - 0 on success; -1 for invalid sizes;
%          positive pivot step if an exact zero pivot is found
%
% At step k, compute the scale of each active row from the absolute
% coefficients in columns k:n. Choose the largest |A(i,k)|/scale,
% swap the active parts of the rows, and eliminate below the pivot.
% These are the same steps used in the Fortran program. A very small
% nonzero pivot may still indicate an ill-conditioned system.
%=======================================================================

n = size(A,1);
x = zeros(n,1);
info = 0;
if n < 1 || size(A,2) ~= n || numel(b) ~= n
    info = -1;
    return
end
b = b(:);

% Forward elimination: select a scaled pivot at each active column.
for k = 1:n-1
    best = 0.0;
    pivot_row = k;
    for i = k:n
        scale = 0.0;
        for j = k:n
            scale = max(scale, abs(A(i,j)));
        end
        if scale == 0.0
            continue
        end
        score = abs(A(i,k))/scale;
        if score > best
            best = score;
            pivot_row = i;
        end
    end

    if best == 0.0
        info = k;
        return
    end

    % Exchange active coefficients and the matching right sides.
    if pivot_row ~= k
        for j = k:n
            store = A(k,j);
            A(k,j) = A(pivot_row,j);
            A(pivot_row,j) = store;
        end
        store = b(k);
        b(k) = b(pivot_row);
        b(pivot_row) = store;
    end

    % Subtract a multiple of the pivot equation from each lower row.
    for i = k+1:n
        factor = A(i,k)/A(k,k);
        A(i,k) = 0.0;
        b(i) = b(i) - factor*b(k);
        for j = k+1:n
            A(i,j) = A(i,j) - factor*A(k,j);
        end
    end
end

% The last pivot was not selected in the loop above.
if A(n,n) == 0.0
    info = n;
    return
end

% Back substitution through the upper triangular system.
for i = n:-1:1
    row_sum = 0.0;
    for j = i+1:n
        row_sum = row_sum + A(i,j)*x(j);
    end
    x(i) = (b(i) - row_sum)/A(i,i);
end
end