% a rendition of the Thomas algorithm to solve a tridiagonal system
% Tx = d without pivoting

% the tridiagonal matrix T is represented as 3 vectors:
% a, b, and c for the sub-diagonal, diagonal, and super-diagonal, respectively

% a comparison of the time required for the direct solution using \
% is also made

clear all, close all

% create a random tridiagonal coefficient matrix T and right-hand side
% vector b of size m
m = 3000;

a = rand(m-1,1);
b = rand(m,1);
c = rand(m-1,1);
d = rand(m,1);

% construct matrix to compare timing
T =  diag(a,-1) + diag(b) + diag(c,1);

% direct solution
disp('Using backslash ...')
tic
T\d;
toc
disp(' ')

% solution using Thomas algorithm

% shift the index of a for better intuition
a = [0; a];

disp('Using the Thomas algorithm ...')
% forward substitution
tic
for i=2:m,
  l = a(i)/b(i-1);
  b(i) = b(i) - l*c(i-1);
  d(i) = d(i) - l*d(i-1);
end

% initialize solution vector
x = zeros(m,1);
% back substitution
x(m) = d(m)/b(m);
for i=m-1:-1:1,
  x(i) = (d(i)-c(i)*x(i+1))/b(i);
end
toc

