function Algorithm_31_BMFS_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.                                                        
%  ************************************************************************
%  ************************************************************************
%
%  Algorithm 3.1 Banded Modified Forward Substitution (BMFS)
%
%  Simple driver function to demonstrate the accuracy of the BMFS algorithm
%  compared to the built-in backslash solver. It demonstrates that the
%  algorithm is able to solve this linear system with accuracy on order
%  of the backslash for very large problems.
%
%  ************************************************************************
%  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, close all

    % generate test linear system
    n = 100; % system dimension, n x n
    
%     p=n; q=n;
%     LHS = rand(n,n);
    
    p = 15;
    q = 20;
    LHS = zeros(n,n);
    for i = 1:n
        for j=1:n
            if i==j+p || j==i+q
                LHS(i,j) = -1;
            elseif i==j
                LHS(i,j) = 2;
            end
        end
    end
    RHS = zeros(n,1);
    RHS(1:n) = 1;
    
    figure, spy(LHS)
        
    % BACKSLASH
    SolutionBackSlash   = callBackSlash(LHS,RHS);
    ErrorBackSlash      = RHS - LHS * SolutionBackSlash;
    L2ErrorBackSlash    = norm(ErrorBackSlash);
    fprintf('Backslash: L2 Error is %e\n',L2ErrorBackSlash)

    % BMFS ALGORITHM
    UseScaling = 1;
    Solution_BMFS        = call_BMFS_Algorithm(LHS,RHS,UseScaling,p,q);
    Error_BMFS           = RHS - LHS * Solution_BMFS;
    L2Error_BMFS         = norm(Error_BMFS);
    fprintf('BMFS Algorithm: L2 Error is %e\n',L2Error_BMFS);
    

    x = linspace(1,n,n);

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

    % plot error
    figure;
    stem(x,ErrorBackSlash,'ko','linew',2), hold all
    stem(x,Error_BMFS,'o','linew',2)
    legend('Backslash','BMFS')
    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_BMFS_Algorithm(a,b,c,p,q)

    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

    % call solver
    x = Algorithm_31_BMFS(a,b,p,q);
    
    if c
        % unscale solution
        x = D1^(-1)*x;
    end




