%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This code demonstrates how to call symtrisolve.m for
% solving symmetric tridiagonal linear systems Tx = y.
% The symtrisolve.m code computes the LBL^T factorization of a 
% symmetric tridiagonal matrix using the modified Bunch-Marcia 
% pivoting strategy.
%
% Author: Roummel F. Marcia
% Date:   July 10, 2014
% 
% Based on the paper, "A simplified pivoting strategy for 
% symmetric tridiagonal matrices", by James R. Bunch and
% Roummel F. Marcia, Numer. Linear Algebra Appl., vol. 13, 
% p. 865-867 (2006).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear
format compact

% Matrix size
n         = 10000;

% Create random symmetric tridiagonal matrix
% TD is a vector corresponding to the diagonal entries
% TS is a vector corresponding to the sub- and super-diagonal entries

TD        = randn(n,1);
TS        = randn(n-1,1);

% Create random solution x 
x         = randn(n,1);

% Create RHS y 
y         = zeros(n,1);
y(1)      = TD(1)*x(1) + TS(1)*x(2);
for j = 2:n-1
    y(j)  = TS(j-1)*x(j-1) + TD(j)*x(j) + TS(j)*x(j+1);
end
y(n)      = TS(n-1)*x(n-1) + TD(n)*x(n);

% Solve using symmetric tridiagonal solver
x0        = symtrisolve(TD, TS, y);

% Verify that the solution is accurate by computing the
% relative residual norm
y0        = zeros(n,1);
y0(1)     = TD(1)*x0(1) + TS(1)*x0(2);
for j = 2:n-1
    y0(j) = TS(j-1)*x0(j-1) + TD(j)*x0(j) + TS(j)*x0(j+1);
end
y0(n)     = TS(n-1)*x0(n-1) + TD(n)*x0(n);

fprintf('\nProblem size           = %12d\n', n)
fprintf('Relative residual norm = %12.4e\n\n', norm(y-y0)/norm(y))



