function x = Algorithm_21_TMFS(a,b,c,d,n)
% ************************************************************************
% ************************************************************************
% THIS IS A U.S. GOVERNMENT WORK.
% NO COPYRIGHT PROTECTION IS AVAILABLE BY 17 U.S.C. 105
% APPROVED FOR PUBLIC RELEASE, DISTRIBUTION IS UNLIMITED
% ************************************************************************
% THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
% WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY
% MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS",
% NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY
% MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF
% USING THE SOFTWARE LIES WITH THE PARTY DOING SO.
% ************************************************************************
% ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE
% ABOVE STATEMENT.
% ************************************************************************
% ************************************************************************
%
% Algorithm 2.1 Tridiagonal Modified Forward Substitution (TMFS)
%
% The TMFS Algorithm solves the tri-diagonal linear system A*x = d
% where the vectors [a,b,c] represent the sub-, main-, and super-diagonals,
% respectively. It is fast because it skips the typical Gaussian elimination
% step and begins a modified forward substitution method immediately.
% It decomposes the original problem into two subproblems. The first subproblem
% uses the actual RHS vector and the second subproblem uses a basis vector.
% Then a scaling constant is determined using the last equation in each subproblem.
% Then superposition is used to construct the final solution x = A^-1*d.
%
% RESTRICTIONS
% All terms in the super-diagonal (c vector) must be non-zero.
% The system matrix must be invertible.
%
% OUTPUT
% x(n) solution vector, x = A^-1*d
%
% INPUTS
% n system dimension
% a(n) sub-diagonal vector
% b(n) main diagonal vector
% c(n) super-diagonal vector
% d(n) RHS vector
%
% Note that a(1) and c(n) is unused by the algorithm.
%
% ***************************************************************************
% Submitted as supplementary material to the publication
% "An Asynchronous Direct Solver for Banded Linear Systems"
% Numerical Algorithms
%
% Michael A. Jandron, Anthony A. Ruffa
% Naval Undersea Warfare Center, Newport, RI 02893, USA
% email: michael.jandron@navy.mil
%
% James Baglama
% Department of Mathematics, University of Rhode Island, Kingston, RI 02881, USA
% ***************************************************************************

x = zeros(n,1);
u = zeros(n,1); % space for internal forward sub null vector
v = zeros(n,1); % space for internal forward sub actual vector

ebar = zeros(n,1); % space for n basis vector

% FIRST TASK
% Compute solution for null RHS vector using forward substitution
[u,alpha1] = TMFS(u,a,b,c,d,n);

% Compute solution for actual RHS vector using forward substitution
ebar(n) = 1;
[v,alpha2] = TMFS(v,a,b,c,ebar,n);

% Solve for extra variable using two solution vectors
alpha = alpha1 / (1.0 - alpha2);

% Solve for x using superposition of two solution vectors
for i = 1:n
x(i) = u(i) + alpha * v(i);
end

function [x,alpha] = TMFS(x,a,b,c,d,n)
x(1) = 1.0;
x(2) = (d(1) - b(1)*x(1)) / c(1);
for i = 3:n
x(i) = (d(i-1) - a(i-1)*x(i-2) - b(i-1)*x(i-1)) / c(i-1);
end
alpha = d(n) - a(n)*x(n-1) - b(n)*x(n);