%=======================================================================
% ch08_thomas.m
%
% Solve a tridiagonal system C*x = b by the Thomas algorithm.
%
% Origin:
%   MATLAB translation of ch08_thomas.f90 by Alexander Godunov, based
%   on Thomas.f90. The elimination multiplier updates both the diagonal
%   and right-hand side, as in the corrected Fortran program.
%
% Example:
%   Eight equations with diagonal 4 and adjacent diagonals -1.
%   All right-hand-side entries are zero except b(8)=16.
%   The example is strictly diagonally dominant, so it does not
%   require row pivoting.
%
% The three stored diagonals, solution, and maximum residual are
% printed. The original C and b are retained for the residual check.
%=======================================================================

clear; clc;

% Row i contains [subdiagonal, diagonal, superdiagonal].
% C(1,1) and C(n,3) are unused.
C = [ 0.0, 4.0, -1.0;
     -1.0, 4.0, -1.0;
     -1.0, 4.0, -1.0;
     -1.0, 4.0, -1.0;
     -1.0, 4.0, -1.0;
     -1.0, 4.0, -1.0;
     -1.0, 4.0, -1.0;
     -1.0, 4.0,  0.0];
n = size(C,1);
b = zeros(n,1);
b(n) = 16.0;

fprintf('Thomas algorithm for a tridiagonal system\n');
fprintf('\nSubdiagonal, diagonal, superdiagonal, and b:\n');
for i = 1:n
    fprintf('%12.6f%12.6f%12.6f%12.6f\n', C(i,:), b(i));
end

% The routine returns changed diagonal and right-hand-side arrays.
[x, C_work, b_work, info] = thomas(C, b); %#ok<ASGLU>
if info ~= 0
    if info < 0
        error('Invalid tridiagonal matrix or right-hand-side size.');
    end
    error('Zero pivot at row %d. This no-pivot method cannot solve it.', info);
end

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

% Check C*x-b row by row with the original coefficients.
max_residual = 0.0;
for i = 1:n
    row_residual = C(i,2)*x(i) - b(i);
    if i > 1
        row_residual = row_residual + C(i,1)*x(i-1);
    end
    if i < n
        row_residual = row_residual + C(i,3)*x(i+1);
    end
    max_residual = max(max_residual, abs(row_residual));
end
fprintf('\nMaximum residual |C*x-b| = %.4e\n', max_residual);


function [x, C, b, info] = thomas(C, b)
%=======================================================================
% Solve a tridiagonal system by elimination without row pivoting.
%
% Input:
%   C - n-by-3 matrix: subdiagonal, diagonal, superdiagonal by row
%       C(1,1) and C(n,3) are unused
%   b - right-hand side
%
% Output:
%   x    - solution if info = 0
%   C    - modified diagonal in column 2
%   b    - modified right-hand side
%   info - 0 on success; -1 for invalid sizes;
%          positive row number of an exact zero pivot otherwise
%
% At row i, coeff = C(i,1)/C(i-1,2) is the elimination multiplier.
% Use coeff for both the diagonal and right-hand-side updates. The
% storage and arithmetic work grow linearly with the number of rows.
% Without pivoting, even a nonsingular matrix may cause a zero pivot
% or an unstable calculation.
%=======================================================================

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

% Eliminate the subdiagonal; retain the superdiagonal unchanged.
for i = 2:n
    coeff = C(i,1)/C(i-1,2);
    C(i,2) = C(i,2) - coeff*C(i-1,3);
    b(i) = b(i) - coeff*b(i-1);
    if C(i,2) == 0.0
        info = i;
        return
    end
end

% Solve the upper triangular equations from bottom to top.
x(n) = b(n)/C(n,2);
for i = n-1:-1:1
    x(i) = (b(i) - C(i,3)*x(i+1))/C(i,2);
end
end