function x = Algorithm_22_TMFBS(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.2 Tridiagonal Modified Forward and Backward Substitution (TMFBS)
%
%  The TMFBS 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 four subproblems.  It is an
%  extension of the Tridiagonal Modified Forward Substitition (TMFS) Algorithm 2.1.
%
%  RESTRICTIONS
%  c(1:m-1) ~= 0 where m = round(n/2)
%  a(m+1:n) ~= 0 where m = round(n/2)
%  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) are 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); % solution vector
    u = zeros(n,1); % space for internal forward sub vector
    v = zeros(n,1); % space for internal backward sub vector 
    e = zeros(n,1); % space for null vector (just to reuse function)
    
    m = round(n/2);
   
    % FIRST TASK
    % Compute solution for actual RHS vector using forward substitution
    [u,alpha1] = TMFS(u,a,b,c,d,m);

    % Compute solution for actual RHS vector using backward substitution
    [u,alpha2] = TMBS(u,a,b,c,d,m,n);
    
    % Compute solution for null RHS vector using forward substitution
    [v,alpha3] = TMFS(v,a,b,c,e,m);

    % Compute solution for null RHS vector using backward substitution
    [v,alpha4] = TMBS(v,a,b,c,e,m,n);
    
    % Note:
    % alpha1 = u(m) from forward substituion
    % alpha2 = u(m) from backward substituion 
    % alpha3 = v(m) from forward substituion
    % alpha4 = v(m) from backward substituion

   
    % SECOND TASK
    % from null forward/backward subsitition
    alpha5 = a(m)*(v(m-1)/alpha3) + b(m) + c(m)*(v(m+1)/alpha4);
    
    % from actual forward/backward substitution
    alpha6 = a(m)*u(m-1) + b(m)*alpha2 + c(m)*u(m+1);
    
    % from mixed cases
	alpha7 = a(m)*(v(m-1)/alpha3*alpha2) + b(m)*alpha1 + c(m)*(v(m+1)/alpha4*alpha1);
    
    % correction factor
    alpha = (d(m) - alpha6 - alpha7) / alpha5;
    
    % new scaling factors
    beta1 = (alpha1 + alpha) / alpha4;
    beta2 = (alpha2 + alpha) / alpha3;
    
    % Compute final solution x by scaling 
    for i = 1:m-1;
        x(i) = u(i) + v(i)*beta1;
    end
    x(m) = alpha2 + alpha4*beta2;
    for i = m+1:n;
        x(i) = u(i) + v(i)*beta2;
    end

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

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




