function Algorithm_22_TMFBS_Driver
%  ************************************************************************
%  ************************************************************************
%  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.                                                        
%  ************************************************************************
%  ************************************************************************
%
%  Driver for Algorithm 2.2 Tridiagonal Modified Forward and Backward Substitution (TMFBS)
%
%  Simple driver function to demonstrate the accuracy of the TMFBS algorithm
%  compared to the built-in backslash solver.
%
%  This version assembles a basic toeplitz problem:
%
%   [  2 -1  0 ] [ x1 ]   [ 1 ]
%   [ -1  2 -1 ] [ x2 ] = [ 1 ]
%   [  0 -1  2 ] [ x3 ]   [ 1 ]
%
%  with an user defined dimension, n x n
% 
%  It demonstrates that each algorithm is able to solve this linear system
%  with accuracy on order of the backslash for very large problems.  If all
%  entries are integers, and the RHS terms are one, then no round off will occur.
%  It is written in sparse format for memory management, but no optimization has
%  been performed on this benchmark routine to eliminate redundant variables.  
%
%  ***************************************************************************
%  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
%  ***************************************************************************

    clc, clear all, close all

    % generate test linear system
    n = 100; % system dimension, n x n
    
    % use sparse storage format to reduce memory
    LHS = gallery('tridiag',n,-1,2,-1);
    RHS = zeros(n,1);
    RHS(1:n) = 1;
    
    % BACKSLASH
    SolutionBackSlash   = callBackSlash(LHS,RHS);
    ErrorBackSlash      = RHS - LHS * SolutionBackSlash;
    L2ErrorBackSlash    = norm(ErrorBackSlash);
    fprintf('Backslash: L2 Error is %e\n',L2ErrorBackSlash)
    
    % TMFBS ALGORITHM
    UseScaling = 0;
    Solution_TMFBS       = call_TMFBS_Algorithm(LHS,RHS,UseScaling);
    Error_TMFBS          = RHS - LHS * Solution_TMFBS;
    L2Error_TMFBS        = norm(Error_TMFBS);
    fprintf('TMFBS Algorithm: L2 Error is %e\n',L2Error_TMFBS);

    x = linspace(1,n,n);

    % plot solution
    figure(1);
    plot(x,SolutionBackSlash,'k-o','linew',2), hold all
    plot(x,Solution_TMFBS,'o','linew',2)
    legend('Backslash','MFBS')
    xlabel('Unknown (k)','fontsize',14)
    ylabel('Solution (x)','fontsize',14)
    set(gca,'fontsize',14)
    grid on
    hold off

    % plot error
    figure(2);
    stem(x,ErrorBackSlash,'ko','linew',2), hold all
    stem(x,Error_TMFBS,'o','linew',2)
    legend('Backslash','MFBS')
    xlabel('Unknown (k)','fontsize',14)
    ylabel('Error (b-Ax)','fontsize',14)
    set(gca,'fontsize',14)
    grid on
    hold off
    
    

function x = callBackSlash(a,b)

    x = a \ b;


function x = call_TMFBS_Algorithm(a,b,c)

    n = length(b);
    
    if c;
        % compute scaling factors
        sc = max(max(abs(a)));
        D1 = sparse(sc*diag(ones(n,1),0));
        a = D1^(-1) * a;
    end

    % prepare arguments
    arg1 = [0;full(diag(a,-1))];
    arg2 = full(diag(a,0));
    arg3 = [full(diag(a,+1));0];
    arg4 = b;
    arg5 = n;

    % call solver
    x = Algorithm_22_TMFBS(arg1,arg2,arg3,arg4,arg5);
    
    if c;
        % unscale solution
        x = D1^(-1)*x;
    end


